Python/Neuron model runs on Windows but does not on Linux.(h.finitialize)

Post Reply
Cbswimmer23
Posts: 1
Joined: Thu Mar 30, 2023 1:40 pm

Python/Neuron model runs on Windows but does not on Linux.(h.finitialize)

Post by Cbswimmer23 »

I have a python/Neuron model which runs on windows but will not run on our server. No error flags are presented and we are using the multiprocessing.pool command. There seems to be an issue with the h.finitialize command where it has an issue with its children processes. Any help would be greatly appreciated.


Code: Select all

import os
from neuron import h
import matplotlib.pyplot as plt
import numpy as np
import multiprocessing
import time
from scipy import interpolate
import h5py

#------------------------------------------------------------------------
#Load Motor Neuron template with incorporated mechanisms. (Kim 2017)
h.load_file('kim.hoc')
#Load standard NEURON run controls.
h.load_file("stdrun.hoc")
#------------------------------------------------------------------------


#Function Definitions----------------------------------------------------
def heav(x):#Needed function to calculate current of stimulation source.
    if x<0:
        return 0
    else:
        return 1

def Soma_Diam_Change(NMU):#Determine Motor Pool Motor Neuron size. (Volk 2021)
    num_MU = NMU
    RR = 50.9
    a =np.log(RR)/num_MU
    RTE=np.empty(num_MU)
    for j in range(1,num_MU+1):
        RTE[j-1]=np.exp(a*j)

    diam_ori =48.8
    diameters = RTE + diam_ori

    return diameters#Return list of soma sizes in order of Motor Units.

def MU_List(N):#Generate Python Object from the Motor Neuron template.
    MU_list = []

    MU=h.Kim()
    MU_list.append(MU)

    return MU_list#Returns a list with a Motor Unit Py Object

def Soma_List(MU_list):#Isolate the Soma Segment of each Motor Unit Object
    soma_list=[]#Initialize list
    dend_list=[]
    for section in h.allsec(MU_list):#For each Soma section on each Motor Neuron Object add them to a list
            if h.issection('.*so.*', sec=section):
                soma_list.append(section)
            if h.issection('.*dend.*', sec=section):
                dend_list.append(section)

    return soma_list, dend_list#Return list of Somas

def Change_Dia(mu_list, soma_list,dend_list,N,NMU):#Change Diameters of the Somas of each Motor Neuron Object
    NMU=NMU #Number of Motor Neurons desired
    diameters = Soma_Diam_Change(NMU) #Calculate Soma Diameters
    
    for i in range(len(mu_list)):#For each Motor Neuron Object in list
        for j in range(soma_list[i].n3d()):#For each section of each Soma of each Motor Neuron Object
            if j == 0:
                mu_list[i].soma.pt3dchange(j, 0, 0, 0, diameters[N])
            else:
                mu_list[i].soma.pt3dchange(j, diameters[N], 0, 0, diameters[N])

    for i in range(len(mu_list)):
        mu_list[i].soma.gnafbar_Naf = ((mu_list[i].soma.gnafbar_Naf)-np.random.normal(0.0,0.1))
        # mu_list[i].soma.ena = (mu_list[i].soma.ena+np.random.normal(0,5))
        # mu_list[i].soma.g_pas = (mu_list[i].soma.g_pas+np.random.normal(0,1e-6))  
        # mu_list[i].soma.ek= (mu_list[i].soma.ek+np.random.normal(0,5))
        mu_list[i].hillock.gnafbar_Naf =  ((mu_list[i].hillock.gnafbar_Naf*np.random.normal(0.5,0.1)))
        mu_list[i].hillock.gkdrbar_KDr =  ((mu_list[i].hillock.gkdrbar_KDr+np.random.normal(0,0.1)))
        # for j in range(len(dend_list)):
        #     mu_list[i].dend[j].g_pas =mu_list[i].dend[j].g_pas+np.random.normal(0,1e-6)
        #     mu_list[i].dend[j].e_pas =mu_list[i].dend[j].e_pas+np.random.normal(0,5)
  
    # for i in range(len(mu_list)):
    #     for j in range(len(mu_list[i].iCaL)):
    #         if mu_list[i].iCaL[j] is None:
    #           ww=0  
    #         else:    
    #             mu_list[i].iCaL[j].gcalbar =mu_list[i].iCaL[j].gcalbar + np.random.normal(0,0.0008)*i




    return mu_list#Return list of updated Motor Neuron Objects

def Record(stimulus, MU_list):#Record simulation characteristics of interest
    #Initialize lists
    t = []
    v = []
    amp = []
    mgi = []
    ina = []
    # h_naf = []

    for i in range(len(MU_list)):#For each Motor Neuron Object record Time, Soma Voltage, Stimulus Current, and Muscle Unit force
        t.append(h.Vector().record(h._ref_t))
        v.append(h.Vector().record(MU_list[i].soma(0.5)._ref_v))
        amp.append(h.Vector().record(stimulus[i]._ref_amp))
        mgi.append(h.Vector().record(MU_list[i].muscle_unit( 0.5 )._ref_mgi))
        cell_t, cell_id = Spike_Detector(MU_list)
        ina.append(h.Vector().record(MU_list[i].soma(0.5)._ref_ina))

    return t, v, amp, mgi, cell_t, cell_id, ina #Return each variable recorded above as a NEURON Vector

def Spike_Detector(mu_list):
    cell_t = h.Vector()
    cell_id = h.Vector()
    for i in range(len(mu_list)):
        nc = h.NetCon(mu_list[i].soma(0.5)._ref_v,None, sec=mu_list[i].soma)
        nc.record(cell_t, cell_id)
    return cell_t, cell_id

def Run_Sim(N):#Simulation parameters

    mu_list = MU_List(N)#Make Motor Neuron object list
    soma_list, dend_list = Soma_List(mu_list)#Make Soma list
    MU_Fin = Change_Dia(mu_list,soma_list,dend_list,N,200)#Updated Motor Neuron object list
    
    #Print Coordinates of cells (Used to check results)

    # print(MU_Fin[0].soma.x3d(0),MU_Fin[0].soma.y3d(1),MU_Fin[0].soma.z3d(0),MU_Fin[0].soma.diam3d(0))
    # print(MU_Fin[1].soma.x3d(0),MU_Fin[1].soma.y3d(0),MU_Fin[1].soma.z3d(0),MU_Fin[1].soma.diam3d(0))
    # print(MU_Fin[2].soma.x3d(0),MU_Fin[2].soma.y3d(0),MU_Fin[2].soma.z3d(0),MU_Fin[2].soma.diam3d(0))

    # print(MU_Fin[0].soma.x3d(1),MU_Fin[0].soma.y3d(1),MU_Fin[0].soma.z3d(1),MU_Fin[0].soma.diam3d(1))
    # print(MU_Fin[1].soma.x3d(1),MU_Fin[1].soma.y3d(1),MU_Fin[1].soma.z3d(1),MU_Fin[1].soma.diam3d(1))
    # print(MU_Fin[2].soma.x3d(1),MU_Fin[2].soma.y3d(1),MU_Fin[2].soma.z3d(1),MU_Fin[2].soma.diam3d(1))

    mul=40
    WU = 5000
    WUf = 5000
    time1=50000 + WU+ WUf#Total time of simulation (ms)
    dela=0
    dur=32500

    dura = (dur+WU)*mul
    ramp_dur = (17500+WU)*mul
    WUd = WU*mul
    WUfd = WUf*mul
    pkamp=18.1
    bias=0
    min = 6
    std=1

    #(nA)
    on = 0*mul #ms
#---------------------------------------------------
    fr_mean = np.genfromtxt('Trap_signal_T3.txt', delimiter=",")
    fr_mean = np.interp(fr_mean, (fr_mean.min(), fr_mean.max()), (min, pkamp))
    
    fr_array = np.ones([len(fr_mean),2])
    for i in range(len(fr_mean)):
        fr_array[i,:] = [i,fr_mean[i]]

    desired_time = 5
    time1= int(desired_time)

    iv=np.arange(0,time1,0.001) #Time vector for current calculation.
    #Make Neuron Vector.
    istim_cur = np.empty(len(iv))
    ival=0
    
    time_fact = desired_time/np.max(fr_array[:,0])
    fr_array[:,0]=fr_array[:,0]*time_fact    

    f = interpolate.interp1d(fr_array[:,0],fr_array[:,1])
    
    istim_cur = f(iv)

    frequency = 1
    amplitude=1
    sampling_rate= len(istim_cur)/time1
    duration=time1
    t=np.arange(0, duration, 1/sampling_rate)
    sine_wave= amplitude * np.sin(2*np.pi*frequency*t)

    istim_cur += sine_wave

    for i in range(len(istim_cur)):
        istim_cur[i] += np.random.normal(0,1)
   
    

    stim = [] #Initialize list
    iv=iv*1000
    ivv = h.Vector().from_python(iv) 
    for i in range(len(MU_Fin)): #Generate and apply Stimulation object to each Soma of each Motor Neuron object.
        stim.append(h.IClamp(MU_Fin[i].soma(0.5)))

    istim_cur_1 = istim_cur
    for i in range(len(stim)): #Generate stimulation parameters for each Stimulation object.
        istim_cur_1 = istim_cur

        for j in range(len(istim_cur)):
            istim_cur_1[j] = istim_cur[j] + np.random.normal(0,1.5)

        i_tv=h.Vector().from_python(istim_cur_1) #Generate Neuron Vector.
        stim[i].dur = time1*1000
        stim[i].delay = dela
        i_tv.play(stim[i]._ref_amp, ivv)


    stimulus = stim

    t, v, amp, mgi, cell_t, cell_id, ina = Record(stimulus, MU_Fin) #Record desired characteristics.

    #Initialize Simulation.
    h.finitialize(-70)
    h.celsius = 36
    h.dt = 0.01
    h.continuerun(time1*1000)

    return t, v, amp, mgi, cell_t, cell_id, ina #Return desired recrodings.
#------------------------------------------------------------------------


if __name__ == "__main__": #Semantic statement for Multiprocessing.
    
    hdf5_file_path = 'simulation_data.h5'
    
    # Delete the existing HDF5 file if it exists
    if os.path.exists(hdf5_file_path):
        os.remove(hdf5_file_path)
        print(f"Deleted existing {hdf5_file_path}")
    
    print("I Started")
    num=200 #Num of desired Neurons
    N = [i for i in range(num)] #Generate range for Multiprocessing command.

    MU_Model = [] #Model List.

    start = time.perf_counter() #For elapsed time counter
    with multiprocessing.Pool(6) as pool: #Run simulation with as many available cpus as possible.        
        MU_Model.append(pool.map(Run_Sim, N))
    
    print('Finished Sim ready to write')
    #MU_model[i][j][k][l]
    # i: access the list
    # j: access desired Motor Neuron
    # k: access desired recorded characteristic (Time, Voltage, Current, Force, mgi)
    # l: acces desired characteristic vector
    holder = np.ones(len(MU_Model[0]))
    for i in range(len(MU_Model[0])):
        holder[i] = len(MU_Model[0][i][4])
    mval = np.argmax(holder)
    t_array = np.empty([len(MU_Model[0][0][0][0]),len(MU_Model[0])])
    v_array = np.empty([len(MU_Model[0][0][0][0]),len(MU_Model[0])])
    amp_array = np.empty([len(MU_Model[0][0][0][0]),len(MU_Model[0])])
    mgi_array = np.empty([len(MU_Model[0][0][3][0]),len(MU_Model[0])])
    cell_t_array = np.ones([len(MU_Model[0][mval][4]),len(MU_Model[0])])
    cell_id_array = np.ones([len(MU_Model[0][0][5]),len(MU_Model[0])])
    ina_array = np.ones([len(MU_Model[0][0][6][0]),len(MU_Model[0])])

    # Create an HDF5 file to store your data
    with h5py.File('simulation_data.h5', 'w') as file:
        # Create datasets for each type of data
        for i in range(len(MU_Model[0])):
            file.create_dataset(f't_array_{i}', data=np.array(MU_Model[0][i][0][0]))
            file.create_dataset(f'v_array_{i}', data=np.array(MU_Model[0][i][1][0]))
            file.create_dataset(f'amp_array_{i}', data=np.array(MU_Model[0][i][2][0]))
            file.create_dataset(f'mgi_array_{i}', data=np.array(MU_Model[0][i][3][:]))
            if len(MU_Model[0][i][4])>0:
                file.create_dataset(f'cell_t_array_{i}', data=np.array(MU_Model[0][i][4]))
            else:
                file.create_dataset(f'cell_t_array_{i}', data=[0])
            file.create_dataset(f'ina_array_{i}', data=np.array(MU_Model[0][i][6][0]))

    stop = time.perf_counter()
    print(f"elapsed time: {stop - start} sec")
    plt.plot(MU_Model[0][0][0][0],MU_Model[0][0][2][0])
    plt.show()
hines
Site Admin
Posts: 1706
Joined: Wed May 18, 2005 3:32 pm

Re: Python/Neuron model runs on Windows but does not on Linux.(h.finitialize)

Post by hines »

I'm afraid the code you copied above is incomplete. After installing the necessary modules in a python virtual environment, I'm stuck at

Code: Select all

$ python -i test.py
NEURON: Couldn't find: kim.hoc
 near line 0
 objref hoc_obj_[2]
                   ^
I Started
multiprocessing.pool.RemoteTraceback: 
"""
Traceback (most recent call last):
  File "/usr/lib/python3.12/multiprocessing/pool.py", line 125, in worker
    result = (True, func(*args, **kwds))
                    ^^^^^^^^^^^^^^^^^^^
  File "/usr/lib/python3.12/multiprocessing/pool.py", line 48, in mapstar
    return list(map(*args))
           ^^^^^^^^^^^^^^^^
  File "/home/hines/tmp/test.py", line 120, in Run_Sim
    mu_list = MU_List(N)#Make Motor Neuron object list
              ^^^^^^^^^^
  File "/home/hines/tmp/test.py", line 41, in MU_List
    MU=h.Kim()
       ^^^^^
AttributeError: 'hoc.HocObject' object has no attribute 'Kim'
"""

Post Reply