Neural Network Simulator for Blender2.5

#!/usr/bin/python
## RPython Neural Network - 0.2
## by Brett Hartshorn goatman.py@gmail.com
## License GNU GPL
## demo1 Python vs PyPy: http://www.youtube.com/watch?v=HwbDG3Rdi_Q
## demo2 Blender Rendering: http://www.youtube.com/watch?v=4sMrG2g7LmA

## user defined options ##
BLENDER_DIR = '../Desktop/blender25/blender'
## number of neurons, be careful with large numbers!
COLUMNS = 24
LAYERS = 3
STEM = 4

'''
Getting Started:

    Tested with Ubuntu Lucid and PyPy 1.3
    (in theory works on Windows, but installing the PyPy toolchain and SDL headers would be a serious pain)
    Why PyPy?  Because Python is too slow to run a simulation with any modest number of neurons.

    1. Download PyPy source code - see http://pypy.org

    2. Install GCC and SDL headers if you haven't already
        apt-get install build-essential libsdl-dev

    3. Make this file executable: "chmod +x rAI.py"

    4. Modify the line above BLENDER_DIR to point to your blender25 installation

    5. Copy this file to your pypy root directory and run "./rAI.py"
        (you can test just pypy compilation with "./rAI.py --pypy --subprocess")

How It Works:
    This program is actually three programs in a single Python script,
    some trickery with command line arguments is used to divide the programs.
    When you run "./rAI.py" the top level begins, opening a GTK window,
    When the drawing-area in the GTK window is realized it calls itself using the
    subprocess module and passing "--pypy --subprocess" on the command line,
    when the script is run this way it imports pypy and translates the simulator
    into C, the SDL_WINDOWID environment variable is used to tell the PyPy
    subprocess to draw into that xwindow ID.
    After the simulation has run, you can dump it to Blender25 by clicking "view model"
    When you do this blender is called with the -P argument pointing to itself (this script)
    The script knows its being run from blender and creates the BlenderExport class
    that handles converting the data.

Unfinished Work:

    For Blender visualization it would be nice to see which dendrites are polarized.

    Realtime Brain Input/Output: a brain sim without any input or output is pretty boring,
        to input spikes into the network provide a function that writes a command to the pipe,
        ( see "spike-all" and "spike-one" )
        for example you could modify this program to open a socket, and then from the
        Blender Game Engine your actuator sends a signal to the socket when something it touched,
        depending on whats touched, the toplevel app will write a command to the pipe that
        the PyPy subprocess is reading.
        Getting output from the network is done the same way, pipe/socket -> client - do something.
        ( note that single neuron spiking gives little information - its better you FFT the spikes into a
        useable pattern that can drive some action )

    Skeleton code for learning/training the network, needs alot of work.



-----------------------------------------------------------
PyPy Tips for Coders:

math.radians    not available
random.uniform  not available

incorrect:
    string.strip()  will not work without an arg
    list.sort()     not available

correct:
    string.strip(' ').strip('\n').strip('\t')       #strips all white space
    list = sortd( list )                        # but the JIT prefers TimSort
    list.pop( index )                       # no list.pop() with end as default

'''


import os, sys, time, pickle
import math             # math.radians is missing in pypy?

degToRad = math.pi / 180.0
def radians(x):
    """radians(x) -> converts angle x from degrees to radians
    """
    return x * degToRad

if '--blender' not in sys.argv:
    from pypy.rlib import streamio
    from pypy.rlib import rpoll
    # apt-get install libsdl-dev
    from pypy.rlib.rsdl import RSDL, RSDL_helper
    from pypy.rlib.rarithmetic import r_uint
    from pypy.rpython.lltypesystem import lltype, rffi
    from pypy.rlib.listsort import TimSort
    from pypy.rlib.jit import hint, we_are_jitted, JitDriver, purefunction_promote

#random.randrange(0, up)        # not a replacement for random.uniform
#Choose a random item from range(start, stop[, step]).
#This fixes the problem with randint() which includes the
#endpoint; in Python this is usually not what you want.


if '--pypy' in sys.argv:        # random.random works in pypy, but random.uniform is missing?
    from pypy.rlib import rrandom
    RAND = rrandom.Random()
    RAND.init_by_array([1, 2, 3, 4])
    def random(): return RAND.random()
    def uniform(start,end): return ( RAND.random() * (end-start) ) - start
else:
    from random import *

def distance( v1, v2 ):
    dx = v1[0] - v2[0]
    dy = v1[1] - v2[1]
    dz = v1[2] - v2[2]
    t = dx*dx + dy*dy + dz*dz
    return math.sqrt( float(t) )

if '--blender' not in sys.argv:
    Njitdriver = JitDriver(
        greens = 'cur train times branes last_spike spikers temporal thresh triggers self abs_refactory abs_refactory_value'.split(),
        reds = 'i t a b c bias neuron'.split()
    )
class RecurrentSpikingModel(object):
    '''
    Model runs in realtime and lossy - stores recent
    spikes in a list for realtime learning.
    Uses spike train, with time delay based on distance,        Spikes will trigger absolute refactory period.
    Membrane has simple linear falloff.

        notes:
            log(1.0) = 0.0
            log(1.5) = 0.40546510810816438
            log(2.0) = 0.69314718055994529
            log(2.7) = 0.99325177301028345
            log(3.0) = 1.0986122886681098

    '''
    def iterate( self ):
        now = float( time.time() )
        self._state_dirty = True; train = self._train
        brane = self._brane; rest = self._rest
        branes = self._branes; temporal = self._temporal
        abs_refactory = self._abs_refactory
        abs_refactory_value = self._abs_refactory_value
        fps = self._fps
        elapsed = now - self._lasttime
        clip = self._clip
        cur = self._lasttime
        last_spike = self._last_spike
        spikers = self._spikers
        thresh = self._thresh
        triggers = self._triggers

        ## times of incomming spikes ##
        times = train.keys()
        ## To do, if seziure, lower all connection strengths
        TimSort(times).sort()       # pypy note - list have no .sort(), and JIT dislikes builtin sorted(list)

        if not times:
            a = now - self._lasttime
            if a > 0.1: brane = rest
            else:
                b = a * 10  # 0.0-1.0
                brane += (rest - brane) * b
            branes.append( (brane,now) )

        i = 0
        t = a = b = c = bias = .0
        neuron = self
        ntrain = len(train)
        while train:        #the JIT does only properly support a while loop as the main dispatch loop
            Njitdriver.can_enter_jit(
                cur=cur, train=train, times=times, branes=branes, last_spike=last_spike, spikers=spikers,
                temporal=temporal, thresh=thresh, triggers=triggers, self=self,
                abs_refactory=abs_refactory, abs_refactory_value=abs_refactory_value,
                i=i, t=t, a=a, b=b, c=c,
                neuron=neuron, bias=bias
            )
            Njitdriver.jit_merge_point(
                cur=cur, train=train, times=times, branes=branes, last_spike=last_spike, spikers=spikers,
                temporal=temporal, thresh=thresh, triggers=triggers, self=self,
                abs_refactory=abs_refactory, abs_refactory_value=abs_refactory_value,
                i=i, t=t, a=a, b=b, c=c,
                neuron=neuron, bias=bias
            )
            ##bias, neuron = train.pop( t ) ## not pypy
            t = times[i]; i += 1
            bias,neuron = train[t]
            del train[t]


            if t <= cur:
                #print 'time error'
                pass
            else:
                a = t - cur
                cur += a
                #print 'delay', a
                if cur - last_spike < abs_refactory:
                    brane = abs_refactory_value
                    branes.append( (brane,cur) )
                else:
                    spikers.append( neuron )
                    if a > 0.1:
                        brane = rest + bias
                        branes.append( (brane,cur) )
                    else:
                        b = a * 10  # 0.0-1.0
                        brane += (rest - brane) * b
                        c = b * temporal
                        bias *= math.log( (temporal-c)+2.8 )
                        brane += bias
                        branes.append( (brane,cur) )

                        if brane > thresh:
                            triggers.append( neuron )
                            self.spike( cur )
                            brane = abs_refactory_value
                            #fps *= 0.25        # was this to play catch up?
                            break # this is efficient!

        #self._train = {}
        self._brane = brane
        #for lst in [branes, spikers, triggers]:    # not allowed in pypy
        while len(branes) > 512: branes.pop(0)      # Adjust this if you run out of memory
        while len(spikers) > clip: spikers.pop(0)
        while len(triggers) > clip: triggers.pop(0)

        self._lasttime = end = float(time.time())
        #print ntrain, end-now


    def detach( self ):
        self._inputs = {}       # time:neuron
        self._outputs = []  # children (neurons)
        self._train = {}
        self._branes = []
        self._spikers = []  # last spikers
        self._triggers = [] # last spikers to cause a spike

    def __init__( self, name='neuron', x=.0,y=.0,z=.0, column=0, layer=0, fps=12, thresh=100.0, dendrite_bias=35.0, dendrite_noise=1.0, temporal=1.0, distance_factor=0.1, red=.0, green=.0, blue=.0 ):
        self._name = name
        self._fps = fps     # active fps 10?
        self._thresh = thresh
        self._column = column
        self._layer = layer
        self._brane = self._rest = -65.0
        self._abs_refactory = 0.05
        self._abs_refactory_value = -200.0
        self._spike_value = 200.0
        self._temporal = temporal   # ranges 1.0 - inf
        self._distance_factor = distance_factor
        self._dendrite_bias = dendrite_bias
        self._dendrite_noise = dendrite_noise
        self._clip = 128
        self._learning = False
        self._learning_rate = 0.1

        ## the spike train of incomming spikes
        self._train = {}

        ## list of tuples, (brane value, abs time)
        ## use this list to render wave timeline ##
        self._branes = []

        self._spikers = []  # last spikers
        self._triggers = [] # last spikers to cause a spike
        self._lasttime = time.time()
        self._last_spike = 0
        self._inputs = {}       # time:neuron
        self._outputs = []  # children (neurons)

        self._color = [ red, green, blue ]
        #x = uniform( 0.2, 0.8 )
        #y = random()   #uniform( 0.2, 0.8 )
        #z = random()   #uniform( 0.2, 0.8 )
        self._pos = [ x,y,z ]
        self._spike_callback = None
        self._state_dirty = False
        self._active = False
        self._cached_distances = {}
        self._clipped_spikes = 0    # for debugging

        self._draw_spike = False
        #self._braneRect = self._spikeRect = None

    def randomize( self ):
        for n in self._inputs:
            bias = self._dendrite_bias + uniform( -self._dendrite_noise, self._dendrite_noise )
            self._inputs[ n ] = bias


    def mutate( self, v ):
        for n in self._spikers:
            self._inputs[ n ] += uniform(-v*2,v*2)
        for n in self._triggers:
            self._inputs[ n ] += uniform(-v,v)

        if not self._spikes or not self._triggers:
            for n in self._inputs:
                self._inputs[ n ] += uniform(-v,v)


    def reward( self, v ):
        if not self._spikers: print( 'no spikers to reward' )
        for n in self._spikers:
            bias = self._inputs[ n ]
            if abs(bias) < 100:
                if bias > 15:
                    self._inputs[ n ] += v
                else:
                    self._inputs[ n ] -= v


    def punish( self, v ):
        for n in self._inputs:
            self._inputs[n] *= 0.9

        for n in self._spikers:
            self._inputs[ n ] += uniform( -v, v )



    def attach_dendrite( self, neuron ):
        bias = self._dendrite_bias + uniform( -self._dendrite_noise, self._dendrite_noise )
        if random()*random() > 0.5: bias = -bias
        # add neuron to input list
        if neuron not in self._inputs:
            self._inputs[ neuron ] = bias
            #print 'attached neuron', neuron, bias
        if self not in neuron._outputs:
            neuron._outputs.append( self )
            neuron.update_distances()
        return bias

    def update_distances( self ):
        for child in self._outputs:
            dist = distance( self._pos, child._pos )
            self._cached_distances[ child ] = dist

    def spike( self, t ):
        #print 'spike', t
        self._draw_spike = True
        self._last_spike = t
        self._branes.append( (self._spike_value,t) )
        self._brane = self._abs_refactory_value
        self._branes.append( (self._brane,t) )
        if self._learning and self._triggers:
            #print 'learning'
            n = self._triggers[-1]
            bias = self._inputs[ n ]
            if bias > 0: self._inputs[n] += self._learning_rate
            else: self._inputs[n] -= self._learning_rate

        if self._spike_callback:
            self._spike_callback( t )

        for child in self._outputs:
            #print 'spike to', child
            bias = child._inputs[ self ]
            #dist = distance( self._pos, child._pos )
            dist = self._cached_distances[ child ]
            dist *= self._distance_factor
            child._train[ t+dist ] = (bias,self)

    def stop( self ):
        self._active = False
        print( 'clipped spikes', self._clipped_spikes )
    def start( self ):
        self._active = True
        self.iterate()

    def setup_draw( self, format, braneRect, groupRect, spikeRect, colors ):
        self._braneRect = braneRect
        self._groupRect = groupRect
        self._spikeRect = spikeRect
        self._sdl_colors = colors
        r,g,b = self._color
        self._group_color = RSDL.MapRGB(format, int(r*255), int(g*255), int(b*255))

    def draw( self, surf ):
        #if self._braneRect:
        fmt = surf.c_format
        b = int(self._brane * 6)
        if b > 255: b = 255
        elif b < 0: b = 0
        color = RSDL.MapRGB(fmt, 0, 0, b)
        RSDL.FillRect(surf, self._braneRect, color)
        RSDL.FillRect(surf, self._groupRect, self._group_color)
        if self._draw_spike:
            RSDL.FillRect(surf, self._spikeRect, self._sdl_colors['white'])
        else:
            RSDL.FillRect(surf, self._spikeRect, color )#self._sdl_colors['black'])
        self._draw_spike = False

    def dump_model( self, brain ):
        #a = ''
        x,y,z = self._pos
        a = '{ "pos" : (%s, %s, %s,), ' %(x,y,z)        # self._pos  not allowed in pypy, unpack tuple first
        b = ''
        for brane,seconds in self._branes: b += '(%s,%s),' %(brane,seconds)
        a += ' "anim" : ( %s ), ' %b
        b = ''
        for n in self._outputs:
            idx = brain.get_neuron_index( n )
            b += '%s,' %idx
        a += ' "outputs" : ( %s ) }' %b
        return a



class Brain( object ):
    def loop( self ):
        self._active = True
        fmt = self.screen.c_format
        RSDL.FillRect(self.screen, lltype.nullptr(RSDL.Rect), self.ColorGrey)
        layers = self._layers
        neurons = self._neurons
        pulse_layers = self._pulse_layers
        screen = self.screen
        loops = 0
        now = start = float( time.time() )
        while self._active:
            #Bjitdriver.can_enter_jit( layers=layers, neurons=neurons, pulse_layers=pulse_layers, loops=loops )
            #Bjitdriver.jit_merge_point( layers=layers, neurons=neurons, pulse_layers=pulse_layers, loops=loops)
            now = float( time.time() )
            self._fps = loops / float(now-start)
            for i,lay in enumerate(self._layers):
                if self._pulse_layers[i] and False:
                    #print 'pulse layer: %s neurons: %s ' %(i, len(lay))
                    for n in lay:
                        if random()*random() > 0.8:
                            n.spike( now )
            for i,col in enumerate(self._columns):
                if self._pulse_columns[i]:
                    for n in col: n.spike(now)
            for n in self._neurons:
                n.iterate()
                n.draw(self.screen)
            #r,w,x = rpoll.select( [self._stdin], [], [], 1 )   # wait
            rl,wl,xl = rpoll.select( [0], [], [], 0.000001 )    # wait
            if rl:
                cmd = self._stdin.readline().strip('\n').strip(' ')
                self.do_command( cmd )
            loops += 1
            self._iterations = loops
            #print loops        # can not always print in mainloop, then select can never read from stdin
            RSDL.Flip(self.screen)
            #self._fps = float(time.time()) - now
        #return loops
        return 0

    def __init__(self):
        start = float(time.time())
        self._neurons = []
        self._columns = []
        self._layers = [ [] ] * LAYERS
        self._pulse_layers = [0] * LAYERS
        self._pulse_layers[ 0 ] = 1

        self._pulse_columns = [0] * COLUMNS
        self._pulse_columns[ 0 ] = 1
        self._pulse_columns[ 1 ] = 1
        self._pulse_columns[ 2 ] = 1
        self._pulse_columns[ 3 ] = 1


        inc = 360.0 / COLUMNS
        scale = float( LAYERS )
        expansion = 1.333
        linc = scale / LAYERS
        for column in range(COLUMNS):
            colNeurons = []
            self._columns.append( colNeurons )
            X = math.sin( radians(column*inc) )
            Y = math.cos( radians(column*inc) )
            expanding = STEM
            width = 1.0 / scale
            for layer in range(LAYERS):
                Z = layer * linc
                r = random() * random()
                g = 0.2
                b = 0.2
                for i in range(int(expanding)):
                    x = uniform( -width, width )
                    rr = random()*random()      # DJ's trick
                    y = uniform( -width*rr, width*rr ) + X
                    z = Z + Y
                    # create 50/50 exitatory/inhibitory
                    n = RecurrentSpikingModel(x=x, y=y, z=z, column=column, layer=layer, red=r, green=g, blue=b )
                    self._neurons.append( n )
                    colNeurons.append( n )
                    self._layers[ layer ].append( n )

                expanding *= expansion
                width *= expansion

        dendrites = 0
        interlayer = 0
        for lay in self._layers:
            for a in lay:
                for b in lay:
                    if a is not b and a._column == b._column:
                        a.attach_dendrite( b )
                        dendrites += 1
                        interlayer += 1

        intercol = 0
        for col in self._columns:
            for a in col:
                for b in col:
                    if a is not b and random()*random() > 0.75:
                        a.attach_dendrite( b )
                        intercol += 1
                        dendrites += 1

        intercore = 0
        core = self._layers[-1]
        for a in core:
            for b in core:
                if a is not b and random()*random() > 0.85:
                    a.attach_dendrite( b )
                    intercore += 1
                    dendrites += 1

        print( 'brain creation time (seconds)', float(time.time())-start )
        print( 'neurons per column', len(self._columns[0]) )
        print( 'inter-layer dendrites', interlayer )
        print( 'inter-column dendrites', intercol )
        print( 'inter-neocoretex dendrites', intercore )
        print( 'total dendrites', dendrites )
        print( 'total neurons', len(self._neurons) )
        for i,lay in enumerate(self._layers):
            print( 'layer: %s   neurons: %s' %(i,len(lay)) )
        for i,col in enumerate(self._columns):
            print( 'column: %s  neurons: %s' %(i,len(col)) )



        self._stdin = streamio.fdopen_as_stream(0, 'r', 1)
        #self._stdout = streamio.fdopen_as_stream(1, 'w', 1)
        #self._stderr = streamio.fdopen_as_stream(2, 'w', 1)

        self._width = 640; self._height = 480
        assert RSDL.Init(RSDL.INIT_VIDEO) >= 0
        self.screen = RSDL.SetVideoMode(self._width, self._height, 32, 0)
        assert self.screen
        fmt = self.screen.c_format
        self.ColorWhite = white = RSDL.MapRGB(fmt, 255, 255, 255)
        self.ColorGrey = grey = RSDL.MapRGB(fmt, 128, 128, 128)
        self.ColorBlack = black = RSDL.MapRGB(fmt, 0, 0, 0)
        self.ColorBlue = blue = RSDL.MapRGB(fmt, 0, 0, 200)

        colors = {'white':white, 'grey':grey, 'black':black, 'blue':blue}

        x = 1; y = 1
        for i,n in enumerate(self._neurons):
            braneRect = RSDL_helper.mallocrect(x, y, 12, 12)
            groupRect = RSDL_helper.mallocrect(x, y, 12, 2)
            spikeRect = RSDL_helper.mallocrect(x+4, y+4, 4, 4)
            n.setup_draw( self.screen.c_format, braneRect, groupRect, spikeRect, colors )
            x += 13
            if x >= self._width-14:
                x = 1
                y += 13

    def do_command( self, cmd ):
        if cmd == 'spike-all':
            t = float(time.time())
            for n in self._neurons: n.spike(t)
        elif cmd == 'spike-one':
            t = float(time.time())
            self._neurons[0].spike(t)
        elif cmd == 'spike-column':
            t = float(time.time())
            for n in self._columns[0]:
                n.spike(t)
        elif cmd == 'info':
            info = self.info()
            #sys.stderr.write( info )
            #sys.stderr.flush()
            print( info )
        elif cmd == 'view-model':
            a = ''
            for n in self._neurons: a += n.dump_model( self ) + ','
            print( '<view_model> ( %s )' %a )

    def info(self):
        r = ' "num-layers": %s,' %len(self._layers)
        r += ' "num-neurons": %s,' %len(self._neurons)
        r += ' "fps" : %s, ' %self._fps
        r += ' "iterations" : %s, ' %self._iterations
        return '<load_info> { %s }' %r

    def get_neuron_index( self, n ): return self._neurons.index( n )




import subprocess, select, time
if '--blender' not in sys.argv: import gtk, glib

class App:
    def load_info( self, arg ): print( arg )
    def view_model( self, arg ):
        pickle.dump( arg, open('/tmp/brain.data','wb'), -1 )
        os.system( '%s -P rAI.py --blender' %BLENDER_DIR )

    def __init__(self):
        self._commands = cmds = []
        self.win = gtk.Window()
        self.win.connect('destroy', lambda w: gtk.main_quit())
        self.root = gtk.VBox(False,10); self.win.add( self.root )
        self.root.set_border_width(20)
        self.header = header = gtk.HBox()
        self.root.pack_start( header, expand=False )
        b = gtk.Button('spike all neurons')
        b.connect('clicked', lambda b,s: s._commands.append('spike-all'), self )
        self.header.pack_start( b, expand=False )

        b = gtk.Button('spike one neuron')
        b.connect('clicked', lambda b,s: s._commands.append('spike-one'), self )
        self.header.pack_start( b, expand=False )

        b = gtk.Button('spike column 1')
        b.connect('clicked', lambda b,s: s._commands.append('spike-column'), self )
        self.header.pack_start( b, expand=False )

        b = gtk.Button('view model')
        b.connect('clicked', lambda b,s: s._commands.append('view-model'), self )
        self.header.pack_start( b, expand=False )

        self.header.pack_start( gtk.SeparatorMenuItem() )

        b = gtk.Button('debug')
        b.connect('clicked', lambda b,s: s._commands.append('info'), self )
        self.header.pack_start( b, expand=False )

        da = gtk.DrawingArea()
        da.set_size_request( 640,480 )
        da.connect('realize', self.realize)
        self.root.pack_start( da )

        self._read = None
        glib.timeout_add( 33, self.loop )
        self.win.show_all()


    def realize(self, da ):
        print( 'gtk drawing area realize' )
        xid = da.window.xid
        self._process = process = subprocess.Popen( 'python rAI.py --pypy --subprocess %s' %xid, stdin=subprocess.PIPE, stdout=subprocess.PIPE, bufsize=32, shell=True )
        self._write = write = process.stdin
        self._read = read = process.stdout
        print( 'read pipe', read )
        print( 'read write', write )

    def loop( self ):
        if self._read:
            rlist,wlist,xlist = select.select( [self._read], [], [], 0.001 )
            while self._commands:
                cmd = self._commands.pop()
                print( 'sending cmd ->', cmd )
                self._write.write( '%s\n'%cmd )
                self._write.flush()
            if rlist:
                a = self._read.readline().strip()
                if a:
                    print( a )
                    if a.startswith('<'):
                        func = a[ 1 : a.index('>') ]
                        arg = a[ a.index('>')+1 : ].strip()
                        func = getattr(self, func)
                        func( eval(arg) )
        return True

class BlenderExport(object):
    #bpy.ops.anim.change_frame( i )     # can not be used from the command line, invalid context
    #mat.diffuse_color.r = 0.1
    #mat.keyframe_insert('diffuse_color')

    def animate_material( self, mat, anim ):
        mat.keyframe_insert('emit')
        start = anim[0][1]
        ecurve = mat.animation_data.action.fcurves[0]
        i = 0
        for brane,seconds in anim:
            t = seconds-start
            t *= 100
            if brane < 99: v = (brane+80) * .01
            else: v = 50
            bez = ecurve.keyframe_points.add( frame=t, value=v )
            i += 1

    def create_taper_curve( self, thickness=20.0, steps=3 ):
        # the Y+ defines the tappering
        data = bpy.data.curves.new(name='dendriteTaper', type='CURVE')
        spline = data.splines.new('BEZIER')
        #spline.bezier_points.add()     # new spline always has one bez
        x = .0; y=0.5 * thickness; z = .0
        bez = spline.bezier_points[-1]
        bez.co.x = x; bez.co.y = y; bez.co.z = z
        bez.handle1.x = x + 0.1
        bez.handle1.y = y / 2
        bez.handle1.z = z
        bez.handle2.x = x + 0.1
        bez.handle2.y = y / 2
        bez.handle2.z = z

        if 0:
            inc = 1.0 / steps
            x = inc
            for i in range(steps):
                spline.bezier_points.add()
                y= (0.1 * thickness) + uniform(.0,0.01)
                bez = spline.bezier_points[-1]
                bez.co.x = x; bez.co.y = y; bez.co.z = z
                bez.handle1.x = x
                bez.handle1.y = y
                bez.handle1.z = z
                bez.handle2.x = x + inc
                bez.handle2.y = y
                bez.handle2.z = z
                x += inc

        spline.bezier_points.add()
        x = 1.0; y=0.05 * thickness
        bez = spline.bezier_points[-1]
        bez.co.x = x; bez.co.y = y; bez.co.z = z
        bez.handle1.x = x - 0.5
        bez.handle1.y = y
        bez.handle1.z = z
        bez.handle2.x = x - 0.5
        bez.handle2.y = y
        bez.handle2.z = z

        obj = bpy.data.objects.new(name='dendriteTraper', object_data=data)
        bpy.context.scene.objects.link(obj)
        return obj

    def __init__(self):
        data = bpy.data.curves.new(name='brain', type='CURVE')
        data.taper_object = self.create_taper_curve()
        self._dendrite_curves = data
        data.dimensions = '3D'
        data.resolution_u = 3
        #print( dir(data) )
        brain = pickle.load( open('/tmp/brain.data','rb') )
        steps = 5
        self._neuron_cell_bodies = ncb = []
        self._neuron_materials = mats = []
        for i,n in enumerate(brain):
            x,y,z = n['pos']
            bpy.ops.mesh.primitive_ico_sphere_add() # does not return mesh
            for ob in bpy.data.objects:
                if ob.name.startswith('Icosphere') and ob not in ncb:
                    ncb.append( ob )
                    ob.location.x = x
                    ob.location.y = y
                    ob.location.z = z
                    ob.scale.x = ob.scale.y = ob.scale.z = .03
                    ob.selected = True
                    bpy.ops.object.shade_smooth()
                    mat = bpy.data.materials.new('cellbody')
                    mats.append( mat )
                    ob.data.add_material( mat )
                    self.animate_material( mat, n['anim'] )
                    for vert in ob.data.verts:
                        if random() > 0.5:
                            p =  uniform(1.0,1.2)
                            vert.co.x *= p
                            vert.co.y *= p
                            vert.co.z *= p
                    break
            for dendrite in n['outputs']:

                spline = data.splines.new('BEZIER')
                #spline.bezier_points.add()     # new spline always has one bez
                bez = spline.bezier_points[-1]
                bez.co.x = x; bez.co.y = y; bez.co.z = z
                bez.handle1.x = x
                bez.handle1.y = y
                bez.handle1.z = z
                bez.handle2.x = x + uniform(-.1,.1)
                bez.handle2.y = y + uniform(-.1,.1)
                bez.handle2.z = z + uniform(-.1,.1)

                dx,dy,dz = brain[ dendrite ]['pos']
                inc = 1.0 / steps
                m = inc
                a = 0.05
                for step in range(steps):
                    mx = x + ((dx-x)*uniform(m-a,m+a))
                    my = y + ((dy-y)*uniform(m-a,m+a))
                    mz = z + ((dz-z)*uniform(m-a,m+a))
                    spline.bezier_points.add()
                    bez = spline.bezier_points[-1]
                    bez.co.x = mx; bez.co.y = my; bez.co.z = mz
                    bez.handle1.x = mx + uniform(-.1,.1)
                    bez.handle1.y = my + uniform(-.1,.1)
                    bez.handle1.z = mz + uniform(-.1,.1)
                    bez.handle2.x = mx + uniform(-.1,.1)
                    bez.handle2.y = my + uniform(-.1,.1)
                    bez.handle2.z = mz + uniform(-.1,.1)
                    m += inc

                spline.bezier_points.add()
                bez = spline.bezier_points[-1]
                bez.co.x = dx; bez.co.y = dy; bez.co.z = dz
                bez.handle1.x = dx
                bez.handle1.y = dy
                bez.handle1.z = dz
                bez.handle2.x = dx + uniform(-.1,.1)
                bez.handle2.y = dy + uniform(-.1,.1)
                bez.handle2.z = dz + uniform(-.1,.1)


        obj = bpy.data.objects.new(name='brain', object_data=data)
        bpy.context.scene.objects.link(obj)

if '--subprocess' in sys.argv:
    os.putenv('SDL_WINDOWID', sys.argv[-1])
    def pypy_entry_point():
        #def jitpolicy(*args):
        #   from pypy.jit.metainterp.policy import JitPolicy
        #   return JitPolicy()

        brain = Brain()
        brain.loop()
    if '--pypy' in sys.argv:
        from pypy.translator.interactive import Translation
        t = Translation( pypy_entry_point )
        ## NotImplementedError: --gcrootfinder=asmgcc requires standalone ##
        #t.config.translation.suggest(jit=True, jit_debug='steps', jit_backend='x86', gc='boehm')
        t.annotate()
        t.rtype()
        f = t.compile_c()
        f()
    else:
        pypy_entry_point()

else:
    if '--blender' in sys.argv:
        import bpy
        BlenderExport()

    else:
        a = App()
        gtk.main()
        print( '-------------------exit toplevel-----------------' )