Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Model review #343

Merged
merged 17 commits into from
Apr 10, 2017
Merged

Model review #343

merged 17 commits into from
Apr 10, 2017

Conversation

PTraeder
Copy link
Contributor

@PTraeder PTraeder commented Jan 30, 2017

Clear warnings in models folder. Some bugfixes.
Some warnings in generation pipeline still exist because Units are not yet printed correctly for generation.

Some cleanup.

Signed-off-by: Philip Traeder <philip.traeder@rwth-aachen.de>
…g plain names.

Signed-off-by: Philip Traeder <philip.traeder@rwth-aachen.de>
Instead of equalizing "ignoreMagnitude" use prettyprinted type names for type comparison
This additionally fixes a bug where a unitless Unit ("[0,0,0,0,0,0,0,0]")
and a real would not be compatible in additons.

Signed-off-by: Philip Traeder <philip.traeder@rwth-aachen.de>
where the differential order was not taken into account

Signed-off-by: Philip Traeder <philip.traeder@rwth-aachen.de>
Signed-off-by: Philip Traeder <philip.traeder@rwth-aachen.de>
Signed-off-by: Philip Traeder <philip.traeder@rwth-aachen.de>
@DimitriPlotnikov
Copy link
Contributor

@PTraeder please, provide a more informative PR description.

@jougs and @Silmathoron are kindly invited to review the PR.

@Silmathoron
Copy link
Member

Hi there, sorry, the email address I'm using for GitHub has been blocked, so I did not receive the notification... I've been quite busy lately, but I hope I'll be able to have a look at this by the end of the week, or next week, at worse.

Copy link
Member

@Silmathoron Silmathoron left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi there, I started reading the PR, but I'm really surprise by what I see, so I'd like to discuss this first before I continue.

  1. I think adding * 1 [ms/s] everywhere is really eavy and should be something we want to avoid at all cost since it is not something we can ask user to remember (it makes no sense from the point of view of an end user).
    If this is really necessary, it must not take place in the frontend.

  2. Same remark applies for I_stim = currents.get_sum() * 1 pA & co: I thought we agreed that the get_sum() function was supposed to return the right unit depending on the buffer (see discussion in Address #324 and Implement compound assignments #328) so why is this 1 pA here?

Signed-off-by: Philip Traeder <philip.traeder@rwth-aachen.de>
Fix models accordingly.

Signed-off-by: Philip Traeder <philip.traeder@rwth-aachen.de>
Fix bug where magnitudes would be dumped in a factor without regard
for the factors exponent.

Signed-off-by: Philip Traeder <philip.traeder@rwth-aachen.de>
new subclass to expressionPrinter: legacyExpressionPrinter = old behaviour, does not print units
legacyPrinter used for solverGenerator since python does not know SIunits.

Signed-off-by: Philip Traeder <philip.traeder@rwth-aachen.de>
Signed-off-by: Philip Traeder <philip.traeder@rwth-aachen.de>
Copy link
Member

@Silmathoron Silmathoron left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

On the whole, this looks much better!
There are still a few things I'm unsure about, especially when it comes to seemingly arbitrary units popping in some places... I'll try to have a look at this in more details during the week.

Act_m' = ( alpha_m - ( alpha_m + beta_m ) * Act_m ) * 1 [1/s]
Act_h' = ( alpha_h - ( alpha_h + beta_h ) * Act_h ) * 1 [1/s]
Inact_n' = ( alpha_n - ( alpha_n + beta_n ) * Inact_n ) * 1 [1/s]
Act_m' = ( alpha_m - ( alpha_m + beta_m ) * Act_m ) / 1 ms
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[this is a memo] I need to check where this unit comes from

@@ -66,10 +66,10 @@ neuron hh_psc_alpha_implicit:
equations:
# synapses: alpha functions
I_syn_ex'' = -I_syn_ex' / tau_syn_ex
I_syn_ex' = I_syn_ex' - ( I_syn_ex / tau_syn_ex ) * 1 [ms/s]
I_syn_ex' = I_syn_ex' - ( I_syn_ex / tau_syn_ex )
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[memo2] this is a highly disturbing way of writing this ODE, there has to be another way!

@@ -99,7 +99,7 @@ neuron ht_neuron_nestml:
recordable function I_h pA = -h_g_peak * Ih_m * ( V_m - h_E_rev )
# The spike current is only activate immediately after a spike.
function I_spike mV = (g_spike) ? -( V_m - E_K ) / Tau_spike : 0
V_m' = ( ( I_Na * 1 pA + I_K * 1 pA + I_syn + I_NaP + I_KNa + I_T + I_h + I_stim ) / Tau_m + I_spike * 1 [pA/(ms * mV)] ) * 1 [s/F]
V_m' = ( ( I_Na * 1 pA + I_K * 1 pA + I_syn + I_NaP + I_KNa + I_T + I_h + I_stim ) / Tau_m + I_spike * 1 [pA/(ms * mV)] ) * 1 [Gs/F]
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

How did this giga second get there? also the extraneous spaces should be removed

Copy link
Contributor Author

@PTraeder PTraeder Feb 27, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I seem to have accidentally deleted my comment:
could have as well been s/nF, I merely dumped the necessary factor into time here.
LHS: mV/ms = V/s
RHS: (pA/ms)*(Gs/F) = (pA/ms) * (s/nF) = A/F = V/s

Changed it to s/nF by now.

Copy link
Member

@Silmathoron Silmathoron Feb 27, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

RHS = mV/ms
I = [pA] = dq/dt = d(CV)/dt = C dV / dt = [pF] * [mV / ms], i.e. your missing factor is just 1pF (at the denominator)

Copy link
Contributor Author

@PTraeder PTraeder Feb 27, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Note that on the RHS we have pA/ms not just pA. Hence the missing factor is ms/pF = s/nF

EDIT: Or am I missing something here?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@heplesser: in your role as the original author of the ht_neuron in NEST, can you please have a look at the use of units in the NESTML version? I have a feeling that some of the /1ms or /1mV that get rid of the units are not actually needed but stem from wrong units for some of the variables in the first place.

@@ -124,8 +124,8 @@ neuron iaf_psc_delta_neuron:
# for decay until end of refractory period
if with_refr_input:
refr_spikes_buffer += spikes * exp(-r * h / tau_m) * 1 mV
else:
spikes.get_sum() # clear buffer entry, ignore spike
#else:
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

is this not required? I think this should be the default behavior...

Copy link
Contributor Author

@PTraeder PTraeder Feb 27, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

With the removal of get_sum(), I am not sure what NESTML construct would be used to clear the buffer. Hence my commenting out the old, dysfunct code, pending review.

Signed-off-by: Philip Traeder <philip.traeder@rwth-aachen.de>
Signed-off-by: Philip Traeder <philip.traeder@rwth-aachen.de>
Copy link

@heplesser heplesser left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@PTraeder At @jougs request, I had a look at some of the nestml examples, see my comments in the code. I have the impression that unit handling is rather unsystematic and that we should have a more high-level look at the approach to units before we work on individual examples.

I_stim = currents.get_sum()
g_ex += spikeExc * 1 nS
g_in += spikeInh * 1 nS
I_stim = currents

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It looks inconsistent that spikeExc and spikeInh are multiplied by a unit, but currents is not. Shouldn't units for the spike input also be handled where the spikes actually arrive, since the connection weights carry weight and the ring buffers in NEST store sums of weights, not spike counts?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Spike buffers are typed as real. g_ex and g_in are both nS.
Current buffers are already typed as pA (same as I_stim in the model) so no factor is necessary.

Copy link
Contributor

@DimitriPlotnikov DimitriPlotnikov Mar 6, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

currents are per construction always pA, since it models an external current which has by definition pA as its dimension. In contrast to current buffers, spike buffers may be either current-based or conductance-based. Since it is not always possible (or it is not transparent for the modeler) we decided to make spike buffers unitless. Thus, the modeler is forced to state the desired type of the input by multiplying it with a corret unit literal (1pA in the case of current based models, 1nS in the case of the conductance based models).

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@heplesser
We will tackle explicit typing of ports in the upcoming release.

function alpha_m_init real = 0.32 * ( 13. - V_m / 1 mV) / ( exp( ( 13. - V_m / 1 mV) / 4. ) - 1. )
function beta_m_init real = 0.28 * ( V_m /1 mV - 40. ) / ( exp( ( V_m / 1 mV - 40. ) / 5. ) - 1. )
function alpha_h_init real = 0.128 * exp( ( 17. - V_m / 1 mV) / 18. )
function beta_h_init real = 4. / ( 1. + exp( ( 40. - V_m / 1 mV ) / 5. ) )

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The units here are the wrong way around. This is mostly the fault of our neurobiologist colleagues who very rarely put proper units into there equations, a tradition going back to at least Hodgkin & Huxley (1952). Assuming that alpha and beta here are rate constants in standard HH-style, they have units of 1/ms. Since V has unit mV, we can deduce all other units and insert them, e.g., for the first equation (slightly simplifying, hope it is possible to write 15 mV instead of 15 * 1 mV)

0.032  / ( 1 ms * 1 mV ) * ( 15 mV - V_m ) / ( exp( ( 15 mV - V_m ) / ( 5 mV ) ) - 1 )

And is it necessary to use decimal points everywhere?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It should not be necessary to points everywhere (but I have to check what happens if it translated to C++). Yes you can use 15mV as an unit literal.

Copy link
Contributor Author

@PTraeder PTraeder Mar 6, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I changed it in both hh_cond_exp_traub models.

Inact_n' = alpha_n - ( alpha_n + beta_n ) * Inact_n
Act_m' = ( alpha_m - ( alpha_m + beta_m ) * Act_m ) / 1 ms
Act_h' = ( alpha_h - ( alpha_h + beta_h ) * Act_h ) / 1 ms
Inact_n' = ( alpha_n - ( alpha_n + beta_n ) * Inact_n ) / 1 ms

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why are these equations divided by 1 ms? The alpha and beta are rates, so they have unit 1/ms and the derivative should have that unit, too, shouldn't it?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Those rates are defined as real, not as 1/ms. I was not aware that this should be the case when I modified the models. Changing the definition of alpha and beta rates from real type to 1/ms would indeed make these extra factors unnecessary.

Inact_n' = alpha_n - ( alpha_n + beta_n ) * Inact_n
Act_m' = ( alpha_m - ( alpha_m + beta_m ) * Act_m ) * 1 [1/ms]
Act_h' = ( alpha_h - ( alpha_h + beta_h ) * Act_h ) * 1 [1/ms]
Inact_n' = ( alpha_n - ( alpha_n + beta_n ) * Inact_n ) * 1 [1/ms]

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Similar to above, but I find the notation here more confusing; / 1 ms would probably be easier, but I think no unit manipulation should happen here at all.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

See above

@@ -99,21 +99,20 @@ neuron ht_neuron_nestml:
recordable function I_h pA = -h_g_peak * Ih_m * ( V_m - h_E_rev )
# The spike current is only activate immediately after a spike.
function I_spike mV = (g_spike) ? -( V_m - E_K ) / Tau_spike : 0
V_m' = ( I_Na + I_K + I_syn + I_NaP + I_KNa + I_T + I_h + I_stim ) / Tau_m + I_spike

V_m' = ( ( I_Na * 1 pA + I_K * 1 pA + I_syn + I_NaP + I_KNa + I_T + I_h + I_stim ) / Tau_m + I_spike * 1 [pA/(ms * mV)] ) * 1 [s/nF]

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why are some currents multiplied with units and others not? This seems to point to inconsistencies and errors elsewhere. All currents should "arrive" here with proper units, including I_spike. I would also suggest that, for easier comparison with the equation in the paper, we should have I_spike/tau_spike here, instead of incorporating the time constant in I_spike; as it is now, I_spike looks like a current but isn't, and that is going to cause misunderstandings later.

Note one challenge with the Hill-Tononi model: In that model, all conductances are dimensionless, which also means that currents are measured in mV. Especially, the 1 s/nF term is not needed.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

They would arrive there with proper units if they were defined with proper units.
Fixed in update

# I_KNa
function D_influx_peak real = 0.025
function tau_D real = 1250.0 # yes, 1.25s
function D_thresh real = -10.0
function D_slope real = 5.0
function D_influx real = 1.0 / ( 1.0 + exp( -( V_m - D_thresh ) / D_slope ) )
function D_influx real = 1.0 / ( 1.0 + exp( -( V_m / 1 mV- D_thresh ) / D_slope ) )

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

D_thresh and D_slope have units mV.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

see above

@@ -85,7 +85,7 @@ neuron iaf_psc_exp_neuron:
shape I_shape_in = exp(-1/tau_syn_in*t)
shape I_shape_ex = exp(-1/tau_syn_ex*t)
function I_syn mV = (curr_sum(I_shape_in, in_spikes) + curr_sum(I_shape_ex, ex_spikes) + I_e + currents)
V_abs' = -1/tau_m * V_abs + 1/C_m * I_syn
V_abs' = -V_abs/tau_m + (I_syn/C_m) * 1 [s**3 * A**2 /(kg * m**2)] * 1 [s/Gs]

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The units here seem very confusing. I would expect I_syn and C_m to have their units defined elsewhere, so that the second term can be written without explicit units just as the first. And if units are needed, they should be given as high-level units with suitable prefactors (e.g. mV, pF) instead of basic MSKA and a G prefactor.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, this one is complete nonsense. It is likely an aggregate of several fixes. This all equates to a factor of 1 nS. Fixed in update.

Signed-off-by: Philip Traeder <philip.traeder@rwth-aachen.de>
Signed-off-by: Philip Traeder <philip.traeder@rwth-aachen.de>
Signed-off-by: Philip Traeder <philip.traeder@rwth-aachen.de>
Signed-off-by: Philip Traeder <philip.traeder@rwth-aachen.de>
@PTraeder
Copy link
Contributor Author

PTraeder commented Mar 14, 2017

@heplesser @Silmathoron @jougs I have adressed several requests regarding this PR. Can you look over the changes/current state of the models and give your Ok/further concerns for the PR?

@Silmathoron
Copy link
Member

I'm sorry, I cannot look at this before 3/23... I will as soon as possible after that date.

@DimitriPlotnikov
Copy link
Contributor

I'm going to merge this PR today evening, if there are no strong concerns regarding aginst this. Generally, this PR enforces type soundness of SI types. There is still lot room for improvement, but it already increases the quality of models significantly.

@DimitriPlotnikov
Copy link
Contributor

I'm merging the PR. Of course, the general discussion is not over.

@DimitriPlotnikov DimitriPlotnikov merged commit 279c571 into nest:master Apr 10, 2017
@@ -73,13 +73,13 @@ neuron iaf_chxk_2008_implicit:
internals:
# Impulse to add to DG_EXC on spike arrival to evoke unit-amplitude
# conductance excursion.
PSConInit_E real = 1.0 * e / tau_syn_ex
PSConInit_E real = 1.0 ms * e / tau_syn_ex
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm very sorry this comes so late (I must say I completely forgot this PR) but I just noticed the combination of lines 76/89 and 120/121 is redundant (therefore less clear) in terms of units: explicit units are used twice because PSConInit_I/E is real, whereas it could directly be nS.
This would lead to:

PSConInit_E nS = 1.0 nS * e / tau_syn_ex
PSConInit_I nS = 1.0 nS * e / tau_syn_ex

then

g_ex' += spikesExc * PSConInit_E
g_in' += spikesInh * PSConInit_I

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@Silmathoron it sounds reasonalbe. @PTraeder what is your opinion?

@@ -65,11 +65,11 @@ neuron iaf_cond_beta_neuron:
internals:
# Impulse to add to GE' on spike arrival to evoke unit-amplitude
# conductance excursion.
PSConInit_E real = 1.0 * e / tau_syn_rise_E #
PSConInit_E real = 1.0 * e / tau_syn_rise_E * 1 ms
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same as before: typing this as nS would make things simpler

@@ -61,11 +61,11 @@ neuron iaf_cond_exp_implicit:
internals:
# Impulse to add to g_ex' on spike arrival to evoke unit-amplitude
# conductance excursion.
PSConInit_E real = 1.0 * e / tau_synE
PSConInit_E real = 1.0 ms * e / tau_synE
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

See above

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In this model, the PSConInit variables are not used. Is this intended?

@@ -150,8 +150,8 @@ neuron terub_neuron_gpe:
end

internals:
PSCurrInit_E real = 1.0 * e / tau_syn_ex
PSCurrInit_I real = 1.0 * e / tau_syn_in
PSCurrInit_E real = 1.0 ms * e / tau_syn_ex
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

See above but with pA

@PTraeder
Copy link
Contributor Author

@Silmathoron I implemented your propositions in #346

DimitriPlotnikov pushed a commit that referenced this pull request Apr 18, 2017
Implement changes proposed in discussion of #343
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

5 participants