# Ivan Raikov

I am a researcher at Stanford University.

You can reach me at ivan (dot) g (dot) raikov (at) gmail (dot) com.

## Software

I have a Github page.

Most of my recent software is written in Scheme for the Chicken Scheme compiler, or in Standard ML for the MLton compiler.

A list of Scheme libraries that I have developed and/or maintain is available on my Chicken user page.

Previous work in which I am no longer involved includes MRCI and RTXI.

## Literature

### Master and Margarita

The Master and Margarita by Mikhail Bulgakov
original in Russian
The Master and Margarita by Mikhail Bulgakov
1967 translation by Michael Glenny, based on the Soviet-censored Russian text.
The Master and Margarita by Mikhail Bulgakov
1997 translation by Richard Pevear and Larissa Volokhonsky. Based on the uncensored original.
Bulgakov's Master and Margarita
An annotation to Bulgakov's novel, created by Kevin Moss at Middlebury College.

### Atanas Dalchev

Collected works of Bulgarian poet Atanas Dalchev

### Yordan Yovkov

The Harvester, a novel by Bulgarian author Yordan Yovkov.

## Notes

Arithmetic Expression Syntax
posted on 2014-04-08

# Definitions

We define an arithmetic expression as a combination of variables, numerical values and operations over these values. For example, the arithmetic expression 2 + 3 uses two numeric values 2 and 3 and an operation +. The syntax of an expression specifies the allowed combinations of expressions. Each expression has a meaning (or value), which is defined by the evaluation of that expression. Evaluation is a process where expressions composed of various components get simplified until eventually we get a value. For example, evaluating 2 + 3resultsin5.

# Syntax

The conceptual structure of an expression is called the abstract syntax. The particular details and rules for writing expressions as strings of characters is called the concrete syntax, e.g. MathML. The abstract syntax for arithmetic expressions is very simple, while the concrete syntax can include additional punctuation and other lexical features.

## Abstract Syntax

The abstract syntax of expressions can be represented by the following EBNF grammar:

Expr ::=   Number   Real
| BinOp    B Expr Expr
| UnOp     U Expr

B ::= + | - | * | / | ^

U ::= neg | abs | atan | asin | acos |  sin | cos | exp | ln |
sqrt | tan | cosh | sinh | tanh | gamma | lgamma | log10 | log2
where Number, BinOp, UnOp indicate the different kinds of expressions: numeric value, binary operation and unary operation, respectively, and the set of binary and unary operations are the the usual arithmetic operators on real numbers.

# Variables and Substitution

Arithmetic expressions contain variables in addition to constants and arithmetic operations. Thus we will need to extend the abstract syntax of expressions to include variables. We can represent variables as an additional kind of expression. The following data definition modifies Expr to include a Variable case:

Expr ::= Number     Real
| BinOp    B Expr Expr
| UnOp     U Expr
| Variable Symbol

(B and U as before)

Further, in order to assign a meaning to a variable in a particular context, we will define the association of a variable x with a value v as a binding, which can be written as x \mapsto v. Bindings can be represented as pair. For example, the binding of x \mapsto 5 can be represented as (symbol\;"x",\;5).

## Substitution

The substitution operation replaces a variable with a value in an expression. Here are some examples of substitution:

1. substitute (x \mapsto 5) in (x + 2) \rightarrow{} (5 + 2)
2. substitute (x \mapsto 5) in (2) \rightarrow{} (2)
3. substitute (x \mapsto 5) in (x * x + x) \rightarrow{} (5 * 5 + 5)
4. substitute (x \mapsto 5) in (x + y) \rightarrow{} (5 + y)

If the variable names don't match, they are not substituted. Given the syntax defined for expressions, the process of substitution can be defined by cases:

simpleSubstitute (var, newVar) exp = subst exp

where

subst (Number r)      = Number r
subst (BinOp B a b)   = BinOp B (subst a) (subst b)
subst (UnOp U a b)    = UnOp U (subst a) (subst b)
subst (Variable name) = if var == name
then Variable newVar
else Variable name

## Environments

There can be multiple variables in a single expression. For example, evaluating (2 * x + y) where x = 3 and y = -2 results in 4.

A collection of bindings is called an environment. Since a binding is represented as a tuple, an environment can be represented as a list of tuples. The environment from the example above would be

env = [ ("x", 3), ("y", -1) ]

An important operation on environments is variable lookup. Variable lookup is an operation that given a variable name and an environment looks up that variable in the environment. For example:

1. lookup (x) in (env)\rightarrow3
2. lookup (y) in (env)\rightarrow-1

The substitution function can then be redefined to use environments rather than single bindings:

substitute env exp = subst exp where

subst (Number r)      = Number r
subst (BinOp B a b)   = BinOp B (subst a) (subst b)
subst (UnOp U a)      = UnOp U (subst a)
subst (Variable name) = case lookup name env of
Some r => Number r
| None  => Variable name

## Local Variables

It is also useful to allow variables to be defined within an expression. We will define local variables with a let expression:

let x = 1 in 2*x + 3

Local variable declarations are themselves expressions, cand can be used inside other expressions:

2 * (let x = 3 in x + 5)

Local variable declarations can be represented in the abstract syntax by adding another clause:

Expr ::= ...
| Let Symbol Expr Expr

When substituting a variable into an expression, it must correctly take into account the scope of the variable. In particular, when substituting for x in an expression, if the expression is of the form let x = e in body then x should be substituted within e but not in body.

simpleSubstitute (var, val) exp = subst exp

subst (Let x exp body)  = Let x (subst exp) body1

where body1 = if x == var
then body
else subst body

In the Let case for subst, the variable is always substituted into the bound expression e. But the substitution is only performed on the body b if the variable var is not the same as x.

TODO

# Summary

Here is the complete definition, substitution and evaluation using an expression language with local variables.

Expr ::=   Number   Real
| BinOp    B Expr Expr
| UnOp     U Expr
| Variable Symbol
| Let      Symbol Expr Expr

simpleSubstitute (var, val) exp = subst exp where

subst (Number r)      = Number r
subst (BinOp B a b)   = BinOp B (subst a) (subst b)
subst (UnOp U a b)    = UnOp U (subst a) (subst b)

subst (Variable name) = if var == name
then Number val
else Variable name

subst (Let x exp body)  = Let x (subst exp) body1

where body1 = if x == var
then body
else subst body

evaluate (Number r)       = r
evaluate (BinOp + a b)    = (evaluate a) + (evaluate b)
...
evaluate (UnOp sqrt a)    = sqrt (evaluate a)
...
evaluate (Let x exp body) = evaluate (simpleSubstitute (x, evaluate exp) body)
evaluate (Variable x)     = Error
An algorithm for simulating neuronal networks
posted on 2013-06-04

The simulation of a neuronal network model with N neurons can be represented as a sequence of transformations. We assume the constant matrix W of size N by N contains the synaptic weights of each connection in the network. We further assume that the global delay of synaptic events is represented by a delay matrix D of dimensions N by T, where T is the global delay as number of time steps.

• The first transformation computes the current state of all neurons given an initial state vector A^{t-1} of size N and initial input vector I^{t-1}, also of size N,
A_{i}^{t} = f_{i} (A_{i}^{t-1}, I_{i}^{t-1})
where f is the function that computes neuronal and synaptic dynamics.
• The next transformation computes the set S of neurons that have spiked during the current time step,
S = \forall{i}(spike(A_{i}^{t}) == true)
where spike is a function that determines whether a spike has ocurred given the current state of a neuron.
• The next transformation computes weighted sums of the inputs,
J_{j} = \sum_{i \in S} W_{i,j} spikecount(A_{j}^{t})
where J is a of the same dimensionality as I, and spikecount is a function that returns the number of spikes that have occurred given the current state of a neuron.
• The final transformation computes the new delay matrix and input vector,
I_{i}^{t} = D_{i,1}^{t-1}
D^{t} = nshift (cat (D, J, 2), 1)
where function cat concatenates J to the right side of D, and function nshift shifts the columns of the resulting matrix to the left by one.

The algorithm outlined above assumes a single type of synaptic connection. This restriction can be lifted by extending W to be of dimensions N by N by K, where K is the number of synaptic connection types, and extending I to be of dimensions N by K. The assumption that W is constant can be lifted by introducing a function that computesthe new value of W for each time step.

An algorithm for solving chemical mechanisms using sparse matrices
posted on 2013-04-26

Sparse Matrix Technique for solving chemical kinetic equations

The solution of any chemical reaction mechanism requires the use of the stoichiometric coefficients to calculate the right-hand sides of the kinetic equations. If the number of species participating in any chemical reaction is always much less than the total number of chemical species involved, then the matrices of stoichiometric coefficients will be sparse. A sparse matrix is one in which there are many more zero elements than nonzero ones. In such a case, one can exploit the prevalence of zero matrix elements with a sparse matrix technique in which only nonzero matrix elements are stored and processed.

This technique involves three one-dimensional arrays: two integer index arrays, IA, JA, and a real one, A. Nonzero elements of the original matrix M are stored row-by-row in A. To identify the nonzero elements in a row, we need to know the column of each entry. The array JA is used to store the column indices which correspond to the nonzero entries of M, i.e., if A(k)=M(i,j), then JA(k)=j.

We also need to know the number of nonzero elements in each row. The index positions in JA and A where each row of M begins are stored in the IA array. If M(i,j) is the leftmost entry of the i-th row and A(k)=M(i,j), then IA(k)=k. Moreover, IA(N+1) is defined as the index in JA and A of the first location following the last (N-th), element in the row. Thus, the number of entries in the i-th row is given by the difference IA(i+1)-IA(i), and the nonzero elements of the i-th row are stored in a sequence

A(IA(i)), A(IA(i)+1), \ldots, A(IA(i+1)-1)
while the corresponding column indices are stored in
JA(IA(i)), JA(IA(i)+1), \ldots, JA(IA(i+1)-1)

An algorithm for computing the compartmental model equations
posted on 2013-03-12

An algorithm for computing the equations of a compartmental neuron model for a general tree structure:

• For each cylinder i with radius and length a_i and L_i in micrometers, let the surface area,
A_i = 2 \pi a_i L_i
and the axial resistance,
Q_i = L_i / (\pi a_i^2)
• Given specific membrane capacitance and resistance c_i and r_{m_i}, respectively, the compartmental capacitance is
C_i = c_i A_i 10^{-8}
and the compartmental resistance is
R_i = (r_{m_i} / A_i) 10^8
• Given longitudinal (intracellular) resistivity r_l, the coupling resistance between compartments i and j is
R_{ij} = 0.5 r_l (Q_i + Q_j) 10^4
• The membrane potential equation for compartment i is then
C_i \frac{dV_i}{dt} = -\frac{V_i}{R_i} + \sum_{j,i} \frac{(V_j - V_i)}{R_{ij}} + I_i A_i {10^{-2}}

The factors of 10^-8, 10^8 and 10^4 are conversion factors from micrometers to centimeters.

The main idea of Brep is to combine experimental data on cell morphology and synthetic cell geometry. Cell geometry is described as sets of 3D points and can be one of the following:

Neurolucida data
cylinders described by start point, end point, diameter
Synthetic geometry
points along a line, polynomial curve, or angular displacement

Cell placement can be defined by rectangular grids, regular tilings (e.g. bricks or hexagons), random point distributions. Connectivity currently is simply projections based on Euclidean distances between cell objects.

### Layers

In Brep terminology, a layer describes the placement of cell objects in 3D space. The placement is determined by a set of points, which can be generated by one of the following methods:

1. sampling from typical distributions such as uniform or exponential
2. point coordinates loaded from a file
3. Regular grid defined by spacing in x, y and z direction
4. Regular 2D tiling, such as hexagons or bricks with a given side length. In this case, position x and y coordinates are determined by the center point of each generated tile, and the z coordinate is set to zero.

A layer is defined by a point set described by one of the methods above, and one or more cell objects that are to be placed at positions determined by the point set. If there are less cell objects than there are points in the point set, then the cell objects are copied until an object is placed at each point.

### Cell objects

A cell object consists of a bounding box and a set of connection points. The connection points can be loaded from a Neurolucida file, or can be sampled from the trajectory of a polynomial curve. Currently the types of polynomial curves supported are lines and angular displacements of the form L * cos (theta) * cos (phi) where theta and phi are chosen randomly from a given range of angles.

The set of points that belong to a cell object can be further divided into compartments, e.g. dendrites and axon, which can be later used when specifying connectivity.

In addition, random perturbations can be applied to the points of a cell or one of its compartments.

Parallel fibers can then be described as lines with small random perturbations, and artificial Golgi cell morphologies can be desribed with angular displacements of their dendritic and axon points.

### Projections

Once layers are defined with the different cell types, connectivity is specified in Brep as projections between layers.

The single type of projection currently supported by BRep is based on Euclidean distances between the points of the different cell object. Given a maximum distance, BRep can determine all cell points in two layers that are close enough to be considered connected.

Older articles