2013/uns13/resources: optsparse.py

File optsparse.py, 3.6 KB (added by drasmuss, 5 years ago)

code for creating sparse connection weights

Line 
1"""
2Implementation of sparse connectivity for NEF/SpiNNaker
3The Spinn constructor can be passed a max_fan_in parameters,
4which limits the number of connections a neuron can receive in the network.
5
6Author: Terry Stewart
7"""
8
9import numeric as np
10import random
11import ca.nengo.util.MU as MU
12#import Jama.Matrix as Matrix
13import subprocess
14import os
15import struct
16from array import array
17from java.lang import System
18
19#uncomment this function to calculate the pseudoinverse directly within
20#nengo's jython
21#def pinv(matrix, limit):
22#    m = Matrix([[x for x in row] for row in matrix])
23#    svd = m.svd()
24#    sInv = svd.getS().inverse()
25#   
26#    i = 0
27#    while i<svd.getS().rowDimension and svd.getS().get(i, i)>limit:
28#        i += 1
29#    while i<len(matrix):
30#        sInv.set(i, i, 0)
31#        i+=1
32#    result = svd.getV().times(sInv).times(svd.getU().transpose()).getArray()
33#    return result
34   
35#this function calculates the pseudoinverse through an external callout
36#to a python file, allowing the use of numpy
37def pinv(matrix, limit):
38    shortfile = "matrix_%f" % random.random()
39    filename = os.path.join("external", shortfile)
40   
41    #write matrix data to file
42    infile = file(filename, "wb")
43    float_array = array('f', [x for row in matrix for x in row])
44    float_array.tofile(infile)
45    infile.close()
46   
47    #create output file
48    outfile = file(filename+".inv", "wb")
49    outfile.close()
50   
51#    print System.getProperty("os.name")
52    if "windows" in System.getProperty("os.name").lower():
53        subprocess.call("cmd /c pseudoInverse.bat "+shortfile+" "+shortfile+".inv"+" "+str(limit)+" "+str(-1),
54                        stdout=subprocess.PIPE, stderr=subprocess.PIPE, cwd="external")
55    else:
56        subprocess.call("external/pseudoInverse "+filename+" "+filename+".inv"+" "+str(limit)+" "+str(-1), shell=True)
57   
58    mfile = file(filename+".inv", "rb")
59    invmatrix = array("f")
60    invmatrix.fromfile(mfile, len(matrix)*len(matrix[0]))
61    return [[invmatrix[i*len(matrix[0])+j] for j in range(len(matrix[0]))] for i in range(len(matrix))]
62       
63def compute_sparse_weights(origin, post, transform, fan_in, noise=0.1, num_samples=100):
64    encoder = post.encoders
65    radius = post.radii[0]
66   
67    if hasattr(transform, 'tolist'): transform=transform.tolist()
68   
69    approx = origin.node.getDecodingApproximator('AXON')   
70   
71    # create X matrix
72    X = approx.evalPoints       
73    X = MU.transpose([f.multiMap(X) for f in origin.functions])
74   
75    # create A matrix
76    A = approx.values
77   
78    S = fan_in
79    N_A = len(A)
80    samples = len(A[0])
81    N_B = len(encoder)
82    w_sparse = np.zeros((N_B, N_A),'f')
83    noise_sd = MU.max(A)*noise
84    num_samples = min(num_samples,N_B)
85    decoder_list = [None for _ in range(num_samples)]
86    for i in range(num_samples):
87        indices = random.sample(range(N_A), S)
88        activity = [A[j] for j in indices]
89        n = [[random.gauss(0, noise_sd) for _ in range(samples)] for j in range(S)]
90        activity = MU.sum(activity, n)
91        activityT = MU.transpose(activity)
92        gamma = MU.prod(activity, activityT)
93       
94        upsilon = MU.prod(activity, X)
95       
96        gamma_inv = pinv(gamma, noise_sd*noise_sd)
97       
98       
99        decoder_list[i] = MU.prod([[x for x in row] for row in gamma_inv], upsilon)
100   
101    for i in range(N_B):
102        ww = MU.prod(random.choice(decoder_list), MU.prod(MU.transpose(transform), encoder[i]))
103       
104        for j, k in enumerate(indices):
105            w_sparse[i,k] = float(ww[j])/radius
106   
107    return list(w_sparse)
108