Merge pull request #403 from Hjorthmedh/region_mesh_redux

Region mesh redux

examples/parallel/KTH_PDC/input_tuning/Dardel_run_input_tuning_cortical_signal_dspn.job

@@ -2,8 +2,8 @@
 #SBATCH --partition=main
 #SBATCH -o log/InputTuning-%j-output.txt
 #SBATCH -e log/InputTuning-%j-error.txt
-#SBATCH -t 8:59:00
-#SBATCH --time-min=7:59:00
+#SBATCH -t 3:59:00
+#SBATCH --time-min=2:59:00
 #SBATCH -J InputTuning
 #SBATCH -A naiss2023-5-231
 #SBATCH --nodes=8-12
@@ -14,10 +14,11 @@
 export SNUDDA_TUNE_NEURON="dspn"
 export SEED_LIST="[10,20,30,40,50]"
 export INPUT_TYPE="cortical"
+export INPUT_FREQ=40
 
 # You need to point this as the directory where you created the network in
 #NETWORK_DIR=/cfs/klemming/home/${USER:0:1}/$USER/Snudda/examples/parallel/KTH_PDC/networks/test_10k
-NETWORK_DIR=networks/input_tuning_${SNUDDA_TUNE_NEURON}_${INPUT_TYPE}_signal
+NETWORK_DIR=networks/input_tuning_${SNUDDA_TUNE_NEURON}_${INPUT_TYPE}_signal_${INPUT_FREQ}Hz
 
 # This is used for NEURON
 export N_WORKERS=$SLURM_NTASKS

examples/parallel/KTH_PDC/input_tuning/Dardel_run_input_tuning_thalamic_signal_dspn.job

@@ -2,8 +2,8 @@
 #SBATCH --partition=main
 #SBATCH -o log/InputTuning-%j-output.txt
 #SBATCH -e log/InputTuning-%j-error.txt
-#SBATCH -t 8:59:00
-#SBATCH --time-min=7:59:00
+#SBATCH -t 3:59:00
+#SBATCH --time-min=1:59:00
 #SBATCH -J InputTuning
 #SBATCH -A naiss2023-5-231
 #SBATCH --nodes=8-12
@@ -14,10 +14,11 @@
 export SNUDDA_TUNE_NEURON="dspn"
 export SEED_LIST="[10,20,30,40,50]"
 export INPUT_TYPE="thalamic"
+export INPUT_FREQ=40
 
 # You need to point this as the directory where you created the network in
 #NETWORK_DIR=/cfs/klemming/home/${USER:0:1}/$USER/Snudda/examples/parallel/KTH_PDC/networks/test_10k
-NETWORK_DIR=networks/input_tuning_${SNUDDA_TUNE_NEURON}_${INPUT_TYPE}_signal
+NETWORK_DIR=networks/input_tuning_${SNUDDA_TUNE_NEURON}_${INPUT_TYPE}_signal_${INPUT_FREQ}Hz
 
 # This is used for NEURON
 export N_WORKERS=$SLURM_NTASKS

examples/parallel/KTH_PDC/input_tuning/PD/Dardel_run_input_tuning_background_dspn.job

@@ -0,0 +1,132 @@
+#!/bin/bash -l
+#SBATCH --partition=main
+#SBATCH -o log/InputTuning-%j-output.txt
+#SBATCH -e log/InputTuning-%j-error.txt
+#SBATCH -t 8:59:00
+#SBATCH --time-min=7:59:00
+#SBATCH -J InputTuning
+#SBATCH -A naiss2023-5-231
+#SBATCH --nodes=6-10
+#SBATCH --tasks-per-node=128
+#SBATCH --mem-per-cpu=930M
+#SBATCH --mail-type=ALL
+
+export SNUDDA_TUNE_NEURON="dspn"
+export SEED_LIST="[10,20,30,40,50]"
+
+# You need to point this as the directory where you created the network in
+NETWORK_DIR=../networks/input_tuning_${SNUDDA_TUNE_NEURON}_background
+
+# This is used for NEURON
+export N_WORKERS=$SLURM_NTASKS
+
+# export N_WORKERS=50
+# export IPNWORKERS=$N_WORKERS
+
+# This is used for ipyparallel and Snudda
+export IPNWORKERS=50
+
+# Clear old ipyparallel
+export IPYTHONDIR="/cfs/klemming/scratch/${USER:0:1}/$USER/.ipython"
+rm -r $IPYTHONDIR
+export IPYTHON_PROFILE=default
+
+module load snic-env
+source $HOME/Snudda/snudda_env/bin/activate
+SNUDDA_DIR=/cfs/klemming/home/"${USER:0:1}"/$USER/Snudda
+
+# WE NEED TO SET THE RIGHT PD SNUDDA_DATA directory
+exit
+
+# If the BasalGangliaData directory exists, then use that for our data
+if [[ -d "/cfs/klemming/home/${USER:0:1}/$USER/BasalGangliaData/data" ]]; then
+    export SNUDDA_DATA="/cfs/klemming/home/${USER:0:1}/$USER/BasalGangliaData/data"
+    echo "Setting SNUDDA_DATA to $SNUDDA_DATA"
+    rm mechanisms
+    ln -s $SNUDDA_DATA/neurons/mechanisms/ mechanisms
+else
+    echo "SNUDDA_DATA environment variable not changed (may be empty): $SNUDDA_DATA"
+    rm mechanisms
+    ln -s ../../../../snudda/data/neurons/mechanisms/
+fi
+
+
+# Start ipyparallel
+#.. Start the ipcontroller
+export FI_CXI_DEFAULT_VNI=$(od -vAn -N4 -tu < /dev/urandom)
+srun -n 1 -N 1 -c 2 --exact --overlap --mem=0 ./../../ipcontroller_new.sh &
+
+
+echo ">>> waiting 60s for controller to start"
+sleep 60 
+
+#.. Read in CONTROLLERIP
+CONTROLLERIP=$(<controller_ip.txt)
+
+##.. Start the engines
+echo ">>> starting ${IPNWORKERS} engines "
+echo "SLURM_JOB_NUM_NODES={SLURM_JOB_NUM_NODES}"
+
+export FI_CXI_DEFAULT_VNI=$(od -vAn -N4 -tu < /dev/urandom)
+srun -n ${IPNWORKERS} -c 2 -N ${SLURM_JOB_NUM_NODES} --exact --overlap --mem=0 ipengine \
+--location=${CONTROLLERIP} --profile=${IPYTHON_PROFILE} --mpi \
+--ipython-dir=${IPYTHONDIR}  --timeout=30.0 c.EngineFactory.max_heartbeat_misses=10  c.MPI.use='mpi4py' \
+1> ipe_${SLURM_JOBID}.out 2> ipe_${SLURM_JOBID}.err &
+
+
+echo ">>> waiting 60s for engines to start"
+sleep 30
+
+
+echo "Network dir: "$NETWORK_DIR
+
+export PATH=$SNUDDA_DIR/snudda_env/bin/:$PATH
+export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$CRAY_LD_LIBRARY_PATH
+export PYTHONPATH=$SNUDDA_DIR/snudda_env/lib/python3.9/
+
+# This will stop NEURON from failing with "can't open DISPLAY"
+unset DISPLAY
+
+# Do we need this magic?
+# export FI_CXI_DEFAULT_VNI=$(od -vAn -N4 -tu < /dev/urandom)
+# srun -n 1 python test_script.py &> test_print.txt
+
+export FI_CXI_DEFAULT_VNI=$(od -vAn -N4 -tu < /dev/urandom)
+
+# srun -n 1 python setup_input_tuning_background.py &> log/setup_network_log.txt
+srun -n 1 -N 1 --exact --overlap --mem=0 python ../setup_input_tuning_background.py &> log/setup_network_log.txt
+
+
+# Quit ipyparallel
+
+python ../ipcontroller_shutdown.py
+
+
+##############
+
+rm -r x86_64
+
+export CXX=CC
+export CC=cc
+export FC=ftn
+export MPICC=cc
+export MPICXX=CC
+
+CC --version
+
+echo "About to run nrnivmodl"
+which nrnivmodl
+
+# srun -n nrnivmodl mechanisms/
+
+# Do we need this magic?
+export FI_CXI_DEFAULT_VNI=$(od -vAn -N4 -tu < /dev/urandom)
+
+srun -n 1 nrnivmodl -incflags "-lltdl=/usr/lib64/libltdl.so.7 -lreadline=/lib64/libreadline.so.7 -lncurses=/lib64/libncurses.so.6.1" -loadflags "-DLTDL_LIBRARY=/usr/lib64/libltdl.so.7 -DREADLINE_LIBRARY=/lib64/libreadline.so.7 -DNCURSES_LIBRARY=/lib64/libncurses.so.6.1" mechanisms/
+
+# Do we need this magic?
+export FI_CXI_DEFAULT_VNI=$(od -vAn -N4 -tu < /dev/urandom)
+
+srun -n $N_WORKERS $SNUDDA_DIR/examples/parallel/KTH_PDC/input_tuning/x86_64/special -mpi -python $SNUDDA_DIR/snudda/input/input_tuning.py simulate $NETWORK_DIR --seed_list $SEED_LIST
+
+

examples/parallel/KTH_PDC/input_tuning/setup_input_tuning_signal.py

@@ -20,14 +20,20 @@
 else:
     input_type = "cortical"
 
+if os.getenv("INPUT_FREQ"):
+    input_freq = int(os.getenv("INPUT_FREQ"))
+else:
+    input_freq = 10
+
+
 if os.getenv("SEED_LIST"):
     seed_list = ast.literal_eval(os.getenv("SEED_LIST"))
 else:
     seed_list = None
 
 print(f"Optimising input for neuron type {neuron_type}")
 
-network_path = os.path.join("networks", f"input_tuning_{neuron_type}_{input_type}_signal")
+network_path = os.path.join("networks", f"input_tuning_{neuron_type}_{input_type}_signal_{input_freq}Hz")
 
 
 from ipyparallel import Client
@@ -56,7 +62,7 @@
                          num_input_min=1,
                          num_input_max=200,
                          input_duration=10.0,
-                         input_frequency_range=[10.0],
+                         input_frequency_range=[input_freq],
                          use_meta_input=True)
 
 

snudda/init/init.py

@@ -451,7 +451,7 @@ def add_neurons(self, name,
                     print(tstr)
                     sys.exit(-1)
 
-                print("Checking boundaries, to make sure P is not too high")
+                # print("Checking boundaries, to make sure P is not too high")
                 x = np.zeros((8, 1))
                 y = np.zeros((8, 1))
                 z = np.zeros((8, 1))
@@ -468,10 +468,10 @@ def add_neurons(self, name,
                 # axon_density function of x,y,z defined above
                 p_corner = numexpr.evaluate(axon_density[1]) * (3e-6 ** 3)
 
-                for P, xx, yy, zz in zip(p_corner, x, y, z):
-                    print(f"{name} axon density P({xx}, {yy}, {zz}) = {P}")
-
                 if (p_corner > 0.01).any():
+                    for P, xx, yy, zz in zip(p_corner, x, y, z):
+                        print(f"{name} axon density P({xx}, {yy}, {zz}) = {P}")
+
                     print("Axon density too high at boundary!!")
                     print("Please increase bounding box")
                     sys.exit(-1)
@@ -854,11 +854,13 @@ def define_striatum(self,
                         f_FS=0.013,
                         f_ChIN=0.011,
                         f_LTS=0.007,
+                        f_NGF=0.0019,
                         num_dSPN=None,
                         num_iSPN=None,
                         num_FS=None,
                         num_ChIN=None,
                         num_LTS=None,
+                        num_NGF=None,
                         volume_type=None,
                         side_len=None,
                         slice_depth=None,
@@ -881,8 +883,9 @@ def define_striatum(self,
             self.num_FS = get_val(num_FS)
             self.num_ChIN = get_val(num_ChIN)
             self.num_LTS = get_val(num_LTS)
+            self.num_NGF = get_val(num_NGF)
 
-            self.num_neurons_total += self.num_FS + self.num_dSPN + self.num_iSPN + self.num_ChIN + self.num_LTS
+            self.num_neurons_total += self.num_FS + self.num_dSPN + self.num_iSPN + self.num_ChIN + self.num_LTS + self.num_NGF
             num_neurons = self.num_neurons_total
 
             if self.num_neurons_total <= 0:
@@ -893,15 +896,16 @@ def define_striatum(self,
                 # No neurons specified, skipping structure
                 return
 
-            f_tot = f_dSPN + f_iSPN + f_FS + f_ChIN + f_LTS
+            f_tot = f_dSPN + f_iSPN + f_FS + f_ChIN + f_LTS + f_NGF
 
             self.num_FS = np.round(f_FS * num_neurons / f_tot)
             self.num_dSPN = np.round(f_dSPN * num_neurons / f_tot)
             self.num_iSPN = np.round(f_iSPN * num_neurons / f_tot)
             self.num_ChIN = np.round(f_ChIN * num_neurons / f_tot)
             self.num_LTS = np.round(f_LTS * num_neurons / f_tot)
+            self.num_NGF = np.round(f_NGF * num_neurons / f_tot)
 
-            n_neurons = int(self.num_FS + self.num_dSPN + self.num_iSPN + self.num_ChIN + self.num_LTS)
+            n_neurons = int(self.num_FS + self.num_dSPN + self.num_iSPN + self.num_ChIN + self.num_LTS + self.num_NGF)
 
             self.num_neurons_total += n_neurons
             if abs(num_neurons - self.num_neurons_total) > 5:
@@ -937,6 +941,7 @@ def define_striatum(self,
             self.add_neuron_density(volume_id="Striatum", neuron_type="FS", density_file=density_file)
             self.add_neuron_density(volume_id="Striatum", neuron_type="LTS", density_file=density_file)
             self.add_neuron_density(volume_id="Striatum", neuron_type="ChIN", density_file=density_file)
+            # Todo: Density missing for NGF
 
         elif volume_type == "slice":
             self.define_structure(struct_name="Striatum",
@@ -997,6 +1002,8 @@ def define_striatum(self,
         iSPN_dir = os.path.join(neurons_dir, "striatum", "ispn")
         ChIN_dir = os.path.join(neurons_dir, "striatum", "chin")
         LTS_dir = os.path.join(neurons_dir, "striatum", "lts")
+        NGF_dir = os.path.join(neurons_dir, "striatum", "ngf")
+
 
         # Add the neurons
 
@@ -1074,6 +1081,15 @@ def define_striatum(self,
                          axon_density=LTS_axon_density,
                          volume_id="Striatum")
 
+        # NGF
+        if os.path.isdir(snudda_parse_path(NGF_dir, self.snudda_data)):
+            self.add_neurons(name="NGF", neuron_dir=NGF_dir,
+                             num_neurons=self.num_NGF,
+                             volume_id="Striatum")
+
+        else:
+            print(f"No directory {NGF_dir}, skipping NGF cells.")
+
         # Define FS targets
 
         # Szydlowski SN, Pollak Dorocic I, Planert H, Carlen M, Meletis K,
@@ -1476,6 +1492,70 @@ def define_striatum(self,
                                mod_file="tmGabaA",
                                channel_param_dictionary=None)
 
+        ####
+
+        if True:
+            # Connections to and from NGF
+            # NGF -> SPN 25/29 connected within 100 micrometers (Ibanez-Sandoval, et al 2011)
+            # NGF -> SPN 11/14 connected (English et al, 2012)
+
+            # TODO: Parameters not tuned yet
+            self.add_neuron_target(neuron_name="NGF",
+                                   target_name="dSPN",
+                                   connection_type="GABA",
+                                   dist_pruning=None,
+                                   f1=0.2308, soft_max=None, mu2=0.5659, a3=1.0,
+                                   conductance=0.5e-9,
+                                   mod_file="ngf_tmGabaA")  # This file does not yet exist
+
+            self.add_neuron_target(neuron_name="NGF",
+                                   target_name="iSPN",
+                                   connection_type="GABA",
+                                   dist_pruning=None,
+                                   f1=0.3527, soft_max=None, mu2=0.2811, a3=0.9997,
+                                   conductance=0.5e-9,
+                                   mod_file="ngf_tmGabaA")  # This file does not yet exist
+
+            # NGF -> ChIN 3/14 (English et al, 2012)
+            self.add_neuron_target(neuron_name="NGF",
+                                   target_name="ChIN",
+                                   connection_type="GABA",
+                                   dist_pruning=None,
+                                   f1=0.088, soft_max=None, mu2=0.17141, a3=1,
+                                   conductance=0.5e-9,
+                                   mod_file="ngf_tmGabaA")  # This file does not yet exist
+
+            # NGF -> FS, Kocaturk et al, 2022 -- 12/20 ???
+
+            # Gap junctions, 1/2 English et al 2012
+            # TODO: Optimise!!
+            self.add_neuron_target(neuron_name="NGF",
+                                   target_name="NGF",
+                                   connection_type="GapJunction",
+                                   dist_pruning=None,
+                                   f1=0.1364, soft_max=None, mu2=0.4625, a3=1.0,
+                                   conductance=0.5e-9)
+
+            # Move these to respective neuron later...
+
+            # FS -> NGF  9/9, Lee et al, 2022
+            self.add_neuron_target(neuron_name="FS",
+                                   target_name="NGF",
+                                   connection_type="GABA",
+                                   dist_pruning=None,
+                                   f1=0.0988, soft_max=None, mu2=0.0624, a3=0.9997,
+                                   conductance=0.5e-9,
+                                   mod_file="tmGabaA")
+
+            # ChIN -> NGF 8/14, English et al, 2012
+            self.add_neuron_target(neuron_name="ChIN",
+                                   target_name="NGF",
+                                   connection_type="GABA",
+                                   dist_pruning=None,
+                                   f1=0.9208, soft_max=10, mu2=0.3393, a3=1.0,
+                                   conductance=0.5e-9,
+                                   mod_file=None)  # This file does not yet exist
+
     ############################################################################
 
     def define_GPe(self, num_neurons, d_min=None, neurons_dir=None):