diff --git a/Dockerfile b/Dockerfile
index fa8a230b55989478700a9c62473290eb7981e901..3d2901a30c188ef8c2bbdea45522a0cf775ffbe7 100644
--- a/Dockerfile
+++ b/Dockerfile
@@ -68,7 +68,7 @@ RUN python -c "import julia; julia.install()"
RUN echo "redo"
RUN /opt/conda/bin/pip install rickpy
RUN /opt/conda/bin/pip install git+https://github.com/scidash/sciunit@dev
-RUN echo "Redo"
+RUN echo "Redo 1"
WORKDIR $HOME/work
RUN git clone -b unittest https://github.com/russelljjarvis/NeuronunitOpt
WORKDIR NeuronunitOpt
diff --git a/examples/lhhneuron.jl b/examples/lhhneuron.jl
index 4b0513ea23850a5fa17ab2e401f579a027228fee..01194f225504243e1e63b60e8e7d8633fec9e810 100644
--- a/examples/lhhneuron.jl
+++ b/examples/lhhneuron.jl
@@ -1,10 +1,46 @@
-using Distributed
-addprocs(1)
-@everywhere include("local_hh_neuron.jl")
+#using Distributed
+#addprocs(1)
+#@everywhere
+include("local_hh_neuron.jl")
using NSGAII
-@everywhere include("../src/units.jl")
-@everywhere using .HHNSGA
-pop = pmap(HHNSGA.init_function,1:3)
+#@everywhere
+include("../src/units.jl")
+#@everywhere
+using .HHNSGA
+#=
+seed = pop[1]
+bincoded = NSGAII.encode(seed, HHNSGA.bc)
+decoded = NSGAII.decode(bincoded, HHNSGA.bc)
+scores = map(HHNSGA.z,pop)
+@show(scores)
+
+pop = []
+for i in 1:3
+ append!(pop,[HHNSGA.init_function()])
+end
+=#
+
+
+function fdecode(x)
+ decoded = NSGAII.decode(x, HHNSGA.bc)
+
+ #y = [ convert(Int64,i) for i in Int(x) ]
+ return decoded
+end
+#@test
+#@show decoded
+#@show seed
+#assert decoded .≈ seed
+
+#lower_list = [v[0] for k,v in ranges.items()]
+#upper_list = [v[1] for k,v in ranges.items()]
+
+#NSGAII.create_indiv(pop[1], fdecode, HHNSGA.z, HHNSGA.init_function2)
+
+
+#@show(pop)
+#pop = map(HHNSGA.init_function,1:3)
+#NSGAII.create_indiv(pop[1],HHNSGA.z)
#=
println(pop)
@@ -16,24 +52,39 @@ include("../src/plot.jl")
#SNN = SpikingNeuralNetworks.SNN
=#
-scores = pmap(HHNSGA.z,pop)
+#repop = NSGAII.nsga(2,2,HHNSGA.z,HHNSGA.init_function())
+#append!(pop,
+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)
+ #show()
+ #sleep(0.1)
+end
+
+repop = NSGAII.nsga(10,10,HHNSGA.z,HHNSGA.init_function, fplot = plot_pop)
+x1 = sort(repop, by = ind -> ind.y[1])[end];
+decoded = fdecode(x1.pheno)
+println("z = $(x1.y)")
#=
-NSGAIII.indiv(copy(pop[1]),HHNSGA.z)
+NSGAII.indiv(copy(pop[1]),HHNSGA.z)
println(scores)
non_dom = NSGAIII.fast_non_dominated_sort!(scores)[1]
pop = setdiff(pop, non_dom)
HHNSGA.plot_pop(scores)
-=#
@show(pop)
@show(scores)
try
pop2 = Any[]
for x in pop
- pop2.append(NSGAII.create_indiv(x, HHNSGA.init_function2, HHNSGA.z))
+ append!(pop,NSGAII.create_indiv(x, HHNSGA.init_function2, HHNSGA.z))
end
catch
@show(pop)
@show(scores)
end
+=#
diff --git a/examples/local_hh_neuron.jl b/examples/local_hh_neuron.jl
index 7cebbdb1d1a97159dcd36fc330c0cc9c6c13590d..080bcd2e0542dcccf75d40b5cb2a242170bbefdd 100644
--- a/examples/local_hh_neuron.jl
+++ b/examples/local_hh_neuron.jl
@@ -10,7 +10,7 @@ module HHNSGA
using PyCall
using NSGAII
export NSGAII
- export nS
+ #export nS
catch
include("install.jl")
@@ -27,6 +27,10 @@ module HHNSGA
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
+ import joblib
+
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
@@ -50,13 +54,16 @@ module HHNSGA
"""
-
simple = py"SimpleModel()"
-
py"""
import pickle
import numpy as np
- #from neuronunit.tests.fi import RheobaseTestP
+ import random
+ from sciunit.scores.collections import ScoreArray
+ from sciunit import TestSuite
+ import sciunit
+
+ from neuronunit.tests.fi import RheobaseTestP
from collections import OrderedDict
cell_tests = pickle.load(open('multicellular_constraints.p','rb'))
for test in cell_tests.values():
@@ -64,7 +71,7 @@ module HHNSGA
temp_test = {k:v for k,v in test.items()}
break
rt = temp_test["Rheobase test"]
-
+ #rtp = RheobaseTestP(rt.observation)
JHH = {
'Vr': -68.9346,
'Cm': 0.0002,
@@ -77,52 +84,62 @@ module HHNSGA
'Vt': -63.0
}
- ranges = OrderedDict({k:[v-0.01*np.abs(v),v+0.01*np.abs(v)] for k,v in copy.copy(JHH).items()})
+ ranges = OrderedDict({k:[v-0.5*np.abs(v),v+0.5*np.abs(v)] for k,v in copy.copy(JHH).items()})
N = len(JHH)
model = SimpleModel()
"""
N = py"N"
-
- #
ranges = py"ranges"
H1=[values(ranges)]
-
current_params = py"rt.params"
- #print(current_params)
simple.attrs = py"JHH"
-
py"""
def evaluate(test,god):
model = SimpleModel()
model.attrs = god
rt = test['Rheobase test']
+ #rtp = RheobaseTestP(rt.observation)
rheo = rt.generate_prediction(model)
-
- from sciunit.scores.collections import ScoreArray
- from sciunit import TestSuite
- import sciunit
scores_ = []
for key, temp_test in test.items():
- if key == 'Rheobase test':
- temp_test.score_type = sciunit.scores.ZScore
+ #if key == 'Rheobase test':
+ 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']:
- rt.params['amplitude'] = rheo['value']
+ 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:
score = temp_test.judge(model)
score = np.abs(score.log_norm_score)
except:
- score = 100#np.inf
+ try:
+ score = np.abs(float(score.raw_score))
+ except:
+ score = 100#np.inf
if isinstance(score, sciunit.scores.incomplete.InsufficientDataScore):
score = 100#np.inf
+ if score == 100:
+ print(temp_test.name)
+ print(temp_test.params['amplitude'])
+ #import pdb
+ #pdb.set_trace()
scores_.append(score)
tt = TestSuite(list(test.values()))
SA = ScoreArray(tt, scores_)
errors = tuple(SA.values,)
+
return errors
"""
@@ -152,26 +169,42 @@ module HHNSGA
py"""
def one(x):
if str('x') in locals():
-
genes_out = x
genes_out_dic = OrderedDict({k:genes_out[i] for i,k in enumerate(ranges.keys()) })
-
else:
- import random
lower_list = [v[0] for k,v in ranges.items()]
upper_list = [v[1] for k,v in ranges.items()]
genes_out = [random.uniform(lower, upper) for lower, upper in zip(lower_list, upper_list)]
genes_out_dic = OrderedDict({k:genes_out[i] for i,k in enumerate(ranges.keys()) })
+
return genes_out, genes_out_dic
"""
- genes_out, genes_out_dic = py"one"(x)
- return genes_out, genes_out_dic#py"genes_out",py"genes_out_dic"
+ decoded = x
+
+ try
+ decoded = NSGAII.decode(x, HHNSGA.bc)
+ catch
+ decoded = x
+ println("no")
+ end
+ genes_out, genes_out_dic = py"one"(decoded)
+ #@show(genes_out)
+
+ bincoded = NSGAII.encode(genes_out, bc)
+ #@show(bincoded)
+ return bincoded, genes_out_dic#py"genes_out",py"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()]
+ #print(len(lower_list))
+ """
+ const bc = BinaryCoding(9, [:Int,:Int,:Int,:Int,:Int,:Int,:Int,:Int,:Int], py"lower_list", py"upper_list")
export init_function
- function init_function(x)
+ function init_function()
py"""
- def one(x):
+ def one():
import random
lower_list = [v[0] for k,v in ranges.items()]
upper_list = [v[1] for k,v in ranges.items()]
@@ -179,8 +212,11 @@ module HHNSGA
genes_out_dic = OrderedDict({k:genes_out[i] for i,k in enumerate(ranges.keys()) })
return genes_out, genes_out_dic
"""
- genes_out, genes_out_dic = py"one"(x)
- return genes_out
+ genes_out, genes_out_dic = py"one"()
+ #@show(genes_out)
+ bincoded = NSGAII.encode(genes_out, bc)
+ #@show(bincoded)
+ return bincoded
end
export z
@@ -189,9 +225,8 @@ module HHNSGA
#x = genes_out
# The slow part can it be done in parallel.
contents = py"z_py"(py"evaluate",god)
- # fCV)
-
- #contents = [ convert(Float64,c) for c in contents ]
+ @show(contents)
+ println("gene fitness calculated")
return contents
end
diff --git a/examples/simple_with_injection.py b/examples/simple_with_injection.py
index 42719687fa2d27992a1f5f35f1257a09a8c5a560..82e8efbfbd4b562455d3bb067ac2c5454a5b6f12 100755
--- a/examples/simple_with_injection.py
+++ b/examples/simple_with_injection.py
@@ -18,6 +18,16 @@ def Id(t,delay,duration,tmax,amplitude):
else:
return 0.0
+import julia
+jl = julia.Julia()
+from julia import Main
+
+#Main.eval('include("../src/SpikingNeuralNetworks.jl")')
+#Main.eval('include("../src/units.jl")')
+#Main.eval('include("../src/plot.jl")')
+#Main.eval('SNN = SpikingNeuralNetworks.SNN')
+
+
# from sciunit.models.runnable import RunnableModel
# probably runnable is the julia crasher
class SimpleModel(sciunit.Model,
@@ -32,10 +42,8 @@ class SimpleModel(sciunit.Model,
"""
self.attrs = attrs
self.vm = None
- import julia
- jl = julia.Julia()
- from julia import Main
- self.Main = Main
+ #self.Main = Main
+
self._backend = self
self.backend = backend
@@ -64,9 +72,9 @@ class SimpleModel(sciunit.Model,
}
if not len(attrs):
attrs = JHH
- Main = self.Main
+ #Main = self.Main
- #Main.eval("SNN.HH(;N = 1)")
+ Main.eval("SNN.HH(;N = 1)")
Main.attrs = attrs
Main.eval('param = SNN.HHParameter(;El = attrs["El"], Ek = attrs["EK"], En = attrs["ENa"], gl = attrs["gl"], gk = attrs["gK"], gn = attrs["gNa"])')#', N = attrs["N"])')
self.attrs.update(attrs)
@@ -77,7 +85,7 @@ class SimpleModel(sciunit.Model,
Main.eval('E2 = SNN.HH(;N = 1)')
Main.eval('E2.param = param')
Main.eval("N = Int32(1)")
- Main.eval("@show(N)")
+ #Main.eval("@show(N)")
Main.eval('E2.v = ones(N).*attrs["Vr"]')
Main.eval("E2.m = zeros(N)")
Main.eval("E2.n = zeros(N)")
@@ -119,9 +127,6 @@ class SimpleModel(sciunit.Model,
Iext_.append(Id(t,delay,duration,tmax,amp))
self.attrs['N'] = len(Iext_)
self.set_attrs(self.attrs)
-
-
- Main = self.Main
Main.eval('pA = 0.001nA')
Main.eval("SNN.monitor(E2, [:v])")
Main.dur = current["duration"]
@@ -129,16 +134,14 @@ class SimpleModel(sciunit.Model,
Main.delay = float(current["delay"])
Main.temp_current = float(amp)
Main.eval("E2.I = [deepcopy(temp_current)*pA]")
- Main.eval("@show(E2.I)")
+ #Main.eval("@show(E2.I)")
Main.eval('SNN.sim!([E2], []; dt ='+str(DT)+'*ms, delay=delay,stimulus_duration=1000,simulation_duration = 1300)')
- #print('gets here')
#Main.eval('SNN.sim!([E2], []; dt = 0.015*ms, delay=current["delay"]*ms,stimulus_duration=1000*ms,simulation_duration = 1300*ms)')
Main.eval("v = SNN.getrecord(E2, :v)")
v = Main.v
#Main.eval('SNN.vecplot(E2, :v) |> display')
self.vM = AnalogSignal(v,units = pq.mV,sampling_period = DT * pq.ms)
- #print("done one.")
return self.vM
def get_membrane_potential(self):
return self.vM
@@ -151,6 +154,11 @@ class SimpleModel(sciunit.Model,
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):