# Intermezzo: permutations

As I haven’t progressed very much with my WebGL apps due to a variety of reasons/excuses (lack of willpower, a crash course on nuclear reactors engineering since the Fukushima nuclear crisis started, etc.), I decided to make a quick post showing a solution to a simple problem. In other words, the main purpose of this post is just to avoid leaving a month without posts 馃榾

##### The problem

Given an arithmetic expression without parentheses, insert them to get a specified result. For example, if the expression is 2+2/2 and the desired result is 2, there is a single solution: (2+2)/2. On the other hand, if the desired result is 0 there is no solution.

It’s easy to see the answers in simple cases like the last one. But it’s not so easy to see if we can get a result of 18 by adding parentheses to 18-6/2+4-1*2… so let’s write code to do the hard work for us.

##### A solution

The parentheses only change the order in which the different operators apply and the effects of any operation order can be obtained with some parenthesization. So, if we search for a permutation of the order in which the operators are applied that gives us the desired result, we are getting a solution for the original problem.

This new problem is not exactly the same as the original one by two reasons that can be exemplified as follows:

• Given the expression 2+2*2+2, the orders of operation 1-3-2 and 3-1-2 are associated to the same parenthesized expression: (2+2)*(2+2).
• The parenthesizations (2+2)+2, 2+(2+2) and 2+2+2 are equivalent due to the associativity of addition.

But, as this only gives us some redundant solutions that can be easily eliminated with some post-processing, it’s not a problem from the correctness point of view. Being able to remove this redundant solutions without needing to enumerate them would be a great thing performance wise, but it doesn’t seem easy and it’s not required for small expressions (10! < 107).

###### Permutations in lexicographic order

To avoid accidentally skipping some permutations, it’s useful to follow a specific order when enumerating them. In this section we will see how we can go through all the permutations of an array [0, 1, 2, …, n – 1] in lexicographic order.

Given some permutation of the array, that we will call a, the next one will probably share some prefix and will differ in the order of the last elements. As we cannot get a lexicographically greater permutation by permuting a suffix that is in reverse order (starting from the greater elements), we need to permute a suffix that is not in reverse order. Then we will make our first step to search for this suffix:

Find the largest index k such that a[k] < a[k + 1].

The specification of the largest index is equivalent to search for the minimal suffix to permute, something required to get a permutation that is lexicographically immediate. If this suffix cannot be found, we are at the last permutation.

Now we need to decide what to permute. It’s clear that we cannot get a greater permutation just by permuting the suffix a[k+1 .. n-1], as they are in reverse order. So we need to replace a[k] by some element and the obvious choice is to search for the immediately greater element in the suffix or, as the suffix a[k+1 .. n-1] is in reverse order:

Find the largest index l such that a[k] < a[l].

The lexicographically smallest permutation starting with the original a[l] can be obtained by appending all the other elements of the original suffix in ascending order. But, as the suffix is in descending order and l is the largest index such that a[k] < a[l], we can get the same result in an easier way:

Swap a[k] with a[l].
Reverse the sequence from a[k + 1] up to and including the final element a[n].

The final code in JS can be seen here.

###### Evaluating an arithmetic expression with different orders of operations

I started using an array whose members were alternately numbers and operators. The main problem is that, while it’s easy to go from an operator to the adjacent numbers in an array, this is not what we want to do. The desired behavior is to go from an operator to the operands, and they can be the result of a sequence of previous operations.

To solve this problem in a moderately efficient way, I adapted the parent node concept from the disjoint-set data structure. Each operation adds the result to the operator node and puts this node as the parent node of both operands, allowing efficient access to the value of the whole operand from any node that is included in it. The code that does this operations can be seen at the SVN repository of this blog.

##### Testing the solution

The whole solution can be tested at the same SVN repository, but let’s see some specific examples:

Input expression: 3+2*2-1*0.5
Input result 1: 3
Output: (3+2*2-1)*0.5
Input result 2: 6
Output: 3+2*(2-1*0.5)

Input expression: 18-6/2+4-1*2
Input result: 18
Output 1: ((18-6)/2+4-1)*2
Output 2: 18-(6/(2+4)-1)*2

# GPGPU with WebGL: solving Laplace’s equation

This is the first post in what will hopefully be a series of posts exploring how to use WebGL to do GPGPU (General-purpose computing on graphics processing units). In this installment we will solve a partial differential equation using WebGL, the Laplace’s equation more specifically.

##### Discretizing the Laplace’s equation

The Laplace’s equation, $\nabla^2 \phi = 0$, is one of the most ubiquitous partial differential equations in physics. It appears in lot of areas, including electrostatics, heat conduction and fluid flow.

To get a numerical solution of a differential equation, the first step is to replace the continuous domain by a lattice and the differential operators with their discrete versions. In our case, we just have to replace the Laplacian by its discrete version:

$\displaystyle \nabla^2 \phi(x) = 0 \rightarrow \frac{1}{h^2}\left(\phi_{i-1\,j} + \phi_{i+1\,j} + \phi_{i\,j-1} + \phi_{i\,j+1} - 4\phi_{i\,j}\right) = 0$,

where $h$ is the grid size.

If we apply this equation at all internal points of the lattice (the external points must retain fixed values if we use Dirichlet boundary conditions) we get a big system of linear equations whose solution will give a numerical approximation to a solution of the Laplace’s equation. Of the various methods to solve big linear systems, the Jacobi relaxation method seems the best fit to shaders, because it applies the same expression at every lattice point and doesn’t have dependencies between computations. Applying this method to our linear system, we get the following expression for the iteration:

$\displaystyle \phi_{i\,j}^{(k+1)} = \frac{1}{4}\left(\phi_{i-1\,j}^{(k)} + \phi_{i+1\,j}^{(k)} + \phi_{i\,j-1}^{(k)} + \phi_{i\,j+1}^{(k)}\right)$,

where $k$ is a step index.

##### Solving the discretized problem using WebGL shaders

If we use a texture to represent the domain and a fragment shader to do the Jacobi relaxation steps, the shader will follow this general pseudocode:

1. Check if this fragment is a boundary point. If it’s one, return the previous value of this point.
2. Get the four nearest neighbors’ values.
3. Return the average of their values.

To flesh out this pseudocode, we need to define a specific representation for the discretized domain. Taking into account that the currently available WebGL versions don’t support floating point textures, we can use 32 bits RGBA fragments and do the following mapping:

R: Higher byte of $\phi$.
G: Lower byte of $\phi$.
B: Unused.
A: 1 if it’s a boundary value, 0 otherwise.

Most of the code is straightforward, but doing the multiprecision arithmetic is tricky, as the quantities we are working with behave as floating point numbers in the shaders but are stored as integers. More specifically, the color numbers in the normal range, [0.0, 1.0], are multiplied by 255 and rounded to the nearest byte value when stored at the target texture.

My first idea was to start by reconstructing the floating point numbers for each input value, do the required operations with the floating numbers and convert the floating point numbers to color components that can be reliably stored (without losing precision). This gives us the following pseudocode for the iteration shader:

// wc is the color to the "west", ec is the color to the "east", ...
float w_val = wc.r + wc.g / 255.0;
float e_val = ec.r + ec.g / 255.0;
// ...
float val = (w_val + e_val + n_val + s_val) / 4.0;
float hi = val - mod(val, 1.0 / 255.0);
float lo = (val - hi) * 255.0;
fragmentColor = vec4(hi, lo, 0.0, 0.0);


The reason why we multiply by 255 in place of 256 is that we need val_lo to keep track of the part of val that will be lost when we store it as a color component. As each byte value of a discrete color component will be associated with a range of size 1/255 in its continuous counterpart, we need to use the “low byte” to store the position of the continuous component within that range.

Simplifying the code to avoid redundant operations, we get:

float val = (wc.r + ec.r + nc.r + sc.r) / 4.0 +
(wc.g + ec.g + nc.g + sc.g) / (4.0 * 255.0);
float hi = val - mod(val, 1.0 / 255.0);
float lo = (val - hi) * 255.0;
fragmentColor = vec4(hi, lo, 0.0, 0.0);


The result of running the full code, implemented in GLSL, is:

Solving the Laplace's equation using a 32x32 grid. Click the picture to see the live solving process (if your browser supports WebGL).

As can be seen, it has quite low resolution but converges fast. But if we just crank up the number of points, the convergence gets slower:

Incompletely converged solution in a 512x512 grid. Click the picture to see a live version.

How can we reconcile these approaches?

##### Multigrid

The basic idea behind multigrid methods is to apply the relaxation method on a hierarchy of increasingly finer discretizations of the problem, using in each step the coarse solution obtained in the previous grid as the “starting guess”. In this mode, the long wavelength parts of the solution (those that converge slowly in the finer grids) are obtained in the first coarse iterations, and the last iterations just add the finer parts of the solution (those that converge relatively easily in the finer grids).

The implementation is quite straightforward, giving us fast convergence and high resolution at the same time:

Multigrid solution using grids from 8x8 to 512x512. Click the picture to see the live version.

##### Conclusions

It’s quite viable to use WebGL to do at least basic GPGPU tasks, though it is, in a certain sense, a step backward in time, as there is no CUDA, floating point textures or any feature that helps when working with non-graphic problems: you are on your own. But with the growing presence of WebGL support in modern browsers, it’s an interesting way of partially accessing the enormous computational power present in modern video cards from any JS application, without requiring the installation of a native application.

In the next posts we will explore other kinds of problem-solving where WebGL can provide a great performance boost.

# Visualizing the Mandelbrot set with JavaScript

#### Fractals and the Mandelbrot set

Even though mathematical structures that are now called fractals were first analysed in the 19th century, the study of self-similar structures got a significant boost when computers able to visualize them started to appear in the 1970s. In fact, the term fractal was coined in 1975 by the late Benoit Mandelbrot to describe structures whose parts are similar to the whole.

The Mandelbrot set is a particular fractal set defined as:

$\displaystyle M = \{ c \in \mathbb{C}: \lim_{n\to\infty} P_c^n(0) \neq \infty \}$,

where $P_c^n(z)$ describes the polynomial resulting of iterating $n$ times the parameterized quadratic polynomial $P_c(z)$:

$\displaystyle P_c(z) = z^2 + c$.

#### More about the Mandelbrot set

It is said that the limit of a complex sequence $z_n$ is infinite if for any positive number $M$ there is an element of the sequence $z_N$ such that its module and the modules of all succeeding sequence elements are greater than $M$. In formal notation:

$\displaystyle \lim_{n\to\infty} z_n = \infty \iff \forall M > 0\,\exists N > 0\,\forall n \geq N: |z_n| \geq N$.

All the points of the complex plane that have a modulus greater than 2 are outside the set. To show that, let’s take a point $c_0$ whose module is strictly greater than 2:

$\displaystyle |c_0| > 2$.

To avoid cluttering the inequalities, let’s make two definitions:

$\displaystyle a_n = |P_{c_0}^n(0)|$ and

$\displaystyle k = |c_0|$.

Then the statement that all points with modulus greater than two are outside the Mandelbrot set is equivalent to

$\displaystyle \lim_{n \to \infty} a_n = \infty$.

It’s clear from the definition that $a_n$ is non-negative. Let’s analyse carefully the relationship between a term and the previous one:

$\displaystyle a_{n+2} = |P_{c_0}^{n+2}(0)| = \left|\left(P_{c_0}^{n+1}(0)\right)^2+c_0\right|$

$\displaystyle a_{n+2} \geq \left|\left|\left(P_{c_0}^{n+1}(0)\right)^2\right|-|c_0|\right|$

$\displaystyle a_{n+2} \geq \left|\left|P_{c_0}^{n+1}(0)\right|^2-k\right|$

$\displaystyle a_{n+2} \geq \left|a_{n+1}^2-k\right|$.

We can now inductively explore the hypothesis $a_{n+2} \geq (1 + \epsilon) a_{n+1}$, with $\epsilon = k - 2$.

Checking the base case:

$\displaystyle a_2 \geq |a_1^2 - k| = |k^2 - k| = |(k-1)k| = (k-1)k = (1+\epsilon)k = (1+\epsilon)a_1$

Applying the strong inductive hypothesis to $a_{n+1}$ we get:

$\displaystyle a_{n+1} \geq (1 + \epsilon)a_n \geq \ldots \geq (1 + \epsilon)^n a_1 = (1 + \epsilon)^n k$

Using this result in the previous inequality:

$\displaystyle a_{n+2} \geq a_{n+1}^2-k \quad (a_{n+1}^2 > a_{n+1} > k)$

$\displaystyle a_{n+2} > a_{n+1}^2-a_{n+1} \quad (a_{n+1} > k)$

$\displaystyle a_{n+2} > a_{n+1}(a_{n+1} - 1)$

$\displaystyle a_{n+2} > a_{n+1}(k - 1) \quad (a_{n+1} > k)$

$\displaystyle a_{n+2} > (1 + \epsilon) a_{n+1} \quad (1 + \epsilon = k - 1)$

This implies that $a_n$ will diverge at least geometrically, with ratio $1+ \epsilon$. Then any point $c$ whose modulus is greater than 2 will be outside the set.

It can be seen by a slight extension of this argument that, if there is some $n$ such that $|P_c^n(0)|$ is greater than 2, $c$ will be outside the Mandelbrot set too.

#### Plotting the Mandelbrot set

The simplest way to plot a Mandelbrot set is to iterate the quadratic polynomial $P_c(z)$ over the points $c$ that are associated with an image in the normal way (via some choice of origin and scale factor). If the iterated point reaches a modulus greater than 2 we know that its module will diverge and that it is outside the set. There isn’t any simple test to know if a point is inside the set, but we can choose a fixed number of iterations or to stop the process when no new point outside the set is detected.

Any point known to be outside the set (reached a modulus greater than 2 when iterated) can be colored white, using black to color the points of the set:

A Mandelbrot set plotted using black for the points inside the set and white for those outside the set. (Click to see a live version.)

We can make prettier visualizations if we use some measure of the “difficulty of divergence” for each point. One natural measure for this quantity is the number of iterations required to reach a given modulus. Combining this number of iterations with a function translating the potentially unbounded number of iterations to a color in the range black – white, we get:

Mandelbrot set in gray scale. The gray level of the points outside the set are related to the number of iterations required to reach a given modulus. (Click to see a live version.)

This visualization is nicer, and shows more clearly the existence of filament-like structures (in fact, the set is connected, but I don’t know how to prove that… 馃榾 ). But we can do better by interpolating the iteration number in a continuous way, giving us a less discrete measure of the “difficulty of divergence”. If we use the algorithm described by Linas Vepstas,

$\mu = n + 1 - \frac{1}{\log 2}\log \log P_c^n(0)$,

being $\mu$ the “continuous” iteration count, the we get this picture:

A representation of the Mandelbrot set using the "continuous iteration count" to colorize the exterior. (Click to see a live version.)

The source code of the visualization software can be seen in the SVN repository of this blog and it can be used to check the self similarity of the fractal (left click zooms in, right click zooms out):

A "copy" of the Mandelbrot set in a filament. Zoom out in the live version to see its location.

# JS recursive puzzle

What happens when we execute the following Javascript code and why does it happen?

(function(f){f(f)})(function(g){g(g)})


In the next post, I will show the answer to this problem and the results of the subset sum collision detection. Stay tuned 馃檪

# Din谩mica de fluidos en Javascript

Inspirado por un ejemplo portado por Oliver Hunt de C a JS (y frustrado porque no se me ocurr铆a como resolver el nivel Ophanim en Manufactoria 馃榾 ), decid铆 portar mi c贸digo de simulaci贸n de fluidos mediante Lattice Boltzmann (LB) a Javascript y agregarle interactividad. La visualizaci贸n es algo primitiva, ya que me concentr茅 en optimizar los c谩lculos sin prestarle mucha atenci贸n a la optimizaci贸n del manejo de canvas; posiblemente mejore ese aspecto en una pr贸xima versi贸n.

El c贸digo presenta algunas limitaciones propias de utilizar una grilla de 64 X 64 como base (el m谩ximo tama帽o con performance “realtime“) y, muy probablemente de mi desconocimiento sobre como optimizar los par谩metros de LB. Es particularmente clara la compresibilidad del fluido y no es muy dif铆cil hacer “explotar” la simulaci贸n (se recupera mediante un reset).

Como WordPress no permite incluir c贸digo Javascript en los posts, para visualizar la simulaci贸n pueden hacer click en la imagen de referencia de los controles.

Referencia de los controles de la simulaci贸n. Hacer click para arrancar la simulaci贸n.

El uso de los controles es muy simple:

• Play/pause: arranca o detiene la simulaci贸n.
• Reset: reinicia la simulaci贸n desde cero.
• Visualization reset: reinicia la visualizaci贸n, creando una nueva serie de part铆culas para marcar el movimiento del fluido (sin alterar el estado del mismo).

Para darle un impulso a una parte del fluido basta con arrastrar sobre el 谩rea de visualizaci贸n.

En un posterior post incluir茅 detalles acerca de como funciona LB y de las optimizaciones realizadas (tengo que leer m谩s, sobre todo respecto a la conexi贸n con los valores f铆sicos 馃榾 ).

# Porqu茅 no se puede ganar siempre

##### Respuesta al post anterior

La t茅cnica mostrada en el post anterior se denomina martingala y es muy antigua. Obviamente, a pesar de lo que puedan decir algunos, no funciona. Pero porqu茅 la simulaci贸n anterior parec铆a mostrar lo contrario?

El motivo principal del fallo de la simulaci贸n es no limitar la cantidad de dinero del jugador disponible para apostar. Como la cantidad que se requiere apostar crece exponencialmente con la cantidad de jugadas, en una cantidad relativamente peque帽a de estas se exceder谩n las reservas del jugador y no podr谩 continuar apostando.

Otro problema de la simulaci贸n es que, incluso una cantidad grande de simulaciones, puede obviar acontecimientos muy improbables. Esto no es problem谩tico cuando los eventos en cuesti贸n no tienen consecuencias excepcionales pero, cuando conducen a un impacto proporcional a su improbabilidad, puede distorsionarse mucho el valor esperado (probar la nueva simulaci贸n Javascript con M = 100000 para N = 1000 y N = 1000000; observar las diferencias en los valores esperados estimados).

Pero puede objetarse (con raz贸n!) que solo demostramos el fallo de una estrategia de apuesta, no de todas las posibles. Por lo tanto, a continuaci贸n, demostraremos que esto se aplica a cualquier estrategia.

##### Porqu茅 los sistemas no funcionan

Como en toda demostraci贸n matem谩tica, necesitaremos hacer ciertas suposiciones sobre el problema en cuesti贸n:

1. El juego consiste en una cierta serie de jugadas, cada una asociada con una variable aleatoria $X_i$. Asumimos que todas las variables $X_i$ tienen la misma distribuci贸n.
2. En la jugada $i$ el jugador apuesta una cantidad $m_i$ a un evento aleatorio $e_i$ (donde tomamos $e_i$ como un subconjunto de lso valores que podr铆a tomar $X_i$).
3. Si el evento ocurre al que el jugador apost贸 ocurre, este recibe el monto apostado multiplicado por una constante $k(e_i)$ dependiente del evento al que apost贸. En caso contrario, pierde lo apostado.
4. $k(e) P(e) \le 1$.
5. Los valores de las variables aleatorias $X_1, X_2, ..., X_{i-1}$ no aportan informaci贸n sobre el valor de $X_i$.

Las primeras 3 suposiciones son simplemente una descripci贸n general de un juego de azar, mientras que 5 nos indica que las anteriores jugadas no nos dan informaci贸n sobre la jugada actual. La suposici贸n 4 parece m谩s compleja, pero esencialmente nos indica que los factores $k$ han sido correctamente elegidos por el casino.

Para ver esto, supongamos que existe un evento $e$ tal que $k(e) P(e) > 1$. Entonces, si elegimos una estrategia de apostar $N$ veces al evento $e$ un monto unitario, nuestro balance esperado ser谩:

$\displaystyle E[b] = E[\sum_{i=1}^N b_i]$

$\displaystyle = \sum_{i=1}^N E[b_i]$

$\displaystyle = \sum_{i=1}^N ((k(e)-1)P(e) - (1-P(e)))$

$\displaystyle = N (k(e)P(e) - 1)$

$\displaystyle > 0$,

donde $b$ es el balance general del jugador y $b_i$ el balance de la jugada $i$.

Eso implicar铆a que podr铆amos ganar simplemente repitiendo una apuesta una y otra vez (puede demostrarse sin mucha dificultad que ni siquiera necesitar铆amos reservas muy grandes, solo de orden $\sqrt(N)$). Como eso claramente no ocurre en los casinos, es obvio que los multiplicadores est谩n elegidos correctamente (de hecho, por ejemplo en la ruleta, est谩n elegidos “a favor” del casino por la existencia del cero).

Busquemos ahora el valor esperado de una estrategia de apuestas $m_i(i, X_1, X_2, ..., X_{i-1}), e_i(i, X_1, X_2, ..., X_{i-1})$:

$\displaystyle E[b] = E[\sum_{i=1}^N b_i]$

$\displaystyle = E[\sum_{i=1}^N b_i]$

$\displaystyle = \sum_{i=1}^N E[b_i]$

$\displaystyle = \sum_{i=1}^N ((k(e_i)-1)P(e_i) - (1-P(e_i)))m_i$

$\displaystyle = \sum_{i=1}^N (k(e_i)P(e_i) - 1)m_i$

siendo $e_i$ y $m_i$ las funciones anteriormente descritas. Pero la suposici贸n 5 nos dice que

$P(e_i|i,X_1=x_1, X_2=x_2, ..., X_{i-1} = x_{i-1}) = P(e_i)$,

es decir que esa informaci贸n no nos modifica la probabilidad de ning煤n evento y, en particular, de $e_i$. Combinando con la suposici贸n 4 llegamos a

$\displaystyle E[b] \le 0$,

siempre que $m_i$ sea positivo, como es razonable.

##### Caminos de escape

Aunque una demostraci贸n fuera perfecta, sus conclusiones nunca ser谩 m谩s fuertes que sus hip贸tesis. En esta secci贸n analizaremos algunas formas en las que las hip贸tesis pueden ser falsas:

mi < 0: Este es esencialmente el m茅todo utilizado por los casinos. Suele ser bastante complejo de seguir, pero muy redituable 馃檪

Fallo en la hip贸tesis 4: Los multiplicadores son f谩ciles de determinar conociendo las probabilidades y generalmente se utilizan probabilidades “nominales” (por ejemplo, suponer equiprobables a los 37 o 38 n煤meros de la ruleta). Pero si alg煤n observador encuentra evidencia que lo mueva a probabilidades distintas de las nominales, puede que existan apuestas con valor esperado positivo.

Fallo en la hip贸tesis 5: Pueden existir correlaciones no esperadas entre las variables aleatorias. Esto es m谩s com煤n en los juegos por computadora, que suelen usar n煤meros pseudoaleatorios y, en algunos casos, tienen errores de implementaci贸n que vuelven predecibles las secuencias generadas.

# Como ganar siempre :-)

Un sistema de apuestas, muy utilizado por lo que puede verse en Google, sigue el siguiente esquema (donde lo aplicamos a una apuesta doble o nada con una moneda “justa”):

2. Apostar a “cara”.
3. En caso de ganar, retirarse.
4. En caso de perder, doblar la apuesta y volver al paso 2.

Analicemos este sistema mediante un mill贸n de simulaciones de la t茅cnica, utilizando 3 lenguajes de programaci贸n: C++ (codepad), Python (codepad, N = 1000) y Javascript.

Como puede observarse, los resultados son esencialmente id茅nticos:

Number of simulations: 1000000
Minimum player balance: 1
Maximum player balance: 1
Average number of rounds before winning: 2.000831

En todas las ejecuciones observamos una ganancia fija de una unidad y que se obtiene en promedio en una peque帽a cantidad de tiradas.

Donde est谩 el problema con este sistema? En el pr贸ximo post la respuesta… 馃榾

# El poder del SAT

SAT es el ejemplo can贸nico de problema NP-completo. Consiste en determinar, dada una f贸rmula booleana en forma de conjunci贸n de disyunciones (donde cada disyunci贸n se llama cla煤sula), si existe una asignaci贸n de valores a sus variables que hagan verdadera la f贸rmula (esto se denomina “satisfacer” la f贸rmula y de ah铆 viene el nombre SATisfiability asignado al problema). Por ejemplo,

$\mathrm{SAT}\left( (x_1 \vee x_2) \wedge (\neg x_1 \vee \neg x_2) \right) = \top$,

ya que la asignaci贸n $x_1 = \top, x_2 = \bot$ hace verdadera a la expresi贸n. Mientras que

$\mathrm{SAT}\left( (x_1 \vee x_2) \wedge (\neg x_1 \vee \neg x_2) \wedge x_1 \wedge x_2\right) = \bot$,

ya que no existe forma de asignarle valores de verdad a $x_1$ y $x_2$ sin hacer falsa a la f贸rmula.

Un algoritmo simple para resolverlo es probar distintas combinaciones de valores para las variables hasta encontrar alguna que haga verdadera a la f贸rmula o agotar las posibles combinaciones de valores. En el peor caso (cuando no se pueda encontrar una asignaci贸n), el algoritmo examinar谩 $2^n$ posibles asignaciones (siendo $n$ el n煤mero de variables). Por lo tanto, este algoritmo tiene una complejidad temporal $\Theta(2^n)$.

Pero, aunque no se cree posible encontrar un algoritmo de complejidad sub-exponencial para este problema, si se conocen algoritmos que mejoran significativamente al algoritmo anteriormente descrito. A continuaci贸n estudiaremos uno de los m谩s simples, conocido como DPLL.

#### DPLL

El algoritmo DPLL, llamado as铆 por las iniciales de sus creadores, mejora al algoritmo anteriormente descrito mediante el uso de dos t茅cnicas:

• Backtracking: la asignaci贸n de variables para satisfacer la f贸rmula se construye variable por variable, volvi茅ndose atr谩s apenas se detecta que la asignaci贸n de valor efectuada a una variable hace imposible satisfacer la f贸rmula.
• Unit propagation: cuando una de las cl谩usulas contiene una sola variable (negada o no), la conjunci贸n obliga a que esta variable tome un valor espec铆fico. El proceso de unit propagation consiste en propagar el efecto de esta asignaci贸n de valor sobre el resto de las cl谩usulas.

Si bien la referencia original utilizaba una t茅cnica adicional de eliminaci贸n de literales puros, esta no proporciona ventajas significativas en performance por lo que la omitiremos del an谩lisis.

#### Implementaci贸n de DPLL en Javascript

A continuaci贸n analizaremos una implementaci贸n de este algoritmo en Javascript. La implementaci贸n enlazada toma una cla煤sula por l铆nea, utilizando - como operador de negaci贸n y sin requerir operadores expl铆citos para la disyunci贸n y la conjunci贸n. Por ejemplo, la f贸rmula $(x_1 \vee x_2 \vee \neg x_3) \wedge (\neg x_2 \vee x_3) \wedge \neg x_1$ se representar铆a como

 x1 x2 -x3 -x2 x3 -x1 

Un script incluido en el archivo HTML se encarga de transformar el problema a un formato can贸nico, en el que las variables son reemplazadas por n煤meros enteros mayores que uno y la negaci贸n l贸gica se transforma en una negaci贸r aritm茅tica. Las cl谩usulas se convierten en arrays Javascript, con lo que la f贸rmula anterior quedar铆a transformada en [[1, 2, -3], [-2, 3], [-1]].

La primera funci贸n que analizaremos es simplifyClause():

function simplifyClause(c, a) {
var nc, i;
nc = [];
for (i = 0; i < c.length; i++) {
if (c[i] in a) {
return null;
} else if (!(-c[i] in a)) {
nc.push(c[i]);
}
}
return nc;
}


Toma dos par谩metros, c y a, que representan una cla煤sula dada como array de valores num茅ricos y una asignaci贸n de valores a las variables dada como un diccionario, respectivamente. El diccionario se utiliza como conjunto, empleando el n煤mero correspondiente a la variable como clave y true como valor; la clave num茅rica se niega para indicar la asignaci贸n de un valor falso a la variable. Por ejemplo, la asignaci贸n $x_1 = \top$ y $x_2 = \bot$ se representar铆a con {1: true, -2: true}.

La funci贸n simplemente copia la cla煤sula c a una nueva cl谩usula nc, examinando las apariciones de variables a la luz de la asignaci贸n de valores que se utiliza en ese momento. Si encuentra una aparici贸n de variable que se corresponda con una asignaci贸n incluida en a, significa que la cl谩usula es verdadera; por lo tanto devuelve null indicando la “desaparici贸n” de la misma. Por el contrario, si encuentra la aparici贸n negada en a, simplemente no la incluye en la nueva cla煤sula (ya que es falsa, el valor neutro de la disyunci贸n).

Continuando con las funciones, podemos examinar la funci贸n applyAssignment():

function applyAssignment(f, a) {
var nf, i, nc;
nf = [];
for (i = 0; i < f.length; i++) {
nc = simplifyClause(f[i], a);
if (nc !== null) {
nf.push(nc);
}
}
return nf;
}


Esta funci贸n hace lo mismo que la anteriormente descrita pero sobre toda una f贸rmula. Como la simplificaci贸n se efect煤a cl谩usula por cl谩usula, el 煤nico procesamiento no delegado a simplifyClause() es la supresi贸n de las cla煤sulas que son devueltas como null (que son las que simplifyClause() determin贸 como tautol贸gicas dada la asignaci贸n a).

Dejando a un lado funciones auxiliares, como cloneAssignment(), solo resta analizar la funci贸n principal: recDPLL().

function recDPLL(f, a) {
var i, na, v, ret;
f = applyAssignment(f, a);
if (f.length === 0) {
return [true, a];
}
for (i = 0; i < f.length; i++) {
if (f[i].length === 0) {
return [false, {}];
} else if (f[i].length === 1) {
na = cloneAssignment(a);
na[f[i][0]] = true;
return recDPLL(f, na);
}
}
na = cloneAssignment(a);
na[f[0][0]] = true;
ret = recDPLL(f, na);
if (ret[0]) {
return ret;
}
delete na[f[0][0]];
na[-f[0][0]] = true;
return recDPLL(f, na);
}


En primer lugar la funci贸n simplifica la f贸rmula f con la asignaci贸n de valores a, empleando las t茅cnicas anteriormente descritas. Si esta simplificaci贸n conduzca a una f贸rmula vac铆a, dado que una conjunci贸n vac铆a es verdadera, devuelve true y la asignaci贸n de valores que llev贸 a este valor de verdad.

Si no se lleg贸 a una f贸rmula vac铆a es obvio que restan cl谩usulas no tautol贸gicas dentro de la misma (todo esto tomando como “dada” la asignaci贸n a de variables). Estas pueden ser de dos tipos: cl谩usulas contradictorias o contingencias. Como las cl谩usulas restantes no pueden contener variables con valores asignados, dado el anterior paso de simplificaci贸n, una cl谩usula contradictoria solo puede ser vac铆a.

El procesamiento restante de la f贸rmula consiste en recorrer las cl谩usulas buscando cl谩usulas vac铆as o con una sola aparici贸n de variable. Si se encuentra un cl谩usula vac铆a, estaremos en presencia de una contradicci贸n que imposibilita aumentar a hasta convertirla en una asignaci贸n que satisfaga a la f贸rmula en cuesti贸n. Por lo tanto se realiza backtracking devolviendo [false, {}].

La aparici贸n de una cl谩usula con una sola variable es m谩s interesante, ya que puede verse como “forzando” una asignaci贸n. En este caso es cuando se realiza el proceso de unit propagation, que consiste simplemente en incorporar la 煤nica asignaci贸n compatible con la cl谩usula unitaria dentro de a y continuar la b煤squeda recursiva. No es necesario simplificar la f贸rmula, ya que esto se realizar谩 apenas comience la ejecuci贸n de la nueva instancia de recDPLL().

Por 煤ltimo, si la f贸rmula no cae en ninguno de estos dos casos se realiza una exploraci贸n recursiva normal.

#### Utilizando SAT para resolver otros problemas

Para mostrar como SAT puede utilizarse para resolver problemas l贸gicos en general, veremos como aplicarlo a la resoluci贸n de un problema de l贸gica conocido como “el problema de Einstein” (aunque es bastante poco probable que Einstein efectivamente lo haya creado). La formulaci贸n que utilizaremos es la siguiente:

1. There are 5 houses (along the street) in 5 different colors: blue, green, red, white and yellow.
2. In each house lives a person of a different nationality: Brit, Dane, German, Norwegian and Swede.
3. These 5 owners drink a certain beverage: beer, coffee, milk, tea and water,
smoke a certain brand of cigar: Blue Master, Dunhill, Pall Mall, Prince and blend,
and keep a certain pet: cat, bird, dog, fish and horse.
4. No owners have the same pet, smoke the same brand of cigar, or drink the same beverage.

Hints

1. The Brit lives in a red house.
2. The Swede keeps dogs as pets.
3. The Dane drinks tea.
4. The green house is on the left of the white house (next to it).
5. The green house owner drinks coffee.
6. The person who smokes Pall Mall rears birds.
7. The owner of the yellow house smokes Dunhill.
8. The man living in the house right in the center drinks milk.
9. The Norwegian lives in the first house.
10. The man who smokes blend lives next to the one who keeps cats.
11. The man who keeps horses lives next to the man who smokes Dunhill.
12. The owner who smokes Blue Master drinks beer.
13. The German smokes Prince.
14. The Norwegian lives next to the blue house.
15. The man who smokes blend has a neighbor who drinks water.

The question is: Who keeps fish?

Si bien probablemente la t茅cnica m谩s simple de aplicar en este caso sea constraint-based programming, podemos aplicar SAT sin problemas. Para ello debemos decidir primero que variables utilizaremos, siendo una elecci贸n razonable numerar las casas, colores, nacionalidades, bebidas, cigarros y mascotas e introducir variables booleanas asociando las casas con los distintos atributos de la persona que vive en ellas. Eso nos da un total de $5^3 = 125$ variables.

La traducci贸n de estas relaciones a cl谩usulas SAT puede ser algo tediosa, por lo que nos conviene emplear un script simple para este proceso. En caso de no disponer de Python, puede consultarse la salida del script en ep_sat.txt. La nomenclatura de las variables es simple: hijk, donde i indica el n煤emro de casa, j indica el atributo en cuesti贸n y k indica un posible valor de este atributo. En el c贸digo Python pueden verse los detalles de como se traducen las implicaciones y restricciones a disyunciones.

Introduciendo las cl谩usulas SAT en la p谩gina de prueba de DPLL, podemos encontrar el siguiente resultado:

 SATISFIABLE h0c4 h1c0 h2c2 h3c1 h4c3 h0n3 h0b4 h0s1 h0p0 h1n1 h1b3 h1s4 h1p4 h2n0 h2b2 h2s2 h2p1 h3n2 h3b1 h3s3 h3p3 h4n4 h4b0 h4s0 h4p2 -h0c3 -h2b0 -h2b1 -h2b3 -h2b4 -h0b2 -h1b2 -h3b2 -h4b2 -h2n1 -h2c1 -h3c3 -h0n0 -h0n1 -h0n2 -h0n4 -h1n3 -h2n3 -h3n3 -h4n3 -h2s0 -h1c1 -h1c2 -h1c3 -h1c4 -h0c0 -h2c0 -h3c0 -h4c0 -h1n0 -h2c3 -h0c1 -h0c2 -h2c4 -h3c2 -h4c2 -h3c4 -h4c1 -h4c4 -h3n0 -h4n0 -h5c3 -h3b0 -h3b3 -h3b4 -h0b1 -h1b1 -h4b1 -h3n1 -h0s0 -h0s2 -h0s3 -h0s4 -h1s1 -h2s1 -h3s1 -h4s1 -h1p0 -h1p1 -h1p2 -h1p3 -h0p4 -h2p4 -h3p4 -h4p4 -h1n4 -h1s2 -h3s0 -h4s4 -h0b0 -h0b3 -h1b4 -h4b4 -h2s4 -h3s4 -h0p1 -h0p2 -h0p3 -h2p0 -h3p0 -h4p0 -h1n2 -h4n1 -h1b0 -h4b3 -h1s0 -h1s3 -h2n2 -h2n4 -h2s3 -h3s2 -h4s2 -h4s3 -h2p2 -h2p3 -h3p1 -h4p1 -h4n2 -h3n4 -h4p3 -h3p2 

Como el pez ser铆a la mascota n煤mero 3, buscamos apariciones de p3 en una variable no negada y encontramos h3p3. Buscando la nacionalidad del habitante de la casa n煤mero 3 encontramos h3n2, por lo que podemos concluir que el alem谩n es quien tiene un pez como mascota (coincidiendo con la soluci贸n dada a este problema).