March 28, 2011 by mchouza

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! < 10⁷).

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

February 21, 2011 by mchouza

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:

Check if this fragment is a boundary point. If it’s one, return the previous value of this point.
Get the four nearest neighbors’ values.
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.

November 8, 2010 by mchouza

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.

August 11, 2010 by mchouza

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 🙂

May 30, 2010 by mchouza

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 😀 ).

May 2, 2010 by mchouza

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:

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.
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$ ).
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.
$k(e) P(e) \le 1$ .
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:

m_i < 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.

April 29, 2010 by mchouza

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”):

Empezar apostando una unidad.
Apostar a “cara”.
En caso de ganar, retirarse.
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… 😀

January 3, 2010 by mchouza

El poder del SAT

Generalidades

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?

(Puede verse también una formulación equivalente en español.)

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).

Spin Foam

Just another WordPress.com weblog

Category Archives: javascript