diff --git a/CSD.py b/CSD.py deleted file mode 100644 index 97bb553c..00000000 --- a/CSD.py +++ /dev/null @@ -1,112 +0,0 @@ -""" -This script is used to generate Current Source Density Estimates, -And serves as base class for the kCSD method Jan et.al (2012). - -This was written by : -Chaitanya Chintaluri, -Laboratory of Neuroinformatics, -Nencki Institute of Exprimental Biology, Warsaw. -""" -import numpy as np - -import utility_functions as utils - -class CSD(object): - """CSD - The base class for CSD methods.""" - def __init__(self, ele_pos, pots): - """Initialize CSD Class. - - Parameters - ---------- - ele_pos : numpy array - positions of electrodes - pots : numpy array - potentials measured by electrodes - - Returns - ------- - None - """ - self.validate(ele_pos, pots) - self.ele_pos = ele_pos - self.pots = pots - self.n_ele = self.ele_pos.shape[0] - self.n_time = self.pots.shape[1] - self.dim = self.ele_pos.shape[1] - self.cv_error = None - return - - def validate(self, ele_pos, pots): - """Basic checks to see if inputs are okay - - Parameters - ---------- - ele_pos : numpy array - positions of electrodes - pots : numpy array - potentials measured by electrodes - - Returns - ------- - None - """ - if ele_pos.shape[0] != pots.shape[0]: - raise Exception("Number of measured potentials is not equal " - "to electrode number!") - if ele_pos.shape[0] < 1+ele_pos.shape[1]: #Dim+1 - raise Exception("Number of electrodes must be at least :", - 1+ele_pos.shape[1]) - if utils.check_for_duplicated_electrodes(ele_pos) is False: - raise Exception("Error! Duplicated electrode!") - return - - def method(self): - """Place holder for the actual method that computes the CSD. - - Parameters - ---------- - None - - Returns - ------- - None - """ - pass - - def values(self): - """Place holder for obtaining CSD at the pos_csd locations, it uses the method - function to obtain the CSD. - - Parameters - ---------- - None - - Returns - ------- - None - """ - pass - - def sanity(self, true_csd, pos_csd): - """Useful for comparing TrueCSD with reconstructed CSD. Computes, the RMS error - between the true_csd and the reconstructed csd at pos_csd using the - method defined. - - Parameters - ---------- - - true_csd : csd values used to generate potentials - pos_csd : csd estimatation from the method - - Returns - ------- - RMSE : root mean squared difference - """ - csd = self.values(pos_csd) - RMSE = np.sqrt(np.mean(np.square(true_csd - csd))) - return RMSE - -if __name__ == '__main__': - ele_pos = np.array([[-0.2, -0.2], [0, 0], [0, 1], [1, 0], [1,1], [0.5, 0.5], [1.2, 1.2]]) - pots = np.array([[-1], [-1], [-1], [0], [0], [1], [-1.5]]) - test_class = CSD(ele_pos, pots) diff --git a/KCSD.py b/KCSD.py index f29baf37..8a12bb77 100644 --- a/KCSD.py +++ b/KCSD.py @@ -1,18 +1,74 @@ -""" -This script is used to generate Current Source Density Estimates, -using the kCSD method Jan et.al (2012). +"""This script is used to generate Current Source Density Estimates, using the +kCSD method Jan et.al (2012). This was written by : -Chaitanya Chintaluri, +[1]Chaitanya Chintaluri, +[2]Michal Czerwinski, Laboratory of Neuroinformatics, Nencki Institute of Exprimental Biology, Warsaw. +KCSD1D[1][2], KCSD2D[1], KCSD3D[1], MoIKCSD[1] + """ +from __future__ import division + import numpy as np -from scipy import integrate, interpolate +from scipy import special, integrate, interpolate from scipy.spatial import distance from numpy.linalg import LinAlgError -from CSD import CSD +skmonaco_available = False + +from . import utility_functions as utils +from . import basis_functions as basis + +class CSD(object): + """CSD - The base class for CSD methods.""" + def __init__(self, ele_pos, pots): + self.validate(ele_pos, pots) + self.ele_pos = ele_pos + self.pots = pots + self.n_ele = self.ele_pos.shape[0] + self.n_time = self.pots.shape[1] + self.dim = self.ele_pos.shape[1] + self.cv_error = None + + def validate(self, ele_pos, pots): + """Basic checks to see if inputs are okay + + Parameters + ---------- + ele_pos : numpy array + positions of electrodes + pots : numpy array + potentials measured by electrodes + """ + if ele_pos.shape[0] != pots.shape[0]: + raise Exception("Number of measured potentials is not equal " + "to electrode number!") + if ele_pos.shape[0] < 1+ele_pos.shape[1]: #Dim+1 + raise Exception("Number of electrodes must be at least :", + 1+ele_pos.shape[1]) + if utils.check_for_duplicated_electrodes(ele_pos) is False: + raise Exception("Error! Duplicated electrode!") + + def sanity(self, true_csd, pos_csd): + """Useful for comparing TrueCSD with reconstructed CSD. Computes, the RMS error + between the true_csd and the reconstructed csd at pos_csd using the + method defined. + + Parameters + ---------- + + true_csd : csd values used to generate potentials + pos_csd : csd estimatation from the method + + Returns + ------- + RMSE : root mean squared difference + """ + csd = self.values(pos_csd) + RMSE = np.sqrt(np.mean(np.square(true_csd - csd))) + return RMSE class KCSD(CSD): """KCSD - The base class for all the KCSD variants. @@ -29,7 +85,6 @@ def __init__(self, ele_pos, pots, **kwargs): self.place_basis() self.create_src_dist_tables() self.method() - return def parameters(self, **kwargs): """Defining the default values of the method passed as kwargs @@ -37,41 +92,29 @@ def parameters(self, **kwargs): ---------- **kwargs Same as those passed to initialize the Class - - Returns - ------- - None """ - self.src_type = kwargs.get('src_type', 'gauss') - self.sigma = kwargs.get('sigma', 1.0) - self.h = kwargs.get('h', 1.0) - self.n_src_init = kwargs.get('n_src_init', 1000) - self.lambd = kwargs.get('lambd', 0.0) - self.R_init = kwargs.get('R_init', 0.23) - self.ext_x = kwargs.get('ext_x', 0.0) - self.xmin = kwargs.get('xmin', np.min(self.ele_pos[:, 0])) - self.xmax = kwargs.get('xmax', np.max(self.ele_pos[:, 0])) - self.gdx = kwargs.get('gdx', 0.01*(self.xmax - self.xmin)) + self.src_type = kwargs.pop('src_type', 'gauss') + self.sigma = kwargs.pop('sigma', 1.0) + self.h = kwargs.pop('h', 1.0) + self.n_src_init = kwargs.pop('n_src_init', 1000) + self.lambd = kwargs.pop('lambd', 0.0) + self.R_init = kwargs.pop('R_init', 0.23) + self.ext_x = kwargs.pop('ext_x', 0.0) + self.xmin = kwargs.pop('xmin', np.min(self.ele_pos[:, 0])) + self.xmax = kwargs.pop('xmax', np.max(self.ele_pos[:, 0])) + self.gdx = kwargs.pop('gdx', 0.01*(self.xmax - self.xmin)) if self.dim >= 2: - self.ext_y = kwargs.get('ext_y', 0.0) - self.ymin = kwargs.get('ymin', np.min(self.ele_pos[:, 1])) - self.ymax = kwargs.get('ymax', np.max(self.ele_pos[:, 1])) - self.gdy = kwargs.get('gdy', 0.01*(self.ymax - self.ymin)) + self.ext_y = kwargs.pop('ext_y', 0.0) + self.ymin = kwargs.pop('ymin', np.min(self.ele_pos[:, 1])) + self.ymax = kwargs.pop('ymax', np.max(self.ele_pos[:, 1])) + self.gdy = kwargs.pop('gdy', 0.01*(self.ymax - self.ymin)) if self.dim == 3: - self.ext_z = kwargs.get('ext_z', 0.0) - self.zmin = kwargs.get('zmin', np.min(self.ele_pos[:, 2])) - self.zmax = kwargs.get('zmax', np.max(self.ele_pos[:, 2])) - self.gdz = kwargs.get('gdz', 0.01*(self.zmax - self.zmin)) - return - - def estimate_at(self): - pass - - def place_basis(self): - pass - - def create_src_dist_tables(self): - pass + self.ext_z = kwargs.pop('ext_z', 0.0) + self.zmin = kwargs.pop('zmin', np.min(self.ele_pos[:, 2])) + self.zmax = kwargs.pop('zmax', np.max(self.ele_pos[:, 2])) + self.gdz = kwargs.pop('gdz', 0.01*(self.zmax - self.zmin)) + if kwargs: + raise TypeError('Invalid keyword arguments:', kwargs.keys()) def method(self): """Actual sequence of methods called for KCSD @@ -81,32 +124,23 @@ def method(self): Parameters ---------- None - - Returns - ------- - None """ self.create_lookup() #Look up table self.update_b_pot() #update kernel self.update_b_src() #update crskernel self.update_b_interp_pot() #update pot interp - return def create_lookup(self, dist_table_density=20): """Creates a table for easy potential estimation from CSD. Updates and Returns the potentials due to a given basis source like a lookup - table whose shape=(dist_table_density,)--> set in KCSD2D_Helpers.py + table whose shape=(dist_table_density,) Parameters ---------- dist_table_density : int number of distance values at which potentials are computed. Default 100 - - Returns - ------- - None """ xs = np.logspace(0., np.log10(self.dist_max+1.), dist_table_density) xs = xs - 1.0 #starting from 0 @@ -118,7 +152,6 @@ def create_lookup(self, dist_table_density=20): self.sigma, self.basis) self.interpolate_pot_at = interpolate.interp1d(xs, dist_table, kind='cubic') - return def update_b_pot(self): """Updates the b_pot - array is (#_basis_sources, #_electrodes) @@ -131,15 +164,10 @@ def update_b_pot(self): Parameters ---------- None - - Returns - ------- - None """ self.b_pot = self.interpolate_pot_at(self.src_ele_dists) self.k_pot = np.dot(self.b_pot.T, self.b_pot) #K(x,x') Eq9,Jan2012 self.k_pot /= self.n_src - return def update_b_src(self): """Updates the b_src in the shape of (#_est_pts, #_basis_sources) @@ -151,15 +179,10 @@ def update_b_src(self): Parameters ---------- None - - Returns - ------- - None """ self.b_src = self.basis(self.src_estm_dists, self.R).T self.k_interp_cross = np.dot(self.b_src, self.b_pot) #K_t(x,y) Eq17 self.k_interp_cross /= self.n_src - return def update_b_interp_pot(self): """Compute the matrix of potentials generated by every source @@ -170,15 +193,10 @@ def update_b_interp_pot(self): Parameters ---------- None - - Returns - ------- - None """ self.b_interp_pot = self.interpolate_pot_at(self.src_estm_dists).T self.k_interp_pot = np.dot(self.b_interp_pot, self.b_pot) self.k_interp_pot /= self.n_src - return def values(self, estimate='CSD'): """Computes the values of the quantity of interest @@ -191,24 +209,37 @@ def values(self, estimate='CSD'): Returns ------- - estimated quantity of shape (ngx, ngy, ngz, nt) + estimation : np.array + estimated quantity of shape (ngx, ngy, ngz, nt) """ if estimate == 'CSD': #Maybe used for estimating the potentials also. estimation_table = self.k_interp_cross elif estimate == 'POT': estimation_table = self.k_interp_pot else: - print 'Invalid quantity to be measured, pass either CSD or POT' + print('Invalid quantity to be measured, pass either CSD or POT') k_inv = np.linalg.inv(self.k_pot + self.lambd * np.identity(self.k_pot.shape[0])) estimation = np.zeros((self.n_estm, self.n_time)) - for t in xrange(self.n_time): + for t in range(self.n_time): beta = np.dot(k_inv, self.pots[:, t]) - for i in xrange(self.n_ele): + for i in range(self.n_ele): estimation[:, t] += estimation_table[:, i] *beta[i] # C*(x) Eq 18 return self.process_estimate(estimation) def process_estimate(self, estimation): + """Function used to rearrange estimation according to dimension, to be + used by the fuctions values + + Parameters + ---------- + estimation : np.array + + Returns + ------- + estimation : np.array + estimated quantity of shape (ngx, ngy, ngz, nt) + """ if self.dim == 1: estimation = estimation.reshape(self.ngx, self.n_time) elif self.dim == 2: @@ -218,35 +249,25 @@ def process_estimate(self, estimation): return estimation def update_R(self, R): - """Used in Cross validation + """Update the width of the basis fuction - Used in Cross validation Parameters ---------- R : float - - Returns - ------- - None """ self.R = R self.dist_max = max(np.max(self.src_ele_dists), np.max(self.src_estm_dists)) + self.R self.method() - return def update_lambda(self, lambd): - """Used in Cross validation + """Update the lambda parameter of regularization, Used in Cross validation Parameters ---------- lambd : float - - Returns - ------- - None """ self.lambd = lambd - return def cross_validate(self, lambdas=None, Rs=None): """Method defines the cross validation. @@ -267,7 +288,7 @@ def cross_validate(self, lambdas=None, Rs=None): Lambda : post cross validation """ if lambdas is None: #when None - print 'No lambda given, using defaults' + print('No lambda given, using defaults') lambdas = np.logspace(-2,-25,25,base=10.) #Default multiple lambda lambdas = np.hstack((lambdas, np.array((0.0)))) elif lambdas.size == 1: #resize when one entry @@ -278,12 +299,12 @@ def cross_validate(self, lambdas=None, Rs=None): index_generator = [] for ii in range(self.n_ele): idx_test = [ii] - idx_train = range(self.n_ele) + idx_train = list(range(self.n_ele)) idx_train.remove(ii) #Leave one out index_generator.append((idx_train, idx_test)) for R_idx,R in enumerate(Rs): #Iterate over R self.update_R(R) - print 'Cross validating R (all lambda) :', R + print('Cross validating R (all lambda) :', R) for lambd_idx,lambd in enumerate(lambdas): #Iterate over lambdas errs[R_idx, lambd_idx] = self.compute_cverror(lambd, index_generator) @@ -293,7 +314,7 @@ def cross_validate(self, lambdas=None, Rs=None): self.cv_error = np.min(errs) #otherwise is None self.update_R(cv_R) #Update solver self.update_lambda(cv_lambda) - print 'R, lambda :', cv_R, cv_lambda + print('R, lambda :', cv_R, cv_lambda) return cv_R, cv_lambda def compute_cverror(self, lambd, index_generator): @@ -325,10 +346,769 @@ def compute_cverror(self, lambd, index_generator): V_est[:, tt] += beta_new[ii, tt] * B_test[:, ii] err += np.linalg.norm(V_est-V_test) except LinAlgError: - print 'Encoutered Singular Matrix Error: try changing ele_pos' - #err = 10000. #singluar matrix errors! - raise + raise LinAlgError('Encoutered Singular Matrix Error: try changing ele_pos slightly') return err +class KCSD1D(KCSD): + """KCSD1D - The 1D variant for the Kernel Current Source Density method. + + This estimates the Current Source Density, for a given configuration of + electrod positions and recorded potentials, in the case of 1D recording + electrodes (laminar probes). The method implented here is based on the + original paper by Jan Potworowski et.al. 2012. + """ + def __init__(self, ele_pos, pots, **kwargs): + """Initialize KCSD1D Class. + + Parameters + ---------- + ele_pos : numpy array + positions of electrodes + pots : numpy array + potentials measured by electrodes + **kwargs + configuration parameters, that may contain the following keys: + src_type : str + basis function type ('gauss', 'step', 'gauss_lim') + Defaults to 'gauss' + sigma : float + space conductance of the tissue in S/m + Defaults to 1 S/m + n_src_init : int + requested number of sources + Defaults to 300 + R_init : float + demanded thickness of the basis element + Defaults to 0.23 + h : float + thickness of analyzed cylindrical slice + Defaults to 1. + xmin, xmax : floats + boundaries for CSD estimation space + Defaults to min(ele_pos(x)), and max(ele_pos(x)) + ext_x : float + length of space extension: x_min-ext_x ... x_max+ext_x + Defaults to 0. + gdx : float + space increments in the estimation space + Defaults to 0.01(xmax-xmin) + lambd : float + regularization parameter for ridge regression + Defaults to 0. + + Raises + ------ + LinAlgException + If the matrix is not numerically invertible. + KeyError + Basis function (src_type) not implemented. See basis_functions.py for available + """ + super(KCSD1D, self).__init__(ele_pos, pots, **kwargs) + + def estimate_at(self): + """Defines locations where the estimation is wanted + Defines: + self.n_estm = self.estm_x.size + self.ngx = self.estm_x.shape + self.estm_x : Locations at which CSD is requested. + + Parameters + ---------- + None + """ + nx = (self.xmax - self.xmin)/self.gdx + self.estm_x = np.mgrid[self.xmin:self.xmax:np.complex(0,nx)] + self.n_estm = self.estm_x.size + self.ngx = self.estm_x.shape[0] + + def place_basis(self): + """Places basis sources of the defined type. + Checks if a given source_type is defined, if so then defines it + self.basis, This function gives locations of the basis sources, + Defines + source_type : basis_fuctions.basis_1D.keys() + self.R based on R_init + self.dist_max as maximum distance between electrode and basis + self.nsx = self.src_x.shape + self.src_x : Locations at which basis sources are placed. + + Parameters + ---------- + None + """ + source_type = self.src_type + try: + self.basis = basis.basis_1D[source_type] + except KeyError: + raise KeyError('Invalid source_type for basis! available are:', + basis.basis_1D.keys()) + (self.src_x, self.R) = utils.distribute_srcs_1D(self.estm_x, + self.n_src_init, + self.ext_x, + self.R_init ) + self.n_src = self.src_x.size + self.nsx = self.src_x.shape + + def create_src_dist_tables(self): + """Creates distance tables between sources, electrode and estm points + + Parameters + ---------- + None + """ + src_loc = np.array((self.src_x.ravel())) + src_loc = src_loc.reshape((len(src_loc), 1)) + est_loc = np.array((self.estm_x.ravel())) + est_loc = est_loc.reshape((len(est_loc), 1)) + self.src_ele_dists = distance.cdist(src_loc, self.ele_pos, 'euclidean') + self.src_estm_dists = distance.cdist(src_loc, est_loc, 'euclidean') + self.dist_max = max(np.max(self.src_ele_dists), np.max(self.src_estm_dists)) + self.R + + def forward_model(self, x, R, h, sigma, src_type): + """FWD model functions + Evaluates potential at point (x,0) by a basis source located at (0,0) + Eq 26 kCSD by Jan,2012 + + Parameters + ---------- + x : float + R : float + h : float + sigma : float + src_type : basis_1D.key + + Returns + ------- + pot : float + value of potential at specified distance from the source + """ + pot, err = integrate.quad(self.int_pot_1D, + -R, R, + args=(x, R, h, src_type)) + pot *= 1./(2.0*sigma) + return pot + + def int_pot_1D(self, xp, x, R, h, basis_func): + """FWD model function. + Returns contribution of a point xp,yp, belonging to a basis source + support centered at (0,0) to the potential measured at (x,0), + integrated over xp,yp gives the potential generated by a + basis source element centered at (0,0) at point (x,0) + Eq 26 kCSD by Jan,2012 + + Parameters + ---------- + xp : floats or np.arrays + point or set of points where function should be calculated + x : float + position at which potential is being measured + R : float + The size of the basis function + h : float + thickness of slice + basis_func : method + Fuction of the basis source + + Returns + ------- + pot : float + """ + m = np.sqrt((x-xp)**2 + h**2) - abs(x-xp) + m *= basis_func(abs(xp), R) #xp is the distance + return m + +class KCSD2D(KCSD): + """KCSD2D - The 2D variant for the Kernel Current Source Density method. + + This estimates the Current Source Density, for a given configuration of + electrod positions and recorded potentials, in the case of 2D recording + electrodes. The method implented here is based on the original paper + by Jan Potworowski et.al. 2012. + """ + def __init__(self, ele_pos, pots, **kwargs): + """Initialize KCSD2D Class. + + Parameters + ---------- + ele_pos : numpy array + positions of electrodes + pots : numpy array + potentials measured by electrodes + **kwargs + configuration parameters, that may contain the following keys: + src_type : str + basis function type ('gauss', 'step', 'gauss_lim') + Defaults to 'gauss' + sigma : float + space conductance of the tissue in S/m + Defaults to 1 S/m + n_src_init : int + requested number of sources + Defaults to 1000 + R_init : float + demanded thickness of the basis element + Defaults to 0.23 + h : float + thickness of analyzed tissue slice + Defaults to 1. + xmin, xmax, ymin, ymax : floats + boundaries for CSD estimation space + Defaults to min(ele_pos(x)), and max(ele_pos(x)) + Defaults to min(ele_pos(y)), and max(ele_pos(y)) + ext_x, ext_y : float + length of space extension: x_min-ext_x ... x_max+ext_x + length of space extension: y_min-ext_y ... y_max+ext_y + Defaults to 0. + gdx, gdy : float + space increments in the estimation space + Defaults to 0.01(xmax-xmin) + Defaults to 0.01(ymax-ymin) + lambd : float + regularization parameter for ridge regression + Defaults to 0. + + Raises + ------ + LinAlgError + Could not invert the matrix, try changing the ele_pos slightly + KeyError + Basis function (src_type) not implemented. See basis_functions.py for available + """ + super(KCSD2D, self).__init__(ele_pos, pots, **kwargs) + + def estimate_at(self): + """Defines locations where the estimation is wanted + Defines: + self.n_estm = self.estm_x.size + self.ngx, self.ngy = self.estm_x.shape + self.estm_x, self.estm_y : Locations at which CSD is requested. + + Parameters + ---------- + None + """ + nx = (self.xmax - self.xmin)/self.gdx + ny = (self.ymax - self.ymin)/self.gdy + self.estm_x, self.estm_y = np.mgrid[self.xmin:self.xmax:np.complex(0,nx), + self.ymin:self.ymax:np.complex(0,ny)] + self.n_estm = self.estm_x.size + self.ngx, self.ngy = self.estm_x.shape + + def place_basis(self): + """Places basis sources of the defined type. + Checks if a given source_type is defined, if so then defines it + self.basis, This function gives locations of the basis sources, + Defines + source_type : basis_fuctions.basis_2D.keys() + self.R based on R_init + self.dist_max as maximum distance between electrode and basis + self.nsx, self.nsy = self.src_x.shape + self.src_x, self.src_y : Locations at which basis sources are placed. + + Parameters + ---------- + None + """ + source_type = self.src_type + try: + self.basis = basis.basis_2D[source_type] + except KeyError: + raise KeyError('Invalid source_type for basis! available are:', + basis.basis_2D.keys()) + (self.src_x, self.src_y, self.R) = utils.distribute_srcs_2D(self.estm_x, + self.estm_y, + self.n_src_init, + self.ext_x, + self.ext_y, + self.R_init ) + self.n_src = self.src_x.size + self.nsx, self.nsy = self.src_x.shape + + def create_src_dist_tables(self): + """Creates distance tables between sources, electrode and estm points + + Parameters + ---------- + None + """ + src_loc = np.array((self.src_x.ravel(), self.src_y.ravel())) + est_loc = np.array((self.estm_x.ravel(), self.estm_y.ravel())) + self.src_ele_dists = distance.cdist(src_loc.T, self.ele_pos, 'euclidean') + self.src_estm_dists = distance.cdist(src_loc.T, est_loc.T, 'euclidean') + self.dist_max = max(np.max(self.src_ele_dists), np.max(self.src_estm_dists)) + self.R + + def forward_model(self, x, R, h, sigma, src_type): + """FWD model functions + Evaluates potential at point (x,0) by a basis source located at (0,0) + Eq 22 kCSD by Jan,2012 + + Parameters + ---------- + x : float + R : float + h : float + sigma : float + src_type : basis_2D.key + + Returns + ------- + pot : float + value of potential at specified distance from the source + """ + pot, err = integrate.dblquad(self.int_pot_2D, + -R, R, + lambda x: -R, + lambda x: R, + args=(x, R, h, src_type)) + pot *= 1./(2.0*np.pi*sigma) #Potential basis functions bi_x_y + return pot + + def int_pot_2D(self, xp, yp, x, R, h, basis_func): + """FWD model function. + Returns contribution of a point xp,yp, belonging to a basis source + support centered at (0,0) to the potential measured at (x,0), + integrated over xp,yp gives the potential generated by a + basis source element centered at (0,0) at point (x,0) + + Parameters + ---------- + xp, yp : floats or np.arrays + point or set of points where function should be calculated + x : float + position at which potential is being measured + R : float + The size of the basis function + h : float + thickness of slice + basis_func : method + Fuction of the basis source + + Returns + ------- + pot : float + """ + y = ((x-xp)**2 + yp**2)**(0.5) + if y < 0.00001: + y = 0.00001 + dist = np.sqrt(xp**2 + yp**2) + pot = np.arcsinh(h/y)*basis_func(dist, R) + return pot + +class MoIKCSD(KCSD2D): + """MoIKCSD - CSD while including the forward modeling effects of saline. + + This estimates the Current Source Density, for a given configuration of + electrod positions and recorded potentials, in the case of 2D recording + electrodes from an MEA electrode plane using the Method of Images. + The method implented here is based on kCSD method by Jan Potworowski + et.al. 2012, which was extended in Ness, Chintaluri 2015 for MEA. + """ + def __init__(self, ele_pos, pots, **kwargs): + """Initialize MoIKCSD Class. + + Parameters + ---------- + ele_pos : numpy array + positions of electrodes + pots : numpy array + potentials measured by electrodes + **kwargs + configuration parameters, that may contain the following keys: + src_type : str + basis function type ('gauss', 'step', 'gauss_lim') + Defaults to 'gauss' + sigma : float + space conductance of the tissue in S/m + Defaults to 1 S/m + sigma_S : float + conductance of the saline (medium) in S/m + Default is 5 S/m (5 times more conductive) + n_src_init : int + requested number of sources + Defaults to 1000 + R_init : float + demanded thickness of the basis element + Defaults to 0.23 + h : float + thickness of analyzed tissue slice + Defaults to 1. + xmin, xmax, ymin, ymax : floats + boundaries for CSD estimation space + Defaults to min(ele_pos(x)), and max(ele_pos(x)) + Defaults to min(ele_pos(y)), and max(ele_pos(y)) + ext_x, ext_y : float + length of space extension: x_min-ext_x ... x_max+ext_x + length of space extension: y_min-ext_y ... y_max+ext_y + Defaults to 0. + gdx, gdy : float + space increments in the estimation space + Defaults to 0.01(xmax-xmin) + Defaults to 0.01(ymax-ymin) + lambd : float + regularization parameter for ridge regression + Defaults to 0. + MoI_iters : int + Number of interations in method of images. + Default is 20 + """ + self.MoI_iters = kwargs.pop('MoI_iters', 20) + self.sigma_S = kwargs.pop('sigma_S', 5.0) + self.sigma = kwargs.pop('sigma', 1.0) + W_TS = (self.sigma - self.sigma_S) / (self.sigma + self.sigma_S) + self.iters = np.arange(self.MoI_iters) + 1 #Eq 6, Ness (2015) + self.iter_factor = W_TS**self.iters + super(MoIKCSD, self).__init__(ele_pos, pots, **kwargs) + + def forward_model(self, x, R, h, sigma, src_type): + """FWD model functions + Evaluates potential at point (x,0) by a basis source located at (0,0) + Eq 22 kCSD by Jan,2012 + + Parameters + ---------- + x : float + R : float + h : float + sigma : float + src_type : basis_2D.key + + Returns + ------- + pot : float + value of potential at specified distance from the source + """ + pot, err = integrate.dblquad(self.int_pot_2D_moi, -R, R, + lambda x: -R, + lambda x: R, + args=(x, R, h, src_type)) + pot *= 1./(2.0*np.pi*sigma) + return pot + + def int_pot_2D_moi(self, xp, yp, x, R, h, basis_func): + """FWD model function. Incorporates the Method of Images. + Returns contribution of a point xp,yp, belonging to a basis source + support centered at (0,0) to the potential measured at (x,0), + integrated over xp,yp gives the potential generated by a + basis source element centered at (0,0) at point (x,0) + #Eq 20, Ness(2015) + + Parameters + ---------- + xp, yp : floats or np.arrays + point or set of points where function should be calculated + x : float + position at which potential is being measured + R : float + The size of the basis function + h : float + thickness of slice + basis_func : method + Fuction of the basis source + + Returns + ------- + pot : float + """ + L = ((x-xp)**2 + yp**2)**(0.5) + if L < 0.00001: + L = 0.00001 + correction = np.arcsinh((h-(2*h*self.iters))/L) + np.arcsinh((h+(2*h*self.iters))/L) + pot = np.arcsinh(h/L) + np.sum(self.iter_factor*correction) + dist = np.sqrt(xp**2 + yp**2) + pot *= basis_func(dist, R) #Eq 20, Ness et.al. + return pot + +class KCSD3D(KCSD): + """KCSD3D - The 3D variant for the Kernel Current Source Density method. + + This estimates the Current Source Density, for a given configuration of + electrod positions and recorded potentials, in the case of 2D recording + electrodes. The method implented here is based on the original paper + by Jan Potworowski et.al. 2012. + """ + def __init__(self, ele_pos, pots, **kwargs): + """Initialize KCSD3D Class. + + Parameters + ---------- + ele_pos : numpy array + positions of electrodes + pots : numpy array + potentials measured by electrodes + **kwargs + configuration parameters, that may contain the following keys: + src_type : str + basis function type ('gauss', 'step', 'gauss_lim') + Defaults to 'gauss' + sigma : float + space conductance of the tissue in S/m + Defaults to 1 S/m + n_src_init : int + requested number of sources + Defaults to 1000 + R_init : float + demanded thickness of the basis element + Defaults to 0.23 + h : float + thickness of analyzed tissue slice + Defaults to 1. + xmin, xmax, ymin, ymax, zmin, zmax : floats + boundaries for CSD estimation space + Defaults to min(ele_pos(x)), and max(ele_pos(x)) + Defaults to min(ele_pos(y)), and max(ele_pos(y)) + Defaults to min(ele_pos(z)), and max(ele_pos(z)) + ext_x, ext_y, ext_z : float + length of space extension: xmin-ext_x ... xmax+ext_x + length of space extension: ymin-ext_y ... ymax+ext_y + length of space extension: zmin-ext_z ... zmax+ext_z + Defaults to 0. + gdx, gdy, gdz : float + space increments in the estimation space + Defaults to 0.01(xmax-xmin) + Defaults to 0.01(ymax-ymin) + Defaults to 0.01(zmax-zmin) + lambd : float + regularization parameter for ridge regression + Defaults to 0. + + Raises + ------ + LinAlgError + Could not invert the matrix, try changing the ele_pos slightly + KeyError + Basis function (src_type) not implemented. See basis_functions.py for available + """ + super(KCSD3D, self).__init__(ele_pos, pots, **kwargs) + + def estimate_at(self): + """Defines locations where the estimation is wanted + Defines: + self.n_estm = self.estm_x.size + self.ngx, self.ngy, self.ngz = self.estm_x.shape + self.estm_x, self.estm_y, self.estm_z : Pts. at which CSD is requested + + Parameters + ---------- + None + """ + nx = (self.xmax - self.xmin)/self.gdx + ny = (self.ymax - self.ymin)/self.gdy + nz = (self.zmax - self.zmin)/self.gdz + self.estm_x, self.estm_y, self.estm_z = np.mgrid[self.xmin:self.xmax:np.complex(0,nx), + self.ymin:self.ymax:np.complex(0,ny), + self.zmin:self.zmax:np.complex(0,nz)] + self.n_estm = self.estm_x.size + self.ngx, self.ngy, self.ngz = self.estm_x.shape + + def place_basis(self): + """Places basis sources of the defined type. + Checks if a given source_type is defined, if so then defines it + self.basis, This function gives locations of the basis sources, + Defines + source_type : basis_fuctions.basis_2D.keys() + self.R based on R_init + self.dist_max as maximum distance between electrode and basis + self.nsx, self.nsy, self.nsz = self.src_x.shape + self.src_x, self.src_y, self.src_z : Locations at which basis sources are placed. + + Parameters + ---------- + None + """ + source_type = self.src_type + try: + self.basis = basis.basis_3D[source_type] + except KeyError: + raise KeyError('Invalid source_type for basis! available are:', + basis.basis_3D.keys()) + (self.src_x, self.src_y, self.src_z, self.R) = utils.distribute_srcs_3D(self.estm_x, + self.estm_y, + self.estm_z, + self.n_src_init, + self.ext_x, + self.ext_y, + self.ext_z, + self.R_init) + + self.n_src = self.src_x.size + self.nsx, self.nsy, self.nsz = self.src_x.shape + + def create_src_dist_tables(self): + """Creates distance tables between sources, electrode and estm points + + Parameters + ---------- + None + """ + src_loc = np.array((self.src_x.ravel(), + self.src_y.ravel(), + self.src_z.ravel())) + est_loc = np.array((self.estm_x.ravel(), + self.estm_y.ravel(), + self.estm_z.ravel())) + self.src_ele_dists = distance.cdist(src_loc.T, self.ele_pos, 'euclidean') + self.src_estm_dists = distance.cdist(src_loc.T, est_loc.T, 'euclidean') + self.dist_max = max(np.max(self.src_ele_dists), np.max(self.src_estm_dists)) + self.R + + def forward_model(self, x, R, h, sigma, src_type): + """FWD model functions + Evaluates potential at point (x,0) by a basis source located at (0,0) + Utlizies sk monaco monte carlo method if available, otherwise defaults + to scipy integrate + + Parameters + ---------- + x : float + R : float + h : float + sigma : float + src_type : basis_3D.key + + Returns + ------- + pot : float + value of potential at specified distance from the source + """ + if src_type.__name__ == "gauss_3D": + if x == 0: x=0.0001 + pot = special.erf(x/(np.sqrt(2)*R/3.0)) / x + elif src_type.__name__ == "gauss_lim_3D": + if x == 0: x=0.0001 + d = R/3. + if x < R: + e = np.exp(-(x/ (np.sqrt(2)*d))**2) + erf = special.erf(x / (np.sqrt(2)*d)) + pot = 4* np.pi * ( (d**2)*(e - np.exp(-4.5)) + + (1/x)*((np.sqrt(np.pi/2)*(d**3)*erf) - x*(d**2)*e)) + else: + pot = 15.28828*(d)**3 / x + pot /= (np.sqrt(2*np.pi)*d)**3 + elif src_type.__name__ == "step_3D": + Q = 4.*np.pi*(R**3)/3. + if x < R: + pot = (Q * (3 - (x/R)**2)) / (2.*R) + else: + pot = Q / x + pot *= 3/(4*np.pi*R**3) + else: + if skmonaco_available: + pot, err = mcmiser(self.int_pot_3D_mc, + npoints=1e5, + xl=[-R, -R, -R], + xu=[R, R, R], + seed=42, + nprocs=num_cores, + args=(x, R, h, src_type)) + else: + pot, err = integrate.tplquad(self.int_pot_3D, + -R, + R, + lambda x: -R, + lambda x: R, + lambda x, y: -R, + lambda x, y: R, + args=(x, R, h, src_type)) + pot *= 1./(4.0*np.pi*sigma) + return pot + + def int_pot_3D(self, xp, yp, zp, x, R, h, basis_func): + """FWD model function. + Returns contribution of a point xp,yp, belonging to a basis source + support centered at (0,0) to the potential measured at (x,0), + integrated over xp,yp gives the potential generated by a + basis source element centered at (0,0) at point (x,0) + + Parameters + ---------- + xp, yp, zp : floats or np.arrays + point or set of points where function should be calculated + x : float + position at which potential is being measured + R : float + The size of the basis function + h : float + thickness of slice + basis_func : method + Fuction of the basis source + + Returns + ------- + pot : float + """ + y = ((x-xp)**2 + yp**2 + zp**2)**0.5 + if y < 0.00001: + y = 0.00001 + dist = np.sqrt(xp**2 + yp**2 + zp**2) + pot = 1.0/y + pot *= basis_func(dist, R) + return pot + + def int_pot_3D_mc(self, xyz, x, R, h, basis_func): + """ + The same as int_pot_3D, just different input: x,y,z <-- xyz (tuple) + FWD model function, using Monte Carlo Method of integration + Returns contribution of a point xp,yp, belonging to a basis source + support centered at (0,0) to the potential measured at (x,0), + integrated over xp,yp gives the potential generated by a + basis source element centered at (0,0) at point (x,0) + + Parameters + ---------- + xp, yp, zp : floats or np.arrays + point or set of points where function should be calculated + x : float + position at which potential is being measured + R : float + The size of the basis function + h : float + thickness of slice + basis_func : method + Fuction of the basis source + + Returns + ------- + pot : float + """ + xp, yp, zp = xyz + return self.int_pot_3D(xp, yp, zp, x, R, h, basis_func) + if __name__ == '__main__': - print 'Invalid usage, use this an inheritable class only' + print('Checking 1D') + ele_pos = np.array(([-0.1],[0], [0.5], [1.], [1.4], [2.], [2.3])) + pots = np.array([[-1], [-1], [-1], [0], [0], [1], [-1.5]]) + k = KCSD1D(ele_pos, pots, + gdx=0.01, n_src_init=300, + ext_x=0.0, src_type='gauss') + k.cross_validate() + print(k.values()) + + print('Checking 2D') + ele_pos = np.array([[-0.2, -0.2],[0, 0], [0, 1], [1, 0], [1,1], [0.5, 0.5], + [1.2, 1.2]]) + pots = np.array([[-1], [-1], [-1], [0], [0], [1], [-1.5]]) + k = KCSD2D(ele_pos, pots, + gdx=0.05, gdy=0.05, + xmin=-2.0, xmax=2.0, + ymin=-2.0, ymax=2.0, + src_type='gauss') + k.cross_validate() + print(k.values()) + + print('Checking MoIKCSD') + k = MoIKCSD(ele_pos, pots, + gdx=0.05, gdy=0.05, + xmin=-2.0, xmax=2.0, + ymin=-2.0, ymax= 2.0) + k.cross_validate() + + print('Checking KCSD3D') + ele_pos = np.array([(0, 0, 0), (0, 0, 1), (0, 1, 0), (1, 0, 0), + (0, 1, 1), (1, 1, 0), (1, 0, 1), (1, 1, 1), + (0.5, 0.5, 0.5)]) + pots = np.array([[-0.5], [0], [-0.5], [0], [0], [0.2], [0], [0], [1]]) + k = KCSD3D(ele_pos, pots, + gdx=0.02, gdy=0.02, gdz=0.02, + n_src_init=1000, src_type='gauss_lim') + k.cross_validate() + print(k.values()) + diff --git a/KCSD1D.py b/KCSD1D.py deleted file mode 100644 index 5fe6f3a4..00000000 --- a/KCSD1D.py +++ /dev/null @@ -1,224 +0,0 @@ -""" -This script is used to generate Current Source Density Estimates, -using the kCSD method Jan et.al (2012) for 1D case - -These scripts are based on Grzegorz Parka's, -Google Summer of Code 2014, INFC/pykCSD - -This was written by : -Michal Czerwinski, Chaitanya Chintaluri -Laboratory of Neuroinformatics, -Nencki Institute of Experimental Biology, Warsaw. -""" -import numpy as np -from scipy import integrate, interpolate -from scipy.spatial import distance - -from KCSD import KCSD -import utility_functions as utils -import basis_functions as basis - -class KCSD1D(KCSD): - """KCSD1D - The 1D variant for the Kernel Current Source Density method. - - This estimates the Current Source Density, for a given configuration of - electrod positions and recorded potentials, in the case of 1D recording - electrodes (laminar probes). The method implented here is based on the - original paper by Jan Potworowski et.al. 2012. - """ - def __init__(self, ele_pos, pots, **kwargs): - """Initialize KCSD1D Class. - - Parameters - ---------- - ele_pos : numpy array - positions of electrodes - pots : numpy array - potentials measured by electrodes - **kwargs - configuration parameters, that may contain the following keys: - src_type : str - basis function type ('gauss', 'step', 'gauss_lim') - Defaults to 'gauss' - sigma : float - space conductance of the medium - Defaults to 1. - n_src_init : int - requested number of sources - Defaults to 300 - R_init : float - demanded thickness of the basis element - Defaults to 0.23 - h : float - thickness of analyzed cylindrical slice - Defaults to 1. - xmin, xmax : floats - boundaries for CSD estimation space - Defaults to min(ele_pos(x)), and max(ele_pos(x)) - ext_x : float - length of space extension: x_min-ext_x ... x_max+ext_x - Defaults to 0. - gdx : float - space increments in the estimation space - Defaults to 0.01(xmax-xmin) - lambd : float - regularization parameter for ridge regression - Defaults to 0. - - Returns - ------- - None - - Raises - ------ - LinAlgException - If the matrix is not numerically invertible. - KeyError - Basis function (src_type) not implemented. See basis_functions.py for available - """ - super(KCSD1D, self).__init__(ele_pos, pots, **kwargs) - return - - def estimate_at(self): - """Defines locations where the estimation is wanted - Defines: - self.n_estm = self.estm_x.size - self.ngx = self.estm_x.shape - self.estm_x : Locations at which CSD is requested. - - Parameters - ---------- - None - - Returns - ------- - None - """ - nx = (self.xmax - self.xmin)/self.gdx - self.estm_x = np.mgrid[self.xmin:self.xmax:np.complex(0,nx)] - self.n_estm = self.estm_x.size - self.ngx = self.estm_x.shape[0] - return - - def place_basis(self): - """Places basis sources of the defined type. - Checks if a given source_type is defined, if so then defines it - self.basis, This function gives locations of the basis sources, - Defines - source_type : basis_fuctions.basis_1D.keys() - self.R based on R_init - self.dist_max as maximum distance between electrode and basis - self.nsx = self.src_x.shape - self.src_x : Locations at which basis sources are placed. - - Parameters - ---------- - None - - Returns - ------- - None - """ - # check If valid basis source type passed: - source_type = self.src_type - try: - self.basis = basis.basis_1D[source_type] - except: - print 'Invalid source_type for basis! available are:', basis.basis_1D.keys() - raise KeyError - #Mesh where the source basis are placed is at self.src_x - (self.src_x, self.R) = utils.distribute_srcs_1D(self.estm_x, - self.n_src_init, - self.ext_x, - self.R_init ) - self.n_src = self.src_x.size - self.nsx = self.src_x.shape - return - - def create_src_dist_tables(self): - """Creates distance tables between sources, electrode and estm points - - Parameters - ---------- - None - - Returns - ------- - None - """ - src_loc = np.array((self.src_x.ravel())) - src_loc = src_loc.reshape((len(src_loc), 1)) - est_loc = np.array((self.estm_x.ravel())) - est_loc = est_loc.reshape((len(est_loc), 1)) - self.src_ele_dists = distance.cdist(src_loc, self.ele_pos, 'euclidean') - self.src_estm_dists = distance.cdist(src_loc, est_loc, 'euclidean') - self.dist_max = max(np.max(self.src_ele_dists), np.max(self.src_estm_dists)) + self.R - return - - def forward_model(self, x, R, h, sigma, src_type): - """FWD model functions - Evaluates potential at point (x,0) by a basis source located at (0,0) - Eq 26 kCSD by Jan,2012 - - Parameters - ---------- - x : float - R : float - h : float - sigma : float - src_type : basis_1D.key - - Returns - ------- - pot : float - value of potential at specified distance from the source - """ - pot, err = integrate.quad(self.int_pot_1D, - -R, R, - args=(x, R, h, src_type)) - pot *= 1./(2.0*sigma) - return pot - - def int_pot_1D(self, xp, x, R, h, basis_func): - """FWD model function. - Returns contribution of a point xp,yp, belonging to a basis source - support centered at (0,0) to the potential measured at (x,0), - integrated over xp,yp gives the potential generated by a - basis source element centered at (0,0) at point (x,0) - Eq 26 kCSD by Jan,2012 - - Parameters - ---------- - xp : floats or np.arrays - point or set of points where function should be calculated - x : float - position at which potential is being measured - R : float - The size of the basis function - h : float - thickness of slice - basis_func : method - Fuction of the basis source - - Returns - ------- - pot : float - """ - #1/2sigma normalization not here - m = np.sqrt((x-xp)**2 + h**2) - abs(x-xp) - m *= basis_func(abs(xp), R) #xp is the distance - return m - - -if __name__ == '__main__': - ele_pos = np.array(([-0.1],[0], [0.5], [1.], [1.4], [2.], [2.3])) - pots = np.array([[-1], [-1], [-1], [0], [0], [1], [-1.5]]) - k = KCSD1D(ele_pos, pots, - gdx=0.01, n_src_init=300, - ext_x=0.0, src_type='gauss') - #k.cross_validate(lambdas=np.array((0.0)), Rs=np.array([0.21, 0.23, 0.24])) - #k.cross_validate() - #print k.values() - #k.cross_validate(Rs=np.array(0.14).reshape(1)) - #k.cross_validate(Rs=np.array((0.01,0.02,0.04))) - #print k.values() diff --git a/KCSD2D.py b/KCSD2D.py deleted file mode 100644 index 18a9cf22..00000000 --- a/KCSD2D.py +++ /dev/null @@ -1,224 +0,0 @@ -""" -This script is used to generate Current Source Density Estimates, -using the kCSD method Jan et.al (2012) for 2D case. - -These scripts are based on Grzegorz Parka's, -Google Summer of Code 2014, INFC/pykCSD - -This was written by : -Chaitanya Chintaluri, -Laboratory of Neuroinformatics, -Nencki Institute of Exprimental Biology, Warsaw. -""" -import numpy as np -from scipy import integrate, interpolate -from scipy.spatial import distance -from numpy.linalg import LinAlgError - -from KCSD import KCSD -import utility_functions as utils -import basis_functions as basis - -class KCSD2D(KCSD): - """KCSD2D - The 2D variant for the Kernel Current Source Density method. - - This estimates the Current Source Density, for a given configuration of - electrod positions and recorded potentials, in the case of 2D recording - electrodes. The method implented here is based on the original paper - by Jan Potworowski et.al. 2012. - """ - def __init__(self, ele_pos, pots, **kwargs): - """Initialize KCSD2D Class. - - Parameters - ---------- - ele_pos : numpy array - positions of electrodes - pots : numpy array - potentials measured by electrodes - **kwargs - configuration parameters, that may contain the following keys: - src_type : str - basis function type ('gauss', 'step', 'gauss_lim') - Defaults to 'gauss' - sigma : float - space conductance of the medium - Defaults to 1. - n_src_init : int - requested number of sources - Defaults to 1000 - R_init : float - demanded thickness of the basis element - Defaults to 0.23 - h : float - thickness of analyzed tissue slice - Defaults to 1. - xmin, xmax, ymin, ymax : floats - boundaries for CSD estimation space - Defaults to min(ele_pos(x)), and max(ele_pos(x)) - Defaults to min(ele_pos(y)), and max(ele_pos(y)) - ext_x, ext_y : float - length of space extension: x_min-ext_x ... x_max+ext_x - length of space extension: y_min-ext_y ... y_max+ext_y - Defaults to 0. - gdx, gdy : float - space increments in the estimation space - Defaults to 0.01(xmax-xmin) - Defaults to 0.01(ymax-ymin) - lambd : float - regularization parameter for ridge regression - Defaults to 0. - - Returns - ------- - None - - Raises - ------ - LinAlgError - Could not invert the matrix, try changing the ele_pos slightly - KeyError - Basis function (src_type) not implemented. See basis_functions.py for available - """ - super(KCSD2D, self).__init__(ele_pos, pots, **kwargs) - return - - def estimate_at(self): - """Defines locations where the estimation is wanted - Defines: - self.n_estm = self.estm_x.size - self.ngx, self.ngy = self.estm_x.shape - self.estm_x, self.estm_y : Locations at which CSD is requested. - - Parameters - ---------- - None - - Returns - ------- - None - """ - #Number of points where estimation is to be made. - nx = (self.xmax - self.xmin)/self.gdx - ny = (self.ymax - self.ymin)/self.gdy - #Making a mesh of points where estimation is to be made. - self.estm_x, self.estm_y = np.mgrid[self.xmin:self.xmax:np.complex(0,nx), - self.ymin:self.ymax:np.complex(0,ny)] - self.n_estm = self.estm_x.size - self.ngx, self.ngy = self.estm_x.shape - return - - def place_basis(self): - """Places basis sources of the defined type. - Checks if a given source_type is defined, if so then defines it - self.basis, This function gives locations of the basis sources, - Defines - source_type : basis_fuctions.basis_2D.keys() - self.R based on R_init - self.dist_max as maximum distance between electrode and basis - self.nsx, self.nsy = self.src_x.shape - self.src_x, self.src_y : Locations at which basis sources are placed. - - Parameters - ---------- - None - - Returns - ------- - None - """ - source_type = self.src_type - try: - self.basis = basis.basis_2D[source_type] - except: - print 'Invalid source_type for basis! available are:', basis.basis_2D.keys() - raise KeyError - #Mesh where the source basis are placed is at self.src_x - (self.src_x, self.src_y, self.R) = utils.distribute_srcs_2D(self.estm_x, - self.estm_y, - self.n_src_init, - self.ext_x, - self.ext_y, - self.R_init ) - self.n_src = self.src_x.size - self.nsx, self.nsy = self.src_x.shape - return - - def create_src_dist_tables(self): - src_loc = np.array((self.src_x.ravel(), self.src_y.ravel())) - est_loc = np.array((self.estm_x.ravel(), self.estm_y.ravel())) - self.src_ele_dists = distance.cdist(src_loc.T, self.ele_pos, 'euclidean') - self.src_estm_dists = distance.cdist(src_loc.T, est_loc.T, 'euclidean') - self.dist_max = max(np.max(self.src_ele_dists), np.max(self.src_estm_dists)) + self.R - return - - def forward_model(self, x, R, h, sigma, src_type): - """FWD model functions - Evaluates potential at point (x,0) by a basis source located at (0,0) - Eq 22 kCSD by Jan,2012 - - Parameters - ---------- - x : float - R : float - h : float - sigma : float - src_type : basis_2D.key - - Returns - ------- - pot : float - value of potential at specified distance from the source - """ - pot, err = integrate.dblquad(self.int_pot_2D, - -R, R, - lambda x: -R, - lambda x: R, - args=(x, R, h, src_type)) - pot *= 1./(2.0*np.pi*sigma) #Potential basis functions bi_x_y - return pot - - def int_pot_2D(self, xp, yp, x, R, h, basis_func): - """FWD model function. - Returns contribution of a point xp,yp, belonging to a basis source - support centered at (0,0) to the potential measured at (x,0), - integrated over xp,yp gives the potential generated by a - basis source element centered at (0,0) at point (x,0) - - Parameters - ---------- - xp, yp : floats or np.arrays - point or set of points where function should be calculated - x : float - position at which potential is being measured - R : float - The size of the basis function - h : float - thickness of slice - basis_func : method - Fuction of the basis source - - Returns - ------- - pot : float - """ - y = ((x-xp)**2 + yp**2)**(0.5) - if y < 0.00001: - y = 0.00001 - dist = np.sqrt(xp**2 + yp**2) - pot = np.arcsinh(h/y)*basis_func(dist, R) - return pot - -if __name__ == '__main__': - #Sample data, do not take this seriously - ele_pos = np.array([[-0.2, -0.2],[0, 0], [0, 1], [1, 0], [1,1], [0.5, 0.5], - [1.2, 1.2]]) - pots = np.array([[-1], [-1], [-1], [0], [0], [1], [-1.5]]) - k = KCSD2D(ele_pos, pots, - gdx=0.05, gdy=0.05, - xmin=-2.0, xmax=2.0, - ymin=-2.0, ymax=2.0, - src_type='gauss') - k.cross_validate() - #print k.values('CSD') - #print k.values('POT') diff --git a/KCSD3D.py b/KCSD3D.py deleted file mode 100644 index 367b19e6..00000000 --- a/KCSD3D.py +++ /dev/null @@ -1,321 +0,0 @@ -""" -This script is used to generate Current Source Density Estimates, -using the kCSD method Jan et.al (2012) for 3D case. - -These scripts are based on Grzegorz Parka's, -Google Summer of Code 2014, INFC/pykCSD - -This was written by : -Chaitanya Chintaluri -Laboratory of Neuroinformatics, -Nencki Institute of Experimental Biology, Warsaw. -""" -import numpy as np -from scipy.spatial import distance -from scipy import special, interpolate, integrate -try: - from skmonaco import mcmiser - skmonaco_available = True - import multiprocessing - num_cores = multiprocessing.cpu_count() -except ImportError: - skmonaco_available = False - -from KCSD import KCSD -import utility_functions as utils -import basis_functions as basis - -class KCSD3D(KCSD): - """KCSD3D - The 3D variant for the Kernel Current Source Density method. - - This estimates the Current Source Density, for a given configuration of - electrod positions and recorded potentials, in the case of 2D recording - electrodes. The method implented here is based on the original paper - by Jan Potworowski et.al. 2012. - """ - def __init__(self, ele_pos, pots, **kwargs): - """Initialize KCSD3D Class. - - Parameters - ---------- - ele_pos : numpy array - positions of electrodes - pots : numpy array - potentials measured by electrodes - **kwargs - configuration parameters, that may contain the following keys: - src_type : str - basis function type ('gauss', 'step', 'gauss_lim') - Defaults to 'gauss' - sigma : float - space conductance of the medium - Defaults to 1. - n_src_init : int - requested number of sources - Defaults to 1000 - R_init : float - demanded thickness of the basis element - Defaults to 0.23 - h : float - thickness of analyzed tissue slice - Defaults to 1. - xmin, xmax, ymin, ymax, zmin, zmax : floats - boundaries for CSD estimation space - Defaults to min(ele_pos(x)), and max(ele_pos(x)) - Defaults to min(ele_pos(y)), and max(ele_pos(y)) - Defaults to min(ele_pos(z)), and max(ele_pos(z)) - ext_x, ext_y, ext_z : float - length of space extension: xmin-ext_x ... xmax+ext_x - length of space extension: ymin-ext_y ... ymax+ext_y - length of space extension: zmin-ext_z ... zmax+ext_z - Defaults to 0. - gdx, gdy, gdz : float - space increments in the estimation space - Defaults to 0.01(xmax-xmin) - Defaults to 0.01(ymax-ymin) - Defaults to 0.01(zmax-zmin) - lambd : float - regularization parameter for ridge regression - Defaults to 0. - - Returns - ------- - None - - Raises - ------ - LinAlgError - Could not invert the matrix, try changing the ele_pos slightly - KeyError - Basis function (src_type) not implemented. See basis_functions.py for available - """ - super(KCSD3D, self).__init__(ele_pos, pots, **kwargs) - return - - def estimate_at(self): - """Defines locations where the estimation is wanted - Defines: - self.n_estm = self.estm_x.size - self.ngx, self.ngy, self.ngz = self.estm_x.shape - self.estm_x, self.estm_y, self.estm_z : Pts. at which CSD is requested - - Parameters - ---------- - None - - Returns - ------- - None - """ - #Number of points where estimation is to be made. - nx = (self.xmax - self.xmin)/self.gdx - ny = (self.ymax - self.ymin)/self.gdy - nz = (self.zmax - self.zmin)/self.gdz - #Making a mesh of points where estimation is to be made. - self.estm_x, self.estm_y, self.estm_z = np.mgrid[self.xmin:self.xmax:np.complex(0,nx), - self.ymin:self.ymax:np.complex(0,ny), - self.zmin:self.zmax:np.complex(0,nz)] - self.n_estm = self.estm_x.size - self.ngx, self.ngy, self.ngz = self.estm_x.shape - return - - def place_basis(self): - """Places basis sources of the defined type. - Checks if a given source_type is defined, if so then defines it - self.basis, This function gives locations of the basis sources, - Defines - source_type : basis_fuctions.basis_2D.keys() - self.R based on R_init - self.dist_max as maximum distance between electrode and basis - self.nsx, self.nsy, self.nsz = self.src_x.shape - self.src_x, self.src_y, self.src_z : Locations at which basis sources are placed. - - Parameters - ---------- - None - - Returns - ------- - None - """ - #If Valid basis source type passed? - source_type = self.src_type - try: - self.basis = basis.basis_3D[source_type] - except: - print 'Invalid source_type for basis! available are:', basis.basis_3D.keys() - raise KeyError - #Mesh where the source basis are placed is at self.src_x - (self.src_x, self.src_y, self.src_z, self.R) = utils.distribute_srcs_3D(self.estm_x, - self.estm_y, - self.estm_z, - self.n_src_init, - self.ext_x, - self.ext_y, - self.ext_z, - self.R_init) - - self.n_src = self.src_x.size - self.nsx, self.nsy, self.nsz = self.src_x.shape - return - - def create_src_dist_tables(self): - """Creates distance tables between sources, electrode and estm points - - Parameters - ---------- - None - - Returns - ------- - None - """ - src_loc = np.array((self.src_x.ravel(), - self.src_y.ravel(), - self.src_z.ravel())) - est_loc = np.array((self.estm_x.ravel(), - self.estm_y.ravel(), - self.estm_z.ravel())) - self.src_ele_dists = distance.cdist(src_loc.T, self.ele_pos, 'euclidean') - self.src_estm_dists = distance.cdist(src_loc.T, est_loc.T, 'euclidean') - self.dist_max = max(np.max(self.src_ele_dists), np.max(self.src_estm_dists)) + self.R - return - - def forward_model(self, x, R, h, sigma, src_type): - """FWD model functions - Evaluates potential at point (x,0) by a basis source located at (0,0) - Utlizies sk monaco monte carlo method if available, otherwise defaults - to scipy integrate - - Parameters - ---------- - x : float - R : float - h : float - sigma : float - src_type : basis_3D.key - - Returns - ------- - pot : float - value of potential at specified distance from the source - """ - if src_type.__name__ == "gauss_3D": - if x == 0: x=0.0001 - pot = special.erf(x/(np.sqrt(2)*R/3.0)) / x - elif src_type.__name__ == "gauss_lim_3D": - if x == 0: x=0.0001 - d = R/3. - if x < R: - #4*pi*((1/a)*(integrate(r**2 * exp(-r**2 / (2*d**2)) *dr ) between 0 and a ) + - #(integrate(r *exp(-r**2 / (2*d**2)) * dr) between a and 3*d)) - e = np.exp(-(x/ (np.sqrt(2)*d))**2) - erf = special.erf(x / (np.sqrt(2)*d)) - pot = 4* np.pi * ( (d**2)*(e - np.exp(-4.5)) + - (1/x)*((np.sqrt(np.pi/2)*(d**3)*erf) - x*(d**2)*e)) - else: - #4*pi*integrate((r**2)*exp(-(r**2 / (2*d**2)))*dr) between 0 and 3*d - pot = 15.28828*(d)**3 / x - pot /= (np.sqrt(2*np.pi)*d)**3 - elif src_type.__name__ == "step_3D": - Q = 4.*np.pi*(R**3)/3. - if x < R: - pot = (Q * (3 - (x/R)**2)) / (2.*R) - else: - pot = Q / x - pot *= 3/(4*np.pi*R**3) - else: - if skmonaco_available: - pot, err = mcmiser(self.int_pot_3D_mc, - npoints=1e5, - xl=[-R, -R, -R], - xu=[R, R, R], - seed=42, - nprocs=num_cores, - args=(x, R, h, src_type)) - else: - pot, err = integrate.tplquad(self.int_pot_3D, - -R, - R, - lambda x: -R, - lambda x: R, - lambda x, y: -R, - lambda x, y: R, - args=(x, R, h, src_type)) - pot *= 1./(4.0*np.pi*sigma) - return pot - - def int_pot_3D(self, xp, yp, zp, x, R, h, basis_func): - """FWD model function. - Returns contribution of a point xp,yp, belonging to a basis source - support centered at (0,0) to the potential measured at (x,0), - integrated over xp,yp gives the potential generated by a - basis source element centered at (0,0) at point (x,0) - - Parameters - ---------- - xp, yp, zp : floats or np.arrays - point or set of points where function should be calculated - x : float - position at which potential is being measured - R : float - The size of the basis function - h : float - thickness of slice - basis_func : method - Fuction of the basis source - - Returns - ------- - pot : float - """ - y = ((x-xp)**2 + yp**2 + zp**2)**0.5 - if y < 0.00001: - y = 0.00001 - dist = np.sqrt(xp**2 + yp**2 + zp**2) - pot = 1.0/y - pot *= basis_func(dist, R) - return pot - - def int_pot_3D_mc(self, xyz, x, R, h, basis_func): - """ - The same as int_pot_3D, just different input: x,y,z <-- xyz (tuple) - FWD model function, using Monte Carlo Method of integration - Returns contribution of a point xp,yp, belonging to a basis source - support centered at (0,0) to the potential measured at (x,0), - integrated over xp,yp gives the potential generated by a - basis source element centered at (0,0) at point (x,0) - - Parameters - ---------- - xp, yp, zp : floats or np.arrays - point or set of points where function should be calculated - x : float - position at which potential is being measured - R : float - The size of the basis function - h : float - thickness of slice - basis_func : method - Fuction of the basis source - - Returns - ------- - pot : float - """ - xp, yp, zp = xyz - return self.int_pot_3D(xp, yp, zp, x, R, h, basis_func) - - -if __name__ == '__main__': - ele_pos = np.array([(0, 0, 0), (0, 0, 1), (0, 1, 0), (1, 0, 0), - (0, 1, 1), (1, 1, 0), (1, 0, 1), (1, 1, 1), - (0.5, 0.5, 0.5)]) - pots = np.array([[-0.5], [0], [-0.5], [0], [0], [0.2], [0], [0], [1]]) - params = {} - k = KCSD3D(ele_pos, pots, - gdx=0.02, gdy=0.02, gdz=0.02, - n_src_init=1000, src_type='gauss_lim') - k.cross_validate() - #k.cross_validate(Rs=np.array(0.14).reshape(1)) - #k.cross_validate(Rs=np.array((0.01,0.02,0.04))) diff --git a/MoIKCSD.py b/MoIKCSD.py deleted file mode 100644 index d97215cc..00000000 --- a/MoIKCSD.py +++ /dev/null @@ -1,156 +0,0 @@ -""" -This script is used to generate Current Source Density Estimates, -using the kCSD method described in Ness et.al (2015) for 2D case. - -This was written by : -Chaitanya Chintaluri, -Laboratory of Neuroinformatics, -Nencki Institute of Exprimental Biology, Warsaw. -""" -import numpy as np -from scipy import integrate - -from KCSD2D import KCSD2D - -class MoIKCSD(KCSD2D): - """MoIKCSD - CSD while including the forward modeling effects of saline. - - This estimates the Current Source Density, for a given configuration of - electrod positions and recorded potentials, in the case of 2D recording - electrodes from an MEA electrode plane using the Method of Images. - The method implented here is based on kCSD method by Jan Potworowski - et.al. 2012, which was extended in Ness, Chintaluri 2015 for MEA. - """ - def __init__(self, ele_pos, pots, **kwargs): - """Initialize MoIKCSD Class. - - Parameters - ---------- - ele_pos : numpy array - positions of electrodes - pots : numpy array - potentials measured by electrodes - **kwargs - configuration parameters, that may contain the following keys: - src_type : str - basis function type ('gauss', 'step', 'gauss_lim') - Defaults to 'gauss' - sigma : float - space conductance of the medium - Defaults to 1. - sigma_S : float - conductance of the saline (medium) - Default is 5 S/m (5 times more conductive) - n_src_init : int - requested number of sources - Defaults to 1000 - R_init : float - demanded thickness of the basis element - Defaults to 0.23 - h : float - thickness of analyzed tissue slice - Defaults to 1. - xmin, xmax, ymin, ymax : floats - boundaries for CSD estimation space - Defaults to min(ele_pos(x)), and max(ele_pos(x)) - Defaults to min(ele_pos(y)), and max(ele_pos(y)) - ext_x, ext_y : float - length of space extension: x_min-ext_x ... x_max+ext_x - length of space extension: y_min-ext_y ... y_max+ext_y - Defaults to 0. - gdx, gdy : float - space increments in the estimation space - Defaults to 0.01(xmax-xmin) - Defaults to 0.01(ymax-ymin) - lambd : float - regularization parameter for ridge regression - Defaults to 0. - MoI_iters : int - Number of interations in method of images. - Default is 20 - Returns - ------- - None - """ - self.MoI_iters = kwargs.get('MoI_iters', 20) - self.sigma_S = kwargs.get('sigma_S', 5.0) - self.sigma = kwargs.get('sigma', 1.0) - #Eq 6, Ness (2015) - W_TS = (self.sigma - self.sigma_S) / (self.sigma + self.sigma_S) - self.iters = np.arange(self.MoI_iters) + 1 - self.iter_factor = W_TS**self.iters - super(MoIKCSD, self).__init__(ele_pos, pots, **kwargs) - return - - def forward_model(self, x, R, h, sigma, src_type): - """FWD model functions - Evaluates potential at point (x,0) by a basis source located at (0,0) - Eq 22 kCSD by Jan,2012 - - Parameters - ---------- - x : float - R : float - h : float - sigma : float - src_type : basis_2D.key - - Returns - ------- - pot : float - value of potential at specified distance from the source - """ - pot, err = integrate.dblquad(self.int_pot_2D_moi, -R, R, - lambda x: -R, - lambda x: R, - args=(x, R, h, src_type)) - pot *= 1./(2.0*np.pi*sigma) - return pot - - def int_pot_2D_moi(self, xp, yp, x, R, h, basis_func): - """FWD model function. Incorporates the Method of Images. - Returns contribution of a point xp,yp, belonging to a basis source - support centered at (0,0) to the potential measured at (x,0), - integrated over xp,yp gives the potential generated by a - basis source element centered at (0,0) at point (x,0) - #Eq 20, Ness(2015) - - Parameters - ---------- - xp, yp : floats or np.arrays - point or set of points where function should be calculated - x : float - position at which potential is being measured - R : float - The size of the basis function - h : float - thickness of slice - basis_func : method - Fuction of the basis source - - Returns - ------- - pot : float - """ - L = ((x-xp)**2 + yp**2)**(0.5) - if L < 0.00001: - L = 0.00001 - correction = np.arcsinh((h-(2*h*self.iters))/L) + np.arcsinh((h+(2*h*self.iters))/L) - pot = np.arcsinh(h/L) + np.sum(self.iter_factor*correction) - dist = np.sqrt(xp**2 + yp**2) - pot *= basis_func(dist, R) #Eq 20, Ness et.al. - return pot - -if __name__ == '__main__': - ele_pos = np.array([[-0.2, -0.2],[0, 0], [0, 1], [1, 0], [1,1], [0.5, 0.5], - [1.2, 1.2]]) - pots = np.array([[-1], [-1], [-1], [0], [0], [1], [-1.5]]) - k = MoIKCSD(ele_pos, pots, - gdx=0.05, gdy=0.05, - xmin=-2.0, xmax=2.0, - ymin=-2.0, ymax= 2.0) - k.cross_validate() - #print k.values() - - - diff --git a/basis_functions.py b/basis_functions.py index bb23a202..c96d5a90 100644 --- a/basis_functions.py +++ b/basis_functions.py @@ -10,6 +10,8 @@ Laboratory of Neuroinformatics, Nencki Institute of Experimental Biology, Warsaw. """ +from __future__ import division + import numpy as np def gauss(d, stdev, dim): diff --git a/tests/csd_profile.py b/tests/csd_profile.py deleted file mode 100644 index 3edcb724..00000000 --- a/tests/csd_profile.py +++ /dev/null @@ -1,290 +0,0 @@ -''' -This script is used to generate dummy CSD sources, -to test the various kCSD methods - -This script is in alpha phase. - -This was written by : -Michal Czerwinski, Chaitanya Chintaluri, -Laboratory of Neuroinformatics, -Nencki Institute of Exprimental Biology, Warsaw. -''' -import numpy as np -import matplotlib.pyplot as plt -import matplotlib.cm as cm - -from matplotlib import gridspec -from mpl_toolkits.mplot3d import axes3d -from numpy import exp - -def get_states_1D(seed, n=1): - """ - Used in the random seed generation - creates a matrix that will generate seeds, here for gaussians: - amplitude (-1., 1.), location (0,1)*ndim, sigma(0,1) - """ - ndim = 1 - if seed == 0: - states = np.array([1.,0.5,0.5], ndmin=2) - rstate = np.random.RandomState(seed) - states = rstate.random_sample(n*(ndim+2)).reshape((n,(ndim+2))) - states[:,0] = (2*states[:,0]) -1. - return states, rstate - -def add_1d_gaussians(x, states): - '''Function used for adding multiple 1D gaussians''' - f = np.zeros(x.shape) - for i in xrange(states.shape[0]): - gauss = states[i,0]*np.exp(-((x-states[i,1])**2)/(2.*states[i,2]))*(2*np.pi*states[i,2])**-0.5 - f += gauss - #f *= np.hanning(gauss.size) #while this is cool in principle, its unpredictable - return f - -def gauss_1d_mono(x, seed=0): - '''Random monopole in 1D''' - states, rstate = get_states_1D(seed, n=1) - f = add_1d_gaussians(x, states) - return f - -def gauss_1d_dipole(x, seed=0): - '''Random dipole source in 1D''' - states, rstate = get_states_1D(seed, n=1) - offset = rstate.random_sample(1) - 0.5 - states = np.tile(states, (2,1)) - states[1,0] *= -1. #A Sink - states[1,1] += offset - f = add_1d_gaussians(x, states) - return f - -def get_states_2D(seed): - """ - Used in the random seed generation for 2d sources - """ - rstate = np.random.RandomState(seed) - states = rstate.random_sample(24) - states[0:12] = 2*states[0:12] -1. - return states - -def gauss_2d_large(x, y, seed=0): - '''random quadpolar'large source' profile in 2012 paper in 2D''' - states = get_states_2D(seed) - z = 0 - zz = states[0:4] - zs = states[4:8] - mag = states[8:12] - loc = states[12:20] - scl = states[20:24] - f1 = mag[0]*exp( (-1*(x-loc[0])**2 - (y-loc[4])**2) /scl[0])* exp(-(z-zz[0])**2 / zs[0]) /exp(-(zz[0])**2/zs[0]) - f2 = mag[1]*exp( (-2*(x-loc[1])**2 - (y-loc[5])**2) /scl[1])* exp(-(z-zz[1])**2 / zs[1]) /exp(-(zz[1])**2/zs[1]); - f3 = mag[2]*exp( (-3*(x-loc[2])**2 - (y-loc[6])**2) /scl[2])* exp(-(z-zz[2])**2 / zs[2]) /exp(-(zz[2])**2/zs[2]); - f4 = mag[3]*exp( (-4*(x-loc[3])**2 - (y-loc[7])**2) /scl[3])* exp(-(z-zz[3])**2 / zs[3]) /exp(-(zz[3])**2/zs[3]); - f = f1+f2+f3+f4 - return f - -def gauss_2d_small(x, y, seed=0): - '''random quadpolar small source in 2D''' - def gauss2d(x,y,p): - """ - p: list of parameters of the Gauss-function - [XCEN,YCEN,SIGMAX,SIGMAY,AMP,ANGLE] - SIGMA = FWHM / (2*sqrt(2*log(2))) - ANGLE = rotation of the X,Y direction of the Gaussian in radians - Returns - ------- - the value of the Gaussian described by the parameters p - at position (x,y) - """ - rcen_x = p[0] * np.cos(p[5]) - p[1] * np.sin(p[5]) - rcen_y = p[0] * np.sin(p[5]) + p[1] * np.cos(p[5]) - xp = x * np.cos(p[5]) - y * np.sin(p[5]) - yp = x * np.sin(p[5]) + y * np.cos(p[5]) - g = p[4]*np.exp(-(((rcen_x-xp)/p[2])**2+ - ((rcen_y-yp)/p[3])**2)/2.) - return g - states = get_states_2D(seed) - angle = states[18]*180. - x_amp = 0.038 - y_amp = 0.056 - f1 = gauss2d(x,y,[states[12],states[14],x_amp,y_amp,0.5,angle]) - f2 = gauss2d(x,y,[states[12],states[15],x_amp,y_amp,-0.5,angle]) - f3 = gauss2d(x,y,[states[13],states[14],x_amp,y_amp,0.5,angle]) - f4 = gauss2d(x,y,[states[13],states[15],x_amp,y_amp,-0.5,angle]) - f = f1+f2+f3+f4 - return f - -def get_states_3D(seed): - """ - Used in the random seed generation for 3D sources - """ - rstate = np.random.RandomState(seed) #seed here! - states = rstate.random_sample(24) - return states - -def gauss_3d_small(x, y, z, seed=0): - '''A random quadpole small souce in 3D''' - states = get_states_3D(seed) - x0, y0, z0 = states[0:3] - x1, y1, z1 = states[3:6] - if states[6] < 0.01: - states[6] *= 25 - sig_2 = states[6] / 75. - p1, p2, p3 = (ii*0.5 for ii in states[8:11]) - A = (2*np.pi*sig_2)**-1 - f1 = A*np.exp( (-(x-x0)**2 -(y-y0)**2 -(z-z0)**2) / (2*sig_2) ) - f2 = -1*A*np.exp( (-(x-x1)**2 -(y-y1)**2 -(z-z1)**2) / (2*sig_2) ) - x2 = np.modf(x0+p1)[0] - y2 = np.modf(y0+p2)[0] - z2 = np.modf(z0+p3)[0] - f3 = A*np.exp( (-(x-x2)**2 -(y-y2)**2 -(z-z2)**2) / (2*sig_2) ) - x3 = np.modf(x1+p1)[0] - y3 = np.modf(y1+p2)[0] - z3 = np.modf(z1+p3)[0] - f4 = -1*A*np.exp( (-(x-x3)**2 -(y-y3)**2 -(z-z3)**2) / (2*sig_2) ) - f = f1+f2+f3+f4 - return f - -def gauss_3d_large(x, y, z, seed=0): - '''A random dipolar Large source in 3D''' - states = get_states_3D(seed) - x0, y0, z0 = states[7:10] - x1, y1, z1 = states[10:13] - if states[1] < 0.01: - states[1] *= 25 - sig_2 = states[1] * 5 - A = (2*np.pi*sig_2)**-1 - f1 = A*np.exp( (-(x-x0)**2 -(y-y0)**2 -(z-z0)**2) / (2*sig_2) ) - f2 = -1*A*np.exp( (-(x-x1)**2 -(y-y1)**2 -(z-z1)**2) / (2*sig_2) ) - f = f1+f2 - return f - -def jan_2d_small_f(x, y): - '''Source from Jan 2012 kCSD paper''' - def gauss2d(x,y,p): - """ - p: list of parameters of the Gauss-function - [XCEN,YCEN,SIGMAX,SIGMAY,AMP,ANGLE] - SIGMA = FWHM / (2*sqrt(2*log(2))) - ANGLE = rotation of the X,Y direction of the Gaussian in radians - Returns - ------- - the value of the Gaussian described by the parameters p - at position (x,y) - """ - rcen_x = p[0] * np.cos(p[5]) - p[1] * np.sin(p[5]) - rcen_y = p[0] * np.sin(p[5]) + p[1] * np.cos(p[5]) - xp = x * np.cos(p[5]) - y * np.sin(p[5]) - yp = x * np.sin(p[5]) + y * np.cos(p[5]) - g = p[4]*np.exp(-(((rcen_x-xp)/p[2])**2+ - ((rcen_y-yp)/p[3])**2)/2.) - return g - f1 = gauss2d(x,y,[0.3,0.7,0.038,0.058,0.5,0.]) - f2 = gauss2d(x,y,[0.3,0.6,0.038,0.058,-0.5,0.]) - f3 = gauss2d(x,y,[0.45,0.7,0.038,0.058,0.5,0.]) - f4 = gauss2d(x,y,[0.45,0.6,0.038,0.058,-0.5,0.]) - f = f1+f2+f3+f4 - return f - -def jan_2d_large_f(x, y): - '''Fixed 'large source' profile in 2012 paper''' - z = 0 - zz = [0.4, -0.3, -0.1, 0.6] - zs = [0.2, 0.3, 0.4, 0.2] - f1 = 0.5965*exp( (-1*(x-0.1350)**2 - (y-0.8628)**2) /0.4464)* exp(-(z-zz[0])**2 / zs[0]) /exp(-(zz[0])**2/zs[0]) - f2 = -0.9269*exp( (-2*(x-0.1848)**2 - (y-0.0897)**2) /0.2046)* exp(-(z-zz[1])**2 / zs[1]) /exp(-(zz[1])**2/zs[1]); - f3 = 0.5910*exp( (-3*(x-1.3189)**2 - (y-0.3522)**2) /0.2129)* exp(-(z-zz[2])**2 / zs[2]) /exp(-(zz[2])**2/zs[2]); - f4 = -0.1963*exp( (-4*(x-1.3386)**2 - (y-0.5297)**2) /0.2507)* exp(-(z-zz[3])**2 / zs[3]) /exp(-(zz[3])**2/zs[3]); - f = f1+f2+f3+f4 - return f - -def gauss_3d_dipole_f(x, y, z): - '''Fixed dipole in 3 dimensions of the volume''' - x0, y0, z0 = 0.3, 0.7, 0.3 - x1, y1, z1 = 0.6, 0.5, 0.7 - sig_2 = 0.023 - A = (2*np.pi*sig_2)**-1 - f1 = A*np.exp( (-(x-x0)**2 -(y-y0)**2 -(z-z0)**2) / (2*sig_2) ) - f2 = -1*A*np.exp( (-(x-x1)**2 -(y-y1)**2 -(z-z1)**2) / (2*sig_2) ) - f = f1+f2 - return f - -def gauss_3d_mono1_f(x, y, z): - '''Fixed monopole in 3D at the center of the volume space''' - x0, y0, z0 = 0.5, 0.5, 0.5 - sig_2 = 0.023 - A = (2*np.pi*sig_2)**-1 - f1 = A*np.exp( (-(x-x0)**2 -(y-y0)**2 -(z-z0)**2) / (2*sig_2) ) - return f1 - -def gauss_3d_mono2_f(x, y, z): - '''Fixed monopole in 3D Offcentered wrt volume''' - x0, y0, z0 = 0.41, 0.41, 0.585 - sig_2 = 0.023 - A = (2*np.pi*sig_2)**-1 - f1 = A*np.exp( (-(x-x0)**2 -(y-y0)**2 -(z-z0)**2) / (2*sig_2) ) - return f1 - -def gauss_3d_mono3_f(x, y, z): - '''Fixed monopole in 3D Offcentered wrt volume''' - x0, y0, z0 = 0.55555556, 0.55555556, 0.55555556 - stdev = 0.3 - h = 1./((2*np.pi)**0.5 * stdev)**3 - c = 0.5*stdev**(-2) - f1 = h*np.exp(-c*((x - x0)**2 + (y - y0)**2 + (z - z0)**2)) - return f1 - -def neat_4d_plot(x, y, z, t, z_steps=5, cmap=cm.bwr_r): - '''Used to show 3D csd profile''' - t_max = np.max(np.abs(t)) - levels = np.linspace(-1*t_max, t_max, 15) - ind_interest = np.mgrid[0:z.shape[2]:np.complex(0,z_steps+2)] - ind_interest = np.array(ind_interest, dtype=np.int)[1:-1] - fig = plt.figure(figsize=(4,12)) - height_ratios = [1 for i in range(z_steps)] - height_ratios.append(0.1) - gs = gridspec.GridSpec(z_steps+1, 1, height_ratios=height_ratios) - for ii, idx in enumerate(ind_interest): - ax = plt.subplot(gs[ii,0]) - im = plt.contourf(chrg_x[:,:,idx], chrg_y[:,:,idx], t[:,:,idx], - levels=levels, cmap=cmap) - ax.get_xaxis().set_visible(False) - ax.get_yaxis().set_visible(False) - title = str(z[:,:,idx][0][0])[:4] - ax.set_title(label=title, fontdict={'x':0.8, 'y':0.7}) - ax.set_aspect('equal') - cax = plt.subplot(gs[z_steps,0]) - cbar = plt.colorbar(im, cax=cax, orientation='horizontal') - cbar.set_ticks(levels[::2]) - cbar.set_ticklabels(np.around(levels[::2], decimals=2)) - gs.tight_layout(fig, rect=[0, 0.03, 1, 0.95]) - #plt.tight_layout() - -if __name__=='__main__': - seed = 3 - - # 1D CASE - csd_profile = gauss_1d_mono - chrg_x = np.arange(0., 1., 1./50.) - f = csd_profile(chrg_x, seed) - plt.plot(chrg_x, f) - - # #2D CASE - # csd_profile = gauss_2d_large - # states[0:12] = 2*states[0:12] -1. #obtain values between -1 and 1 - # chrg_x, chrg_y = np.mgrid[0.:1.:50j, - # 0.:1.:50j] - # f = csd_profile(chrg_x, chrg_y, seed=seed) - # fig = plt.figure(1) - # ax1 = plt.subplot(111, aspect='equal') - # im = ax1.contourf(chrg_x, chrg_y, f, 15, cmap=cm.bwr_r) - # cbar = plt.colorbar(im, shrink=0.5) - # plt.show() - - # #3D CASE - # csd_profile = gauss_3d_small - # chrg_x, chrg_y, chrg_z = np.mgrid[0.:1.:50j, - # 0.:1.:50j, - # 0.:1.:50j] - # f = csd_profile(chrg_x, chrg_y, chrg_z, seed=seed) - # neat_4d_plot(chrg_x, chrg_y, chrg_z, f) - - plt.show() diff --git a/tests/test_kCSD1D.py b/tests/test_kCSD1D.py deleted file mode 100644 index f877c6ea..00000000 --- a/tests/test_kCSD1D.py +++ /dev/null @@ -1,219 +0,0 @@ -''' -This script is used to test the Current Source Density Estimates, -using the kCSD method Potworowski et.al (2012) for 1D case - -This script is in alpha phase. - -This was written by: -Michal Czerwinski, Chaitanya Chintaluri, -Laboratory of Neuroinformatics, -Nencki Institute of Exprimental Biology, Warsaw. -''' -import time -import os -import sys -sys.path.append('..') - -import numpy as np -from scipy.interpolate import interp1d -from scipy.integrate import simps -import matplotlib.pyplot as plt - -import csd_profile as CSD -from KCSD1D import KCSD1D - -def generate_csd_1D(csd_profile, csd_seed, start_x=0., end_x=1., res_x=50): - """ - Gives CSD profile at the requested spatial location, at 'res' resolution - """ - chrg_pos_x = np.linspace(start_x, end_x, res_x) - f = csd_profile(chrg_pos_x, seed=csd_seed) - return chrg_pos_x, f - -def grid(x, z, resX=100): - """ - Convert 2 column data to two (2) matplotlib axes to plot, interpolates to given resolution using scipy.interp1d - """ - z = z.flatten() - xi = np.linspace(min(x), max(x), resX) - X = x.flatten() - f = interp1d(X,z) - Y = f(xi) - return xi, Y - -def generate_electrodes(xlims=[0.1,0.9], res=5): - """ - Places electrodes lineary - """ - ele_x= np.linspace(xlims[0], xlims[1], res) - ele_x = ele_x.flatten() - return ele_x - -def generate_csd_space(xlims=[0.0,1.0], gdX=0.01): - """ - Area of interest for the CSD reconstruction - """ - nx = (xlims[1] - xlims[0])/gdX - space_X = np.linspace(xlims[0], xlims[1], nx) - return space_X - -def make_plots(title, - chrg_x, csd, - ele_x, pots, - csd_x, est_csd, est_pot, - true_pot=None): - """ - Shows 2 subplots - 1_a true CSD generated based on the random seed given - 1_b kernel CSD the output of this method - 2_a pots - true potentials, generated by simpsons rule integration - 2_b Kernel interpolated LFP, - """ - fig = plt.figure(figsize=(7,10)) - #CSDs - ax1 = plt.subplot(211) - if np.array(csd_x).any() != False: - im1b = ax1.plot(csd_x, est_csd[:,0], 'r', label='kCSD', linewidth=3) - im1a = ax1.plot(chrg_x, csd, 'g', label = 'CSD', linestyle='--', linewidth=3) - ax1.plot(ele_x, np.zeros_like(ele_x), 'ko',markersize=2.) - plt.legend() - ax1.set_xlim([0.,1.]) - #ax1.set_ylim(ax1.get_ylim()[::-1]) #Zero on the top --ASK?! - max_csd = np.maximum(max(np.abs(csd)), max(np.abs(est_csd[:,0]))) - max_csd += max_csd*0.2 - ax1.set_ylim([-max_csd, max_csd]) - ax1.set_xlabel('Depth mm') - ax1.set_ylabel('CSD mA/mm') - #Potentails - ax2 = plt.subplot(212) - ax2.plot( ele_x, np.zeros_like(ele_x),'ko',markersize=2.) - im2b = ax2.plot(csd_x, est_pot, 'b', label = 'kPOT', linewidth=3) - im2a = ax2.plot(chrg_x,true_pot, color = 'orange', - linestyle='--', label='TruePOT', linewidth=3) - ax2.set_xlim([0.,1.]) - #ax2.set_ylim(ax2.get_ylim()[::-1]) #Zero on the top --ASK?! - ax2.plot(ele_x, pots, 'kx', markersize=8.) - max_pots = np.maximum(max(np.abs(true_pot)), max(np.abs(est_pot))) - max_pots += max_pots*0.2 - ax2.set_xlabel('Depth mm') - ax2.set_ylabel('Potential mV') - ax2.set_ylim([-max_pots, max_pots]) - fig.suptitle(title) - plt.legend() - plt.show() - return - -def integrate_1D(x0, csd_x, csd, h): - m = np.sqrt((csd_x-x0)**2 + h**2) - abs(csd_x-x0) - y = csd * m - I = simps(y, csd_x) - return I - -def calculate_potential_1D(csd, measure_locations, csd_space_x, h): - sigma = 1.0 - pots = np.zeros(len(measure_locations)) - for ii in range(len(measure_locations)): - pots[ii] = integrate_1D(measure_locations[ii], csd_space_x, csd, h) - pots *= 1/(2.*sigma) #eq.: 26 from Potworowski et al - return pots - -def electrode_config(ele_lims, ele_res, true_csd, t_csd_x, h): - """ - creates electrodes positions, and potentials on them - electrode lims, electrode resolution, profile, states - """ - ele_x = generate_electrodes(ele_lims, ele_res) - pots = calculate_potential_1D(true_csd, ele_x, t_csd_x, h) - ele_pos = ele_x.reshape((len(ele_x), 1)) - return ele_pos, pots - -def do_kcsd(ele_pos, pots, **params): - """ - Function that calls the KCSD2D module - """ - num_ele = len(ele_pos) - pots = pots.reshape(num_ele, 1) - k = KCSD1D(ele_pos, pots, **params) - #k.cross_validate(Rs=np.arange(0.01,0.2,0.01), lambdas= np.logspace(15,-25,25)) - k.cross_validate(Rs=np.array([0.275]), lambdas=np.logspace(15,-25, 35)) - est_csd = k.values() - est_pot = k.values('POT') - return k, est_csd, est_pot - -def main_loop(csd_profile, csd_seed, total_ele): - """ - Loop that decides the random number seed for the CSD profile, - electrode configurations and etc. - """ - csd_name = csd_profile.func_name - print 'Using sources %s - Seed: %d ' % (csd_name, csd_seed) - h = 10. - - #TrueCSD - start_x, end_x, csd_res = [0.,1.,100] - t_csd_x, true_csd = generate_csd_1D(csd_profile, csd_seed, - start_x=start_x, - end_x=end_x, - res_x=csd_res) - - #Electrodes - ele_res = int(total_ele) - ele_lims = [0.10, 0.9] - ele_pos, pots = electrode_config(ele_lims, ele_res, true_csd, t_csd_x, h) - num_ele = ele_pos.shape[0] - print 'Number of electrodes:', num_ele - x_array_pots, true_pots = electrode_config(ele_lims, 100, true_csd, t_csd_x, h) - - #kCSD estimation - gdX = 0.01 - x_lims = [0.,1.] #CSD estimation place - tic = time.time() #time it - k, est_csd, est_pot = do_kcsd(ele_pos, pots, h=h, gdx=gdX, - xmin=x_lims[0], xmax=x_lims[1], n_src_init=300) - toc = time.time() - tic - - #RMS of estimation - gives estimate of how good the reconstruction was - chr_x, test_csd = generate_csd_1D(csd_profile, csd_seed, - start_x=x_lims[0], end_x=x_lims[1], - res_x=int((x_lims[1]-x_lims[0])/gdX)) - rms = np.linalg.norm(abs(test_csd - est_csd[:,0])) - rms /= np.linalg.norm(test_csd) - - #Plots - title ="Lambda: %0.2E; R: %0.2f; CV_Error: %0.2E; RMS_Error: %0.2E; Time: %0.2f" %(k.lambd, k.R, k.cv_error, rms, toc) - make_plots(title, t_csd_x, true_csd, ele_pos, pots, k.estm_x, est_csd, est_pot, true_pots) - return - -def test_calculating_potentials(csd_seed): - csd_profile = CSD.gauss_1d_mono - plt.figure('csd') - _csdx = np.linspace(0,1,100) - plt.plot(_csdx, CSD.gauss_1d_mono(_csdx, csd_seed)) - #constant measure_locations, constant csd_space_lims, different csd resolution - plt.figure('changing csd res') - for i in xrange(5): - csdres = 10+i*10 - (x, V) = electrode_config([0.,1.], 1000, csd_profile, [0.,1.], csd_seed, CSDres=csdres) - #ans = calculate_potential_1D(csd_profile, measure_locations, states, csd_space_lims=[0.,1.], CSDres=csdres) - plt.plot(x, V, label=str(csdres)) - plt.legend() - #changing olny measure resolution - #measure_locations = np.arange(0,1,1000) - plt.figure('changing measure res') - for i in xrange(5): - measure_res= 20+i*50 - (x, V) = electrode_config([0.,1.], measure_res, csd_profile, [0.,1.], csd_seed, CSDres=1000) - plt.plot(x,V, ms = 0.5, label =str(measure_res)) - plt.legend() - plt.show() - return - -if __name__=='__main__': - total_ele = 30 - csd_seed = 11 - csd_profile = CSD.gauss_1d_mono - #test_calculating_potentials(csd_seed) - a = main_loop(csd_profile, csd_seed, total_ele) - #fig.savefig(os.path.join(plots_folder, save_as+'.png')) - #plt.clf() - #plt.close() diff --git a/tests/test_kCSD2D.py b/tests/test_kCSD2D.py deleted file mode 100644 index 5f4d7cf3..00000000 --- a/tests/test_kCSD2D.py +++ /dev/null @@ -1,241 +0,0 @@ -''' -This script is used to test Current Source Density Estimates, -using the kCSD method Jan et.al (2012) - -This script is in alpha phase. - -This was written by : -Chaitanya Chintaluri, -Laboratory of Neuroinformatics, -Nencki Institute of Exprimental Biology, Warsaw. -''' -import os -import time -import sys -sys.path.append('..') - -import numpy as np -from scipy.integrate import simps -from numpy import exp, linspace -import matplotlib.pyplot as plt -import matplotlib.cm as cm -from matplotlib.mlab import griddata - -import csd_profile as CSD -from KCSD2D import KCSD2D - -def generate_csd_2D(csd_profile, csd_seed, - start_x=0., end_x=1., - start_y=0., end_y=1., - res_x=50, res_y=50): - """ - Gives CSD profile at the requested spatial location, at 'res' resolution - """ - csd_x, csd_y = np.mgrid[start_x:end_x:np.complex(0,res_x), - start_y:end_y:np.complex(0,res_y)] - f = csd_profile(csd_x, csd_y, seed=csd_seed) - return csd_x, csd_y, f - -def grid(x, y, z, resX=100, resY=100): - """ - Convert 3 column data to matplotlib grid - """ - z = z.flatten() - xi = linspace(min(x), max(x), resX) - yi = linspace(min(y), max(y), resY) - zi = griddata(x, y, z, xi, yi, interp='linear') - return xi, yi, zi - -def generate_electrodes(xlims=[0.1,0.9], ylims=[0.1,0.9], res=5): - """ - Places electrodes in a square grid - """ - ele_x, ele_y = np.mgrid[xlims[0]:xlims[1]:np.complex(0,res), - ylims[0]:ylims[1]:np.complex(0,res)] - ele_x = ele_x.flatten() - ele_y = ele_y.flatten() - return ele_x, ele_y - -def make_test_plot(csd, est_csd): - print csd.shape, est_csd.shape, 'Shapes of csd' - fig = plt.figure(1) - ax1 = plt.subplot(111, aspect='equal') - yy1 = csd.flatten() - yy1.sort(axis=0) - im = ax1.plot(np.arange(len(yy1)), yy1, 'r') - ax1.hold(True) - yy2 = est_csd.flatten() - yy2.sort(axis=0) - im1 = ax1.plot(np.arange(len(yy2)), yy2, c='b') - plt.show() - return - -def make_plots(title, - t_csd_x, t_csd_y, true_csd, - ele_x, ele_y, pots, - k_csd_x, k_csd_y, est_csd): - """ - Shows 3 plots - 1_ true CSD generated based on the random seed given - 2_ interpolated LFT (NOT kCSD pot though), generated by simpsons rule integration - 3_ results from the kCSD 2D for the default values - """ - #True CSD - fig = plt.figure(figsize=(10,7)) - ax1 = plt.subplot(131, aspect='equal') - t_max = np.max(np.abs(true_csd)) - levels = np.linspace(-1*t_max, t_max, 16) - im = ax1.contourf(t_csd_x, t_csd_y, true_csd, levels=levels, cmap=cm.bwr_r) - ax1.set_title('TrueCSD') - cbar = plt.colorbar(im, orientation='horizontal') - #Potentials - v_max = np.max(np.abs(pots)) - levels_pot = np.linspace(-1*v_max, v_max, 16) - X,Y,Z = grid(ele_x, ele_y, pots) - ax2 = plt.subplot(132, aspect='equal') - im2 = plt.contourf(X, Y, Z, levels=levels_pot, cmap=cm.PRGn) - ax2.hold(True) - #im3 = plt.scatter(ele_x, ele_y, 30, pots, cmap=cm.PRGn) - im3 = plt.scatter(ele_x, ele_y, 5) - ax2.set_xlim([0.,1.]) - ax2.set_ylim([0.,1.]) - ax2.get_xaxis().set_visible(False) - ax2.get_yaxis().set_visible(False) - ax2.set_title('Pots, Ele_pos') - #cbar2 = plt.colorbar(im3, shrink=0.5) - cbar2 = plt.colorbar(im2, orientation='horizontal') - #KCSD - t_max = np.max(np.abs(est_csd[:,:,0])) - levels_kcsd = np.linspace(-1*t_max, t_max, 16) - ax3 = plt.subplot(133, aspect='equal') - im3 = ax3.contourf(k_csd_x, k_csd_y, est_csd[:,:,0], - levels=levels_kcsd, cmap=cm.bwr_r) - ax3.set_xlim([0.,1.]) - ax3.set_ylim([0.,1.]) - ax3.set_title('kCSD') - ax3.get_xaxis().set_visible(True) - ax3.get_yaxis().set_visible(False) - cbar = plt.colorbar(im3, orientation='horizontal') - fig.suptitle("Lambda,R,CV_Error,RMS_Error,Time = "+title) - #Showing/saving - plt.show() - #fig.savefig(os.path.join(plots_folder, save_as+'.png')) - #plt.clf() - #plt.close() - return - -def integrate_2D(x, y, xlim, ylim, csd, h, xlin, ylin, X, Y): - """ - X,Y - parts of meshgrid - Mihav's implementation - """ - Ny = ylin.shape[0] - m = np.sqrt((x - X)**2 + (y - Y)**2) # construct 2-D integrand - m[m < 0.0000001] = 0.0000001 # I increased acuracy - y = np.arcsinh(2*h / m) * csd # corrected - I = np.zeros(Ny) # do a 1-D integral over every row - for i in xrange(Ny): - I[i] = simps(y[:, i], ylin) # I changed the integral - F = simps(I, xlin) # then an integral over the result - return F - -def calculate_potential_2D(true_csd, ele_xx, ele_yy, csd_x, csd_y): - """ - For Mihav's implementation to compute the LFP generated - """ - xlin = csd_x[:,0] - ylin = csd_y[0,:] - xlims = [xlin[0], xlin[-1]] - ylims = [ylin[0], ylin[-1]] - sigma = 1.0 - h = 50. - pots = np.zeros(len(ele_xx)) - for ii in range(len(ele_xx)): - pots[ii] = integrate_2D(ele_xx[ii], ele_yy[ii], - xlims, ylims, true_csd, h, - xlin, ylin, csd_x, csd_y) - pots /= 2*np.pi*sigma - return pots - -def electrode_config(ele_lims, ele_res, true_csd, csd_x, csd_y): - """ - What is the configuration of electrode positions, between what and what positions - """ - #Potentials - ele_x_lims = ele_y_lims = ele_lims - ele_x, ele_y = generate_electrodes(ele_x_lims, ele_y_lims, ele_res) - pots = calculate_potential_2D(true_csd, ele_x, ele_y, csd_x, csd_y) - ele_pos = np.vstack((ele_x, ele_y)).T #Electrode configs - num_ele = ele_pos.shape[0] - print 'Number of electrodes:', num_ele - return ele_pos, pots - -def do_kcsd(ele_pos, pots, **params): - """ - Function that calls the KCSD2D module - """ - num_ele = len(ele_pos) - pots = pots.reshape(num_ele, 1) - k = KCSD2D(ele_pos, pots, **params) - #k.cross_validate(Rs=np.arange(0.01,0.10,0.02)) - k.cross_validate(Rs=np.array(0.08).reshape(1)) - est_csd = k.values('CSD') - return k, est_csd - -def main_loop(csd_profile, csd_seed, total_ele): - """ - Loop that decides the random number seed for the CSD profile, - electrode configurations and etc. - """ - csd_name = csd_profile.func_name - print 'Using sources %s - Seed: %d ' % (csd_name, csd_seed) - #TrueCSD - t_csd_x, t_csd_y, true_csd = generate_csd_2D(csd_profile, csd_seed, - start_x=0., end_x=1., - start_y=0., end_y=1., - res_x=100, res_y=100) - #Electrodes - ele_lims = [0.05, 0.95] #square grid, xy min,max limits - ele_res = int(np.sqrt(total_ele)) #resolution of electrode grid - ele_pos, pots = electrode_config(ele_lims, ele_res, true_csd, t_csd_x, t_csd_y) - ele_x = ele_pos[:, 0] - ele_y = ele_pos[:, 1] - - #kCSD estimation - gdX = 0.01 - gdY = 0.01 #resolution of CSD space - x_lims = [.0,1.] #CSD estimation place - y_lims = [.0,1.] - tic = time.time() #time it - k, est_csd = do_kcsd(ele_pos, pots, h=50., gdx=gdX, gdy=gdY, - xmin=x_lims[0], xmax=x_lims[1], - ymin=y_lims[0], ymax=y_lims[1]) - toc = time.time() - tic - - #RMS of estimation - gives estimate of how good the reconstruction was - chr_x, chr_y, test_csd = generate_csd_2D(csd_profile, csd_seed, - start_x=x_lims[0], end_x=x_lims[1], - start_y=y_lims[0], end_y=y_lims[1], - res_x=int((x_lims[1]-x_lims[0])/gdX), - res_y=int((y_lims[1]-y_lims[0])/gdY)) - rms = np.linalg.norm(abs(test_csd - est_csd[:,:,0])) - rms /= np.linalg.norm(test_csd) - #Plots - title = str(k.lambd)+', '+str(k.R)+', '+str(k.cv_error)+', '+str(rms)+', '+str(toc) - #save_as = csd_name+'_'+str(csd_seed)+'of'+str(total_ele) - make_plots(title, - t_csd_x, t_csd_y, true_csd, - ele_x, ele_y, pots, - k.estm_x, k.estm_y, est_csd) - return - -if __name__=='__main__': - total_ele = 81 - #Normal run - csd_seed = 15 #0-49 are small sources, 50-99 are large sources - if csd_seed < 50: - csd_profile = CSD.gauss_2d_small - else: - csd_seed %= 50 #Modulus of - csd_profile = CSD.gauss_2d_large - main_loop(csd_profile, csd_seed, total_ele) - #fig.savefig(os.path.join(plots_folder, save_as+'.png')) diff --git a/tests/test_kCSD3D.py b/tests/test_kCSD3D.py deleted file mode 100644 index 3b065b37..00000000 --- a/tests/test_kCSD3D.py +++ /dev/null @@ -1,359 +0,0 @@ -''' -This script is used to test Current Source Density Estimates, -using the kCSD method Jan et.al (2012) for the 3D case - -This script is in alpha phase. - -This was written by : -Chaitanya Chintaluri, -Laboratory of Neuroinformatics, -Nencki Institute of Exprimental Biology, Warsaw. -''' -import os -import time -import sys -sys.path.append('..') - -import numpy as np -from scipy.integrate import simps -from numpy import exp, linspace -import matplotlib.pyplot as plt -import matplotlib.cm as cm -from matplotlib import gridspec -from matplotlib.mlab import griddata -from scipy.spatial import distance - -try: - from joblib import Parallel, delayed - import multiprocessing - num_cores = multiprocessing.cpu_count()-1 - parallel_available = True -except ImportError: - parallel_available = False -import csd_profile as CSD -from KCSD3D import KCSD3D - -def generate_csd_3D(csd_profile, csd_seed, - start_x=0., end_x=1., - start_y=0., end_y=1., - start_z=0., end_z=1., - res_x=50, res_y=50, - res_z=50): - """ - Gives CSD profile at the requested spatial location, at 'res' resolution - """ - csd_x, csd_y, csd_z = np.mgrid[start_x:end_x:np.complex(0,res_x), - start_y:end_y:np.complex(0,res_y), - start_z:end_z:np.complex(0,res_z)] - f = csd_profile(csd_x, csd_y, csd_z, seed=csd_seed) - return csd_x, csd_y, csd_z, f - -def grid(x, y, z, resX=100, resY=100): - """ - Convert 3 column data to matplotlib grid - """ - x = x.flatten() - y = y.flatten() - z = z.flatten() - xi = linspace(min(x), max(x), resX) - yi = linspace(min(y), max(y), resY) - zi = griddata(x, y, z, xi, yi, interp='linear') - return xi, yi, zi - -def generate_electrodes(xlims=[0.1,0.9], ylims=[0.1,0.9], zlims=[0.1,0.9], res=5): - """ - Places electrodes in a square grid - """ - ele_x, ele_y, ele_z = np.mgrid[xlims[0]:xlims[1]:np.complex(0,res), - ylims[0]:ylims[1]:np.complex(0,res), - zlims[0]:zlims[1]:np.complex(0,res)] - ele_x = ele_x.flatten() - ele_y = ele_y.flatten() - ele_z = ele_z.flatten() - return ele_x, ele_y, ele_z - -def make_plots(fig_title, - t_csd_x, t_csd_y, t_csd_z, true_csd, - ele_x, ele_y, ele_z, pots, - k_csd_x, k_csd_y, k_csd_z, est_csd): - """ - Shows 3 plots - 1_ true CSD generated based on the random seed given - 2_ interpolated LFT (NOT kCSD pot though), generated by simpsons rule integration - 3_ results from the kCSD 2D for the default values - """ - fig = plt.figure(figsize=(10,16)) - #True CSD - z_steps = 5 - height_ratios = [1 for i in range(z_steps)] - height_ratios.append(0.1) - gs = gridspec.GridSpec(z_steps+1, 3, height_ratios=height_ratios) - t_max = np.max(np.abs(true_csd)) - levels = np.linspace(-1*t_max, t_max, 16) - ind_interest = np.mgrid[0:t_csd_z.shape[2]:np.complex(0,z_steps+2)] - ind_interest = np.array(ind_interest, dtype=np.int)[1:-1] - for ii, idx in enumerate(ind_interest): - ax = plt.subplot(gs[ii, 0]) - im = plt.contourf(t_csd_x[:,:,idx], t_csd_y[:,:,idx], true_csd[:,:,idx], - levels=levels, cmap=cm.bwr_r) - ax.get_xaxis().set_visible(False) - ax.get_yaxis().set_visible(False) - title = str(t_csd_z[:,:,idx][0][0])[:4] - ax.set_title(label=title, fontdict={'x':0.8, 'y':0.8}) - ax.set_aspect('equal') - cax = plt.subplot(gs[z_steps,0]) - cbar = plt.colorbar(im, cax=cax, orientation='horizontal') - cbar.set_ticks(levels[::2]) - cbar.set_ticklabels(np.around(levels[::2], decimals=2)) - #Potentials - v_max = np.max(np.abs(pots)) - levels_pot = np.linspace(-1*v_max, v_max, 16) - ele_res = int(np.ceil(len(pots)**(3**-1))) - ele_x = ele_x.reshape(ele_res, ele_res, ele_res) - ele_y = ele_y.reshape(ele_res, ele_res, ele_res) - ele_z = ele_z.reshape(ele_res, ele_res, ele_res) - pots = pots.reshape(ele_res, ele_res, ele_res) - for idx in range(min(5,ele_res)): - X,Y,Z = grid(ele_x[:,:,idx], ele_y[:,:,idx], pots[:,:,idx]) - ax = plt.subplot(gs[idx, 1]) - im = plt.contourf(X, Y, Z, levels=levels_pot, cmap=cm.PRGn) - ax.hold(True) - plt.scatter(ele_x[:,:,idx], ele_y[:,:,idx], 5) - ax.get_xaxis().set_visible(False) - ax.get_yaxis().set_visible(False) - title = str(ele_z[:,:,idx][0][0])[:4] - ax.set_title(label=title, fontdict={'x':0.8, 'y':0.8}) - ax.set_aspect('equal') - ax.set_xlim([0.,1.]) - ax.set_ylim([0.,1.]) - cax = plt.subplot(gs[z_steps,1]) - cbar2 = plt.colorbar(im, cax=cax, orientation='horizontal') - cbar2.set_ticks(levels_pot[::2]) - cbar2.set_ticklabels(np.around(levels_pot[::2], decimals=2)) - # #KCSD - t_max = np.max(np.abs(est_csd[:,:,:,0])) - levels_kcsd = np.linspace(-1*t_max, t_max, 16) - ind_interest = np.mgrid[0:k_csd_z.shape[2]:np.complex(0,z_steps+2)] - ind_interest = np.array(ind_interest, dtype=np.int)[1:-1] - for ii, idx in enumerate(ind_interest): - ax = plt.subplot(gs[ii, 2]) - im = plt.contourf(k_csd_x[:,:,idx], k_csd_y[:,:,idx], est_csd[:,:,idx,0], - levels=levels_kcsd, cmap=cm.bwr_r) - #im = plt.contourf(k_csd_x[:,:,idx], k_csd_y[:,:,idx], est_csd[:,:,idx,0], - # levels=levels, cmap=cm.bwr_r) - ax.get_xaxis().set_visible(False) - ax.get_yaxis().set_visible(False) - title = str(k_csd_z[:,:,idx][0][0])[:4] - ax.set_title(label=title, fontdict={'x':0.8, 'y':0.8}) - ax.set_aspect('equal') - cax = plt.subplot(gs[z_steps,2]) - cbar3 = plt.colorbar(im, cax=cax, orientation='horizontal') - cbar3.set_ticks(levels_kcsd[::2]) - #cbar3.set_ticks(levels[::2]) - cbar3.set_ticklabels(np.around(levels_kcsd[::2], decimals=2)) - #cbar3.set_ticklabels(np.around(levels[::2], decimals=2)) - fig.suptitle("Lambda,R,CV_Error,RMS_Error,Time = "+fig_title) - gs.tight_layout(fig, rect=[0, 0.03, 1, 0.95]) - # #Showing - #plt.tight_layout() - plt.show() - return - -def integrate_3D(x, y, z, xlim, ylim, zlim, csd, xlin, ylin, zlin, X, Y, Z): - """ - X,Y - parts of meshgrid - Mihav's implementation - """ - Nz = zlin.shape[0] - Ny = ylin.shape[0] - m = np.sqrt((x - X)**2 + (y - Y)**2 + (z - Z)**2) - m[m < 0.0000001] = 0.0000001 - z = csd / m - Iy = np.zeros(Ny) - for j in xrange(Ny): - Iz = np.zeros(Nz) - for i in xrange(Nz): - Iz[i] = simps(z[:,j,i], zlin) - Iy[j] = simps(Iz, ylin) - F = simps(Iy, xlin) - return F - -def calculate_potential_3D(true_csd, ele_xx, ele_yy, ele_zz, - csd_x, csd_y, csd_z): - """ - For Mihav's implementation to compute the LFP generated - """ - xlin = csd_x[:,0,0] - ylin = csd_y[0,:,0] - zlin = csd_z[0,0,:] - xlims = [xlin[0], xlin[-1]] - ylims = [ylin[0], ylin[-1]] - zlims = [zlin[0], zlin[-1]] - sigma = 1.0 - pots = np.zeros(len(ele_xx)) - tic = time.time() - for ii in range(len(ele_xx)): - pots[ii] = integrate_3D(ele_xx[ii], ele_yy[ii], ele_zz[ii], - xlims, ylims, zlims, true_csd, - xlin, ylin, zlin, - csd_x, csd_y, csd_z) - print 'Electrode:', ii - pots /= 4*np.pi*sigma - toc = time.time() - tic - print toc, 'Total time taken - series, sims' - return pots - -def calculate_potential_3D_parallel(true_csd, ele_xx, ele_yy, ele_zz, - csd_x, csd_y, csd_z): - """ - For Mihav's implementation to compute the LFP generated - """ - - xlin = csd_x[:,0,0] - ylin = csd_y[0,:,0] - zlin = csd_z[0,0,:] - xlims = [xlin[0], xlin[-1]] - ylims = [ylin[0], ylin[-1]] - zlims = [zlin[0], zlin[-1]] - sigma = 1.0 - #tic = time.time() - pots = Parallel(n_jobs=num_cores)(delayed(integrate_3D)(ele_xx[ii],ele_yy[ii],ele_zz[ii], - xlims, ylims, zlims, true_csd, - xlin, ylin, zlin, - csd_x, csd_y, csd_z) for ii in range(len(ele_xx))) - pots = np.array(pots) - pots /= 4*np.pi*sigma - #toc = time.time() - tic - #print toc, 'Total time taken - parallel, sims ' - return pots -1 -def electrode_config(ele_lims, ele_res, true_csd, csd_x, csd_y, csd_z): - """ - What is the configuration of electrode positions, between what and what positions - """ - #Potentials - ele_x_lims = ele_y_lims = ele_z_lims = ele_lims - ele_x, ele_y, ele_z = generate_electrodes(ele_x_lims, ele_y_lims, ele_z_lims, ele_res) - if parallel_available: - pots = calculate_potential_3D_parallel(true_csd, - ele_x, ele_y, ele_z, - csd_x, csd_y, csd_z) - else: - pots = calculate_potential_3D(true_csd, - ele_x, ele_y, ele_z, - csd_x, csd_y, csd_z) - ele_pos = np.vstack((ele_x, ele_y, ele_z)).T #Electrode configs - num_ele = ele_pos.shape[0] - print 'Number of electrodes:', num_ele - return ele_pos, pots - -def do_kcsd(ele_pos, pots, **params): - """ - Function that calls the KCSD3D module - """ - num_ele = len(ele_pos) - pots = pots.reshape(num_ele, 1) - k = KCSD3D(ele_pos, pots, **params) - #k.cross_validate(Rs=np.arange(0.2,0.4,0.02)) - #k.cross_validate(Rs=np.arange(0.02,0.27,0.01)) - k.cross_validate(Rs=np.array(0.31).reshape(1)) - est_csd = k.values('CSD') - return k, est_csd - -def main_loop(csd_profile, csd_seed, total_ele, num_init_srcs=1000): - """ - Loop that decides the random number seed for the CSD profile, - electrode configurations and etc. - """ - csd_name = csd_profile.func_name - print 'Using sources %s - Seed: %d ' % (csd_name, csd_seed) - - #TrueCSD - t_csd_x, t_csd_y, t_csd_z, true_csd = generate_csd_3D(csd_profile, csd_seed, - start_x=0., end_x=1., - start_y=0., end_y=1., - start_z=0., end_z=1., - res_x=100, res_y=100, - res_z=100) - - #Electrodes - ele_lims = [0.15, 0.85] #square grid, xy min,max limits - ele_res = int(np.ceil(total_ele**(3**-1))) #resolution of electrode grid - ele_pos, pots = electrode_config(ele_lims, ele_res, true_csd, t_csd_x, t_csd_y, t_csd_z) - ele_x = ele_pos[:, 0] - ele_y = ele_pos[:, 1] - ele_z = ele_pos[:, 2] - - #kCSD estimation - gdX = 0.05 - gdY = 0.05 - gdZ = 0.05 - x_lims = [.0,1.] #CSD estimation place - y_lims = [.0,1.] - z_lims = [.0,1.] - params = {'h':50., - 'gdX': gdX, 'gdY': gdY, 'gdZ': gdZ, - 'xmin': x_lims[0], 'xmax': x_lims[1], - 'ymin': y_lims[0], 'ymax': y_lims[1], - 'zmin': y_lims[0], 'zmax': y_lims[1], - 'ext': 0.0, 'n_srcs_init': num_init_srcs} - tic = time.time() #time it - k, est_csd = do_kcsd(ele_pos, pots, h=50., - gdx=gdX, gdy= gdY, gdz=gdZ, - xmin=x_lims[0], xmax=x_lims[1], - ymin=y_lims[0], ymax=y_lims[1], - zmin=z_lims[0], zmax=z_lims[1], - n_src_init=num_init_srcs, src_type='step') - toc = time.time() - tic - - #RMS of estimation - gives estimate of how good the reconstruction was - chr_x, chr_y, chr_z, test_csd = generate_csd_3D(csd_profile, csd_seed, - start_x=x_lims[0], end_x=x_lims[1], - start_y=y_lims[0], end_y=y_lims[1], - start_z=z_lims[0], end_z=z_lims[1], - res_x=int((x_lims[1]-x_lims[0])/gdX), - res_y=int((y_lims[1]-y_lims[0])/gdY), - res_z=int((z_lims[1]-z_lims[0])/gdZ)) - rms = np.linalg.norm(abs(test_csd - est_csd[:,:,:,0])) - rms /= np.linalg.norm(test_csd) - - #Plots - title = str(k.lambd)+','+str(k.R)+', '+str(k.cv_error)+', '+str(rms)+', '+str(toc) - save_as = csd_name+'_'+str(csd_seed)+'of'+str(total_ele) - #save_as = csd_name+'_'+str(num_init_srcs)+'_'+str(total_ele) - make_plots(title, - chr_x, chr_y, chr_z, test_csd, - ele_x, ele_y, ele_z, pots, - k.estm_x, k.estm_y, k.estm_z, est_csd) - #save - result_kcsd = [k.lambd, k.R, k.cv_error, rms, toc] - return est_csd, result_kcsd - -if __name__=='__main__': - total_ele = 125 - #Normal run - csd_seed = 54 #0-49 are small sources, 50-99 are large sources - csd_profile = CSD.gauss_3d_small - main_loop(csd_profile, csd_seed, total_ele, 1000) - #fig.savefig(os.path.join(plots_folder, save_as+'.png')) - - #Loop! - #for ii in range(100): - # main_loop(ii, total_ele) - - # import h5py - # for n_srcs in [18,20]: #range(12,18,2): - # num_init_srcs = n_srcs**3 - # for ele in range(3,8): - # total_ele = ele**3 - # h = h5py.File('Error_map_3D.h5', 'a') - # print 'Starting KCSD for src, ele:', num_init_srcs, total_ele - # est_csd, result = main_loop(csd_seed, total_ele, num_init_srcs) - # h[str(num_init_srcs)+'/'+str(total_ele)+'/est_csd'] = est_csd - # h[str(num_init_srcs)+'/'+str(total_ele)+'/kcsd_para'] = result - # print 'Finished KCSD for src, ele:', num_init_srcs, total_ele - # h.close() - - - - - diff --git a/utility_functions.py b/utility_functions.py index 66af3865..a1c4bbf2 100644 --- a/utility_functions.py +++ b/utility_functions.py @@ -10,6 +10,8 @@ Laboratory of Neuroinformatics, Nencki Institute of Experimental Biology, Warsaw. """ +from __future__ import division + import numpy as np from scipy import interpolate @@ -87,7 +89,8 @@ def distribute_srcs_2D(X, Y, n_src, ext_x, ext_y, R_init): X_src, Y_src = np.mgrid[(np.min(X) - ext_x_n):(np.max(X) + ext_x_n):np.complex(0,nx), (np.min(Y) - ext_y_n):(np.max(Y) + ext_y_n):np.complex(0,ny)] d = round(R_init/ds) - R = d * ds + #R = d * ds + R = R_init return X_src, Y_src, R def get_src_params_2D(Lx, Ly, n_src): @@ -166,7 +169,7 @@ def distribute_srcs_3D(X, Y, Z, n_src, ext_x, ext_y, ext_z, R_init): (np.min(Y) - ext_y_n):(np.max(Y) + ext_y_n):np.complex(0,ny), (np.min(Z) - ext_z_n):(np.max(Z) + ext_z_n):np.complex(0,nz)] d = np.round(R_init/ds) - R = d * ds + R = R_init return (X_src, Y_src, Z_src, R) def get_src_params_3D(Lx, Ly, Lz, n_src):