Mapa Logístico con PyLab – NumPy – PyOpenCL

En este caso vamos a meternos en PyOpenCL con la excusa de los mapas logísticos .

ACA tienen el código fuente de un excelente tutorial, en el cual me base.

El mapa logístico puede expresarse matematicamente como:

X_{n+1} = \alpha X_n (1-X_n)

Habiendo dicho esto veremos el código para nuestro kernel (.cl):

/*kernel for implementing the Verhulst model.
Fractal and Chaos - Chapter 5 - Iterative Feedback Processes and Chaos. */

__kernel void verhulst(global float* alpha, global float* xd, global float* output)
{
    float temp;
    int idx = get_local_id(0)+get_local_size(0)*get_group_id(0);
    temp = xd[0];
    int i;
    for(i=0;i<=idx;i++)
    {
        temp = alpha[0]*temp*(1-temp);
    }
    output[idx] = temp;
}

Y nuestro programa encargado de crear el contexto, crear buffers, etc.

# http://documen.tician.de/pyopencl/
# Ian Johnson.
#
#
# Celita changes buffers, data and parameters
# for implementation of the logistic map
#

import pyopencl as cl
import numpy
import pylab
import argparse

class CL:
    def __init__(self):
        self.ctx = cl.create_some_context()
        self.queue = cl.CommandQueue(self.ctx)

    def loadProgram(self, filename):
        #read in the OpenCL source file as a string
        f = open(filename, 'r')
        fstr = "".join(f.readlines())
        print fstr
        #create the program
        self.program = cl.Program(self.ctx, fstr).build()

    def logisticMap(self, alpha=0.9, xd=0.2):
        mf = cl.mem_flags

        #initialize client side (CPU) arrays
        self.alpha = numpy.array([alpha], dtype=numpy.float32)
        self.xd = numpy.array([xd], dtype=numpy.float32)
        self.output = numpy.zeros(50, numpy.float32)

        #create OpenCL buffers
        self.alpha_buf = cl.Buffer(self.ctx, mf.READ_ONLY | mf.COPY_HOST_PTR,
                                   hostbuf=self.alpha)
        self.xd_buf = cl.Buffer(self.ctx, mf.READ_ONLY | mf.COPY_HOST_PTR,
                                hostbuf=self.xd)
        self.dest_buf = cl.Buffer(self.ctx, mf.WRITE_ONLY, self.output.nbytes)

    def execute(self):
        self.program.verhulst(self.queue, self.alpha.shape, None,
                              self.alpha_buf, self.xd_buf, self.dest_buf)
        cl.enqueue_copy(self.queue, self.output, self.dest_buf)
        print "alpha: ", self.alpha
        print "xn: ", self.xd
        print "results: ", self.output

if __name__ == "__main__":
    parser = argparse.ArgumentParser(description=
                                     'Logistic Map with control parameter')
    parser.add_argument("--alpha", dest="alpha", default=0.9,
                        help='This option is used to set the alpha')
    args = parser.parse_args()
    alphaIn = args.alpha

    myLogisticMap = CL()
    myLogisticMap.loadProgram("verhulst.cl")
    myLogisticMap.logisticMap(alpha=alphaIn)
    myLogisticMap.execute()

    pylab.plot(myLogisticMap.output, 'ro-')
    pylab.ylim(0, 1)
    pylab.text(40, 0.2, r'$\alpha = %f$' % (myLogisticMap.alpha))
    pylab.title("Successive Iterations of the logistic Map")
    pylab.show()

Como se ve en el codigo se puede pasar el \alpha por parametros o dejar que se tome el por defecto …
Y como resultado veríamos algo asi …

Para mayor seguridad de código actualizado verlo ACA

🙂

Conjunto de Cantor con Python/Mayavi

Reutilizando el código anterior con mínimas modificaciones tenemos como resultado el Conjunto de Cantor (http://es.wikipedia.org/wiki/Conjunto_de_Cantor).

Solo colocaré la parte modifcada del código anterior, no se olviden de importar mayavi!


def generar(self, nivel, cursor, distancia):
    if nivel == 0:
       a = (cursor.x, cursor.x, cursor.x)
       b = (cursor.x + distancia*np.cos(cursor.ang),
            cursor.y + distancia*np.sin(cursor.ang))
       mlab.points3d(a, a, a, colormap="copper", scale_factor=.25)
       
       cursor.x = b[0]
       cursor.y = b[1]
    else:
       print "nivel = ", nivel
       for i in self.gramatica[self.inicial]:
           if i == 'A':
              self.generar(nivel-1, cursor, distancia/3.0)
           elif i == 'B':
              cursor.x += 3*distancia

if __name__ == "__main__":
     miSistema = SistemaL(gramatica = {'A':"ABA", 'B':"BBB"},
                          inicial = 'A',distancia = 1, nivel = 5)
     miSistema.iniciar()
     miSistema.generar(miSistema.nivel, miSistema.cursor, miSistema.distancia)
     print "termine con todos los niveles"
     mlab.show()

Y así quedaría …

El pase de diapositivas requiere JavaScript.

btw: puede haber quedado algo residual del anterior problema 🙂