--- a/ghost/nodes2/Node.py	Mon Jun 16 12:55:29 2008 -0400
+++ b/ghost/nodes2/Node.py	Wed Jun 18 13:43:42 2008 -0400
@@ -1,6 +1,17 @@ from numpy.linalg import norm
 from numpy.linalg import norm
 from numpy import array,concatenate,zeros,prod,exp,meshgrid,matrix,dot,abs
 import sys
+
+
+def get_output_vector(nodes,input_vector):
+    
+    infer=[]
+    for col in nodes:
+        infer.append(col.inference(input_vector))
+    x=array(concatenate(infer,axis=1),order='C')
+    
+    return x
+    
 
 def divide_input_2D(node_list,input_size,node_size,equal=(False,False),overlap=(0,0)):
     

--- a/ghost/nodes2/__init__.py	Mon Jun 16 12:55:29 2008 -0400
+++ b/ghost/nodes2/__init__.py	Wed Jun 18 13:43:42 2008 -0400
@@ -1,4 +1,5 @@ from Node import *
 from Node import *
+from BCL import *
 from LeaderFollower import *
 from TemporalAdjacency import *
-from IncrementalSequence import *
\ No newline at end of file
+from IncrementalSequence import *

--- a/ghost/nodes2/compile	Mon Jun 16 12:55:29 2008 -0400
+++ b/ghost/nodes2/compile	Wed Jun 18 13:43:42 2008 -0400
@@ -2,3 +2,4 @@
 ./cython_compile.py _Node.pyx
 ./cython_compile.py _LeaderFollower.pyx
 ./cython_compile.py _IncrementalSequence.pyx
+./cython_compile.py _BCL.pyx

--- a/ghost/test_simple.py	Mon Jun 16 12:55:29 2008 -0400
+++ b/ghost/test_simple.py	Wed Jun 18 13:43:42 2008 -0400
@@ -1,30 +1,89 @@ from nodes import *
-from nodes import *
+from nodes2 import *
 from inputs import *
 from pylab import *
 
-def test_RBCL_Node():
-    rbcl = RBCL_Node()
-    sig = 0.2
-    for i in range(200):
-        for v in input_vectors3(int(rand()*50),sig):
-            rbcl.update(v)
-            plot([v[0]],[v[1]],'g.')
-        for v in input_vectors1(int(rand()*50),sig):
-            rbcl.update(v)
-            plot([v[0]],[v[1]],'b.')
-    for c in rbcl.clusters:
-        plot([c[0]],[c[1]],'ro')    
-        
-    show()
-
-    return rbcl
-    
 
 def test_BCL_Layers():
     
     # each node must have an output
-    layer1=[BCL_Node() for i in range(5)]
-    layer2=[BCL_Node() for i in range(3)]
+    layer1=[BCL() for i in range(5)]
+    layer2=[BCL() for i in range(3)]
+
+    for col in layer2:
+        col.add_children(layer1)
+
+    for v in input_vectors(300,0.25):
+        for col in layer1:
+            col.update(v)
+
+    for col in layer1:
+        col.learning=False
+
+    for v in input_vectors(300,0.25):
+
+        x=get_output_vector(layer1,v)  # they are all the same
+
+        for col in layer2:
+#            x=get_output_vector(col.children,v)  # they are not all the same
+
+            col.update(x)
+
+    
+def test_BCL_Node():
+
+    col = BCL()
+    for v in input_vectors(300,0.25):
+        col.update(v)
+        plot([v[0]],[v[1]],'b.')
+    for c in col.clusters:
+        plot([c[0]],[c[1]],'ro')    
+    show()
+    
+    return col
+
+def test_LeaderFollower1():
+    col=LeaderFollower()
+    clf()
+    for v in input_vectors(500,0.25):
+        plot([v[0]],[v[1]],'b.')
+        col.update(v)
+    print 'Number of clusters:',len(col.clusters)
+    Y,X=mgrid[-3.:3.:0.1,-3:3:.1]
+    Y=-Y
+    Z=zeros(Y.shape)
+    count=0
+    for x,y in zip(X.ravel(),Y.ravel()):
+        p=col.inference(array([x,y]))
+        Z.ravel()[count]=argmax(p)-0.5+max(p)
+        count+=1
+
+    pcolor(X,Y,Z)
+    colorbar()
+    draw()
+    show()
+
+    return col
+
+def test_2D_laplace(node=BCL()):
+    """Tests the clustering on rotated laplace numbers."""
+    
+    clf()
+    for v in input_vectors2(1000):
+        node.update(v)
+        plot([v[0]],[v[1]],'b.')
+    for c in node.clusters:
+        plot([c[0]],[c[1]],'ro')    
+    axis('equal')
+    show()
+
+
+
+
+if __name__=="__main__":
+    
+    # each node must have an output
+    layer1=[BCL() for i in range(5)]
+    layer2=[BCL() for i in range(3)]
     
     for col in layer2:
         col.add_children(layer1)
@@ -37,64 +96,11 @@ def test_BCL_Layers():
         col.learning=False
                          
     for v in input_vectors(300,0.25):
-        for col in layer1:
-            col.update(v)
+        
+        x=get_output_vector(layer1,v)  # they are all the same
+            
         for col in layer2:
-            col.update()
+#            x=get_output_vector(col.children,v)  # they are not all the same
+            
+            col.update(x)
 
-    return col
-    
-def test_BCL_Node():
-
-    col = BCL_Node()
-    for v in input_vectors(300,0.25):
-        col.update(v)
-        plot([v[0]],[v[1]],'b.')
-    for c in col.clusters:
-        plot([c[0]],[c[1]],'ro')    
-    show()
-    
-    return col
-
-def test_LeaderFollower1():
-    col=LeaderFollower_Node()
-    clf()
-    for v in input_vectors(500,0.25):
-        plot([v[0]],[v[1]],'b.')
-        col.update(v)
-    print 'Number of clusters:',len(col.clusters)
-    Y,X=mgrid[-3.:3.:0.1,-3:3:.1]
-    Y=-Y
-    Z=zeros(Y.shape)
-    count=0
-    for x,y in zip(X.ravel(),Y.ravel()):
-        p=col.inference([x,y])
-        Z.ravel()[count]=argmax(p)-0.5+max(p)
-        count+=1
-
-    pcolor(X,Y,Z)
-    colorbar()
-    draw()
-    show()
-
-    return col
-
-def test_2D_laplace(node=BCL_Node()):
-    """Tests the clustering on rotated laplace numbers."""
-    
-    clf()
-    for v in input_vectors2(1000):
-        node.update(v)
-        plot([v[0]],[v[1]],'b.')
-    for c in node.clusters:
-        plot([c[0]],[c[1]],'ro')    
-    axis('equal')
-    show()
-
-
-
-
-if __name__=="__main__":
-    
-    node=test_RBCL_Node()
-    print node
\ No newline at end of file

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/ghost/nodes2/BCL.py	Wed Jun 18 13:43:42 2008 -0400
@@ -0,0 +1,205 @@
+import numpy
+from numpy import dot
+from pylab import norm,exp,arccos,pi
+from Node import Node
+
+try:
+    from _BCL import get_winner
+except ImportError:
+    pass
+
+
+class BCL(Node):
+    def __init__(self,alpha=1.0,
+                      eta=0.05,
+                      lambd=None,
+                      dist_thresh=0.1,
+                      angle_thresh=90,
+                      activation_thresh=10.0,
+                      delta=1.0,
+                      beta=0.01,
+                      max_N=50):
+
+        super(BCL,self).__init__()
+
+        self.alpha = alpha # time decay factor
+        self.eta = eta # learning rate
+        self.dist_thresh = dist_thresh
+        self.angle_thresh = angle_thresh
+        self.activation_thresh = activation_thresh
+        self.delta = delta # constant activation increment
+        self.beta = beta # minimum frequency for activation
+        self.activations = []
+        self.count = [] # count number of activations
+        self.history = [0] # number of clusters in past (2 timesteps)
+        self.dyn_count = 0 # count of sequential timesteps that len(clusters) has been constant
+        self.y = []
+        self.y_prev = []
+        self.max_N=max_N
+
+        if lambd:
+            self.lambd = lambd # forgetting factor 0 < lambd < 1
+            self.adaptive_lambd=False
+        else:
+            self.adaptive_lambd=True
+        self.sigma=dist_thresh
+        self.N=0
+        
+    def update_feedforward(self,input_vector):
+
+        self.N+=1
+        self.lambd=1.0/self.N
+        if len(self.clusters) == 0:
+            self.LOG('Initializing clusters.')
+            
+            self.clusters = [input_vector]
+            self.activations = [0]
+            self.count = [0]
+            self.y = [0]
+            self.y_prev = [0]
+            
+        self.y_prev=self.y[:]
+        try:
+            min_clust,min_dist,y=get_winner(input_vector,
+                                             array(self.y),
+                                             self.clusters,
+                                             self.alpha)
+            self.y=list(y)
+            
+        except NameError:  # not using cython            
+            for i,cluster in enumerate(self.clusters):
+                try:
+                    self.y[i] = (1-self.alpha)*self.y[i] + self.alpha*(input_vector - cluster)
+                except ValueError:
+                    print input_vector,cluster
+                    raise ValueError
+            
+                distance = norm(self.y[i])
+                if i==0:
+                    min_dist = distance
+                    min_clust = i
+                elif distance < min_dist:
+                    min_dist = distance
+                    min_clust = i
+
+        self.LOG('Input vector',input_vector)
+        self.LOG('Clusters',self.clusters)
+        self.LOG('Min distance',min_dist,"at cluster",min_clust)
+
+        if self.learning: #or True:            
+            # update activation level of winner
+            diffc = self.y[min_clust]
+            diffl = self.y_prev[min_clust]
+            
+            nc=norm(diffc)
+            nl=norm(diffl)
+            d=min(nc,nl)
+            
+            if nc==0.0:
+                dotty=0.0
+            elif nl==0.0:
+                dotty=0.0
+            else:
+                dotty=dot(diffc,diffl)/norm(diffc)/norm(diffl)
+                
+            if dotty>1:
+                if dotty-1e-10<1:
+                    dotty=1.0
+            elif dotty<-1:
+                if dotty+1e-10>-1:
+                    dotty=-1.0
+                
+            if dotty>1 or dotty<-1:
+                self.LOG('Invalid vectors for arccos')
+                self.LOG(diffc)
+                self.LOG(diffl)
+                self.LOG(d)
+                pdb.set_trace()
+                raise ValueError,"Invalid arccos"
+                
+            phi=arccos(dotty)*180/pi
+            
+            if self.DEBUG:
+                self.LOG('Valid Activation Test...')
+                if phi < 0:
+                    self.LOG('  Angle...passed.',phi)
+                else:
+                    self.LOG('  Angle...failed.',phi)
+            
+                d=min(norm(diffc),norm(diffl))
+                if d>self.dist_thresh:
+                    self.LOG('  Distance...passed.',d)
+                else:
+                    self.LOG('  Distance...failed.',d)
+            
+            if phi >= self.angle_thresh and d >= self.dist_thresh:
+                self.activations[min_clust] += self.delta
+                self.LOG('  Adding')
+            else:
+                self.activations[min_clust] -= self.lambd*self.activations[min_clust]
+                self.LOG('  Decay')
+            
+            self.LOG('  Activations',self.activations)
+
+            # add new weight if activation level of winner exceeds threshold
+            if self.activations[min_clust] > self.activation_thresh:
+                self.activations[min_clust] = 0
+                self.clusters += [self.clusters[min_clust] + self.eta*self.y[min_clust]]
+                self.activations += [0]
+                self.count += [1]
+                self.y += [0]
+                self.y_prev += [0]
+
+            else: # else update weight
+                self.clusters[min_clust] += self.eta*self.y[min_clust]
+                self.count[min_clust] += 1
+
+
+            # decrease activation levels of all other weights
+            for i in xrange(len(self.clusters)):
+                if i != min_clust:
+                        self.activations[i] = (1-self.lambd)*self.activations[i]
+
+            # delete weights with activation frequency below threshold beta
+            i = 0
+            N = sum(self.count)
+            while True:
+                if i==len(self.clusters): break
+                elif float(self.count[i])/N < self.beta:
+                    self.clusters.pop(i)
+                    self.activations.pop(i)
+                    self.count.pop(i)
+                    self.y.pop(i)
+                    self.y_prev.pop(i)
+                else: i+=1
+
+            # check to see if network has reached dynamic equilibrium
+            if len(self.history)==2:
+                self.history.pop(0)
+            self.history.append(len(self.clusters))
+            if self.history[1] == self.history[0]:
+                self.dyn_count+=1
+            else: self.dyn_count = 0
+            if self.dyn_count >= self.max_N: 
+                self.learning=False
+#                print 'Dynamic equilibrium reached.'
+
+        else: # dynamic equilibrium has been reached
+            self.clusters[min_clust] += self.eta*self.y[min_clust]
+
+
+    def copy(self):
+        newnode = BCL(self,alpha=self.alpha,
+                          eta=self.eta,
+                          lambd=self.lambd,
+                          dist_thresh=self.dist_thresh,
+                          angle_thresh=self.angle_thresh,
+                          activation_thresh=self.activation_thresh,
+                          delta=self.delta,
+                          beta=self.beta,
+                          max_N=self.max_N)
+        newnode.clusters=self.clusters[:]
+        newnode.learning=self.learning
+        
+
+

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/ghost/nodes2/_BCL.pyx	Wed Jun 18 13:43:42 2008 -0400
@@ -0,0 +1,57 @@
+cimport c_python
+cimport c_numpy
+# Numpy must be initialized
+c_numpy.import_array()
+
+from numpy import zeros
+
+cdef here(c):
+    import sys
+    print c
+    sys.stdout.flush()
+
+cdef extern from "math.h":
+    double floor(double)
+    double exp(double)
+    double log(double)
+    double tanh(double)
+    double sqrt(double)
+    double fabs(double)
+
+cpdef get_winner(c_numpy.ndarray input_vector, c_numpy.ndarray y, clusters, double alpha):
+
+    cdef int n=len(clusters)
+    cdef int i,j,min_clust
+    cdef int inp_len=input_vector.dimensions[0]
+    cdef c_numpy.ndarray cluster,y_i
+    cdef double *cluster_p
+    cdef double *y_p
+    cdef double *input_vector_p=<double *>input_vector.data
+    cdef double distance_squared,min_distance_squared
+    cdef double min_distance
+
+    min_distance_squared=1e500
+    min_clust=0
+
+    for i from 0<=i<n:
+        cluster=clusters[i]
+        y_i=y[i]
+        cluster_p=<double *> cluster.data
+        y_p=<double *>y_i.data
+
+
+        distance_squared=0.0
+        for j from 0<=j<inp_len: # dot product
+            y_i[j]=(1-alpha)*y_i[j] + alpha*(input_vector_p[j]-cluster_p[j])
+            distance_squared+=y_i[j]*y_i[j]
+
+        if distance_squared < min_distance_squared:
+            min_distance_squared = distance_squared
+            min_clust = i
+
+        y[i]=y_i
+
+    min_distance=sqrt(min_distance_squared)
+
+    return min_clust,min_distance,y
+