python 3

.gitignore

@@ -0,0 +1,4 @@
+*.c
+*.o
+*.pyc
+*.dll

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())
 # <headingcell level=1>
 
 # 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())
 
 # <headingcell level=1>
 
@@ -104,17 +117,19 @@
 
 # <codecell>
 
-%pylab inline
+# %pylab inline
 
 # <markdowncell>
 
 # Use matplotlib to display the distance versus threshold results.
 
 # <codecell>
+# 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.
 
 # <codecell>
+# print("5: " + check())
 
 from matplotlib import pyplot
 plot_func = lambda sec:sec.v

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,15 +156,15 @@ 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)
                 d=h.diam3d(ii,sec=sec)
                 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,15 +1079,16 @@ 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:
                         data[ii].append(x.strip())
             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

functions.py

@@ -1,21 +1,27 @@
+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))
     return array([array(x),array(y),array(z)])
 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

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)