Notebook

In [67]:

%pylab inline
rcParams['figure.figsize'] = 14, 7

Populating the interactive namespace from numpy and matplotlib

In [68]:

from scipy.integrate import simps, quad


x = np.linspace(-3.0, 12.0, 300)
R = 0.1
sigma = 1.0  #conductance
true_csd = 1.0*np.exp(-(x-2.)**2 / (2*np.pi*0.5)) + 0.5*np.exp(-(x-7)**2 / (2*np.pi*1.0))

#print true_csd

fig, (ax1, ax2) = plt.subplots(1, 2, sharex=True)
ax1.plot(x, true_csd)
ax1.set_title('True CSD')

def calculate_pot(csd, z, z0):
    return 1.0/(2 *sigma) * np.trapz( (np.sqrt((z-z0)**2 + R**2 ) - np.abs(z-z0))*csd, z)

pots = [calculate_pot(true_csd, x, x0) for x0 in x]
ax2.plot(x, pots)
ax2.set_title('Potentials')

Out[68]:

<matplotlib.text.Text at 0xd656290>

In [69]:

#elec_pos = np.linspace(-3.0, 13.0, 30)
elec_pos = [ -2.0, -1.0, 0.0, 2.2, 4.3, 7.3, 9.0]

meas_pot = np.array([calculate_pot(true_csd, x, x0) for x0 in elec_pos])

plot(x, pots)
scatter(elec_pos, meas_pot)

Out[69]:

<matplotlib.collections.PathCollection at 0xdbd59d0>

In [70]:

simple_csd = np.array([ -(meas_pot[i-1]+meas_pot[i+1] -2*meas_pot[i])/(elec_pos[i]-elec_pos[i-1])**2 for i in xrange(1, len(elec_pos)-1)])

fig = plt.subplot(1,1,1)
plot(elec_pos[1:len(elec_pos)-1], simple_csd, label="approximated csd")
plot(x, true_csd, label="true csd")
fig.legend(loc=(0.01,0.7))

Out[70]:

<matplotlib.legend.Legend at 0xdae07d0>

kCSD workflow¶

In [71]:

n_sources = 100

source_space = linspace(min(elec_pos), max(elec_pos), n_sources)
#source_space = linspace(-4.4, 4.4, n_sources)

In [72]:

#define the positions of the basis functions
xmax=15.0
xmin=-3.0
estimation_area = linspace(xmin, xmax, n_sources)
max_dist = xmax-xmin

dist_density = 300
h = 1.3*(elec_pos[1]-elec_pos[0])
R = 0.15

In [73]:

def pot_intarg(src, arg, current_pos, h, R, sigma, src_type):
    """
    """
    if src_type == "gaussian":
        y = (1./(2*sigma))*(sqrt((arg-current_pos)**2 + R**2) - abs(arg - current_pos)) * gauss_rescale(src, current_pos, h)
        # for this formula look at formula (8) from Pettersen et al., 2006
    if src_type == "step":
        #TODO re-evaluate
        y = (1./(2*sigma))*(sqrt((arg-current_pos)**2+R**2) - abs(arg-current_pos)) * (abs(src-current_pos)<0.5*h)
    if src_type == "wavelet":
        pass
            
    return y

In [74]:

def gauss_rescale(x, mu, three_stdev):
    """
    returns normalized gaussian scale function centered in mu, cut off farther than three_stdev from the center
    """
    variance = (three_stdev/3.0)**2
    g = 1./sqrt(2. * pi * variance) * exp( -(1./(2.*variance)) * (x-mu)**2 ) #* (abs(x-mu)<three_stdev)
    return g

In [75]:

def b_pot_quad(src, arg, h, R, sigma, src_type):
    """
    """
    if src_type in ["gaussian", "step"]:
        x = linspace(src - 4*h, src + 4*h, 150) #manipulate the resolution
        y = [pot_intarg(src, arg, current_pos, h, R, sigma, src_type) for current_pos in x]
        plot(x,y)
        z = np.trapz(y, x)
    else:
        raise error("Source type not implemented!")

    return z

In [76]:

def create_dist_table(dist_dens, dist_max, h, R, sigma, src_type):
    """
    Creates table of source contribution to overall potential as a function of distance
    """
    dist_table = zeros(dist_dens+1)

    for i in xrange(0, dist_dens+1):
        dist_table[i] = b_pot_quad(0,(float(i)/dist_dens) * dist_max, h, R, sigma, src_type)
    
    return dist_table

dist_table = create_dist_table(dist_density, max_dist, h, R, sigma, src_type = "gaussian")
#print dist_table
figure()
plot(dist_table)

Out[76]:

[<matplotlib.lines.Line2D at 0xeae4310>]

In [77]:

def calculate_b_pot_matrix(x, src, elec_pos, dist_table, dist_dens):
    """ 
    compute the matrix of potentials generated by every source basis function
    at every electrode position. This procedure can be optimized
    bPotMatrix(j,k) = b_j(x_k) 
    with x_k electrode position, b_j j-th basis source from the paper
    """
    n_elec = len(elec_pos)
    
    b_pot_matrix = zeros((n_sources, n_elec))
    
    for src_ind, current_src in enumerate(src):
        for arg_ind, arg in enumerate(elec_pos):
            r = abs(current_src - arg)
            
            #dist_table_ind = uint16(0.85 * r * float(dist_dens)/max_dist)
            dist_table_ind = uint16(r * float(dist_dens)/max_dist)
            #print dist_table_ind
            b_pot_matrix[src_ind, arg_ind] = dist_table[dist_table_ind]
            
    return b_pot_matrix
  
b_pot_matrix = calculate_b_pot_matrix( estimation_area, source_space, elec_pos, dist_table, dist_density)
#print bpm
imshow(b_pot_matrix, interpolation="none")
colorbar()

Out[77]:

<matplotlib.colorbar.Colorbar instance at 0xee99248>

In [78]:

k_pot = dot(transpose(b_pot_matrix), b_pot_matrix) # Kernel matrix K from (2.14) from the paper

In [79]:

imshow(k_pot, interpolation="none")
colorbar()

Out[79]:

<matplotlib.colorbar.Colorbar instance at 0xee66c68>

In [80]:

def calculate_b_src_matrix(x, src, h):
    n_gx = len(x)

    b_src_matrix = zeros((n_gx, n_sources))

    for src_ind, curr_src in enumerate(src):
        b_src_matrix[:, src_ind] = gauss_rescale(x, curr_src, h)

    return b_src_matrix
    
b_src_matrix = calculate_b_src_matrix( estimation_area, source_space, h)

imshow(b_src_matrix, interpolation="none")
colorbar()

Out[80]:

<matplotlib.colorbar.Colorbar instance at 0xe959a70>

In [81]:

interp_cross = dot(b_src_matrix, b_pot_matrix) # cross-kernel matrix from (2.17) in the paper
imshow(interp_cross, interpolation="none")
colorbar()

Out[81]:

<matplotlib.colorbar.Colorbar instance at 0xed6d560>

In [82]:

def calculate_b_interp_pot_matrix(x, src, dist_table, dist_dens):
    """
    compute the matrix of potentials generated by every source basis function
    at every electrode position. This procedure can be optimized
    bPotMatrix(j,k) = b_j(x_k) 
    with x_k electrode position, b_j j-th basis source from the paper
    """
    n_points = len(x)
    dist_max = max(x) - min(x)

    b_interp_pot_matrix = np.zeros((n_sources, n_points))
    
    for src_ind, current_src in enumerate(src):
        for arg_ind, arg in enumerate(x):
            r = abs(current_src - arg);
        
            #dist_table_ind = uint16(0.85*r*float(dist_dens)/dist_max)   
            dist_table_ind = uint16(r*float(dist_dens)/dist_max)
            b_interp_pot_matrix[src_ind, arg_ind] = dist_table[dist_table_ind]
    
    return b_interp_pot_matrix

b_interp_pot_matrix = calculate_b_interp_pot_matrix( estimation_area, source_space, dist_table, dist_density)

imshow(b_interp_pot_matrix, interpolation="none")
colorbar()

Out[82]:

<matplotlib.colorbar.Colorbar instance at 0xe6d2560>

In [83]:

interp_pot = dot(transpose(b_interp_pot_matrix), b_pot_matrix)
imshow(interp_pot, interpolation="none")
colorbar()

Out[83]:

<matplotlib.colorbar.Colorbar instance at 0xf4536c8>

In [84]:

from sklearn.cross_validation import KFold, LeaveOneOut,ShuffleSplit

lambd=0.0

def calc_err(lambd, pot, k_pot, ind_test, ind_train):
    k_train = k_pot[ind_train, ind_train]
    
    #print "indtest",ind_test
    #print "indtrain", ind_train
    pot_train = pot[ind_train]
    pot_test = pot[ind_test]
    
    beta = np.dot(np.linalg.inv(k_train + lambd*identity(k_train.shape[0])) , pot_train)
    #print "beta.sh", beta.shape
    
    #plot(beta)
    
    #k_cross = np.array([[k_pot[i,j] for i in ind_test] for j in ind_train])
    #k_cross = k_pot[np.array([ind_test, ind_train])]
    k_cross = k_pot[ind_test][:,ind_train]
    #imshow(k_cross, interpolation="none")
    
    pot_est = np.zeros(size(transpose(ind_test)))
    
    pot_est = dot(k_cross, beta)
    
    #plot(pot_est, pot_test)
    print "pottest:", pot_test, "pot_est:", pot_est
    err = linalg.norm(pot_test - pot_est)
    
    return err

#print calc_err(0.0, estimate_pots(), k_pot, np.array([range(0,len(elec_pos)-2)]), np.array([len(elec_pos)-2]))


def cross_validation(lambd, pot, k_pot, n_folds):
    errors = []
    n = len(elec_pos)
    #print "n:",n
    kf = KFold(n, n_folds=n_folds, indices=True)
    #kf = LeaveOneOut(n, indices=True)
    #kf = ShuffleSplit(5, n_iter=15, test_size=0.25, indices=True)
    for ind_train, ind_test in kf:
        #ind_test, ind_train = np.random.permutation(ind_test), np.random.permutation(ind_train)
        err = calc_err(lambd, pot, k_pot, ind_test, ind_train)
        #print err
        errors.append(err)
        
        #print("%s %s" % (ind_train, ind_test))
        
    error = mean(errors)
    print "l=", lambd, ", err=", error
    return error

def choose_lambda(lambdas, n_folds, n_iter):
    """finds the optimal regularization parameter lambda for Tikhonov regularization"""
    n = len(lambdas)
    errors = np.zeros(n)
    errors_iter = np.zeros(n_iter)
    for i, lambd in enumerate(lambdas):
        for j in xrange(n_iter):
            errors_iter[j] = cross_validation(lambd, meas_pot, k_pot, n_folds)
            #errors_iter[j] = calc_cv_error(lambd, n_folds)
        errors[i] = mean(errors_iter)
    fig = pyplot.figure()
    ax = fig.add_subplot(1,1,1)
    plot(lambdas, errors)
    ax.set_xscale('log')
    ax.set_yscale('log')
    return lambdas[errors == min(errors)]

lambdas = np.array([10./j**2 for j in xrange(1,60) ])
#lambdas = np.array([0,1])
lambd = choose_lambda(lambdas, n_folds = 5,  n_iter=1)[0]
#lambd = 0.0
print "lambda =", lambd

pottest: [ 0.00296697  0.00468031] pot_est: [  5.13806253e-06   8.68440536e-06]
pottest: [ 0.00903358  0.0202294 ] pot_est: [  4.93689518e-06   5.03892972e-06]
pottest: [ 0.01113273] pot_est: [  8.63698763e-06]
pottest: [ 0.01176185] pot_est: [  6.44005818e-06]
pottest: [ 0.00767375] pot_est: [  5.52222352e-06]
l= 10.0 , err= 0.0116454618352
pottest: [ 0.00296697  0.00468031] pot_est: [  2.03703207e-05   3.44171111e-05]
pottest: [ 0.00903358  0.0202294 ] pot_est: [  1.95372037e-05   1.99989385e-05]
pottest: [ 0.01113273] pot_est: [  3.42617946e-05]
pottest: [ 0.01176185] pot_est: [  2.55108432e-05]
pottest: [ 0.00767375] pot_est: [  2.18738925e-05]
l= 2.5 , err= 0.0116233523446
pottest: [ 0.00296697  0.00468031] pot_est: [  4.51664472e-05   7.62638239e-05]
pottest: [ 0.00903358  0.0202294 ] pot_est: [  4.31858417e-05   4.44216441e-05]
pottest: [ 0.01113273] pot_est: [  7.60429023e-05]
pottest: [ 0.01176185] pot_est: [  5.64877419e-05]
pottest: [ 0.00767375] pot_est: [  4.84320279e-05]
l= 1.11111111111 , err= 0.011587378367
pottest: [ 0.00296697  0.00468031] pot_est: [  7.86916946e-05   1.32753945e-04]
pottest: [ 0.00903358  0.0202294 ] pot_est: [  7.49089149e-05   7.75812661e-05]
pottest: [ 0.01113273] pot_est: [ 0.00013268]
pottest: [ 0.01176185] pot_est: [  9.82375196e-05]
pottest: [ 0.00767375] pot_est: [  8.42276376e-05]
l= 0.625 , err= 0.0115387715222
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.00011987  0.000202  ]
pottest: [ 0.00903358  0.0202294 ] pot_est: [ 0.00011344  0.00011854]
pottest: [ 0.01113273] pot_est: [ 0.00020252]
pottest: [ 0.01176185] pot_est: [ 0.00014932]
pottest: [ 0.00767375] pot_est: [ 0.00012804]
l= 0.4 , err= 0.0114791073001
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.00016748  0.00028182]
pottest: [ 0.00903358  0.0202294 ] pot_est: [ 0.00015732  0.00016619]
pottest: [ 0.01113273] pot_est: [ 0.00028367]
pottest: [ 0.01176185] pot_est: [ 0.00020809]
pottest: [ 0.00767375] pot_est: [ 0.00017848]
l= 0.277777777778 , err= 0.0114101842978
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.00022019  0.00036991]
pottest: [ 0.00903358  0.0202294 ] pot_est: [ 0.00020496  0.00021937]
pottest: [ 0.01113273] pot_est: [ 0.00037416]
pottest: [ 0.01176185] pot_est: [ 0.00027283]
pottest: [ 0.00767375] pot_est: [ 0.00023415]
l= 0.204081632653 , err= 0.0113339012848
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.0002767   0.00046395]
pottest: [ 0.00903358  0.0202294 ] pot_est: [ 0.00025474  0.00027686]
pottest: [ 0.01113273] pot_est: [ 0.00047206]
pottest: [ 0.01176185] pot_est: [ 0.00034186]
pottest: [ 0.00767375] pot_est: [ 0.00029364]
l= 0.15625 , err= 0.0112521457687
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.00033574  0.00056172]
pottest: [ 0.00903358  0.0202294 ] pot_est: [ 0.00030508  0.0003375 ]
pottest: [ 0.01113273] pot_est: [ 0.00057554]
pottest: [ 0.01176185] pot_est: [ 0.00041359]
pottest: [ 0.00767375] pot_est: [ 0.0003557]
l= 0.123456790123 , err= 0.011166703392
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.00039617  0.00066119]
pottest: [ 0.00903358  0.0202294 ] pot_est: [ 0.00035453  0.00040019]
pottest: [ 0.01113273] pot_est: [ 0.00068296]
pottest: [ 0.01176185] pot_est: [ 0.0004866]
pottest: [ 0.00767375] pot_est: [ 0.00041921]
l= 0.1 , err= 0.0110791925941
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.00045696  0.00076057]
pottest: [ 0.00903358  0.0202294 ] pot_est: [ 0.00040175  0.00046393]
pottest: [ 0.01113273] pot_est: [ 0.00079292]
pottest: [ 0.01176185] pot_est: [ 0.00055966]
pottest: [ 0.00767375] pot_est: [ 0.00048322]
l= 0.0826446280992 , err= 0.0109910246641
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.00051725  0.00085834]
pottest: [ 0.00903358  0.0202294 ] pot_est: [ 0.00044557  0.00052787]
pottest: [ 0.01113273] pot_est: [ 0.00090426]
pottest: [ 0.01176185] pot_est: [ 0.00063179]
pottest: [ 0.00767375] pot_est: [ 0.00054699]
l= 0.0694444444444 , err= 0.0109033862498
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.00057633  0.00095324]
pottest: [ 0.00903358  0.0202294 ] pot_est: [ 0.00048502  0.00059126]
pottest: [ 0.01113273] pot_est: [ 0.00101604]
pottest: [ 0.01176185] pot_est: [ 0.00070216]
pottest: [ 0.00767375] pot_est: [ 0.00060991]
l= 0.0591715976331 , err= 0.0108172397124
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.00063364  0.00104427]
pottest: [ 0.00903358  0.0202294 ] pot_est: [ 0.00051928  0.0006535 ]
pottest: [ 0.01113273] pot_est: [ 0.00112758]
pottest: [ 0.01176185] pot_est: [ 0.00077016]
pottest: [ 0.00767375] pot_est: [ 0.00067158]
l= 0.0510204081633 , err= 0.0107333362349
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.00068873  0.00113068]
pottest: [ 0.00903358  0.0202294 ] pot_est: [ 0.0005477   0.00071412]
pottest: [ 0.01113273] pot_est: [ 0.00123835]
pottest: [ 0.01176185] pot_est: [ 0.00083537]
pottest: [ 0.00767375] pot_est: [ 0.00073169]
l= 0.0444444444444 , err= 0.0106522369432
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.00074131  0.00121192]
pottest: [ 0.00903358  0.0202294 ] pot_est: [ 0.00056981  0.00077275]
pottest: [ 0.01113273] pot_est: [ 0.00134801]
pottest: [ 0.01176185] pot_est: [ 0.00089748]
pottest: [ 0.00767375] pot_est: [ 0.00079006]
l= 0.0390625 , err= 0.0105743381043
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.00079116  0.00128762]
pottest: [ 0.00903358  0.0202294 ] pot_est: [ 0.00058524  0.00082912]
pottest: [ 0.01113273] pot_est: [ 0.00145636]
pottest: [ 0.01176185] pot_est: [ 0.00095632]
pottest: [ 0.00767375] pot_est: [ 0.0008466]
l= 0.0346020761246 , err= 0.0104998974362
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.00083816  0.00135756]
pottest: [ 0.00903358  0.0202294 ] pot_est: [ 0.00059375  0.00088305]
pottest: [ 0.01113273] pot_est: [ 0.0015633]
pottest: [ 0.01176185] pot_est: [ 0.0010118]
pottest: [ 0.00767375] pot_est: [ 0.00090129]
l= 0.0308641975309 , err= 0.0104290594943
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.00088225  0.00142163]
pottest: [ 0.00903358  0.0202294 ] pot_est: [ 0.00059519  0.00093444]
pottest: [ 0.01113273] pot_est: [ 0.00166882]
pottest: [ 0.01176185] pot_est: [ 0.00106391]
pottest: [ 0.00767375] pot_est: [ 0.00095416]
l= 0.0277008310249 , err= 0.010361878885
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.00092343  0.00147984]
pottest: [ 0.00903358  0.0202294 ] pot_est: [ 0.00058949  0.00098322]
pottest: [ 0.01113273] pot_est: [ 0.00177297]
pottest: [ 0.01176185] pot_est: [ 0.00111272]
pottest: [ 0.00767375] pot_est: [ 0.00100527]
l= 0.025 , err= 0.0102983406665
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.00096173  0.00153225]
pottest: [ 0.00903358  0.0202294 ] pot_est: [ 0.00057665  0.0010294 ]
pottest: [ 0.01113273] pot_est: [ 0.00187587]
pottest: [ 0.01176185] pot_est: [ 0.0011583]
pottest: [ 0.00767375] pot_est: [ 0.00105473]
l= 0.0226757369615 , err= 0.0102383777243
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.00099722  0.00157897]
pottest: [ 0.00903358  0.0202294 ] pot_est: [ 0.0005567  0.001073 ]
pottest: [ 0.01113273] pot_est: [ 0.00197764]
pottest: [ 0.01176185] pot_est: [ 0.00120078]
pottest: [ 0.00767375] pot_est: [ 0.00110264]
l= 0.0206611570248 , err= 0.0101818851949
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.00102999  0.00162017]
pottest: [ 0.00903358  0.0202294 ] pot_est: [ 0.00052972  0.00111409]
pottest: [ 0.01113273] pot_est: [ 0.00207844]
pottest: [ 0.01176185] pot_est: [ 0.00124029]
pottest: [ 0.00767375] pot_est: [ 0.00114913]
l= 0.0189035916824 , err= 0.0101287321729
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.00106014  0.00165602]
pottest: [ 0.00903358  0.0202294 ] pot_est: [ 0.00049582  0.00115274]
pottest: [ 0.01113273] pot_est: [ 0.00217845]
pottest: [ 0.01176185] pot_est: [ 0.00127699]
pottest: [ 0.00767375] pot_est: [ 0.00119433]
l= 0.0173611111111 , err= 0.0100787710243
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.0010878   0.00168673]
pottest: [ 0.00903358  0.0202294 ] pot_est: [ 0.00045513  0.00118904]
pottest: [ 0.01113273] pot_est: [ 0.00227784]
pottest: [ 0.01176185] pot_est: [ 0.00131103]
pottest: [ 0.00767375] pot_est: [ 0.00123836]
l= 0.016 , err= 0.0100318446511
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.00111307  0.00171251]
pottest: [ 0.00903358  0.0202294 ] pot_est: [ 0.00040777  0.00122309]
pottest: [ 0.01113273] pot_est: [ 0.0023768]
pottest: [ 0.01176185] pot_est: [ 0.00134257]
pottest: [ 0.00767375] pot_est: [ 0.00128135]
l= 0.0147928994083 , err= 0.00998779204739
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.00113607  0.00173357]
pottest: [ 0.00903358  0.0202294 ] pot_est: [ 0.00035389  0.00125501]
pottest: [ 0.01113273] pot_est: [ 0.00247549]
pottest: [ 0.01176185] pot_est: [ 0.00137175]
pottest: [ 0.00767375] pot_est: [ 0.00132343]
l= 0.0137174211248 , err= 0.00994645245792
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.00115693  0.00175011]
pottest: [ 0.00903358  0.0202294 ] pot_est: [ 0.00029363  0.00128489]
pottest: [ 0.01113273] pot_est: [ 0.0025741]
pottest: [ 0.01176185] pot_est: [ 0.00139873]
pottest: [ 0.00767375] pot_est: [ 0.00136472]
l= 0.0127551020408 , err= 0.00990766841169
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.00117577  0.00176234]
pottest: [ 0.00903358  0.0202294 ] pot_est: [ 0.00022715  0.00131285]
pottest: [ 0.01113273] pot_est: [ 0.00267278]
pottest: [ 0.01176185] pot_est: [ 0.00142366]
pottest: [ 0.00767375] pot_est: [ 0.00140531]
l= 0.0118906064209 , err= 0.00987128786569
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.00119268  0.00177046]
pottest: [ 0.00903358  0.0202294 ] pot_est: [ 0.00015458  0.001339  ]
pottest: [ 0.01113273] pot_est: [ 0.00277169]
pottest: [ 0.01176185] pot_est: [ 0.00144665]
pottest: [ 0.00767375] pot_est: [ 0.00144532]
l= 0.0111111111111 , err= 0.00983716565287
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.00120779  0.00177466]
pottest: [ 0.00903358  0.0202294 ] pot_est: [  7.60643815e-05   1.36343547e-03]
pottest: [ 0.01113273] pot_est: [ 0.00287097]
pottest: [ 0.01176185] pot_est: [ 0.00146786]
pottest: [ 0.00767375] pot_est: [ 0.00148485]
l= 0.0104058272633 , err= 0.00980516439465
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.00122119  0.00177512]
pottest: [ 0.00903358  0.0202294 ] pot_est: [ -8.26025112e-06   1.38625998e-03]
pottest: [ 0.01113273] pot_est: [ 0.00297077]
pottest: [ 0.01176185] pot_est: [ 0.00148738]
pottest: [ 0.00767375] pot_est: [ 0.00152399]
l= 0.009765625 , err= 0.00977515500684
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.00123298  0.00177201]
pottest: [ 0.00903358  0.0202294 ] pot_est: [ -9.82663115e-05   1.40757005e-03]
pottest: [ 0.01113273] pot_est: [ 0.00307122]
pottest: [ 0.01176185] pot_est: [ 0.00150534]
pottest: [ 0.00767375] pot_est: [ 0.00156281]
l= 0.00918273645546 , err= 0.00974701690203
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.00124325  0.00176548]
pottest: [ 0.00903358  0.0202294 ] pot_est: [-0.00019383  0.00142746]
pottest: [ 0.01113273] pot_est: [ 0.00317243]
pottest: [ 0.01176185] pot_est: [ 0.00152184]
pottest: [ 0.00767375] pot_est: [ 0.0016014]
l= 0.00865051903114 , err= 0.00972063796931
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.00125208  0.00175568]
pottest: [ 0.00903358  0.0202294 ] pot_est: [-0.00029484  0.00144601]
pottest: [ 0.01113273] pot_est: [ 0.00327451]
pottest: [ 0.01176185] pot_est: [ 0.00153698]
pottest: [ 0.00767375] pot_est: [ 0.00163984]
l= 0.00816326530612 , err= 0.00969591439487
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.00125956  0.00174276]
pottest: [ 0.00903358  0.0202294 ] pot_est: [-0.00040118  0.0014633 ]
pottest: [ 0.01113273] pot_est: [ 0.00337758]
pottest: [ 0.01176185] pot_est: [ 0.00155084]
pottest: [ 0.00767375] pot_est: [ 0.00167818]
l= 0.00771604938272 , err= 0.00967275037207
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.00126575  0.00172685]
pottest: [ 0.00903358  0.0202294 ] pot_est: [-0.00051275  0.00147942]
pottest: [ 0.01113273] pot_est: [ 0.00348172]
pottest: [ 0.01176185] pot_est: [ 0.00156352]
pottest: [ 0.00767375] pot_est: [ 0.00171648]
l= 0.0073046018992 , err= 0.00965105773836
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.00127073  0.00170806]
pottest: [ 0.00903358  0.0202294 ] pot_est: [-0.00062946  0.00149443]
pottest: [ 0.01113273] pot_est: [ 0.00358703]
pottest: [ 0.01176185] pot_est: [ 0.00157508]
pottest: [ 0.00767375] pot_est: [ 0.00175481]
l= 0.00692520775623 , err= 0.00963075556727
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.00127456  0.00168651]
pottest: [ 0.00903358  0.0202294 ] pot_est: [-0.0007512  0.0015084]
pottest: [ 0.01113273] pot_est: [ 0.00369359]
pottest: [ 0.01176185] pot_est: [ 0.0015856]
pottest: [ 0.00767375] pot_est: [ 0.00179322]
l= 0.00657462195924 , err= 0.00961176973645
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.0012773  0.0016623]
pottest: [ 0.00903358  0.0202294 ] pot_est: [-0.00087791  0.0015214 ]
pottest: [ 0.01113273] pot_est: [ 0.00380147]
pottest: [ 0.01176185] pot_est: [ 0.00159514]
pottest: [ 0.00767375] pot_est: [ 0.00183175]
l= 0.00625 , err= 0.00959403248723
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.001279    0.00163553]
pottest: [ 0.00903358  0.0202294 ] pot_est: [-0.0010095   0.00153348]
pottest: [ 0.01113273] pot_est: [ 0.00391074]
pottest: [ 0.01176185] pot_est: [ 0.00160376]
pottest: [ 0.00767375] pot_est: [ 0.00187044]
l= 0.0059488399762 , err= 0.00957748198679
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.00127972  0.0016063 ]
pottest: [ 0.00903358  0.0202294 ] pot_est: [-0.0011459   0.00154469]
pottest: [ 0.01113273] pot_est: [ 0.00402146]
pottest: [ 0.01176185] pot_est: [ 0.00161153]
pottest: [ 0.00767375] pot_est: [ 0.00190935]
l= 0.00566893424036 , err= 0.00956206190094
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.00127949  0.00157467]
pottest: [ 0.00903358  0.0202294 ] pot_est: [-0.00128703  0.00155509]
pottest: [ 0.01113273] pot_est: [ 0.00413371]
pottest: [ 0.01176185] pot_est: [ 0.00161848]
pottest: [ 0.00767375] pot_est: [ 0.00194849]
l= 0.00540832882639 , err= 0.0095477209826
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.00127838  0.00154074]
pottest: [ 0.00903358  0.0202294 ] pot_est: [-0.00143285  0.00156473]
pottest: [ 0.01113273] pot_est: [ 0.00424752]
pottest: [ 0.01176185] pot_est: [ 0.00162467]
pottest: [ 0.00767375] pot_est: [ 0.00198791]
l= 0.0051652892562 , err= 0.00953441267966
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.00127641  0.00150456]
pottest: [ 0.00903358  0.0202294 ] pot_est: [-0.00158327  0.00157364]
pottest: [ 0.01113273] pot_est: [ 0.00436295]
pottest: [ 0.01176185] pot_est: [ 0.00163015]
pottest: [ 0.00767375] pot_est: [ 0.00202765]
l= 0.00493827160494 , err= 0.0095220947641
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.00127362  0.00146621]
pottest: [ 0.00903358  0.0202294 ] pot_est: [-0.00173826  0.00158188]
pottest: [ 0.01113273] pot_est: [ 0.00448006]
pottest: [ 0.01176185] pot_est: [ 0.00163494]
pottest: [ 0.00767375] pot_est: [ 0.00206772]
l= 0.0047258979206 , err= 0.00951072898341
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.00127006  0.00142575]
pottest: [ 0.00903358  0.0202294 ] pot_est: [-0.00189776  0.00158947]
pottest: [ 0.01113273] pot_est: [ 0.00459887]
pottest: [ 0.01176185] pot_est: [ 0.0016391]
pottest: [ 0.00767375] pot_est: [ 0.00210815]
l= 0.00452693526483 , err= 0.00950028073453
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.00126574  0.00138324]
pottest: [ 0.00903358  0.0202294 ] pot_est: [-0.00206173  0.00159646]
pottest: [ 0.01113273] pot_est: [ 0.00471944]
pottest: [ 0.01176185] pot_est: [ 0.00164265]
pottest: [ 0.00767375] pot_est: [ 0.00214897]
l= 0.00434027777778 , err= 0.00949071876011
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.0012607   0.00133873]
pottest: [ 0.00903358  0.0202294 ] pot_est: [-0.00223011  0.00160288]
pottest: [ 0.01113273] pot_est: [ 0.00484179]
pottest: [ 0.01176185] pot_est: [ 0.00164563]
pottest: [ 0.00767375] pot_est: [ 0.00219021]
l= 0.00416493127863 , err= 0.00948201486622
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.00125498  0.00129226]
pottest: [ 0.00903358  0.0202294 ] pot_est: [-0.00240287  0.00160875]
pottest: [ 0.01113273] pot_est: [ 0.00496597]
pottest: [ 0.01176185] pot_est: [ 0.00164806]
pottest: [ 0.00767375] pot_est: [ 0.00223187]
l= 0.004 , err= 0.00947414366086
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.00124858  0.00124388]
pottest: [ 0.00903358  0.0202294 ] pot_est: [-0.00257997  0.00161411]
pottest: [ 0.01113273] pot_est: [ 0.005092]
pottest: [ 0.01176185] pot_est: [ 0.00164997]
pottest: [ 0.00767375] pot_est: [ 0.00227399]
l= 0.00384467512495 , err= 0.00946708231197
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.00124155  0.00119363]
pottest: [ 0.00903358  0.0202294 ] pot_est: [-0.00276137  0.00161898]
pottest: [ 0.01113273] pot_est: [ 0.00521991]
pottest: [ 0.01176185] pot_est: [ 0.0016514]
pottest: [ 0.00767375] pot_est: [ 0.00231657]
l= 0.00369822485207 , err= 0.00946081032396
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.00123389  0.00114156]
pottest: [ 0.00903358  0.0202294 ] pot_est: [-0.00294704  0.00162339]
pottest: [ 0.01113273] pot_est: [ 0.00534973]
pottest: [ 0.01176185] pot_est: [ 0.00165235]
pottest: [ 0.00767375] pot_est: [ 0.00235964]
l= 0.00355998576006 , err= 0.00945530933156
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.00122564  0.0010877 ]
pottest: [ 0.00903358  0.0202294 ] pot_est: [-0.00313695  0.00162735]
pottest: [ 0.01113273] pot_est: [ 0.00548149]
pottest: [ 0.01176185] pot_est: [ 0.00165286]
pottest: [ 0.00767375] pot_est: [ 0.00240321]
l= 0.00342935528121 , err= 0.00945056290973
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.00121681  0.00103209]
pottest: [ 0.00903358  0.0202294 ] pot_est: [-0.00333107  0.0016309 ]
pottest: [ 0.01113273] pot_est: [ 0.0056152]
pottest: [ 0.01176185] pot_est: [ 0.00165294]
pottest: [ 0.00767375] pot_est: [ 0.0024473]
l= 0.00330578512397 , err= 0.00944655639861
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.00120741  0.00097474]
pottest: [ 0.00903358  0.0202294 ] pot_est: [-0.00352938  0.00163406]
pottest: [ 0.01113273] pot_est: [ 0.00575089]
pottest: [ 0.01176185] pot_est: [ 0.00165261]
pottest: [ 0.00767375] pot_est: [ 0.00249191]
l= 0.0031887755102 , err= 0.00944327674231
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.00119747  0.0009157 ]
pottest: [ 0.00903358  0.0202294 ] pot_est: [-0.00373184  0.00163683]
pottest: [ 0.01113273] pot_est: [ 0.00588859]
pottest: [ 0.01176185] pot_est: [ 0.00165189]
pottest: [ 0.00767375] pot_est: [ 0.00253707]
l= 0.00307787011388 , err= 0.00944071234059
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.001187  0.000855]
pottest: [ 0.00903358  0.0202294 ] pot_est: [-0.00393844  0.00163924]
pottest: [ 0.01113273] pot_est: [ 0.0060283]
pottest: [ 0.01176185] pot_est: [ 0.0016508]
pottest: [ 0.00767375] pot_est: [ 0.00258277]
l= 0.00297265160523 , err= 0.00943885291235
pottest: [ 0.00296697  0.00468031] pot_est: [ 0.00117601  0.00079264]
pottest: [ 0.00903358  0.0202294 ] pot_est: [-0.00414915  0.00164131]
pottest: [ 0.01113273] pot_est: [ 0.00617005]
pottest: [ 0.01176185] pot_est: [ 0.00164935]
pottest: [ 0.00767375] pot_est: [ 0.00262904]
l= 0.00287273771905 , err= 0.00943768937015
lambda = 0.00287273771905

In [85]:

estimation_table = interp_pot
#lambd=0.0

k_inv = linalg.inv(k_pot + lambd * identity(k_pot.shape[0]))
imshow(k_inv, interpolation="none")
colorbar()

Out[85]:

<matplotlib.colorbar.Colorbar instance at 0xfd5c518>

In [86]:

def estimate_pots():
    estimation_table = interp_pot
    
    k_inv = linalg.inv(k_pot + lambd *identity(k_pot.shape[0]))
    beta = dot(k_inv, meas_pot)                 
    output = dot(estimation_table, beta)
    return output

output = estimate_pots()
fig, ax = plt.subplots()


output_area = linspace(min(estimation_area), max(estimation_area), len(output))
ax.plot(x, pots, label = "true pot")
ax.plot(elec_pos, meas_pot, label = "interpolated measurements")
ax.plot(output_area, output*max(pots)/max(output), label = "kCSD estimated pot")
#ax.set_ylim([min(pots) ,max(pots)])
ax.legend(loc=(0.01,0.7))

Out[86]:

<matplotlib.legend.Legend at 0xffa9450>

In [87]:

estimation_table = interp_cross

k_inv = linalg.inv(k_pot + lambd *identity(k_pot.shape[0]))
imshow(k_inv, interpolation="none")
colorbar()

Out[87]:

<matplotlib.colorbar.Colorbar instance at 0xfd557e8>

In [88]:

def estimate_csd():
    estimation_table = interp_cross

    k_inv = linalg.inv(k_pot + lambd *identity(k_pot.shape[0]))
    beta = dot(k_inv, meas_pot)                 
    csd = dot(estimation_table, beta)
    return csd

output = estimate_csd()

plot(linspace(min(estimation_area), max(estimation_area), len(output)), output/np.linalg.norm(output), label = "kCSD estimated CSD") #csd calculated with kCSD
plot(x, true_csd/np.linalg.norm(true_csd), label = "true CSD")


legend(loc = (0.01, 0.75))

Out[88]:

<matplotlib.legend.Legend at 0xfec75d0>