diff --git a/examples/install.jl b/examples/install.jl
index cdbd7300df830ea6dfa3ba68b9b2aae2e6bfb4dc..06fcf47ca0760cff199304ae4d81d85e910a2dd6 100644
--- a/examples/install.jl
+++ b/examples/install.jl
@@ -1,6 +1,4 @@
-import Pkg; Pkg.add("Plots")
-import Pkg; Pkg.add("UnicodePlots")
-Pkg.add("OrderedCollections")
+
include("../src/SpikingNeuralNetworks.jl")
include("../src/units.jl")
@@ -9,8 +7,15 @@ try
using UnicodePlots
using PyCall
using NSGAII
+ using OrderedCollections
+ using Plotly
catch
import Pkg
+ Pkg.add("Plotly")
+ import Pkg; Pkg.add("Plots")
+ import Pkg; Pkg.add("UnicodePlots")
+ Pkg.add("OrderedCollections")
+
Pkg.clone("https://github.com/gsoleilhac/NSGAII.jl")
Pkg.add("ProgressMeter")
Pkg.add("UnicodePlots")
@@ -19,7 +24,7 @@ catch
Pkg.add("PyCall")
Pkg.build("PyCall")
#Pkg.clone("https://github.com/gsoleilhac/NSGAIII.jl")
- Pkg.clone("https://github.com/gsoleilhac/NSGAII.jl")
+ #Pkg.clone("https://github.com/JuliaCN/Py2Jl.jl")
using PyCall
using UnicodePlots
end
diff --git a/examples/iz_neuron.jl b/examples/iz_neuron.jl
index 4ee546261c8166ce419de49cd7b779dc21ecc570..540d69d9d7e311f540bbfa5e68b5bf420d9e50aa 100644
--- a/examples/iz_neuron.jl
+++ b/examples/iz_neuron.jl
@@ -1,4 +1,9 @@
-using Plots, SNN
+include("../src/SpikingNeuralNetworks.jl")
+include("../src/units.jl")
+
+SNN = SpikingNeuralNetworks
+
+using Plots
RS = SNN.IZ(;N = 1, param = SNN.IZParameter(;a = 0.02, b = 0.2, c = -65, d = 8))
IB = SNN.IZ(;N = 1, param = SNN.IZParameter(;a = 0.02, b = 0.2, c = -55, d = 4))
diff --git a/examples/izhi.py b/examples/izhi.py
new file mode 100755
index 0000000000000000000000000000000000000000..dc1db3764636a797d05f89c55de1e24ed01c82c6
--- /dev/null
+++ b/examples/izhi.py
@@ -0,0 +1,133 @@
+import pickle
+from neo.core import AnalogSignal
+import sciunit
+import sciunit.capabilities as scap
+import neuronunit.capabilities as cap
+import neuronunit.capabilities.spike_functions as sf
+import quantities as pq
+from elephant.spike_train_generation import threshold_detection
+def Id(t,delay,duration,tmax,amplitude):
+ if 0.0 < t < delay:
+ return 0.0
+ elif delay < t < delay+duration:
+
+ return amplitude#(100.0)
+
+ elif delay+duration < t < tmax:
+ return 0.0
+ else:
+ return 0.0
+
+import julia
+jl = julia.Julia()
+from julia import Main
+# from sciunit.models.runnable import RunnableModel
+# probably runnable is the julia crasher
+class IZModel(sciunit.Model,
+ cap.ReceivesSquareCurrent,
+ cap.ProducesActionPotentials,
+ cap.ProducesMembranePotential,
+ scap.Runnable):
+ """A model which produces a frozen membrane potential waveform."""
+
+ def __init__(self,attrs={},backend="JULIA_IZ"):
+ """Create an instace of a model that produces a static waveform.
+ """
+ self.attrs = attrs
+ self.vm = None
+ #self.Main = Main
+
+
+ self._backend = self
+ self.backend = backend
+
+ def set_run_params(self,t_stop=None):
+ #if 'tmax' in self.params.keys():
+ # self.t_stop=self.params['tmax']
+ #else:
+ self.t_stop = t_stop
+ def get_membrane_potential(self, **kwargs):
+ """Return the Vm passed into the class constructor."""
+ return self.vm
+ def set_stop_time(self,tmax=1300*pq.ms):
+ self.tmax = tmax
+ def set_attrs(self, attrs):
+ JUIZI = {
+ 'a': 0.02,
+ 'b': 0.2,
+ 'c': -65,
+ 'd': 8,
+ }
+ if not len(attrs):
+ attrs = JUIZI
+
+ Main.attrs = attrs
+ Main.eval('param = SNN.IZParameter(;a = attrs["a"], b = attrs["b"], c = attrs["c"], d = attrs["d"])')
+ self.attrs.update(attrs)
+ Main.N_ = attrs["N"]
+ Main.eval('E2 = SNN.IZ(;N = 1, param = param)')
+ Main.eval("N = Int32(1)")
+
+
+
+ def inject_square_current(self, current):#, section = None, debug=False):
+ """Inputs: current : a dictionary with exactly three items, whose keys are: 'amplitude', 'delay', 'duration'
+ Example: current = {'amplitude':float*pq.pA, 'delay':float*pq.ms, 'duration':float*pq.ms}}
+ where \'pq\' is a physical unit representation, implemented by casting float values to the quanitities \'type\'.
+ Description: A parameterized means of applying current injection into defined
+ Currently only single section neuronal models are supported, the neurite section is understood to be simply the soma.
+
+ """
+ import numpy as np
+ if 'injected_square_current' in current.keys():
+ c = current['injected_square_current']
+ else:
+ c = current
+ duration = float(c['duration'])#.rescale('ms'))
+ delay = float(c['delay'])#.rescale('ms'))
+ amp = 1000000.0*float(c['amplitude'])#.rescale('uA')
+ tmax = 1.3#00
+ tmin = 0.0
+ DT = 0.25
+ T = np.linspace(tmin, tmax, int(tmax/DT))
+ Iext_ = []
+ for t in T:
+ Iext_.append(Id(t,delay,duration,tmax,amp))
+ self.attrs['N'] = len(Iext_)
+ self.set_attrs(self.attrs)
+ Main.eval('pA = 0.001nA')
+ Main.eval("SNN.monitor(E2, [:v])")
+ Main.dur = current["duration"]
+ Main.current = current
+ Main.delay = float(current["delay"])
+ Main.temp_current = float(amp)
+ Main.eval("E2.I = [deepcopy(temp_current)*pA]")
+ Main.eval('SNN.sim!([E2], []; dt ='+str(DT)+'*ms, delay=delay,stimulus_duration=1000,simulation_duration = 1300)')
+ Main.eval("v = SNN.getrecord(E2, :v)")
+ #Main.eval("SNN.vecplot(E2, :v) |> display")
+ v = Main.v
+ self.vM = AnalogSignal(v,units = pq.mV,sampling_period = DT * pq.ms)
+ return self.vM
+ def get_membrane_potential(self):
+ return self.vM
+
+ def get_spike_count(self):
+ thresh = threshold_detection(self.vM,threshold=0*pq.mV)
+ return len(thresh)
+ def get_backend(self):
+ return self
+ def get_spike_train(self):
+ thresh = threshold_detection(self.vM)
+ return thresh
+ def _backend_run(self):
+ results['vm'] = self.vM.magnitude
+ results['t'] = self.vM.times
+
+ return results
+
+
+ def get_APs(self, **run_params):
+ """Return the APs, if any, contained in the static waveform."""
+ vm = self.get_membrane_potential(**run_params)
+ waveforms = sf.get_spike_waveforms(vm)
+ return waveforms
diff --git a/examples/lizneuron.jl b/examples/lizneuron.jl
new file mode 100644
index 0000000000000000000000000000000000000000..984dc63698a703cc5ef8d97f41474bfe664a88af
--- /dev/null
+++ b/examples/lizneuron.jl
@@ -0,0 +1,62 @@
+try
+ using NSGAII
+catch
+ include("install.jl")
+end
+#using Pkg
+#Pkg.clone("https://github.com/JuliaCN/Py2Jl.jl")
+
+#using Py2Jl
+
+include("local_iz_neuron.jl")
+using NSGAII
+include("../src/units.jl")
+using .HHNSGA
+#=
+function fdecode(x)
+ decoded = NSGAII.decode(x, HHNSGA.bc)
+ return decoded
+end
+=#
+using Plots
+unicodeplots()
+
+function plot_pop(P)
+ P = filter(x -> x.rank == 1, P)
+ plot(map(x -> x.y[1], P), map(x -> x.y[2], P))#, "bo", markersize = 1)
+end
+#@bp
+
+#HHNSGA.init_function()
+#@bp
+@time repop = NSGAII.nsga(4,4,HHNSGA.z,HHNSGA.init_function)#, fplot = plot_pop)
+@time best = sort(repop, by = ind -> ind.y[1])[end];
+@time worst = sort(repop, by = ind -> ind.y[1])[end];
+
+@time decoded = fdecode(best.pheno)
+println("z = $(x1.y)")
+py"""
+from izhi import IZModel
+model = IZModel()
+"""
+py"model.attrs" = decoded
+
+py"""
+rt = test['Rheobase test']
+rheo = rt.generate_prediction(model)
+"""
+
+vm = py"model.get_membrane_potential()"
+Pkg.add("Plotly")
+using Plotly
+
+
+trace1 = [
+ "x" => vm,
+ "y" => vm.times,
+ "mode" => "lines+markers",
+ "type" => "scatter"
+]
+data = [trace1]
+response = Plotly.plot(data, ["filename" => "line-scatter", "fileopt" => "overwrite"])
+plot_url = response["url"]
diff --git a/examples/local_hh_neuron.jl b/examples/local_hh_neuron.jl
index 080bcd2e0542dcccf75d40b5cb2a242170bbefdd..3241b1bf1d0ac961325b6b093072d25aebaee8f5 100644
--- a/examples/local_hh_neuron.jl
+++ b/examples/local_hh_neuron.jl
@@ -93,7 +93,9 @@ module HHNSGA
H1=[values(ranges)]
current_params = py"rt.params"
simple.attrs = py"JHH"
+ #using Py2Jl
py"""
+
def evaluate(test,god):
model = SimpleModel()
model.attrs = god
@@ -101,21 +103,31 @@ module HHNSGA
#rtp = RheobaseTestP(rt.observation)
rheo = rt.generate_prediction(model)
+ # errors = tuple([100.0,100.0,100.0,100.0,100.0,100.0,100.0,100.0])
+ # return errors
scores_ = []
+ #assert 'value' in rheo.keys()
+
+ # for temp_test in test.values():
+ # #print(temp_test.name,temp_test.params)
+ # #print(temp_test.params['injected_square_current']['amplitude'])
+ # #if 'InjectedCurrentAPWidthTest' not in temp_test.params.keys():
+ # # temp_test.params['injected_square_current'] = {}
+ # if temp_test.name in "InjectedCurrentAPWidthTest":
+ # temp_test.params['injected_square_current']['amplitude'] = rheo['value']
+ # if temp_test.name in "InjectedCurrentAPThresholdTest":
+ # temp_test.params['injected_square_current']['amplitude'] = rheo['value']
+ # if temp_test.name in "InjectedCurrentAPThresholdTest":
+ # temp_test.params['injected_square_current']['amplitude'] = rheo['value']
+ # if temp_test.name in "RheobaseTest":
+ # temp_test.params['injected_square_current']['amplitude'] = rheo['value']
for key, temp_test in test.items():
#if key == 'Rheobase test':
+ # def to_map(key_val):
temp_test.score_type = sciunit.scores.ZScore
if key in ['Injected current AP width test', 'Injected current AP threshold test', 'Injected current AP amplitude test', 'Rheobase test']:
temp_test.params['amplitude'] = rheo['value']
- if temp_test.name == "InjectedCurrentAPWidthTest":
- temp_test.params['amplitude'] = rheo['value']
- if temp_test.name == "InjectedCurrentAPThresholdTest":
- temp_test.params['amplitude'] = rheo['value']
- if temp_test.name == "InjectedCurrentAPThresholdTest":
- temp_test.params['amplitude'] = rheo['value']
- if temp_test.name == "RheobaseTest":
- temp_test.params['amplitude'] = rheo['value']
try:
@@ -129,12 +141,15 @@ module HHNSGA
score = 100#np.inf
if isinstance(score, sciunit.scores.incomplete.InsufficientDataScore):
score = 100#np.inf
- if score == 100:
+ #if score == 100:
print(temp_test.name)
- print(temp_test.params['amplitude'])
- #import pdb
- #pdb.set_trace()
+ print(temp_test.params)
+ if not rheo['value'] is None:
+
+ import pdb
+ pdb.set_trace()
scores_.append(score)
+
tt = TestSuite(list(test.values()))
SA = ScoreArray(tt, scores_)
errors = tuple(SA.values,)
@@ -143,12 +158,12 @@ module HHNSGA
return errors
"""
+
+
py"""
def z_py(evaluate,god):
-
cell_tests = pickle.load(open('multicellular_constraints.p','rb'))
tests = cell_tests[list(cell_tests.keys())[0]]
-
SA = evaluate(tests,god)
return SA
"""
@@ -164,6 +179,7 @@ module HHNSGA
genes = py"genes"
return genes
end
+
export init_function2
function init_function2(x)
py"""
diff --git a/examples/local_iz_neuron.jl b/examples/local_iz_neuron.jl
new file mode 100644
index 0000000000000000000000000000000000000000..6211764c7d500d1d99a2ad2e14a591f32fa0018f
--- /dev/null
+++ b/examples/local_iz_neuron.jl
@@ -0,0 +1,233 @@
+
+module HHNSGA
+ using Pkg
+ using Plots
+ using UnicodePlots
+ using OrderedCollections
+ using LinearAlgebra
+ using UnicodePlots
+ using PyCall
+ using NSGAII
+ export NSGAII
+
+ export SNN
+
+
+ include("../src/SpikingNeuralNetworks.jl")
+ include("../src/units.jl")
+ include("../src/plot.jl")
+ SNN = SpikingNeuralNetworks.SNN
+ using Random
+ py"""
+ import copy
+ from neuronunit.optimisation import algorithms
+ from neuronunit.optimisation import optimisations
+ from neuronunit.optimisation.optimization_management import TSD
+ from neuronunit.optimisation.optimization_management import OptMan
+ from neuronunit.optimisation import make_sim_tests
+ import matplotlib
+ import matplotlib.pyplot as plt
+ import numpy as np
+ matplotlib.get_backend()
+ import sys
+ import os
+ sys.path.append(os.getcwd())
+ """
+ py"""
+ from neuronunit import tests
+ from neuronunit.tests import fi
+ from izhi import IZModel
+ from neo import AnalogSignal
+ import quantities as pqu
+ try:
+ from julia import Main
+ except:
+ import julia
+ julia.install()
+ from julia import Main
+ Main.include("../src/SpikingNeuralNetworks.jl")
+ Main.eval("SNN = SpikingNeuralNetworks.SNN")
+
+
+ """
+ ##model = py"IZModel()"
+ simple = py"IZModel()"
+
+ py"""
+ import pickle
+ import numpy as np
+ import random
+ from sciunit.scores.collections import ScoreArray
+ from sciunit import TestSuite
+ import sciunit
+ import random
+
+ from neuronunit.tests.fi import RheobaseTestP
+ from collections import OrderedDict
+ from neuronunit import get_neab
+ import pdb
+ pdb.set_trace()
+ cell_tests = pickle.load(open('multicellular_constraints.p','rb'))
+ for test in cell_tests.values():
+ if "Rheobase test" in test.keys():
+ temp_test = {k:v for k,v in test.items()}
+ break
+ rt = temp_test["Rheobase test"]
+ #rtp = RheobaseTestP(rt.observation)
+ JUIZI = {
+ 'a': 0.02,
+ 'b': 0.2,
+ 'c': -65,
+ 'd': 8,
+ }
+
+ #ranges = OrderedDict({k:[v-0.5*np.abs(v),v+0.5*np.abs(v)] for k,v in copy.copy(JUIZI).items()})
+ temp = {'a':[0.02,0.1],'b':[0.2,0.26],'c':[-65,-50],'d':[0.05,8]}
+ ranges = OrderedDict(temp)
+
+ N = len(JUIZI)
+ """
+ N = py"N"
+ ranges = py"ranges"
+ #H1=[values(ranges)]
+ #current_params = py"rt.params"
+ #simple.attrs = py"JUIZI"
+ py"""
+ from izhi import IZModel
+
+ def evaluate(test,god):
+ model = IZModel()
+ model.attrs = god
+ rt = test['Rheobase test']
+ rheo = rt.generate_prediction(model)
+ scores_ = []
+ test = {k:v for k,v in test.items() if v.observation is not None}
+ for temp_test in test.values():
+ if 'injected_square_current' not in temp_test.params.keys():
+ temp_test.params['injected_square_current'] = {}
+ temp_test.params['amplitude'] = rheo['value']
+
+ if temp_test.name == "InjectedCurrentAPWidthTest":
+ temp_test.params['injected_square_current']['amplitude'] = rheo['value']
+ temp_test.params['amplitude'] = rheo['value']
+
+ if temp_test.name == "Injected CurrentAPThresholdTest":
+ temp_test.params['injected_square_current']['amplitude'] = rheo['value']
+ temp_test.params['amplitude'] = rheo['value']
+
+ if temp_test.name == "InjectedCurrentAPAmplitudeTest":
+ temp_test.params['injected_square_current']['amplitude'] = rheo['value']
+ temp_test.params['amplitude'] = rheo['value']
+
+ if temp_test.name == "RheobaseTest":
+ temp_test.params['injected_square_current']['amplitude'] = rheo['value']
+ temp_test.params['amplitude'] = rheo['value']
+
+ to_map = [(tt,model) for tt in test.values() ]
+ def map_score(content):
+ temp_test = content[0]
+ model = content[1]
+ temp_test.score_type = sciunit.scores.ZScore
+ try:
+ score = temp_test.judge(model)
+ score = np.abs(score.log_norm_score)
+ except:
+ try:
+ score = np.abs(float(score.raw_score))
+ except:
+ score = 100
+ if isinstance(score, sciunit.scores.incomplete.InsufficientDataScore):
+ score = 100
+ return score
+
+ scores_ = list(map(map_score,to_map))
+ tt = TestSuite(list(test.values()))
+ SA = ScoreArray(tt, scores_)
+ errors = tuple(SA.values,)
+ return errors
+ """
+
+ py"""
+ def z_py(evaluate,god):
+
+ cell_tests = pickle.load(open('multicellular_constraints.p','rb'))
+ tests = cell_tests['Neocortex pyramidal cell layer 5-6']
+
+
+ model_type="RAW"
+ from neuronunit.optimisation import make_sim_tests
+ fps = ['a','b','c','d']
+ #sim_tests, OM, target = make_sim_tests.test_all_objective_test(fps,model_type=model_type)
+ SA = evaluate(tests,god)
+ return SA
+ """
+
+
+ export transdict
+ function transdict(x)
+ py"""
+ def map_dict(x):
+ genes_out = x
+ #genes_out_dic = OrderedDict({k:genes_out[i] for i,k in enumerate(ranges.keys()) })
+ return genes_out#, genes_out_dic
+ """
+ decoded = x
+
+ try
+ decoded = NSGAII.decode(x, HHNSGA.bc)
+ catch
+ decoded = x
+ end
+ genes_out = py"map_dict"(decoded)
+ bincoded = NSGAII.encode(genes_out, bc)
+ return bincoded, genes_out_dic
+ end
+ py"""
+ lower_list = [v[0] for k,v in ranges.items()]
+ upper_list = [v[1] for k,v in ranges.items()]
+ """
+
+ const bc = BinaryCoding(4, [:Int,:Int,:Int,:Int], py"lower_list", py"upper_list")
+ export init_function
+ function init_function()
+ py"""
+ def pheno_map():
+ lower_list = [v[0] for k,v in ranges.items()]
+ upper_list = [v[1] for k,v in ranges.items()]
+ gene_out = [random.uniform(lower, upper) for lower, upper in zip(lower_list, upper_list)]
+ #gene_out_dic = OrderedDict({k:gene_out[i] for i,k in enumerate(ranges.keys()) })
+ return gene_out#, gene_out_dic
+ """
+ gene_out = py"pheno_map"()
+ @time bincoded = NSGAII.encode(gene_out, bc)
+ return bincoded
+ end
+
+ function references_for_H()
+ py"""
+ from sklearn.model_selection import ParameterGrid
+ grid = ParameterGrid(ranges)
+ genes = []
+ for g in grid:
+ genes.append(list(g.values()))
+ genes_out_dic = OrderedDict({k:genes[i] for i,k in enumerate(ranges.keys()) })
+
+ """
+ genes_out = py"genes"
+ bcd = []
+ for g in genes_out
+ bincoded = NSGAII.encode(g, bc)
+ append!(bcd,bincoded)
+ end
+ return bcd
+ end
+
+
+ export z
+ function z(x)
+ genes_out,god = transdict(x)
+ contents = py"z_py"(py"evaluate",god)
+ return contents
+ end
+
+end
diff --git a/examples/simple_with_injection.py b/examples/simple_with_injection.py
index 82e8efbfbd4b562455d3bb067ac2c5454a5b6f12..a5501270aeb4061e0f903ae89ef3dc21413e1840 100755
--- a/examples/simple_with_injection.py
+++ b/examples/simple_with_injection.py
@@ -37,7 +37,7 @@ class SimpleModel(sciunit.Model,
scap.Runnable):
"""A model which produces a frozen membrane potential waveform."""
- def __init__(self,attrs={},backend="JULIA"):
+ def __init__(self,attrs={},backend="JULIA_HH"):
"""Create an instace of a model that produces a static waveform.
"""
self.attrs = attrs