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

🙂

PyEEG & Análisis de Series Temporales

En este caso haremos una breve presentación de las funcionalidades que nos otorga PyEEG para el trabajo con electroencefalogramas, fusionando con matplotlib para mostrar los resultados.

PyEGG tiene muchas funcionalidades:

  • bin_power(), Power Spectral Density (PSD), spectrum power in a set of frequency bins, and, Relative Intensity Ratio (RIR), PSD normalized by total power in all frequency bins
  • pfd(), Petrosian Fractal Dimension (PFD)
  • hfd(), Higuchi Fractal Dimension (HFD)
  • hjorth(), Hjorth mobility and complexity
  • spectral_entropy(), spectral entropy, the entropy of RIRs
  • svd_entropy(), Singular Value Decomposition (SVD) entropy
  • fisher_info(), Fisher Information, the entropy of SVD values, i.e., SVD spectrum
  • ap_entropy(), Approximate Entropy (ApEn)
  • samp_entropy(), Sample Entropy (SampEn)
  • dfa(), Detrended Fluctuation Analysis (DFA)
  • hurst(), Hurst Exponent (Hurst)

En nuestro caso hicimos uso de DFA , en pocas palabras trata de cuantificar la autosemejanza de series temporales (En el link DFA esta genialmente explicado).

Los datos con los que trabajamos fueron tomados de  Universidad de Bonn.

Ahora vamos al código:

**Aclaración, solo por cuestiones de ploteo se le agrego a la función dfa que devuelva otros valores mas …

#!/usr/bin/env python2
# -*- coding: utf-8 -*-

#The data analyzed in our study is available on this page. 
# http://epileptologie-bonn.de/cms/front_content.php?idcat=193&lang=3
#The sampling rate of the data was 173.61 Hz. For a more 
#detailed description of the data please refer to the manuscript
#Please note, however, that the time series have the spectral 
#bandwith of the aquisition system, which is 0.5 Hz to 85 Hz.

from pyeeg import dfa
import pylab 
from scipy import log10
from numpy import arange, loadtxt, random
from random import randrange

def randomize(data):
	randomData = []
	for i in range(len(data)-1):
		randomData.append(data[randrange(0, len(data)-1)])
	
	return randomData

if __name__ == "__main__":
	
	originalData = loadtxt('F/F003.txt').T
	randomData = randomize(originalData)
	#get the dfa for both set of values
	alphaO, forplotO, forplotNO = dfa(originalData)
	alphaR, forplotR, forplotNR = dfa(randomData)

	#plot the results
	fig = pylab.figure()
	ax = fig.add_subplot(3, 1, 1)
	ori = ax.plot(log10(forplotNO), forplotO, 'b-', label = "Original Data")
	ran = ax.plot(log10(forplotNR), forplotR, 'g-', label = "Random Data")

	pylab.title("DFA")
	pylab.legend([ori[0], ran[0]], ['Original Data','Random Data'])
	pylab.text(2,3, r'$\alpha Original = %f $'%(alphaO), multialignment = 'center')
	pylab.text(2,2, r'$\alpha Random  = %f $'%(alphaR), multialignment = 'center')

	#plot the original and randomize data
	bx = fig.add_subplot(3, 1, 2)
	bx.plot(randomData, 'g-')
	pylab.title("Random Data")
	cx = fig.add_subplot(3, 1, 3)
	cx.plot(originalData, 'b-')
	pylab.title("Original Data")

	pylab.show()

Para observar la variación del \alpha lo que se hizo fue tomar la muestra original y transformarla en ruido blanco (el cual debe dar simil a 0.5)

Por si alguien le interesa ACA van a encontrar el pyeeg que retorna las funciones que se encuentran ploteadas, un setup.py para los herejes que lo necesiten corriendo en W$.

🙂

Curva de Koch con Numpy y Pylab

Esta vez venimos con código para generar la famosa Curva de Koch 
http://es.wikipedia.org/wiki/Copo_de_nieve_de_Koch , este es solo la curva..no el copo ...
pero el sistema es el mismo. 
Como verán el código es genérico porq tiene toda la intención de usarse 
con otras gramáticas de Lin. .. pero para eso falta un poco de pulido.
Los puntos clave pasan por la recursión, el cursor y los niveles :).
Espero que les sea útil.
#! /usr/bin/python2
# -*- coding: utf-8 -*-

import numpy as np
import pylab as pl

class Cursor:
    def __init__(self,x,y, angulo):
        self.x = x
        self.y = y
        self.ang = angulo

class SistemaL:
     def __init__(self,gramatica, inicial, distancia, nivel):
        self.gramatica = gramatica
        self.inicial = inicial
        self.distancia = distancia
        self.nivel = nivel

     def iniciar(self):
        self.cursor = Cursor(0,0,0)
        pl.plot([self.cursor.x, self.cursor.x+ self.distancia],
                [self.cursor.y, self.cursor.y])

     def generar(self, nivel, cursor, distancia):
        if nivel == 0:
           a = (cursor.x, cursor.y)
           b = (cursor.x + distancia*np.cos(cursor.ang),
                cursor.y + distancia*np.sin(cursor.ang))
           print a, "-->", b
           pl.plot([a[0],b[0]],[a[1],b[1]])
           cursor.x = b[0]
           cursor.y = b[1]
        else:
           print "nivel = ", nivel
           for i in self.gramatica[self.inicial]:
               if i == 'F':
                    self.generar(nivel-1, cursor, distancia/3)
               elif i == '+':
                    cursor.ang = cursor.ang - np.pi/3
               elif i == '-':
                    cursor.ang = cursor.ang + np.pi/3

if __name__ == "__main__":
     miSistema = SistemaL(gramatica = {'F':"F-F++F-F"},inicial = 'F',
                          distancia = 10000, nivel = 5)
     miSistema.iniciar()
     miSistema.generar(miSistema.nivel, miSistema.cursor,
                       miSistema.distancia)
     print "termine con todos los niveles"
     pl.show()

Yyyy unas imágenes de como sería nuestro resultado ...

El pase de diapositivas requiere JavaScript.