diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..4ca3ddd --- /dev/null +++ b/.gitignore @@ -0,0 +1,4 @@ +*.c +*.o +*.pyc +*.dll diff --git a/Optical Stimulation.py b/Optical Stimulation.py index 71b149e..20b3091 100644 --- a/Optical Stimulation.py +++ b/Optical Stimulation.py @@ -59,8 +59,20 @@ from neuron import h h.load_file("stdrun.hoc") -h.cvode_active(1) - +h.load_file("nrngui.hoc") +h.cvode_active(0) + +def check(): + mt = h.MechanismType(1) + mname = h.ref('') + for i in range(mt.count()): + mt.select(i) + mt.selected(mname) + if mname[0] == 'ostim': + return "Good" + return "Nope" + +# print("1st: " + check()) # # Instantiate cell, stimulator and simulation classes @@ -75,6 +87,7 @@ cell = Hu() optrode = Optrode(h.soma) sim = Sim(cell,optrode,output_filename='csv/distance_threshold.csv') +# print("2: " + check()) # @@ -104,17 +117,19 @@ # -%pylab inline +# %pylab inline # # Use matplotlib to display the distance versus threshold results. # +# print("3: " + check()) from matplotlib import pyplot from classes import Data from functions import make_legend +# print("4: " + check()) # Simulation ds = Data('csv/distance_threshold.orig.csv') @@ -126,7 +141,7 @@ for f in [0.1,0.2,0.4]: ds.set_slice(ds.data['Fiber Optic Diameter (mm)']==f) pyplot.semilogy(ds.slice['Distance (um)'],ds.slice['Threshold (W/cm2)'], - styles.next(), + styles.__next__(), label='%.1f' % f) pyplot.plot([1400], [38.0], '*', color='0.5',label='Aravanis (0.2)') # Aravanis data point pyplot.xlabel('Distance (um)') @@ -142,6 +157,7 @@ # 2D visualization of the neuron geometry can be done in matplotlib using the cell's *plot* method. The *plot* method generates the cell's trajectory with a color overlay determined by an arbitrary input function that is run for each NEURON section. In this example, the input function returns the section's intracellular voltage. # +# print("5: " + check()) from matplotlib import pyplot plot_func = lambda sec:sec.v diff --git a/classes.py b/classes.py index d688ae2..795fae9 100644 --- a/classes.py +++ b/classes.py @@ -1,4 +1,6 @@ from neuron import h +h.load_file("stdrun.hoc") +h.load_file("nrngui.hoc") from functions import * class Hu(object): @@ -61,7 +63,7 @@ def _ap_counters(self): try: sec = h.node[h.axonnodes.__int__()-2] except AttributeError: - print "No node compartments!" + print("No node compartments!") return 0 apc = h.APCount(0.5, sec=sec) apc.thresh = 0 # mV @@ -79,7 +81,7 @@ def _ap_counters(self): self.apc_times.append(apc_time) else: if h.number_of_apc>2: - raise ValueError,"Too many apc counters; only 1 or 2 allowed" + raise ValueError("Too many apc counters; only 1 or 2 allowed") def intensity(self,sec): return sec.irradiance_chanrhod * sec.Tx_chanrhod def photon_flux(self, sec): @@ -95,13 +97,13 @@ def build_tree(self, func,segfunc=False): """ func must act on a neuron section """ from numpy import array - print "-"*100 + print("-"*100) def append_data(sec, xyzdv, parent_id, connections,func,segfunc): """ Append data to xyzdv """ if not segfunc: v=func(sec) n = int(h.n3d(sec=sec)) - for ii in xrange(1, n): + for ii in range(1, n): x = h.x3d(ii,sec=sec) y = h.y3d(ii,sec=sec) z = h.z3d(ii,sec=sec) @@ -154,7 +156,7 @@ def move(self, xyz, move_mlab=False): tree = h.SectionList() tree.wholetree(sec=self.root) for sec in tree: - for ii in xrange(h.n3d(sec=sec).__int__()): + for ii in range(h.n3d(sec=sec).__int__()): x=h.x3d(ii,sec=sec) y=h.y3d(ii,sec=sec) z=h.z3d(ii,sec=sec) @@ -162,7 +164,7 @@ def move(self, xyz, move_mlab=False): h.pt3dchange(ii,x+float(xyz[0]),y+float(xyz[1]),z+float(xyz[2]),d) def retrieve_coordinates(self, sec): xyzds = [] - for ii in xrange(int(h.n3d(sec=sec))): + for ii in range(int(h.n3d(sec=sec))): xyzds.append([h.x3d(ii,sec=sec), h.y3d(ii,sec=sec), h.z3d(ii,sec=sec), @@ -242,7 +244,7 @@ def plot(self, func, scaling = 1, segfunc=False, clim=None,cmap=None): edges = self.connections diam = self.xyzdv[:,3] data = self.xyzdv[:,4] - print "DATA RANGE: ",data.min(),data.max() + print("DATA RANGE: ",data.min(),data.max()) # Define colors if not cmap: from matplotlib.cm import jet as cmap @@ -327,7 +329,7 @@ def get_apical_tuft(self): at the branch point """ secs=[] - for ii in xrange(23,len(h.dend11)): + for ii in range(23,len(h.dend11)): secs.append(h.dend11[ii]) return secs def get_apical_shaft(self): @@ -409,7 +411,7 @@ def set_required_aps(self,stimulator,additional_aps=1): assert self.apc_times,'No action potential counters' for apct in self.apc_times: h.required_aps=max((h.required_aps,len(apct))) - print "** NO STIM APS: %d; Goal APS: %d **" % (h.required_aps, + print("** NO STIM APS: %d; Goal APS: %d **" % h.required_aps, h.required_aps+additional_aps) h.required_aps += additional_aps # usually require one additional ap stimulator.amplitude=initial_amplitude @@ -509,7 +511,7 @@ def interpxyz(self): yy = h.Vector(nn) zz = h.Vector(nn) length = h.Vector(nn) - for ii in xrange(nn): + for ii in range(nn): xx.x[ii] = h.x3d(ii,sec=sec) yy.x[ii] = h.y3d(ii,sec=sec) zz.x[ii] = h.z3d(ii,sec=sec) @@ -651,7 +653,7 @@ def rotate_optrode(self, angle, distance, sec, axis_defining_plane='z'): elif axis_defining_plane=='z': uvw = cross(xyz1, array([0,0,1])) else: - raise ValueError,'No such plane: %s' % axis_defining_plane + raise ValueError('No such plane: %s' % axis_defining_plane) uvw /= sqrt(sum(uvw**2)) # Normalize # Rotate the optrode around the vector @@ -858,10 +860,10 @@ def get_distance(self, sec): optrode_output) elif isinstance(sec,nrn.Segment): seg = sec - print find_section_coordinates(seg.sec),seg.x - raise StandardError,"Not yet implemented" + print(find_section_coordinates(seg.sec),seg.x) + raise StandardError("Not yet implemented") else: - raise TypeError, "Wrong type: %s" % type(sec) + raise TypeError( "Wrong type: %s" % type(sec)) def set_distance(self, sec, z_distance): """ Set optrode a certain distance in z-direction below the given section directed upwards """ @@ -910,13 +912,13 @@ def get_xyz(self): def set_xyz(self,val): from numpy import array val = array(val).astype(float) - #print "SHAPE: ",val.shape + #print("SHAPE: ",val.shape) if val.shape == (3,2): self.set_position(val[0,:],val[1,:],val[2,:]) elif val.shape == (2,3): self.set_position(val[:,0],val[:,1],val[:,2]) else: - raise ValueError, "Trying to set xyz with inappropriate array %s" % str(val) + raise ValueError( "Trying to set xyz with inappropriate array %s" % str(val)) def get_length(self): return self.sec.L def set_length(self,val): @@ -1039,12 +1041,12 @@ def find_threshold(self,upper_limit=1e4,error_threshold=1e-3,verbose=True,additi self.cell.set_tstop(h.tstop,self.stimulator,additional_aps) supra=upper_limit/0.9 sub=0 - if verbose:print "+ _SubT____, Amplitude, _SupraT__, _Error___ +" + if verbose:print("+ _SubT____, Amplitude, _SupraT__, _Error___ +") while (supra-sub)/self.stimulator.amplitude > error_threshold: # while error larger than threshold h.run() if self.cell.response:response='+' else:response='-' - if verbose:print "%s %0.3e, %0.3e, %0.3e, %0.3e : %s" % (response, + if verbose:print("%s %0.3e, %0.3e, %0.3e, %0.3e : %s" % response, sub, self.stimulator.amplitude, supra, @@ -1056,7 +1058,7 @@ def find_threshold(self,upper_limit=1e4,error_threshold=1e-3,verbose=True,additi else: if self.stimulator.amplitude >= upper_limit: # Probably not going to reach upper - if verbose:print "** Upper threshold reached **" + if verbose:print("** Upper threshold reached **") self.stimulator.amplitude=0 h.run() return float('nan') @@ -1066,7 +1068,7 @@ def find_threshold(self,upper_limit=1e4,error_threshold=1e-3,verbose=True,additi if not self.cell.response: self.stimulator.amplitude=supra h.run() - if verbose:print "** THRESHOLD: %g **" % self.stimulator.amplitude + if verbose:print("** THRESHOLD: %g **" % self.stimulator.amplitude) return self.stimulator.amplitude class Data(object): def __init__(self, filename, seperator=','): @@ -1077,7 +1079,8 @@ def __init__(self, filename, seperator=','): line_type = self.data_or_header(data_list) if line_type =='data': for ii,x in enumerate(data_list): - if not data.has_key(ii):data[ii]=[] + if not ii in data.keys(): + data[ii]=[] try: data[ii].append(float(x)) except: @@ -1085,7 +1088,7 @@ def __init__(self, filename, seperator=','): elif line_type == 'header': headers=data_list if not data: - raise ValueError,"Empty data file, or not seperated by %s" % seperator + raise ValueError("Empty data file, or not seperated by %s" % seperator) for k in data.keys(): data[k] = array(data[k]) self.filename=filename diff --git a/csv/distance_threshold.csv b/csv/distance_threshold.csv new file mode 100644 index 0000000..e69de29 diff --git a/functions.py b/functions.py index 9ee1c76..64756f4 100644 --- a/functions.py +++ b/functions.py @@ -1,10 +1,14 @@ +import sys + def find_section_coordinates(sec): """ Determine xyz coordinates for a given section """ from numpy import array from neuron import h + h.load_file("stdrun.hoc") + h.load_file("nrngui.hoc") x=[];y=[];z=[]; n=h.n3d(sec=sec).__int__() - for ii in xrange(n): + for ii in range(n): x.append(h.x3d(ii,sec=sec)) y.append(h.y3d(ii,sec=sec)) z.append(h.z3d(ii,sec=sec)) @@ -12,10 +16,12 @@ def find_section_coordinates(sec): def find_mean_section_coordinates(sec): """ Determine average coordinate for a given section """ from neuron import h + h.load_file("stdrun.hoc") + h.load_file("nrngui.hoc") from numpy import array, mean x=[];y=[];z=[]; n=h.n3d(sec=sec).__int__() - for ii in xrange(n): + for ii in range(n): x.append(h.x3d(ii,sec=sec)) y.append(h.y3d(ii,sec=sec)) z.append(h.z3d(ii,sec=sec)) @@ -108,7 +114,7 @@ def rotate_coordinates(X,Y,Z,v,axis='z'): v = array(v) return xyz[:,0].reshape(X.shape),xyz[:,1].reshape(Y.shape),xyz[:,2].reshape(Z.shape),v else: - raise StandardError,'Not implemented' + sys.exit('Not implemented') def find_cylindrical_coords(X,Y,Z,v): " Align v with positive z axis " from numpy import sqrt diff --git a/mwe.py b/mwe.py new file mode 100644 index 0000000..7fa8a1c --- /dev/null +++ b/mwe.py @@ -0,0 +1,23 @@ +# -*- coding: utf-8 -*- +""" +Created on Wed Jul 27 09:36:33 2022 + +@author: maria +""" +from neuron import h +h.load_file("stdlib.hoc") +h.load_file("nrngui.hoc") + +def check(): + mt = h.MechanismType(1) + mname = h.ref('') + for i in range(mt.count()): + mt.select(i) + mt.selected(mname) + if mname[0] == 'ostim': + return "Good" + return "Nope" + +soma = h.Section(name='soma') + +testpp = h.ostim(0.5) \ No newline at end of file