rxd.cpp

Go to the documentation of this file.

#include <../../nrnconf.h>
#include <stdio.h>
#include <assert.h>
#include <string.h>
#include "grids.h"
#include "rxd.h"
#include <../nrnoc/section.h>
#include <../nrnoc/nrn_ansi.h>
#include <../nrnoc/multicore.h>
#include <nrnwrap_Python.h>
#include <nrnpython.h>

static void ode_solve(double, double*, double*);
extern PyTypeObject* hocobject_type;
extern int structure_change_cnt;
extern int states_cvode_offset;
extern int _nrnunit_use_legacy_;
int prev_structure_change_cnt = 0;
int prev_nrnunit_use_legacy = _nrnunit_use_legacy_;
unsigned char initialized = FALSE;

/*
    Globals
*/
extern NrnThread* nrn_threads;
int NUM_THREADS = 1;
pthread_t* Threads = NULL;
TaskQueue* AllTasks = NULL;


extern double *dt_ptr;
extern double *t_ptr;


fptr _setup, _initialize, _setup_matrices, _setup_units;
extern NrnThread* nrn_threads;

/*intracellular diffusion*/
unsigned char diffusion = FALSE;
int _rxd_euler_nrow=0, _rxd_euler_nnonzero=0, _rxd_num_zvi=0;
long* _rxd_euler_nonzero_i = NULL;
long* _rxd_euler_nonzero_j = NULL;
double* _rxd_euler_nonzero_values = NULL;
long* _rxd_zero_volume_indices = NULL;
double* _rxd_a = NULL;
double* _rxd_b = NULL;
double* _rxd_c = NULL;
double* _rxd_d = NULL;
long* _rxd_p = NULL; 
unsigned int*  _rxd_zvi_child_count = NULL;
long** _rxd_zvi_child = NULL;
static int _cvode_offset;
static int _ecs_count;

/*intracellular reactions*/
ICSReactions* _reactions = NULL;

/*Indices used to set atol scale*/
SpeciesIndexList* species_indices = NULL;

/*intracellular reactions*/
double* states;
unsigned int num_states=0;
int _num_reactions = 0;
int _curr_count;
int* _curr_indices = NULL;
double* _curr_scales = NULL; 
double** _curr_ptrs = NULL;
int  _conc_count;
int* _conc_indices = NULL;
double** _conc_ptrs = NULL;

/*membrane fluxes*/
int _memb_curr_total = 0;   /*number of membrane currents (sum of 
                              _memb_species_count)*/
int _memb_curr_nodes = 0;   /*corresponding number of nodes
                              equal to _memb_curr_total if Extracellular is not used*/
int _memb_count = 0;        /*number of membrane currents*/
double* _rxd_flux_scale;    /*array length _memb_count to scale fluxes*/
double* _rxd_induced_currents_scale;    /*array of Extracellular current scales*/
int* _cur_node_indices;      /*array length _memb_count into nodes index*/


int* _memb_species_count;   /*array of length _memb_count
                             number of species involved in each membrane
                             current*/

/*arrays of size _memb_count by _memb_species_count*/
double*** _memb_cur_ptrs; /*hoc pointers TODO: replace with index for _curr_ptrs*/
int** _memb_cur_charges;    
int*** _memb_cur_mapped;      /*array of pairs of indices*/
int*** _memb_cur_mapped_ecs;  /*array of pointer into ECS grids*/

double* _rxd_induced_currents = NULL;       /*set when calculating reactions*/
ECS_Grid_node** _rxd_induced_currents_grid = NULL;

unsigned char _membrane_flux = FALSE;   /*TRUE if any membrane fluxes are in the model*/
int* _membrane_lookup; /*states index -> position in _rxd_induced_currents*/

int _node_flux_count = 0;
long* _node_flux_idx = NULL;
double* _node_flux_scale = NULL;
PyObject** _node_flux_src = NULL;

static void free_zvi_child()
{
    int i;
    /*Clear previous _rxd_zvi_child*/
    if(_rxd_zvi_child != NULL && _rxd_zvi_child_count != NULL)
    {
        for(i = 0; i < _rxd_num_zvi; i++)
            if(_rxd_zvi_child_count[i]>0) free(_rxd_zvi_child[i]);
        free(_rxd_zvi_child);
        free(_rxd_zvi_child_count);
        _rxd_zvi_child_count = NULL;
        _rxd_zvi_child = NULL;
    }
}

static void transfer_to_legacy()
{
    /*TODO: support 3D*/
    int i;
    for(i = 0; i < _conc_count; i++)
    {
        *(double*)_conc_ptrs[i] = states[_conc_indices[i]];
    }
}

static inline void* allocopy(void* src, size_t size)
{
    void* dst = malloc(size);
    memcpy(dst, src, size);
    return dst;
}

extern "C" void rxd_set_no_diffusion()
{
    int i;
    diffusion = FALSE;
    if(_rxd_a != NULL)
    {
        free(_rxd_a);
        free(_rxd_b);
        free(_rxd_c);
        free(_rxd_d);
        free(_rxd_p);
        free(_rxd_euler_nonzero_i);
        free(_rxd_euler_nonzero_j);
        free(_rxd_euler_nonzero_values);
        _rxd_a = NULL;
    }
}

extern "C" void free_curr_ptrs()
{
    _curr_count = 0;
    if(_curr_indices != NULL)
        free(_curr_indices);
    _curr_indices = NULL;
    if(_curr_scales != NULL)
        free(_curr_scales);
    _curr_scales = NULL;
    if(_curr_ptrs != NULL)
        free(_curr_ptrs);
    _curr_ptrs = NULL;
}

extern "C" void free_conc_ptrs()
{
    _conc_count = 0;
    if(_conc_indices != NULL)
        free(_conc_indices);
    _conc_indices = NULL;
    if(_conc_ptrs != NULL)
        free(_conc_ptrs);
    _conc_ptrs = NULL;
}


extern "C" void rxd_setup_curr_ptrs(int num_currents, int* curr_index, double* curr_scale,
                          PyHocObject** curr_ptrs)
{
    int i;
    free_curr_ptrs();
    /* info for NEURON currents - to update states */
    _curr_count = num_currents;
    _curr_indices = (int*)malloc(sizeof(int)*num_currents);
    memcpy(_curr_indices,curr_index, sizeof(int)*num_currents); 

    _curr_scales = (double*)malloc(sizeof(double)*num_currents);
    memcpy(_curr_scales, curr_scale, sizeof(double)*num_currents); 

    _curr_ptrs = (double**)malloc(sizeof(double*)*num_currents);
    for(i = 0; i < num_currents; i++)
        _curr_ptrs[i] = (double*)curr_ptrs[i]->u.px_;
}

extern "C" void rxd_setup_conc_ptrs(int conc_count, int* conc_index, 
                         PyHocObject** conc_ptrs)
{
    /* info for NEURON concentration - to transfer to legacy */
    int i;
    free_conc_ptrs();
    _conc_count = conc_count;
    _conc_indices = (int*)malloc(sizeof(int)*conc_count);
    memcpy(_conc_indices, conc_index, sizeof(int)*conc_count); 

    _conc_ptrs = (double**)malloc(sizeof(double*)*conc_count);
    for(i = 0; i < conc_count; i++)
        _conc_ptrs[i] = (double*)conc_ptrs[i]->u.px_;
}

extern "C" void rxd_include_node_flux3D(int grid_count, int* grid_counts,
                                        int* grids, long* index, double* scales,
                                        PyObject** sources)
{
    Grid_node* g;
    int i = 0, j, k, n, grid_id;
    int offset = 0;

    for(g = Parallel_grids[0]; g != NULL; g = g->next)
    {
        if(g->node_flux_count > 0)
        {
            g->node_flux_count = 0;
            free(g->node_flux_scale);
            free(g->node_flux_idx);
            free(g->node_flux_src);
        }
    }
    if(grid_count == 0)
        return;
    for(grid_id=0, g = Parallel_grids[0]; g != NULL; grid_id++, g = g->next)
    {
#if NRNMPI
        if(nrnmpi_use && dynamic_cast<ECS_Grid_node*>(g))
        {
            if (grid_id == grids[i])
                n = grid_counts[i++];
            else
                n = 0;

            g->proc_num_fluxes[nrnmpi_myid] = n;
            nrnmpi_int_allgather_inplace(g->proc_num_fluxes, 1);

            g->proc_flux_offsets[0] = 0;
            for(j = 1; j < nrnmpi_numprocs; j++)
                g->proc_flux_offsets[j] =  g->proc_flux_offsets[j-1] +
                                      g->proc_num_fluxes[j-1];
            g->node_flux_count = g->proc_flux_offsets[j-1] + g->proc_num_fluxes[j-1];
            /*Copy array of the indexes and scales -- sources are evaluated at runtime*/
            if(n > 0)
            {
                g->node_flux_idx = (long*)malloc(g->node_flux_count*sizeof(long));
                g->node_flux_scale = (double*)malloc(g->node_flux_count*sizeof(double));
                g->node_flux_src = (PyObject**)allocopy(&sources[offset], n*sizeof(PyObject*));
            }

            for(j = 0, k = g->proc_flux_offsets[nrnmpi_myid]; j < n; j++, k++)
            {
                g->node_flux_idx[k] = index[offset + j];
                g->node_flux_scale[k] = scales[offset + j];
            }
            nrnmpi_long_allgatherv_inplace(g->node_flux_idx, g->proc_num_fluxes, g->proc_flux_offsets);
            nrnmpi_dbl_allgatherv_inplace(g->node_flux_scale, g->proc_num_fluxes, g->proc_flux_offsets);

            offset += n;
        }
        else
        {
            if(grid_id == grids[i])
            {
                g->node_flux_count = grid_counts[i];
                if(grid_counts[i] > 0)
                {
                    g->node_flux_idx = (long*)allocopy(&index[offset], grid_counts[i]*sizeof(long));
                    g->node_flux_scale = (double*)allocopy(&scales[offset], grid_counts[i]*sizeof(double));
                    g->node_flux_src = (PyObject**)allocopy(&sources[offset], grid_counts[i]*sizeof(PyObject*));
                }
                offset += grid_counts[i++];
            }
        }
#else
        if(grid_id == grids[i])
        {
            g->node_flux_count = grid_counts[i];
            if(grid_counts[i] > 0)
            {
                g->node_flux_idx = (long*)allocopy(&index[offset], grid_counts[i]*sizeof(long));
                g->node_flux_scale = (double*)allocopy(&scales[offset], grid_counts[i]*sizeof(double));
                g->node_flux_src = (PyObject**)allocopy(&sources[offset], grid_counts[i]*sizeof(PyObject*));
            }
            offset += grid_counts[i++];

        }
#endif
    }
}

extern "C" void rxd_include_node_flux1D(int n, long* index, double* scales,
                                        PyObject** sources)
{
    if (_node_flux_count != 0)
    {
        free(_node_flux_idx);
        free(_node_flux_scale);
        free(_node_flux_src);
    }
    _node_flux_count = n;
    if (n > 0)
    {
        _node_flux_idx = (long*)allocopy(index, n*sizeof(long)*n);
        _node_flux_scale = (double*)allocopy(scales, n*sizeof(double));
        _node_flux_src = (PyObject**)allocopy(sources, n*sizeof(PyObject*));
    }
}



void apply_node_flux(int n, long* index, double* scale, PyObject** source, double dt, double* states)
{
    long i, j;
    PyObject *result;
    PyHocObject *src;

    for(i = 0; i < n; i++)
    {
        if(index == NULL)
            j = i;
        else
            j = index[i];
        if(PyFloat_Check(source[i]))
        {
            states[j] += dt * PyFloat_AsDouble(source[i]) / scale[i];
        }
        else if(PyCallable_Check(source[i]))
        {
            /* It is a Python function or a PyHocObject*/
            if (PyObject_TypeCheck(source[i], hocobject_type))
            {
                src = (PyHocObject*)source[i];
                /*TODO: check it is a reference */
                if(src->type_ == PyHoc::HocRefNum)
                {
                    states[j] +=  dt * (src->u.x_) / scale[i];
                }
                else
                {
                    states[j] +=  dt * *(src->u.px_) / scale[i];
                }
            }
            else
            {
                result = PyEval_CallObject(source[i], NULL);
                if(PyFloat_Check(result))
                {
                    states[j] += dt * PyFloat_AsDouble(result) / scale[i];
                }
                else if(PyLong_Check(result))
                {
                    states[j] += dt * (double)PyLong_AsLong(result) / scale[i];
                }
                else if(PyInt_Check(result))
                {
                    states[j] += dt * (double)PyInt_AsLong(result) / scale[i];
                }

                else
                {
                    PyErr_SetString(PyExc_Exception, "node._include_flux callback did not return a number.\n");
                }
            }
        }
        else
        {
            PyErr_SetString(PyExc_Exception, "node._include_flux unrecognised source term.\n");
        }
    }
}


static void apply_node_flux1D(double dt, double *states)
{
    apply_node_flux(_node_flux_count, _node_flux_idx, _node_flux_scale, _node_flux_src, dt, states);
}

extern "C" void rxd_set_euler_matrix(int nrow, int nnonzero, long* nonzero_i,
                          long* nonzero_j, double* nonzero_values,
                          double* c_diagonal)
{
    long i, j, idx;
    double val;
    unsigned int k, ps;
    unsigned int* parent_count;
    /*free old data*/
    if(_rxd_a != NULL)
    {
        free(_rxd_a);
        free(_rxd_b);
        free(_rxd_c);
        free(_rxd_d);
        free(_rxd_p);
        free(_rxd_euler_nonzero_i);
        free(_rxd_euler_nonzero_j);
        free(_rxd_euler_nonzero_values);
        _rxd_a = NULL;
    }
    diffusion = TRUE;

    _rxd_euler_nrow = nrow;
    _rxd_euler_nnonzero = nnonzero;
    _rxd_euler_nonzero_i = (long*)allocopy(nonzero_i, sizeof(long) * nnonzero);
    _rxd_euler_nonzero_j = (long*)allocopy(nonzero_j, sizeof(long) * nnonzero);
    _rxd_euler_nonzero_values = (double*)allocopy(nonzero_values, sizeof(double) * nnonzero);

    _rxd_a = (double*)calloc(nrow,sizeof(double));
    _rxd_b = (double*)calloc(nrow,sizeof(double));
    _rxd_c = (double*)calloc(nrow,sizeof(double));
    _rxd_d = (double*)calloc(nrow,sizeof(double));
    _rxd_p = (long*)malloc(nrow*sizeof(long));
    parent_count = (unsigned int*)calloc(nrow,sizeof(unsigned int));    

    for(idx = 0; idx < nrow; idx++)
        _rxd_p[idx] = -1;

    for(idx = 0; idx < nnonzero; idx++)
    {
        i = nonzero_i[idx];
        j = nonzero_j[idx];
        val = nonzero_values[idx];
        if(i < j)
        {
            _rxd_p[j] = i;
            parent_count[i]++;
            _rxd_a[j] = val;
        }
        else if(i == j)
        {
            _rxd_d[i] = val;
        }
        else
        {
            _rxd_b[i] = val;
        }
        
    }

    for(idx = 0; idx < nrow; idx++)
    {
        _rxd_c[idx] =  _rxd_d[idx] > 0 ? c_diagonal[idx] : 1.0;
    }
    
    if(_rxd_num_zvi > 0)
    {
        _rxd_zvi_child_count = (unsigned int*)malloc(_rxd_num_zvi*sizeof(unsigned int));
        _rxd_zvi_child = (long**)malloc(_rxd_num_zvi*sizeof(long*));

        /* find children of zero-volume-indices */
        for(i = 0; i < _rxd_num_zvi; i++)
        {
            ps = parent_count[_rxd_zero_volume_indices[i]];
            if(ps == 0) 
            {
                _rxd_zvi_child_count[i] = 0;
                _rxd_zvi_child[i] = NULL;

                continue;
            }
            _rxd_zvi_child[i] = (long*)malloc(ps*sizeof(long));
            _rxd_zvi_child_count[i] = ps;
            for(j = 0, k = 0; k < ps; j++)
            {
                if(_rxd_zero_volume_indices[i] == _rxd_p[j])
                {
                    _rxd_zvi_child[i][k] = j;
                    k++;
                }
            }
        }
    }
    free(parent_count);
}

static void add_currents(double * result)
{
    long i, j, k, idx;
    short side;
    /*Add currents to the result*/
        for (k = 0; k < _curr_count; k++)
            result[_curr_indices[k]] += _curr_scales[k] * (*_curr_ptrs[k]);

    /*Subtract rxd induced currents*/
    if(_membrane_flux)
    {
        for(i = 0, k = 0; i < _memb_count; i++)
        {
            for(j = 0; j < _memb_species_count[i]; j++, k++)
            {
                for(side = 0; side < 2; side++)
                {            
                    idx = _memb_cur_mapped[i][j][side];
                    if (idx != SPECIES_ABSENT)
                    {
                        result[_curr_indices[idx]] -= _curr_scales[idx] * _rxd_induced_currents[k];
                    }
                }
            }
        }
    }
}
static void mul(int nnonzero, long* nonzero_i, long* nonzero_j, const double* nonzero_values, const double* v, double* result) {
    long i, j, k;
    /* now loop through all the nonzero locations */
    /* NOTE: this would be more efficient if not repeatedly doing the result[i] lookup */
    for (k = 0; k < nnonzero; k++) {
        i = *nonzero_i++;
        j = *nonzero_j++;
        result[i] -= (*nonzero_values++) * v[j];
    }
}

extern "C" void set_setup(const fptr setup_fn) {
    _setup = setup_fn;
}

extern "C" void set_initialize(const fptr initialize_fn) {
    _initialize = initialize_fn;
    set_num_threads(NUM_THREADS);
}

extern "C" void set_setup_matrices(fptr setup_matrices) {
    _setup_matrices = setup_matrices;
}

extern "C" void set_setup_units(fptr setup_units) {
    _setup_units = setup_units;
    _setup_units();
}

/* nrn_tree_solve modified from nrnoc/ldifus.c */
static void nrn_tree_solve(double* a, double* b, double* c, double* dbase, double* rhs, long* pindex, long n, double dt) {
    /*
        treesolver
        
        a      - above the diagonal 
        b      - below the diagonal
        c      - used to define diagonal: diag = c + dt * dbase
        dbase  - together with c, used to define diagonal
        rhs    - right hand side, which is changed to the result
        pindex - parent indices
        n      - number of states
    */
    long i, pin;
    double* d = (double*)malloc(sizeof(double) * n);
    double* myd;
    double* myc;
    double* mydbase;

    /* TODO: check that d non-null */
    
    /* construct diagonal */
    myd = d;
    mydbase = dbase;
    myc = c;
    for (i = 0; i < n; i++) {
        *myd++ = *myc++ + dt * (*mydbase++);
    }

    /* triang */
    for (i = n - 1; i > 0; --i) {
        pin = pindex[i];
        if (pin > -1) {
            double p;
            p = dt * a[i] / d[i];
            d[pin] -= dt * p * b[i];
            rhs[pin] -= p * rhs[i];
        }
    }
    /* bksub */
    for (i = 0; i < n; ++i) {
        pin = pindex[i];
        if (pin > -1) {
            rhs[i] -= dt * b[i] * rhs[pin];
        }
        rhs[i] /= d[i];
    }
    
    free(d);
}



static void ode_solve(double dt, double* p1, double* p2)
{
    long i, j;
    double* b = p1 + _cvode_offset;
    double* y = p2 + _cvode_offset;
    double* full_b, *full_y;
    long* zvi = _rxd_zero_volume_indices;
   
    if(_rxd_num_zvi > 0)
    {
        full_b = (double*)calloc(sizeof(double), num_states);
        full_y = (double*)calloc(sizeof(double), num_states);
        for(i = 0, j = 0; i < num_states; i++)
        {
            if(i == zvi[j])
            {
                j++;
            }
            else
            {
                full_b[i] = b[i-j];
                full_y[i] = y[i-j];
            }
        }
    }
    else
    {
        full_b = b;
        full_y = y;
    }
    if(diffusion)
        nrn_tree_solve(_rxd_a, _rxd_b, _rxd_c, _rxd_d, full_b, _rxd_p, _rxd_euler_nrow, dt);

    do_ics_reactions(full_y, full_b, y, b);
   
    if(_rxd_num_zvi > 0)
    {
        for(i = 0, j = 0; i < num_states; i++)
        {
            if(i == zvi[j])
                j++;
            else
                b[i-j] = full_b[i]; 
        }
        free(full_b);
        free(full_y);
    }
    
}

static void ode_abs_tol(double* p1)
{
    int i, j;
    double* y = p1 + _cvode_offset;
    if(species_indices != NULL)
    {
        SpeciesIndexList* list;
        for(list = species_indices; list->next != NULL; list = list->next)
        {
            for(i=0, j=0; i<list->length; i++)
            {
                for(; j < _rxd_num_zvi  && _rxd_zero_volume_indices[j] <= list->indices[i]; j++);
                y[list->indices[i] - j] *= list->atolscale; 
            }

        }
    }
}

static void free_currents()
{
    int i, j;
    if(!_membrane_flux)
        return;
    for(i = 0; i < _memb_count; i++)
    {
        for(j = 0; j < _memb_species_count[i]; j++)
        {
            free(_memb_cur_mapped[i][j]);
        }
        free(_memb_cur_mapped[i]);
        free(_memb_cur_ptrs[i]);
    }
    free(_memb_cur_ptrs);
    free(_memb_cur_mapped);
    free(_memb_species_count);
    free(_cur_node_indices);
    free(_rxd_induced_currents);
    free(_rxd_flux_scale);
    free(_membrane_lookup);
    free(_memb_cur_mapped_ecs);
    free(_rxd_induced_currents_grid);
    free(_rxd_induced_currents_scale);
    _membrane_flux = FALSE; 
}

extern "C" void setup_currents(int num_currents, int num_fluxes,
        int* num_species, int* node_idxs, double* scales,
        PyHocObject** ptrs, int* mapped, int* mapped_ecs)
{
    int i, j, k, id, side, count;
    int* induced_currents_ecs_idx;
    int* induced_currents_grid_id;
    int* ecs_indices;
    double* current_scales;
    PyHocObject** ecs_ptrs;

    Current_Triple* c;
    Grid_node* g;
    ECS_Grid_node* grid;

    
    free_currents();
    
    _memb_count = num_currents;
    _memb_curr_total = num_fluxes;
    _memb_species_count = (int*)malloc(sizeof(int)*num_currents);
    memcpy(_memb_species_count,num_species,sizeof(int)*num_currents);

    _rxd_flux_scale = (double*)calloc(sizeof(double),num_fluxes);
    //memcpy(_rxd_flux_scale,scales,sizeof(double)*num_currents);

    /*TODO: if memory is an issue - replace with a map*/
    _membrane_lookup = (int*)malloc(sizeof(int)*num_states);
     memset(_membrane_lookup, SPECIES_ABSENT, sizeof(int)*num_states);

    _memb_cur_ptrs = (double***)malloc(sizeof(double**)*num_currents);
    _memb_cur_mapped_ecs = (int***)malloc(sizeof(int*)*num_currents);        
    _memb_cur_mapped = (int***)malloc(sizeof(int**)*num_currents);
    induced_currents_ecs_idx = (int*)malloc(sizeof(int)*_memb_curr_total);
    induced_currents_grid_id = (int*)malloc(sizeof(int)*_memb_curr_total);
    // initialize memory here to allow currents from an intracellular species
    // with no corresponding nrn_region='o' or Extracellular species 
    memset(induced_currents_ecs_idx, SPECIES_ABSENT,
           sizeof(int)*_memb_curr_total);

    for(i = 0, k = 0; i < num_currents; i++)
    {
        _memb_cur_ptrs[i] = (double**)malloc(sizeof(double*)*num_species[i]);
        //memcpy(_memb_cur_ptrs[i], &ptrs[k], sizeof(PyHocObject*)*num_species[i]);
        _memb_cur_mapped_ecs[i] = (int**)malloc(sizeof(int*)*num_species[i]);
        _memb_cur_mapped[i] =  (int**)malloc(sizeof(int*)*num_species[i]);


        for(j = 0; j < num_species[i]; j++, k++)
        {
            _memb_cur_ptrs[i][j] = (double*)ptrs[k]->u.px_;
            _memb_cur_mapped[i][j] = (int*)malloc(2*sizeof(int));
            _memb_cur_mapped_ecs[i][j] = (int*)malloc(2*sizeof(int));


            for(side = 0; side < 2; side++)
            {
                _memb_cur_mapped[i][j][side] = mapped[2*k + side];
                _memb_cur_mapped_ecs[i][j][side] = mapped_ecs[2*k + side];
            }
            for(side = 0; side < 2; side++)
            {
                if(_memb_cur_mapped[i][j][side] != SPECIES_ABSENT)
                {
                    _membrane_lookup[_curr_indices[_memb_cur_mapped[i][j][side]]] = k;
                    _rxd_flux_scale[k] = scales[i];
                    if(_memb_cur_mapped[i][j][(side+1)%2] == SPECIES_ABSENT)
                    {
                        induced_currents_grid_id[k] = _memb_cur_mapped_ecs[i][j][0];
                        induced_currents_ecs_idx[k] = _memb_cur_mapped_ecs[i][j][1];
                    }
                }
            }
        }
    }
    _rxd_induced_currents_grid = (ECS_Grid_node**)calloc(_memb_curr_total,sizeof(ECS_Grid_node*));
    _rxd_induced_currents_scale = (double*)calloc(_memb_curr_total,sizeof(double));
    for (id = 0, g = Parallel_grids[0]; g != NULL; g = g -> next, id++)
    {
        grid = dynamic_cast<ECS_Grid_node*>(g);
        if(grid == NULL) continue; // ignore ICS grids
  
        for(count = 0, k = 0; k < _memb_curr_total; k++)
        {
            if(induced_currents_grid_id[k] == id)
            {
                _rxd_induced_currents_grid[k] = grid;
                count++;
            }
        }
        if(count > 0)
        {
            ecs_indices = (int*)malloc(count * sizeof(int));
            ecs_ptrs = (PyHocObject**)malloc(count * sizeof(PyHocObject*));
            for(i = 0, k = 0; k < _memb_curr_total; k++)
            {
                if(induced_currents_grid_id[k] == id)
                {
                    ecs_indices[i] = induced_currents_ecs_idx[k];
                    ecs_ptrs[i++] = ptrs[k];
                }
            }
            current_scales = grid->set_rxd_currents(count, ecs_indices, ecs_ptrs);
            free(ecs_ptrs);

            for(i = 0, k = 0; k < _memb_curr_total; k++)
            {
                if(induced_currents_grid_id[k] == id)
                    _rxd_induced_currents_scale[k] = current_scales[i];
            }
        }
    }
    /*index into arrays of nodes/states*/
    _cur_node_indices = (int*)malloc(sizeof(int)*num_currents);
    memcpy(_cur_node_indices, node_idxs, sizeof(int)*num_currents);
    _membrane_flux = TRUE;
    _rxd_induced_currents = (double*)malloc(sizeof(double)*_memb_curr_total);
    free(induced_currents_ecs_idx);
    free(induced_currents_grid_id);

}

static void _currents(double* rhs)
{
    int i, j, k, idx, side;
    Grid_node* g;
    ECS_Grid_node* grid;
    if(!_membrane_flux)
        return;
    get_all_reaction_rates(states, NULL, NULL);
    for (g = Parallel_grids[0]; g != NULL; g = g -> next)
    {
        grid = dynamic_cast<ECS_Grid_node*>(g);
        if(grid)
            grid->induced_idx = 0;
    }
    for(i = 0, k = 0; i < _memb_count; i++)
    {
        idx = _cur_node_indices[i];

        for(j = 0; j < _memb_species_count[i]; j++, k++)
        {
            rhs[idx] -= _rxd_induced_currents[k];
            *(_memb_cur_ptrs[i][j]) += _rxd_induced_currents[k];
 
            for(side = 0; side < 2; side++)
            {
                if(_memb_cur_mapped[i][j][side] == SPECIES_ABSENT)
                {
                        /*Extracellular region is within the ECS grid*/
                        grid = _rxd_induced_currents_grid[k];
                        if(grid != NULL &&_memb_cur_mapped[i][j][(side+1)%2] != SPECIES_ABSENT)
                            grid->local_induced_currents[grid->induced_idx++] = _rxd_induced_currents[k];
                }
            }
        }
    }
    
}

extern "C" int rxd_nonvint_block(int method, int size, double* p1, double* p2, int) {
        if(initialized)
        {
            if(structure_change_cnt != prev_structure_change_cnt)
            {
                /*TODO: Exclude irrelevant (non-rxd) structural changes*/
                /*Needed for node.include_flux*/
                _setup_matrices();
            }
            if (prev_nrnunit_use_legacy != _nrnunit_use_legacy_)
            {
                _setup_units();
                prev_nrnunit_use_legacy = _nrnunit_use_legacy_;
            }
        }
        switch (method) {
        case 0:
            _setup();
            break;
        case 1:
            _initialize();
            //TODO: is there a better place to initialize multicompartment reactions
            Grid_node* grid;
            ECS_Grid_node* g;
            for (grid = Parallel_grids[0]; grid != NULL; grid = grid -> next)
            {
                g = dynamic_cast<ECS_Grid_node*>(grid);
                if(g) g->initialize_multicompartment_reaction();
            }
            break;
        case 2:
            /* compute outward current to be subtracted from rhs */
            _currents(p1);
            break;
        case 3:
            /* conductance to be added to d */
            break;
        case 4:
            /* fixed step solve */
            _fadvance();
            _fadvance_fixed_step_3D();
            break;
        case 5:
            /* ode_count */
            _cvode_offset = size;
            _ecs_count = ode_count(size + num_states - _rxd_num_zvi);
            return _ecs_count + num_states - _rxd_num_zvi;
        case 6:
            /* ode_reinit(y) */
            _ode_reinit(p1); 
            _ecs_ode_reinit(p1);
            break;
        case 7:
            /* ode_fun(t, y, ydot); from t and y determine ydot */
            _rhs_variable_step(p1, p2);
            _rhs_variable_step_ecs(p1, p2);
            break;
        case 8:
            ode_solve(*dt_ptr, p1, p2); /*solve mx=b replace b with x */
            /* TODO: we can probably reuse the dgadi code here... for now, we do nothing, which implicitly approximates the Jacobian as the identity matrix */
            //y= p1 = states and b = p2 = RHS for x direction
            ics_ode_solve(*dt_ptr, p1, p2);
            break;
        case 9:
            //ode_jacobian(*dt_ptr, p1, p2); /* optionally prepare jacobian for fast ode_solve */
            break;
        case 10:
            ode_abs_tol(p1);
            ecs_atolscale(p1);
            /* ode_abs_tol(y_abs_tolerance); fill with cvode.atol() * scalefactor */
            break;
        default:
            printf("Unknown rxd_nonvint_block call: %d\n", method);
            break;
    }
    return 0;
}


 
/*****************************************************************************
*
* Begin intracellular code
*
*****************************************************************************/


extern "C" void register_rate(int nspecies, int nparam, int nregions, int nseg,
                              int* sidx, int necs, int necsparam, int* ecs_ids,
                              int* ecsidx, int nmult, double* mult,
                              PyHocObject** vptrs, ReactionRate f)
{
    int i,j,k,idx, ecs_id, ecs_index, ecs_offset;
    unsigned char counted;
    Grid_node* g;
    ECS_Grid_node* grid;
    ICSReactions* react = (ICSReactions*)malloc(sizeof(ICSReactions));
    react->reaction = f;
    react->num_species = nspecies;
    react->num_regions = nregions;
    react->num_params = nparam;
    react->num_segments = nseg;
    react->num_ecs_species = necs;
    react->num_ecs_params = necsparam; 
    react->num_mult = nmult;
    react->icsN = 0;
    react->ecsN = 0;
    if (vptrs != NULL)
    {
        react->vptrs = (double**)malloc(nseg*sizeof(double*));
        for(i = 0; i < nseg; i++)
           react->vptrs[i] = vptrs[i]->u.px_; 
    }
    else
    {
        react->vptrs = NULL;
    }
    react->state_idx = (int***)malloc(nseg*sizeof(double**));
    for(i = 0, idx = 0; i < nseg; i++)
    {
        react->state_idx[i] = (int**)malloc((nspecies + nparam)*sizeof(int*));
        for(j = 0; j < nspecies + nparam; j++)
        {
            react->state_idx[i][j] = (int*)malloc(nregions*sizeof(int));
            for(k = 0; k < nregions; k++, idx++)
            {
                if(sidx[idx] < 0)
                {
                    react->state_idx[i][j][k] = SPECIES_ABSENT;
                }
                else
                {
                    react->state_idx[i][j][k] = sidx[idx];
                    if(i==0 && j < nspecies)
                        react->icsN++;
                }
            }
        } 
    }
    if( nmult > 0 )
    {
        react->mc_multiplier = (double**)malloc(nmult*sizeof(double*));
        for(i = 0; i < nmult; i++)
        {
            react->mc_multiplier[i] = (double*)malloc(nseg*sizeof(double));
            memcpy(react->mc_multiplier[i], (mult+i*nseg), nseg*sizeof(double));
        }
    }

    if(react->num_ecs_species + react->num_ecs_params > 0)
    {
        react->ecs_grid = (ECS_Grid_node**)malloc(react->num_ecs_species*sizeof(Grid_node*));
        react->ecs_state = (double***)malloc(nseg*sizeof(double**));
        react->ecs_index = (int**)malloc(nseg*sizeof(int*));
        react->ecs_offset_index = (int*)malloc(react->num_ecs_species*sizeof(int));
        for(i = 0; i < nseg; i++)
        {
            react->ecs_state[i] = (double**)malloc((necs+necsparam)*sizeof(double*));
            react->ecs_index[i] = (int*)malloc((necs+necsparam)*sizeof(int));

        }
        for(j = 0; j < necs + necsparam; j++)
        {
            ecs_offset = num_states - _rxd_num_zvi;

            for(ecs_id = 0, g = Parallel_grids[0]; g != NULL; g = g -> next, ecs_id++)
            {
            
                if (ecs_id == ecs_ids[j])
                {
                    grid = dynamic_cast<ECS_Grid_node*>(g);
                    assert(grid != NULL);
                    if (j < necs)
                    {
                        react->ecs_grid[j] = grid;
                        react->ecs_offset_index[j] = grid->add_multicompartment_reaction(nseg, &ecsidx[j], necs + necsparam);
                    }
                    
                    for(i = 0, counted=FALSE; i < nseg; i++)
                    {
                        //nseg x nregion x nspecies
                        ecs_index = ecsidx[i*(necs + necsparam) + j];

                        if(ecs_index >= 0)
                        {
                            react->ecs_state[i][j] = &(grid->states[ecs_index]);
                            react->ecs_index[i][j] = ecs_offset + ecs_index;
                            if(j < necs && !counted)
                            {
                                react->ecsN++;
                                counted=TRUE;
                            }
                        }
                        else
                        {
                            react->ecs_state[i][j] = NULL;
                        }
                    }
                    ecs_offset += grid->size_x*grid->size_y*grid->size_z;
                }
            }
        }
    }
    else
    {
        react->ecs_state = NULL;
    }
    if(_reactions == NULL)
    {
        _reactions = react;
        react -> next = NULL;

    }
    else
    {
        react -> next = _reactions;
        _reactions = react;
    }

    for (g = Parallel_grids[0]; g != NULL; g = g -> next)
    {
        grid = dynamic_cast<ECS_Grid_node*>(g);
        if(grid)
            grid->initialize_multicompartment_reaction();
    }
}

extern "C" void clear_rates()
{
    ICSReactions *react, *prev;
    int i, j;
    for(react = _reactions; react != NULL;)
    {
        if(react->vptrs != NULL)
            free(react->vptrs);
        for(i = 0; i < react->num_segments; i++)
        {
            for(j = 0; j < react->num_species; j++)
            {
                free(react->state_idx[i][j]);
            }
            free(react->state_idx[i]);
           
            if(react->num_ecs_species + react->num_ecs_params > 0)
            {
                free(react->ecs_state[i]);
            }
        }
        if(react->num_mult > 0)
        {
            for(i = 0; i < react->num_mult; i++)
                free(react->mc_multiplier[i]);
            free(react->mc_multiplier);
        }

        free(react->state_idx);
        SAFE_FREE(react->ecs_state);
        prev = react;
        react = react->next;
        SAFE_FREE(prev);
    }
    _reactions = NULL;
    /*clear extracellular reactions*/
    clear_rates_ecs();
}


extern "C" void species_atolscale(int id, double scale, int len, int* idx)
{
    SpeciesIndexList* list;
    SpeciesIndexList* prev;
    if(species_indices != NULL)
    {
        for(list = species_indices, prev = NULL;
            list != NULL; list = list->next)
        {
            if(list->id == id)
            {
                list->atolscale = scale;
                return;
            }
            prev = list;
        }
        prev->next = (SpeciesIndexList*)malloc(sizeof(SpeciesIndexList));
        list = prev->next;
    }
    else
    {
        species_indices = (SpeciesIndexList*)malloc(sizeof(SpeciesIndexList));
        list = species_indices;
    }
    list->id = id;
    list->indices = (int*)malloc(sizeof(int)*len);
    memcpy(list->indices,idx,sizeof(int)*len);
    list->length = len;
    list->atolscale = scale;
    list->next = NULL;
}

extern "C" void remove_species_atolscale(int id)
{
    SpeciesIndexList* list;
    SpeciesIndexList* prev;
    for(list = species_indices, prev = NULL; 
        list != NULL; prev = list, list = list->next)
    {
        if(list->id == id)
        {
            if(prev == NULL)
                species_indices = list->next;
            else
                prev->next = list->next;
            free(list->indices);
            free(list);
            break;
        }
    }
}

extern "C" void setup_solver(double* my_states, int my_num_states, long* zvi, int num_zvi) {
    free_currents();
    states = my_states;
    num_states = my_num_states;
    free_zvi_child();
    _rxd_num_zvi = num_zvi;
    if (_rxd_zero_volume_indices != NULL)
        free(_rxd_zero_volume_indices);
    if (num_zvi == 0)
        _rxd_zero_volume_indices = NULL;
    else
        _rxd_zero_volume_indices = (long*)allocopy(zvi, num_zvi * sizeof(long));
    dt_ptr = &nrn_threads->_dt;
    t_ptr = &nrn_threads->_t;
    set_num_threads(NUM_THREADS);
    initialized = TRUE;
    prev_structure_change_cnt = structure_change_cnt;
}

void start_threads(const int n)
{
    int i;
    if(Threads == NULL)
    {
        AllTasks = (TaskQueue*)calloc(1,sizeof(TaskQueue));
        Threads = (pthread_t*)malloc(sizeof(pthread_t)*(n-1));
        AllTasks->task_mutex = (pthread_mutex_t*)malloc(sizeof(pthread_mutex_t));
        AllTasks->waiting_mutex = (pthread_mutex_t*)malloc(sizeof(pthread_mutex_t));
        AllTasks->task_cond = (pthread_cond_t*)malloc(sizeof(pthread_cond_t));
        AllTasks->waiting_cond = (pthread_cond_t*)malloc(sizeof(pthread_cond_t));
        pthread_mutex_init(AllTasks->task_mutex, NULL);
        pthread_cond_init(AllTasks->task_cond, NULL);
        pthread_mutex_init(AllTasks->waiting_mutex, NULL);
        pthread_cond_init(AllTasks->waiting_cond, NULL);
        AllTasks->length = 0;
        for(i = 0; i < n-1; i++)
            pthread_create(&Threads[i], NULL, TaskQueue_exe_tasks, AllTasks);
    }

}

void TaskQueue_add_task(TaskQueue* q, void* (*task)(void*), void* args, void* result)
{
    TaskList *t;
    t = (TaskList*)malloc(sizeof(TaskList));
    t->task = task;
    t->args = args;
    t->result = result;
    t->next = NULL;
    
    //Add task to the queue
    pthread_mutex_lock(q->task_mutex);
    if(q->first == NULL)  //empty queue
    {
        q->first=t;
        q->last=t;
    }
    else                    //non-empty
    {
        q->last->next = t;
        q->last = t;
    }

    pthread_mutex_lock(q->waiting_mutex);
    q->length++;
    pthread_mutex_unlock(q->waiting_mutex);
    pthread_mutex_unlock(q->task_mutex);
    
    //signal waiting threads
    pthread_cond_signal(q->task_cond);
 
}

void* TaskQueue_exe_tasks(void* dat)
{
    TaskList* job;
    TaskQueue* q = (TaskQueue*)dat;
    /*int id;
    for(id=0;id<NUM_THREADS;id++)
    {
        if (pthread_equal(Threads[id],pthread_self()))
        {
            break;
        }
    }
    fprintf(stderr,"%i] ready\n",id);
    */
    while(1)    //loop until thread is killed
    {
        pthread_mutex_lock(q->task_mutex);
        while(q->first == NULL) //no tasks
        {
            //Wait for new tasks
            pthread_cond_wait(q->task_cond, q->task_mutex);
        }
        //fprintf(stderr,"%i] running\n",id); 
        job = q->first;
        q->first = job->next;
        pthread_mutex_unlock(q->task_mutex); 
    
        //execute
        job->result = job->task(job->args);
        free(job);

        //fprintf(stderr,"%i] updating\n",id); 
        pthread_mutex_lock(q->waiting_mutex);
        if(--(q->length) == 0)  //all finished
        {
            pthread_cond_broadcast(q->waiting_cond);
        }
        pthread_mutex_unlock(q->waiting_mutex);
        //fprintf(stderr,"%i] done\n",id);
    } 
    
    return NULL;
}


void set_num_threads(const int n)
{
    int k, old_num = NUM_THREADS;
    if(Threads == NULL)
    {
        start_threads(n);
    }
    else
    {
        if(n<old_num)
        {
            //Kill some threads
            for(k=old_num-1; k>=n; k--)
            {
                TaskQueue_sync(AllTasks);
                pthread_cancel(Threads[k]);
            } 
            Threads = (pthread_t*)realloc(Threads,sizeof(pthread_t) * n);
            assert(Threads);
        }
        else if(n>old_num)
        {
            //Create some threads
            Threads = (pthread_t*)realloc(Threads,sizeof(pthread_t) * n);
            assert(Threads);
            
            for (k = old_num-1; k < n; k++) 
            {
                pthread_create(&Threads[k], NULL, TaskQueue_exe_tasks, AllTasks);
            }
        }
    }
    set_num_threads_3D(n);
    NUM_THREADS = n;

}

void TaskQueue_sync(TaskQueue *q)
{
    //Wait till the queue is empty
    pthread_mutex_lock(q->waiting_mutex);
    while(q->length > 0)
        pthread_cond_wait(q->waiting_cond,q->waiting_mutex);
    pthread_mutex_unlock(q->waiting_mutex);
}

int get_num_threads(void) {
    return NUM_THREADS;
}


void _fadvance(void) {
    double dt = *dt_ptr;
    long i;
    /*variables for diffusion*/
    double *rhs; 
    long* zvi = _rxd_zero_volume_indices;
    
    rhs = (double*)calloc(num_states,sizeof(double));
    /*diffusion*/
    if(diffusion)
        mul(_rxd_euler_nnonzero, _rxd_euler_nonzero_i, _rxd_euler_nonzero_j, _rxd_euler_nonzero_values, states, rhs);

    add_currents(rhs);

    /* multiply rhs vector by dt */
    for (i = 0; i < num_states; i++) {
        rhs[i] *= dt;
    }

    if(diffusion)
        nrn_tree_solve(_rxd_a, _rxd_b, _rxd_c, _rxd_d, rhs, _rxd_p, _rxd_euler_nrow, dt);
   
    /* increment states by rhs which is now really deltas */
    for (i = 0; i < num_states; i++) {
        states[i] += rhs[i];
    }

    /* clear zero volume indices (conservation nodes) */
    for (i = 0; i < _rxd_num_zvi; i++) {
        states[zvi[i]] = 0;
    }
    
    free(rhs);

    /*reactions*/
    do_ics_reactions(states, NULL, NULL, NULL);

    /*node fluxes*/
    apply_node_flux1D(dt, states);

    transfer_to_legacy();
}


void _ode_reinit(double* y)
{
    long i, j;
    long* zvi = _rxd_zero_volume_indices;
    y += _cvode_offset;
    /*Copy states to CVode*/ 
    if(_rxd_num_zvi > 0)
    {
        for(i=0, j=0; i < num_states; i++)
        {
            if(zvi[j] == i)
                j++;
            else
            {
                y[i-j] = states[i];
            }
        }
    }
    else
    {
        memcpy(y,states,sizeof(double)*num_states);
    }
}


void _rhs_variable_step(const double* p1, double* p2) 
{
    Grid_node *grid;
    long i, j, p, c;
    unsigned int k;
    const unsigned char calculate_rhs = p2 == NULL ? 0 : 1;
    const double* my_states = p1 + _cvode_offset;
    double* ydot = p2 + _cvode_offset;
    /*variables for diffusion*/
    double *rhs;
    long* zvi = _rxd_zero_volume_indices;

    double const * const orig_states3d = p1 + states_cvode_offset;
    double* const orig_ydot3d = p2 + states_cvode_offset;

    /*Copy states from CVode*/ 
    if(_rxd_num_zvi > 0)
    {
        for(i=0, j=0; i < num_states; i++)
        {
            if(zvi[j] == i)
                j++;
            else{
                states[i] =  my_states[i-j];
           }
        }
    }
    else
    {
        memcpy(states,my_states,sizeof(double)*num_states);
    }

    if(diffusion)
    {
        for(i = 0; i < _rxd_num_zvi; i++)
        {
            j = zvi[i];
            p = _rxd_p[j];
            states[j] = (p>0?-(_rxd_b[j]/_rxd_d[j])*states[p]:0);
            for(k = 0; k < _rxd_zvi_child_count[i]; k++)
            {
                c = _rxd_zvi_child[i][k];
                states[j] -= (_rxd_a[c]/_rxd_d[j])*states[c]; 
            }
        }
    }

    transfer_to_legacy();

    if(!calculate_rhs)
    {
        for(i = 0; i < _rxd_num_zvi; i++)
            states[zvi[i]] = 0;
        return;
    }

    /*diffusion*/
    rhs = (double*)calloc(num_states,sizeof(double));
    
    if(diffusion)
        mul(_rxd_euler_nnonzero, _rxd_euler_nonzero_i, _rxd_euler_nonzero_j, _rxd_euler_nonzero_values, states, rhs);

    /*reactions*/
    MEM_ZERO(&ydot[num_states - _rxd_num_zvi], sizeof(double)*_ecs_count);
    get_all_reaction_rates(states, rhs, ydot);


    const double* states3d = orig_states3d;
    double* ydot3d = orig_ydot3d;
    int grid_size;
    for (grid = Parallel_grids[0]; grid != NULL; grid = grid -> next) {
        grid_size = grid->size_x * grid->size_y * grid->size_z;
        if(grid->hybrid)
        {
            grid->variable_step_hybrid_connections(states3d, ydot3d, states, rhs);
        }
        ydot3d += grid_size;
        states3d += grid_size;        
    }

    /*Add currents to the result*/
    add_currents(rhs);

    /*Add node fluxes to the result*/
    apply_node_flux1D(1.0, rhs);

    /* increment states by rhs which is now really deltas */
    if(_rxd_num_zvi > 0)
    {
        for (i = 0, j = 0; i < num_states; i++)
        {
            if(zvi[j] == i)
            {
                states[i] = 0;
                j++;
            }
            else
            {
                ydot[i-j] = rhs[i];
            }
        }
    }
    else
    {
        memcpy(ydot,rhs,sizeof(double)*num_states);
    }
    free(rhs);
}

void get_reaction_rates(ICSReactions* react, double* states, double* rates, double* ydot)
{
    int segment;
    int i, j, k, idx;
    double** states_for_reaction = (double**)malloc(react->num_species*sizeof(double*));
    double** params_for_reaction = (double**)malloc(react->num_params*sizeof(double*));
    double** result_array = (double**)malloc(react->num_species*sizeof(double*));
    double* mc_mult = NULL;
    double** flux = NULL;
    if(react->num_mult > 0)
        mc_mult = (double*)malloc(react->num_mult*sizeof(double));

    double* ecs_states_for_reaction = NULL;
    double* ecs_params_for_reaction = NULL;
    double* ecs_result = NULL;
    int * ecsindex = NULL;
    double v = 0;
    if(react->num_ecs_species > 0)
    {
        ecs_states_for_reaction = (double*)malloc(react->num_ecs_species*sizeof(double));
        ecs_result = (double*)malloc(react->num_ecs_species*sizeof(double));
    }
    if(react->num_ecs_params > 0)
        ecs_params_for_reaction = (double*)calloc(react->num_ecs_params,sizeof(double));

    if(_membrane_flux)
    {
        flux = (double**)malloc(react->icsN*sizeof(double*));
        for(i = 0; i < react->icsN; i++)
            flux[i] = (double*)calloc(react->num_regions, sizeof(double));
    }
    
    for(i = 0; i < react->num_species; i++)
    {
        states_for_reaction[i] = (double*)calloc(react->num_regions,sizeof(double));
        result_array[i] = (double*)malloc(react->num_regions*sizeof(double));
        
    }
    for(i = 0; i < react->num_params; i++)
        params_for_reaction[i] = (double*)calloc(react->num_regions,sizeof(double));
    ecsindex = (int*)malloc(react->num_ecs_species*sizeof(int));
    for(i = 0; i < react->num_ecs_species; i++)
        ecsindex[i] = react->ecs_grid[i]->react_offsets[react->ecs_offset_index[i]];

    for(segment = 0; segment < react->num_segments; segment++)
    {
        for(i = 0; i < react->num_species; i++)
        {
            for(j = 0; j < react->num_regions; j++)
            {
                if(react->state_idx[segment][i][j] != SPECIES_ABSENT)
                {
                    states_for_reaction[i][j] = states[react->state_idx[segment][i][j]];
                }
                else
                {
                    states_for_reaction[i][j] = NAN;
                }
            }
            MEM_ZERO(result_array[i],react->num_regions*sizeof(double));
        }
        for(k = 0; i < react->num_species + react->num_params; i++, k++)
        {
            for(j = 0; j < react->num_regions; j++)
            {
                if(react->state_idx[segment][i][j] != SPECIES_ABSENT)
                {
                    params_for_reaction[k][j] = states[react->state_idx[segment][i][j]];
                }
                else
                {
                    params_for_reaction[k][j] = NAN;
                }
            }
        }

        for(i = 0; i < react->num_ecs_species; i++)
        {
            if(react->ecs_state[segment][i] != NULL)
            {
                ecs_states_for_reaction[i] = *(react->ecs_state[segment][i]);
            }
            else
            {
                ecs_states_for_reaction[i] = NAN;

            }
        }
        for(k = 0; i < react->num_ecs_species + react->num_ecs_params; i++, k++)
        {
            if(react->ecs_state[segment][i] != NULL)
            {
                ecs_params_for_reaction[k] = *(react->ecs_state[segment][i]);
            }
            else
            {
                ecs_params_for_reaction[k] = NAN;
            }
        }
        MEM_ZERO(ecs_result,react->num_ecs_species*sizeof(double));

        for(i = 0; i < react->num_mult; i++)
        {
            mc_mult[i] = react->mc_multiplier[i][segment];
        }

        if(react->vptrs != NULL)
        {
            v = *(react->vptrs[segment]);
        }

        react->reaction(states_for_reaction, params_for_reaction, result_array, mc_mult, ecs_states_for_reaction, ecs_params_for_reaction, ecs_result, flux, v);
        
        for(i = 0; i < react->num_species; i++)
        {
            for(j = 0; j < react->num_regions; j++)
            {
                idx = react->state_idx[segment][i][j];
                if(idx != SPECIES_ABSENT)
                {
                    if(_membrane_flux && _membrane_lookup[idx] != SPECIES_ABSENT)
                    {
                        _rxd_induced_currents[_membrane_lookup[idx]] -= _rxd_flux_scale[_membrane_lookup[idx]] * flux[i][j];
                    }
                    if(rates != NULL)
                    {
                        rates[idx] += result_array[i][j];
                    }
                }
            }
        }
        if(ydot != NULL)
        {
            for(i = 0;  i < react->num_ecs_species; i++)
            {
                if(react->ecs_state[segment][i] != NULL)
                    react->ecs_grid[i]->all_reaction_states[ecsindex[i]++] = ecs_result[i];
                    //ydot[react->ecs_index[segment][i]] += ecs_result[i];
            }
        }
    }
    /* free allocated memory */
    if(react->num_mult > 0)
        free(mc_mult);
    if(_membrane_flux)
    {
        for(i = 0; i < react->icsN; i++)
            free(flux[i]);
        free(flux);
    }
    if(react->num_ecs_species>0)
    {
        free(ecs_states_for_reaction);
        free(ecs_result);
    }
    for(i = 0; i < react->num_species; i++)
    {
        free(states_for_reaction[i]);
        free(result_array[i]);
    }
    free(states_for_reaction);
    free(result_array);
    for(i = 0; i < react->num_params; i++)
    {
        free(params_for_reaction[i]);
    }
    free(params_for_reaction);
    if(react->num_ecs_params>0)
    {
        free(ecs_params_for_reaction);
    }


}

void solve_reaction(ICSReactions* react, double* states, double *bval, double* cvode_states, double* cvode_b)
{
    int segment;
    int i, j, k, idx, jac_i, jac_j, jac_idx;
    int N = react->icsN + react->ecsN;  /*size of Jacobian (number species*regions for a segments)*/
    double pd;
    double dt = *dt_ptr;
    double dx = FLT_EPSILON;
    MAT* jacobian = m_get(N,N);
    VEC* b = v_get(N);
    VEC* x = v_get(N);
    PERM* pivot = px_get(N);
    
    double** states_for_reaction = (double**)malloc(react->num_species*sizeof(double*));
    double** states_for_reaction_dx = (double**)malloc(react->num_species*sizeof(double*));
    double** params_for_reaction = (double**)malloc(react->num_params*sizeof(double*));
    double** result_array = (double**)malloc(react->num_species*sizeof(double*));
    double** result_array_dx = (double**)malloc(react->num_species*sizeof(double*));
    double* mc_mult = NULL;
    if(react->num_mult > 0)
        mc_mult = (double*)malloc(react->num_mult*sizeof(double));

    double* ecs_states_for_reaction = NULL;
    double* ecs_states_for_reaction_dx = NULL;
    double* ecs_params_for_reaction = NULL;
    double* ecs_result = NULL;
    double* ecs_result_dx = NULL;
    double v = 0;
    int* ecsindex = NULL;
    
    if(react->num_ecs_species > 0)
    {
        ecsindex = (int*)malloc(react->num_ecs_species*sizeof(int));
        for(i = 0; i < react->num_ecs_species; i++)
            ecsindex[i] = react->ecs_grid[i]->react_offsets[react->ecs_offset_index[i]];
        ecs_states_for_reaction = (double*)malloc(react->num_ecs_species*sizeof(double));
        ecs_states_for_reaction_dx = (double*)malloc(react->num_ecs_species*sizeof(double));
        ecs_result = (double*)malloc(react->num_ecs_species*sizeof(double));
        ecs_result_dx = (double*)malloc(react->num_ecs_species*sizeof(double));
    }

    if(react->num_ecs_params > 0)
        ecs_params_for_reaction = (double*)malloc(react->num_ecs_params*sizeof(double));


    for(i = 0; i < react->num_species; i++)
    {
        states_for_reaction[i] = (double*)malloc(react->num_regions*sizeof(double));
        states_for_reaction_dx[i] = (double*)malloc(react->num_regions*sizeof(double));
        result_array[i] = (double*)malloc(react->num_regions*sizeof(double));
        result_array_dx[i] = (double*)malloc(react->num_regions*sizeof(double));
    }
    for(i = 0; i < react->num_params; i++)
        params_for_reaction[i] = (double*)malloc(react->num_regions*sizeof(double));

    for(segment = 0; segment < react->num_segments; segment++)
    {
        if(react->vptrs != NULL)
            v = *(react->vptrs[segment]);

        for(i = 0; i < react->num_species; i++)
        {
            for(j = 0; j < react->num_regions; j++)
            {
                if(react->state_idx[segment][i][j] != SPECIES_ABSENT)
                {
                    states_for_reaction[i][j] = states[react->state_idx[segment][i][j]];
                    states_for_reaction_dx[i][j] = states_for_reaction[i][j];
                }
                else
                {
                    states_for_reaction[i][j] = SPECIES_ABSENT;
                    states_for_reaction_dx[i][j] = states_for_reaction[i][j];
                }

            }
            MEM_ZERO(result_array[i],react->num_regions*sizeof(double));
            MEM_ZERO(result_array_dx[i],react->num_regions*sizeof(double));
        }
        for(k = 0; i < react->num_species + react->num_params; i++, k++)
        {
            for(j = 0; j < react->num_regions; j++)
            {
                if(react->state_idx[segment][i][j] != SPECIES_ABSENT)
                {
                    params_for_reaction[k][j] = states[react->state_idx[segment][i][j]];
                }
                else
                {
                    params_for_reaction[k][j] = SPECIES_ABSENT;
                }

            }
        }


        for(i = 0; i < react->num_ecs_species; i++)
        {
            if(react->ecs_state[segment][i] != NULL)
            {
                if(cvode_states != NULL)
                {
                    ecs_states_for_reaction[i] = cvode_states[react->ecs_index[segment][i]];
                }
                else
                {
                    ecs_states_for_reaction[i] = *(react->ecs_state[segment][i]);
                }
                ecs_states_for_reaction_dx[i] = ecs_states_for_reaction[i];
            }
        }
        for(k = 0; i < react->num_ecs_species + react->num_ecs_params; i++, k++)
        {
            if(react->ecs_state[segment][i] != NULL)
            {
                    ecs_params_for_reaction[k] = *(react->ecs_state[segment][i]);
            }
        }

        if(react->num_ecs_species > 0)
        {
            MEM_ZERO(ecs_result,react->num_ecs_species*sizeof(double));
            MEM_ZERO(ecs_result_dx,react->num_ecs_species*sizeof(double));
        }
        
        for(i = 0; i < react->num_mult; i++)
        {
            mc_mult[i] = react->mc_multiplier[i][segment];
        }

        react->reaction(states_for_reaction, params_for_reaction, result_array,
                        mc_mult, ecs_states_for_reaction,
                        ecs_params_for_reaction, ecs_result, NULL, v);

        /*Calculate I - Jacobian for ICS reactions*/
        for(i = 0, idx = 0; i < react->num_species; i++)
        {
            for(j = 0; j < react->num_regions; j++)
            {
                if(react->state_idx[segment][i][j] != SPECIES_ABSENT)
                {
                    if(bval == NULL)
                        v_set_val(b, idx, dt*result_array[i][j]);
                    else
                        v_set_val(b, idx, bval[react->state_idx[segment][i][j]]);

    
                    // set up the changed states array
                    states_for_reaction_dx[i][j] += dx;
    
                    /* TODO: Handle approximating the Jacobian at a function upper
                    * limit, e.g. acos(1)
                    */
                    react->reaction(states_for_reaction_dx, params_for_reaction,
                                    result_array_dx, mc_mult,
                                    ecs_states_for_reaction,
                                    ecs_params_for_reaction, ecs_result_dx,
                                    NULL, v);
    
                    for (jac_i = 0, jac_idx = 0; jac_i < react->num_species; jac_i++)
                    {
                        for (jac_j = 0; jac_j < react->num_regions; jac_j++)
                        {
                            // pd is our Jacobian approximated
                            if(react->state_idx[segment][jac_i][jac_j] != SPECIES_ABSENT)
                            {
                                pd = (result_array_dx[jac_i][jac_j] - result_array[jac_i][jac_j])/dx;
                                m_set_val(jacobian, jac_idx, idx, (idx==jac_idx) - dt*pd);
                                jac_idx += 1;
                            }
                            result_array_dx[jac_i][jac_j] = 0;
                        }
                    }
                    for (jac_i = 0; jac_i < react->num_ecs_species; jac_i++)
                    {
                        // pd is our Jacobian approximated
                        if(react->ecs_state[segment][jac_i] != NULL)
                        {
                            pd = (ecs_result_dx[jac_i] - ecs_result[jac_i])/dx;
                            m_set_val(jacobian, jac_idx, idx, -dt*pd);
                            jac_idx += 1;
                        }
                        ecs_result_dx[jac_i] = 0;
                    }
                    // reset dx array
                    states_for_reaction_dx[i][j] -= dx;
                    idx++;
                }
            }
        }
    
        /*Calculate I - Jacobian for MultiCompartment ECS reactions*/
        for(i = 0; i < react->num_ecs_species; i++)
        {
            if(react->ecs_state[segment][i] != NULL)
            {
                if(bval == NULL)
                    v_set_val(b, idx, dt*ecs_result[i]);
                else
                    v_set_val(b, idx, cvode_b[react->ecs_index[segment][i]]);

                    
                // set up the changed states array
                ecs_states_for_reaction_dx[i] += dx;

                /* TODO: Handle approximating the Jacobian at a function upper
                 * limit, e.g. acos(1)
                 */
                react->reaction(states_for_reaction, params_for_reaction, result_array_dx,
                                mc_mult, ecs_states_for_reaction_dx,
                                ecs_params_for_reaction, ecs_result_dx, NULL, v);
    
                for (jac_i = 0, jac_idx = 0; jac_i < react->num_species; jac_i++)
                {
                    for (jac_j = 0; jac_j < react->num_regions; jac_j++)
                    {
                        // pd is our Jacobian approximated
                        if(react->state_idx[segment][jac_i][jac_j] != SPECIES_ABSENT)
                        {
                            pd = (result_array_dx[jac_i][jac_j] - result_array[jac_i][jac_j])/dx;
                            m_set_val(jacobian, jac_idx, idx, - dt*pd);
                            jac_idx += 1;
                        }
                    }
                }
                for (jac_i = 0; jac_i < react->num_ecs_species; jac_i++)
                {
                    // pd is our Jacobian approximated
                    if(react->ecs_state[segment][jac_i] != NULL)
                    {
                        pd = (ecs_result_dx[jac_i] - ecs_result[jac_i])/dx;
                        m_set_val(jacobian, jac_idx, idx, (idx==jac_idx) - dt*pd);
                        jac_idx += 1;
                    }
                    else
                    {
                        m_set_val(jacobian, idx, idx, 1.0);
                    }
                    // reset dx array
                    ecs_states_for_reaction_dx[i] -= dx;
                }
                idx++;
            }
        }
        // solve for x, destructively
        LUfactor(jacobian, pivot);
        LUsolve(jacobian, pivot, b, x);

        if(bval != NULL) //variable-step
        {
            for(i = 0, jac_idx=0; i < react->num_species; i++)
            {
                for(j = 0; j < react->num_regions; j++)
                {
                    idx = react->state_idx[segment][i][j];
                    if(idx != SPECIES_ABSENT)
                    {
                        bval[idx] = v_get_val(x, jac_idx++);
                    }
                }
            }
            for(i = 0; i < react->num_ecs_species; i++)
            {
                if(react->ecs_state[segment][i] != NULL)
                    react->ecs_grid[i]->all_reaction_states[ecsindex[i]++] = v_get_val(x,jac_idx++);
                    //cvode_b[react->ecs_index[segment][i]] = v_get_val(x, jac_idx++);
            }
        }
        else  //fixed-step
        {
            for(i = 0, jac_idx=0; i < react->num_species; i++)
            {
                for(j = 0; j < react->num_regions; j++)
                {
                    idx = react->state_idx[segment][i][j];
                    if(idx != SPECIES_ABSENT)
                        states[idx] += v_get_val(x, jac_idx++);
                }
            }
            for(i = 0; i < react->num_ecs_species; i++)
            {
                if(react->ecs_state[segment][i] != NULL)
                    react->ecs_grid[i]->all_reaction_states[ecsindex[i]++] = v_get_val(x,jac_idx++);
            }
        }
    }
    free(ecsindex);
    m_free(jacobian);
    v_free(b);
    v_free(x);
    px_free(pivot);
    for(i = 0; i < react->num_species; i++)
    {
        free(states_for_reaction[i]);
        free(states_for_reaction_dx[i]);
        free(result_array[i]);
        free(result_array_dx[i]);
    }
    for(i = 0; i < react->num_params; i++)
        free(params_for_reaction[i]);
    if(react->num_mult > 0)
        free(mc_mult);
    free(states_for_reaction_dx);
    free(states_for_reaction);
    free(params_for_reaction);
    free(result_array);
    free(result_array_dx);
    if(react->num_ecs_species > 0)
    {
        free(ecs_states_for_reaction);
        free(ecs_states_for_reaction_dx);
        free(ecs_result);
        free(ecs_result_dx);
    }
    if(react->num_ecs_params > 0)
        free(ecs_params_for_reaction);
}

void do_ics_reactions(double* states, double* b, double* cvode_states, double* cvode_b)
{
    ICSReactions* react;
    for(react = _reactions; react != NULL; react = react->next)
    {
        if(react->icsN + react->ecsN > 0)
            solve_reaction(react, states, b, cvode_states, cvode_b);
    }
}

void get_all_reaction_rates(double* states, double* rates, double* ydot)
{
    ICSReactions* react;
    if(_membrane_flux)
        MEM_ZERO(_rxd_induced_currents, sizeof(double)*_memb_curr_total);
    for(react = _reactions; react != NULL; react = react->next)
    {
        if(react->icsN + react->ecsN > 0)
            get_reaction_rates(react, states, rates, ydot);
    }
}

/*****************************************************************************
*
* End intracellular code
*
*****************************************************************************/