russelljjarvis
commented
4 years ago backend/model/test interaction has built in assumptions now about model.set_stop_time and `test.protocol['sim_length']``
Consequences of this bug:
- A waveform that should last a whole 1second lasts for 500 ms.
- The NU tests cannot tell that waves are contracted in time.
- Scores for Capacitance Test Time Constant Test and Spike Width test are artificially improved due to the time contraction distortion.
The notebook contents need translating to be compatible with the newer versions of NU. Code below reproduces minimal version of chapter3 containing the bug
# Chapter 3
Back to [Chapter 2](chapter2.ipynb)
%matplotlib inline
import os,sys
import numpy as np
import matplotlib.pyplot as plt
import quantities as pq
import sciunit
import neuronunit
from neuronunit import aibs
from neuronunit.models.reduced import ReducedModel
import quantities as pq
from neuronunit import tests as nu_tests, neuroelectro
neuron = {'nlex_id': 'nifext_50'} # Layer V pyramidal cell
tests = []
dataset_id = 354190013 # Internal ID that AIBS uses for a particular Scnn1a-Tg2-Cre
# Primary visual area, layer 5 neuron.
try:
assert suite is not None
assert len(suite) == 8
except:
import quantities as pq
from neuronunit import tests as nu_tests, neuroelectro
neuron = {'nlex_id': 'nifext_50'} # Layer V pyramidal cell
tests = []
dataset_id = 354190013 # Internal ID that AIBS uses for a particular Scnn1a-Tg2-Cre
# Primary visual area, layer 5 neuron.
# Obtain the empirical rheobase current from the Allen Brain Insitute Cell Types database.
observation = aibs.get_observation(dataset_id,'rheobase')
rheobase_test = nu_tests.RheobaseTest(observation=observation)
tests += [rheobase_test]
test_class_params = [(nu_tests.InputResistanceTest,None),
(nu_tests.TimeConstantTest,None),
(nu_tests.CapacitanceTest,None),
(nu_tests.RestingPotentialTest,None),
(nu_tests.InjectedCurrentAPWidthTest,None),
(nu_tests.InjectedCurrentAPAmplitudeTest,None),
(nu_tests.InjectedCurrentAPThresholdTest,None)
]
# Obtain all other parameters from neuroelectro.org.
for cls,params in test_class_params:
observation = cls.neuroelectro_summary_observation(neuron)
tests += [cls(observation,params=params)]
# A hook to update all tests after the RheobaseTest to use the rheobase current (to produce exactly one AP)
def update_amplitude(test,tests,score):
rheobase = score.prediction['value']
for test in tests[1:]:
if 'Injected' in test.name:
# Set current injection to just suprathreshold
test.params['injected_square_current']['amplitude'] = rheobase*1.01
hooks = {tests[0]:{'f':update_amplitude}}
suite = sciunit.TestSuite(tests,name="vm_suite",hooks=hooks)
suite[0].observation['mean'] = suite[0].observation['value']
suite[0].observation['std'] = suite[0].observation['value']
with open('suite.p','wb') as f:
pickle.dump(suite,f)
# This example is from https://github.com/OpenSourceBrain/IzhikevichModel.
for rel_path in ['.','models','../models']:
# Check to see if this is the neuronunit/docs directory.
DOCS_PATH = os.path.abspath(os.path.join(os.getcwd(),rel_path))
if DOCS_PATH.endswith('neuronunit/models'):
break
assert DOCS_PATH.endswith('neuronunit/models'), ("Could not find the path to neuronunit/docs. "
"Change the DOCS_PATH variable to the path to "
"your neuronunit/docs directory.")
LEMS_MODEL_PATH = os.path.join(DOCS_PATH,'NeuroML2/LEMS_2007One.xml')
model = ReducedModel(LEMS_MODEL_PATH,name='vanilla',backend='jNeuroML')
suite[0].params
for t in suite:
t.params2 = {}
t.params2['injected_sqaure_current'] = t.params
t.params = t.params2
suite[0].params
model = None
model = ReducedModel(LEMS_MODEL_PATH,name='vanilla',backend='jNeuroML')
try:
model.inject_square_current(suite[0].params['injected_sqaure_current'])
except:
model.inject_square_current(suite[0].params)
vm = model.get_membrane_potential()
len(vm)
vm.sampling_rate
vm.times[-1]/len(vm)
print(len(vm))
import matplotlib.pyplot as plt
plt.plot(vm.times,vm.magnitude)
plt.show()
suite[0].params = suite[0].params['injected_sqaure_current']
print(suite[0].params['delay'])
print(suite[0].params['duration'])This last part reveals the bug:
sim_length = suite[0].params['delay'] + suite[0].params['duration']
print(vm.times[-1] == sim_length)
print(vm.times[-1],sim_length)
print(vm.times[-1] == sim_length/2.0)
print(vm.times[-1],sim_length/2.0)For contrast you can compare with a newer version of Izhikevich model which has accurate Neo sampling period, and scores less well on those NU tests.
from neuronunit.optimisation.data_transport_container import DataTC
from neuronunit.optimisation.optimization_management import inject_and_plot_model, inject_and_not_plot_model
dtc = DataTC()
dtc.backend = "RAW"
fast_model = dtc.dtc_to_model()
dtc.attrs = fast_model.default_attrs
results = inject_and_not_plot_model(dtc)
print(len(results))
print(len(vm))
inject_and_plot_model(dtc)
rgerkin
@rgerkin I am not sure if the significance of this result was understood. From my perspective it meant that jNeuroML models could not be regarded as a ground truth until this issue was resolved. Perhaps the issue was resolved, since I could not run the final notebooks it was not possible to evaluate the status of this issue. @ChihweiLHBird, since your working installation of NU is likely to be a lot closer to Rick's I was wondering if you would be able to externally verify Rick's fixes at the very bottom of this thread?