#!/usr/bin/python
# cellpars2manyhocs.pl by Ted Ballou, Version 2. 2/7/2011.
# The purpose of this script is to create a set of cell models that
# incorporate mechanism and geometric parameters and global values that are
# extrapolated between two extremes; for convenience these two extremes are
# identified as "low threshold" and "high threshold." Each desired parameter
# or global along with their associated extreme values must be represented by
# a line in the "csv" input file according to the format described below. The
# number of cells to be generated is specified by the user, with a default
# value of 11. For each generated cell a directory is created and a "hoc" file
# that contains the hoc statements to assign the interpolated values
# appropriate to that cell is placed in its directory. The directory and the
# hoc file names are based on the csv file name, with the cell number
# appended.
#
# There are two applications for this at present:
#
# 1. Generate a set of cells for sensitivity analysis, to illustrate the
# effects of gradually changing various parameters.
#
# 2. Generate sets of cells with gradually changing properties in
# order to simulate populations of neurons with incrementally varying
# properties; for example, spinal motoneurons with continuously varying
# properties representing slow and fatiguing muscle motor units.
# This code takes a comma-delimited file (with extension "csv" standing for
# "commma-separated-values") with four required fields and one optional field:
#
# 1. NEURON segment name, OR the string "Global", OR "*" (see below)
# 2. NEURON mechanism, OR the global parameter name
# 3. low-threshold value ("lt") for the parameter or global
# 4. high-threshold value ("ht") for the parameter or global
# 5. (optional) nonlinear "indicator"
#
# If "*" is in the first field then what ever is in the second field is
# transferred directly to the output cell files so you can place hoc
# statements in the csv file to be transferred into the cell files.
#
# On execution the code will prompt for the "Number of cells" to generate,
# which will be assigned to the variable $Ncells. The default value for
# $Ncells is 11, resulting in cell numbers 0, 1, 2, ..., 10.
#
# The code will generate hoc files for $Ncells cells, containing an assignment
# line for each parameter or global, with assignment values extrapolated
# between the two extremes. The low threshold value will be assigned to cell
# number 0, and the high threshold value will be assigned to cell number
# ($Ncells-1); in the default case $Ncells is ll and the ht value is assigned
# to cell number 10. For cell numbers 1 through $Ncells-2, an extrapolated
# value will be assigned based on the nonlinear indicator; if this 5th field
# is absent, the default nonlinear indicator of 0 will be used, which assigns
# a linear extrapolation between the lt and ht values.
#
# For positive nonlinear indicator values, the population will be skewed with
# most of the cells having extrapolated values less than the mean of lt and
# ht; this "$skew" is caluclated as the fraction of the range from lt to ht,
# taken to the power of (1+$nonlin_ind). For negative nonlinear indicator
# values, the population will be skewed with most of the cells having
# extrapolated values more than the mean of lt and ht; this "$skew" is
# calculated as (1 minus the fraction of the range from lt to ht), taken to
# the power of (1-$nonlin_ind), subtracted from 1. These formulas result in a
# nicely symmetric chart for positive, negative, and fractional values of
# $nonlin_ind. See subroutine intrpolate() for the details. Also see the
# documentation and examples in the excel file doc/NonlinearityScaling.xls.
#
# EXAMPLE: say the parameter soma.diam has lt=45 and ht=65. The following
# table shows sample interpolations for nonlinear indicator values of 0 (the
# default), 1, and 2. Note that most of the values for the last two columns
# are LESS than the mean of 45 and 65 (55).
#
# Indicator 0 1 2
# Cell 0 45.00 45.00 45.00
# Cell 1 47.00 45.20 45.02
# Cell 2 49.00 45.80 45.16
# Cell 3 51.00 46.80 45.54
# Cell 4 53.00 48.20 46.28
# Cell 5 55.00 50.00 47.50
# Cell 6 57.00 52.20 49.32
# Cell 7 59.00 54.80 51.86
# Cell 8 61.00 57.80 55.24
# Cell 9 63.00 61.20 59.58
# Cell 10 65.00 65.00 65.00
#
# I find it convenient to create the input csv file as a table using excel or
# Open Office calc "OOcalc" and then "SaveAs" a comma-delimited csv file, but
# any text editor can be used to create the csv file.
# File structure:
# First line: ,,lt header, ht header (this line is not used by the script)
# The remaining lines have four OR five fields:
# section, mechanism, lt value, ht value, (nonlin_ind)
# OR
# designation as "Global", global parameter name, lt value, ht value, (nonlin_ind)
###############
# CODE BEGINS
###############
# Invoke perl module to get current working directory
import os # has os.getcwd()
import sys # for args passed on command line
Ncells = 11 # the total number of cells, including lt and ht
cur_cell_no = 0 # globalize this variable, "current cell #"
skew = 0.0 # a global variable
# optionally put csv file on command line, or as reply to prompt when program runs
# [Note: it's easier on command line since <TAB> provides file expansion]
global nonlin_ind
nonlin_ind=0.0 # global variable
fields=[]
print "%%%%%%%% to abort type <ctrl-c> %%%%%%%%\n"
if (len(sys.argv)!=2):
myf=raw_input("Enter the name of the ht/lt parameters \".csv\" file: ")
else:
myf = sys.argv[1]
if (len(myf)>1):
try:
fid=open(myf,"r")
except:
print("Cannot open file " + myf + " for reading" )
print("Invoke with csv file in current directory")
sys.exit(1)
else:
print "No txt file supplied as input."
exit(1)
print("Default number of cells is "+str(Ncells))
test=raw_input("If this is ok hit <return>; otherwise type the desired number: ")
if len(test):
test_num=eval(test)
if test_num>0:
Ncells = float(test_num)
else:
print "non-numeric or zero input; using default Ncells value = ", Ncells
#Generate a base name for generated sensitivity directories & files
# mydir array is the list of directories in the full path to current dir
absolute_dir=os.getcwd()
mydir=absolute_dir.split('/')
# get the name of the parent directory: last member of the array split by '/'
basen = mydir[-1].strip()
# strip it down to only the characters preceding the first underscore
basen=basen.split('_')[0]
# Read the data into an array
datalist = fid.readlines()
fid.close()
# scrub the quotes generated by the conversion of spreadsheet to cvs
dl=[]
for line in datalist:
dl.append(line.replace('"', '').strip())
# The first line of the data file contains the low threshold and high
# threshold headers
# [these variable are not currently used]
fields = dl[0].split(',') # first line of the csv file
ltidx = 2
lthdr = fields[ltidx]
htidx = 3
hthdr = fields[htidx]
print("lthdr = "+lthdr+", hthdr = " + hthdr)
# scrub the first line
del dl[0]
# The remaining lines in this file have four fields:
# section, mechanism, lt value, ht value
# OR
# designation as "Global", global parameter name, lt value, ht value
# now treat each segment mechanism, generating a cell array between the
# lt and ht values
# Use values linearly extrapolated between the lt and ht cells, and
# generate hoc files for Ncells cells, each in its own directory
##########################################
# functions
#######################
#################
# isnumeric() subroutine
#
# arg 1: string to be tested for being a numeric value
#
# Utility to validate argument as numeric
# from
# http://stackoverflow.com/questions/354038/how-do-i-check-if-a-string-is-a-number-float-in-python
def isnumeric(s):
try:
float(s)
return True
except ValueError:
return False
########################################
# getintrp subroutine
#
# arg 1: lt variable value
# arg 2: ht variable value
#
# interpolates lt and ht values; uses global variables Ncells, skew, and
# cur_cell_no. lt and ht non-ranging parameter values are passed as
# arguments
def getintrp(start, end):
global skew
#print "start="+repr(start)+",end="+repr(end)+", skew="+repr(skew)
if ((not isnumeric(start)) or (not isnumeric(end))):
print("getintrp() bad args: "+ start +" and "+end)
sys.exit("getintrp() requires numeric arguments")
# The following line calculates the interpolation for the cur_cell_no'th
# interval between the start and end values. This is a nonlinear mechanism
# for the nonlin_ind being non-zero: skewed towards lt for positive values,
# and skewed towards ht for negative values. See calculation of skew above.
tmp=eval(start) + skew*(eval(end)-eval(start))
#print "returning "+repr(tmp)
return tmp
##########################
# extract_range subroutine
#
# arg 1: colon-delimited lt ranging variable
# arg 2: colon-delimited ht ranging variable
#
# Returns: interpolated colon-delimited ranging variable
#
# This function extracts the start and end range values for lt and for ht,
# uses them to generate interpolated start and end range values, and returns
# the interpolated ranging variable
def extract_range(ltarg, htarg):
# extract values from colon-delimited fields
ltarg = ltarg.split(":")
htarg = htarg.split(":")
# and construct the interpolated ranging variable that is returned
#print "debugging:"+repr(ltarg)+"|"+repr(htarg)
#print len(ltarg)
#print len(htarg)
if len(ltarg)==1:
return repr(getintrp(ltarg[0], htarg[0]))
else:
return repr(getintrp(ltarg[0], htarg[0]))+":"+repr(getintrp(ltarg[1], htarg[1]))
# intrpolate subroutine
#
# No arguments for this function; the cell number is accessed via the global
# variable cur_cell_no, and the low threshold (lt) and high threshold (ht)
# values are accessed via the glabal array elements fields[2] and fields[3].
#
# The "nonlinear indicator" nonlin_ind value is used to skew the interpolation,
# where skew is a value between 0 and 1, according to the formula below.
#
# Single values are interpolated between lt and ht.
# For colon-delimited ranges ("ranging variables") the range limits are
# extracted for low and high thresholds, and extrapolations between both
# elements are used, to generate the desired colon-delimited result, in the
# extract_range() subroutine.
#
# The cell number is accessed via the global variable cur_cell_no
def intrpolate():
global fields
global skew
global nonlin_ind
# Generate the "skew" for nonlinear interpolations. This converts the
# nonlinear indicator nonlin_ind from a single number to a list of values,
# one for each of the generated cells.
fraction = float(cur_cell_no)/(Ncells-1) # fraction of the range from start to
#print "fraction = "+repr(fraction)+" for cell num " +repr(cur_cell_no) # end
if (nonlin_ind >= 0):
skew = fraction**(1+nonlin_ind)
else:
skew = 1-(1-fraction)**(1-nonlin_ind)
#print "*** skew = ", skew;
# Generate linearly interpolated value between lt and ht
#print "len(fields)="+repr(len(fields))
#print "fields = "+repr(fields)
#print "fields[2]="+repr(fields[2])
if (":" not in fields[2]): # here if lt is NOT a ranging variable
if (":" in fields[3]):
sys.exit("high threshold (ht) is ranging but low threshold (lt) is not in "+line)
# just a number, linearly interpolated based on the distance cur_cell_no
# is between 0 and Ncells-1
getintrp(fields[2], fields[3]);
else:
# Here if lt IS a ranging variable
if ":" not in fields[3]:
sys.exit("low threshold (lt) is ranging but high threshold (ht) is not")
# generate the interpolated ranging value (two colon-delimited numbers)
return extract_range(fields[2], fields[3])
# The make_cell() subroutine: create the interpolated hoc files.
# make_cell() argument is the cell # being generated
def make_cell(i):
global fields
global cur_cell_no
global nonlin_ind
#print "in make_cell("+str(i)+")"
# Assign this global variable to the argument passed by caller
cur_cell_no = i # shift; # global variable used throughout the subroutines
# Generate the output directory and file for the cell
myfile = basen+"_"+str(cur_cell_no)
# Continue even if the directory was previously present
if not os.path.exists(myfile):
os.makedirs(myfile)
# Create the new hoc file; overwrite if already present
MYF=open(myfile+"/"+myfile+".hoc","w")
# calculate and write the values interpolated between lt and ht
for line in dl:
line=line.strip()
# Populate the fields[] array with the current line from the csv file
fields = line.split(',')
#print "Processing line = "+repr(line)+" with fields = "+repr(fields)
if fields[0]=='*': # asterisk transfers next field into program
hocstatement=fields[1]
MYF.write(hocstatement+"\n")
continue
else:
if len(fields)>3:
# Discard lines that have the lt or ht value field blank
if (not len(fields[2])) or (not len(fields[3])):
continue
if len(fields)<4:
# Discard lines that have less than 4 fields
if not cur_cell_no: # limit error print to one per line
print "Less than four fields in discarded line="+repr(line)
continue
# fields 2 and 3 are low and high threshold values; if they are BOTH
# colon-delimited ranging variables then field 1 (mechanism) MUST ALSO be
# a colon-delimited ranging variable, ie Diam(10:4) intrp will then be
# the interpolated colon-delimited ranging values
# The nonlinear indicator field is set to 0 (NO nonlinearity) by default,
# if 5th field is absent.
if len(fields)>=5:
if not len(fields[4]):
nonlin_ind = 0.0
else:
# any rational number is ok
nonlin_ind = float(fields[4])
# print "setting nonlin_ind ="+repr(nonlin_ind)
else:
nonlin_ind = 0.0
if fields[0]=='':
hocstatement=""
else:
intrp = intrpolate()
#print "intrp = "+intrp
# generate hoc statement for:
# Global variables
if ("Global" in fields[0]):
hocstatement = fields[1]+" = "+intrp
# OR "forall" directive
elif ("forall" in fields[0]):
hocstatement = fields[0]+"{"+fields[1]+" = "+intrp+"}"
# OR (default) section.mechanism=value
elif len(fields)>1:
hocstatement = fields[0]+"."+fields[1]+" = "+intrp
# save in the hoc file
MYF.write(hocstatement+"\n")
MYF.close()
# end of make_cell
for i in range(int(Ncells)):
make_cell(i) # Create each hoc file in its own directory
print("SUCCESS: created directories and hoc files for "+str(Ncells)+" cells")
############################
# END of MAIN
############################
