## Development and release log

**Febrary 11 2008**: bug in the stochastic calcluation voltage interpolation.
A variable name error in the preprocessing stage for the stochastic calculation meant that the
interpolation for transition rates was always using the the rate for the lower end of a voltage interval
rather than the linearly interpolated pointe between the two ends. This is unlikely to have had a
significant effect on spiking models, since the statistics smooth out as the voltage changes, but
does affect voltage clamped models where the applied rates could correspond to a membrane potential
that is up to 1mV away from the actual potential. The problem only applied to the stochastic calculation
since this does not use the transition matrices themselves, but rather sorted cumulative transition
probabilities for the two ends of each interval (which have to be sorted the same way). The upper end
array was getting the same values as the lower end.

**January 11 2008 3D visualization**
Java 3D visualiztion option for cells and channels: details to come in icing documentation.

**January 2008 Channel positioning**
The way channels are allocated has changed. This is a fairly large change that has affected al
range of other components. Previously the cell was discretized, the metrics were computed on the
discretized version and and channels were allocated to compartments accordingly.
Now, the channel allocation is independent of the discretization. Each channel has an exact position
on the cell. For the computation, the discretization is performed after the channel allocation and
channels are allocated to compartments according to their positions. The key benefit is that now the
cell specification is complete in itself, and has no relation to any eventual disretization. It also
allows for more intuitive visualization of the model and permits the use of distribution rules with
structure on scales smaller than the discretization. For example, you can can treat a whole cell as
three compartments to see what effect it has and still be sure that the number of channels is correct.

The previous method is still available via the "quickChannels" attribute on PSICSRun. Along with this change are a number of modifications and extensions to the ways channels are distributed as documented in the Channel Distribution section of the formats reference.

**October 15 Temperature dependence**
All KSTransition objects now support temperture dependence of the rates via
attributes q10 and baseTemperature. Rates are multiplied by q10^{x}
where x = (temperature - baseTemperature) / 10. Computed transitions can use this
form or can implement their own temperature dependence internally as before.

**October 15 Sign conventions**
Changed the convention for the sign of the scale attribute ion ExpLinearTransition
and SigmoidTransition. Formerly, the expression for the first was
rate = A x / (exp(x) - 1) where x = (v - vhalf) / scale.
Now it is rate = A x / (1 - exp(-x)) which has an extra minus sign, but means that a
positive scale corresponds to a rate increases with deploarization. Likewise for the sigmoid
expression. The HH models have been updated accordingly.

**October 10 Alternative morphology formats** The main PSICSRun element now
accepts MorphologySource sub-elements to specify external morphology sources.
The first supported format is the .swc format as avaiable on neuromorpho.org. This provides
an alternative to the morphology attribute in PSICSRun: supplying a MorphologySiource element
pointing to the swc file is equivalent to directly pointing to a morphology format in
PSICS native format.

**October 6 Time series data for stimulation** There is a new element TimeSeries
that can be included inside a clamp definition to specify an file from which the stimulation
data should be red. The file should just contain two columns: the time and the values. In the
Time Series element you must specify what the base units are for each column so PSICS knows
how to rescale them for its calculations. It is not necessary for the times to be at
uniform increments so, for example, long stretches at the same value can be left out of the
file completely. If two successive time values are the same it is interpreted as a sudden
step from one value to the other. Currently, the contents of the file are read during the
preprocessing stage, rescaled and written to the input file for the core calculation. This
could need extending to allow the core to read the file itself if large input
sequences (more than a few tens of thousands of points) are used but the latter case is probably
best handled later within a general framework for external input.

**October 4 Discretization of channel densities**
The problem to be addressed is that channel distribution rules give a floating point number of
for each compartment, but the calculation requires a whole number of channels. Worse,
balance conditions that adjust channel densities to acheive a stable equilibrium potential
may generate negative densities on some compartments. The current solution provides a deterministic
method for adjusting densities into the most appropriate non-negative values on each
compartment. It starts by propagating negative densities inward from terminal points until
they are cancelled out by positive channel densities. Then it propagates in from the terminals
again integerizing the densities and propagating the remainders up the tree. The remainder
arriving at the root point is rounded to the nearest whole number of channels and a warning
issued if it is negative. This allows large conductance channels to be used as leaks in the
channel balance equtions when studying the effects of stochasticity.

**September 28 Step transitions** There are now two ways of handling step transitions to
reduce possible artifacts when timestep boundaries do not fall exactly on step boundaries.
The choices are either to evaluate the command profile at the middle of a timestep and use that
value for the whole step, or to apply the average value of the command profile over a timestep.
These can be specified with the stepType attribute on voltage and current clamps which should
take the value "midpoint" or "average". In general, the average approach is best for current
clamps (it gives the right charge delivery, irrespective of the timestep) and the midpoint
method should be applied to voltage clamps since a step to an imtermediate potential is generally
not equivalent to a shorter step to the full step potential. The cost is some uncertainty about
effective onset times (plus or minus half a timestep) but this is almost inevitable with
fixed step methods. There are examples of the two cases and further explanation in the
two models:
step-onset and
step-onset-avg.

**September 13 Initial equilibrium** A new feature in channel distribution specifications
allows the densities of one or two channel types to be adjusted so as to achieve a steady resting
potential. This is typically best done with sodium and potassium leaks, and provides a more phisiologically
relevant alternative to the approach of simply adjusting hte reversal potential of the leak on a compartment
by compartment basis.

**August - Formal documentation** Most features are now documented in the user guide
and refernce pages, including recent developments on core efficieny and model specification.
We have reversed the previous decision about supporting arbitrary channel models: it is now possible to
embed code fragments in CDATA sections of the channel specification to implement ad-hoc expressions for
transition rates. The main motivation here is to allow direct comparison with models from other systems.
It uses an embedded compiler and compiles on the fly from generated source code to make classes to evaluate
the rates. These are tabulated as before prior to being passed to the Fortran core.

**25 July - Porting models** The choice not to support arbitrary expressions for state transition
rates means that to re-implement existing models requires fitting a Boltzmann type transition to the legacy models.
The image shows the type of problem that can arise. In black are the forward and reverse transition rates as a
function of membrane potential (log rate in ms against mV). In red and blue are the best fits to these rates.
It remains to be seen if there are any functional differences in the behavior of the model.

**12 July - Efficient calculation of HH Models**
So far, HH models (multiple independent two-state gating complexes) have been computed either
as the single equivalent scheme, or as a collection of serial two-state gates. This is a little inefficient
in the continuous regime, since normally several gates have the same kinetics. This is now exploited in the
usual way with an "instances" attribute on a complex to give the power of the gating variable. In the stochastic
regime it makes no difference since it is necessary to follow each complex (or the single equivalent complex)
in detail in any case.

**4 July - Noise and Conductance inputs**
TBD

**6 June 2007 - Stochastic simulation and random numbers**. How many random numbers should it take to
advance a population of n channels? It can certainly be done with n random numbers. If, as is often the case,
for a given start state there is an overridingly popular destination state (probability p), then you only
need to use a few times n(1-p) random numbers - just those channels that have a chance of not being in the
default destination (often, but not always, the current state). But for some cases it should still be lower
than this: eg, first generate the number that change at all, then the number that aren't in the second
most popular destination, then you may well have none left at all. But the above two steps could cost more
than quite a lot of random numbers, so when is it worth it?

**28 May 2007 - KS Scheme equivalents of HH Channels**. Hodgkin-Huxley style models - multiple serial
independent two state complexes, frequently with several sharing the same properties - are a simple extension
of single-complex schemes for the continuous case, but not in the stochastic case. For the former, the complexes
can be treated independently, just multiplying up the relative conductances at the end. But for a population of
stochastic channels it is not enough to know how many of each type of complex are in each state: you need
to know which complex of one type goes with which complex of the other type.
So the slow stochastic algorithm is fine, where the state of every complex
is recorded, but population based methods do not work.

The current solution is to **automatically** construct the minimal equivalent single scheme specification
from a HH-type model. As yet the procedure is specific to models of
type a b^{n} (one type a complex, n type b complexes) but that should cover all interesting
cases: for more sophisticated multi-complex models it is probably worth constructing a single-complex
equivalent by hand.

Predictable, the performance of single-scheme equivalents is not as good in the continuous regime as the multi-complex version (one 8x8 matrix by vector multiplication costs more than two or even four 2*2 matrix by vector multiplications), so the multi-complex option has been left in for continuous calculations.

**May 2007 - Fortran core**. Now running, and producing the same results as the Java reference at
something between 15 and 50 times faster (though the comparison is not just on the language since they
use different methods). The main observation, performance-wise, is that the elegant new (well, newish) features
of Fortran 95 are not much use in inner loops: data parallel statements and pointer arrays are substantially
slower (using a single core) than their old fashioned equivalents both with GFortran and Absoft. Presumably
this is from all the additional copying where "A = A + b" leads to a copy of A being incremented and assigned
back to A, instead of the element-wise increment in a do loop. Likewise "forall" is much more costly
in the PSICS inner loops than explicit "do" loops. Rather disappointing really, since the data parallel
versions are much more readable! And if only Fortran would introduce a "+=" construct to avoid writing
complicated index lists twice
("pop%gcs(igc)%states_stoch(p, icpt) = pop%gcs(igc)%states_stoch(p,icpt) + 1" is just something you
shouldn't have to write). Pointers would solve the problem, but at the cost of setting the target
attribute on the arrays used in inner loops which is probably not smart. On the other hand, automatic
allocation (calling a subroutine and declaring that new arrays should be available inside it with sizes
that depend in its arguments) is pretty nice.

**30 April 07 - Multiple runs and parameter access**. The RunSet block now has access to the
whole model and allows embedded CommandConfig blocks to do two-dimensional parameter variation.
The default behavior is to process a family of commands, such as a sequence of voltage steps,
through a single model. Steps are specified in the same way as the RunSet itself by identifying
the parameter to vary and setting the values it should take in successive cases. For the RunSet,
the target parameter can now be anywhere in the model. It is sufficient to set the id on the component
containing the parameter and then to target it from the RunSet with a target of the form
id:field where id is the id of the target component and field is the name of the quantity
to be varied. This is mainly implemented in ModelMap.java, which could easily be extended to
provide unix directory style access to the model - if a need arises.

**26 Aripl 2007 - JNI and Fortran** There is now an embryonic version of FPSICS (Fortran PSICS)
that provides an entry point form Java and has access to a fortran-style serialization of the
pre-processed task specification. The outer-core generates the internal structures ready to start
the calculation (discretizes the cell, tabulates channel transition rates etc) and writes them to
a text file in a "sequentially readable" style - that is, such that a line-based reader
can reinstantiate the whole lot without stacking or backtracking (which has a lot going for it in fact).
No lines mixing doubles and ints, dimensions declared in advance rather than deduced from
the context etc. The latter could be fragile for manual editing, but these things shouldn't be edited by
hand.

On an implementation level, Photran with CDT (the Fortran and C development systems for Eclipse respectively) does the trick nicely for calling Fortran from Java. The process is to declare the native method in Java, generate a c header file, write a c wrapper function that calls fortran, write the fortran and then join them together, putting in the name-mangling and string mappings that are needed between the two until there are no linkage errors. With a bit of configuration, the Fortran build can also compile the c and make a shared library with the c as entry point, which is all that java needs. Not much fun, but once it works that should be it.

**16 April 07 - Numerical Methods**. There is now a "method" attribute to the main run configuration
that can take the values "implicit_euler", "crank_nicolson", "weighted_crank_nicolson" (the default) and
"forward_euler". These only vary in the values of the temporal difference weighting factor
(1.0, 0.5, 0.51, and 0.0 respectively). This factor can also be set explicitly with the tdWeighting
attribute, which is mainly intended for exploring the dependence of the results on the numerical
method, rather than for normal modeling applications. However, it may also allow better
accuracy with large timesteps for quick approximate calculations without the normal Crank Nicolson
oscillations by picking values in the range 0.5 to 1.

**11 April 07 - Validation site** The results of rallpack calculations and any
other models in the validation source tree are now published to validation.psics.org.
The SiteBuilder class in org.psics.validation regenerates all the html.

**April 07 - Rallpack**
Creating equivalents of the rallpack models required a couple of
extensions to the preprocessing step. The first is completion of the
implementation of the "minor" flag for points in the structure. This
can now be applied to any point and specifies that the segment from
its parent is a cylinder of constant radius, rather than tapered as is the
default. The discretization also pays attention to minor segments and
treats them as the beginning of a new compartment rather than splitting
the parent compartment some distance along the segment as it would for a
tapered segment. This allows the branching "tin-can" structure of Rallpack 2
to be implemented by declaring all the points in the structure except the first,
to be minor.

The Other extension was to add two more possible transition types to cover the expressions used in Rallpack 3 for the HH equations. These are the Exponential, A exp ((V - V0) / B) and Sigmoid, A / (1 + exp ((V - V0) / B) transitions.

**March 30 07 - Mesh generation 2**
The comments below (February 07) turn out to be misguided. The "remeshing problem" is
only a problem if you require compartments to be straight cylinders. But there
is no reason why they should be. Removing this constraint and returning to the
numerical objectives yields a simple solution. Just use the structure exactly as it
is and discretizes into subsections according to the needs of the numerical method
and accuracy constraints (just as one would for any other pde problem...). The result
is that you get compartments with bends or even branches in them, so there is a slight
increase in computational bookkeeping to compute areas and axial resistances, and
a bit more work to draw them. But these are all straightforward processes and
implementing them has the double benefits of not changing the structure for the convenience of the
calculation, and not tying the hands of the calculation by the details of the
structure.

The two images illustrate the results of using different discretization criteria on a section of dendrite. In each case, some compartments comprise multiple points from the original structure, but the shape is preserved.

**March 07 - Java inner core**
First complete runs using channel models and extended cells. The stochastic calculation is
only channel-by-channel, but gives quantitative agreement with the ensemble method.
Channel distribution is handled with an external expression parser - jep - which is wrapped
in an Evaluator class so alternatives can be easily switched in if necessary (a graph reduction
parser might be easier to modify for mixing extra morphology conditions in with the expression).
Output is currently only ad-hoc text for viewing with CCViz

**2 February 07 - Mesh generation**
Probably the trickiest part of preprocessing morphologies is the compartmentalization since it presents a
number of potentially incompatible requirements. Between two successive branch points one would like
to change then number of segments, but preserve: end point positions; total surface area;
total volume; axial resistance; radius profile with path-length; position profile. In general they can't
all be satisfied (except for refining-only discretizations). The simplest case is for a wiggly section
of dendrite that should be approximated by a few longer segments. Just evening out the wiggles will give
too low a total length, capacitance, etc. The optimum for most criteria is likely to be a zig-zag between
the end points.

The current strategy has two steps:

- decide on the number of segments and interpolate between the two branch points to position them on the profile.
- iterative refinement to adjust point positions and optimize some function of the possible metrics to achieve the best fit under given criteria.

**28 January 07 - Neuron Imports** Two new packages, org.psics.models.neuron and ...neuron.lc (Lower Case)
for importing morphologies from neuron. They contain mostly small java classes that are instantiated
by reflection on reading a neuron file. The "getFinal" methods generate standard pSICS morphologies, enabling
the formats to be read without any specific handling code other than to recognize which packages should be
searched when parsing these formats.

**January 07** Format specifications and website.