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

🙂

Anuncios

Responder

Introduce tus datos o haz clic en un icono para iniciar sesión:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s