#include "model.h"

#include <math.h>
#include <stdlib.h>

const char VERSION[] = "0.5.0";
const char LIBCELLML_VERSION[] = "0.5.0";

/* The counts would only be for exposed variables in each type. State variables should always be
   exposed i.e., it would be an error to set them to not exposed. Could have the counts having
   an EXPOSED_ prefix but then the generated code would not be the same as previously generated
   code and thus not be backwardsly compatibile. */
const size_t STATE_COUNT = 4;
const size_t CONSTANT_COUNT = 2;
const size_t COMPUTED_CONSTANT_COUNT = 0;
const size_t ALGEBRAIC_COUNT = 1;

/* The not exposed counts are also listed although I would not have cause for them outside of this code. */
#define NOT_EXPOSED_STATE_COUNT 0
#define NOT_EXPOSED_CONSTANT_COUNT 3
#define NOT_EXPOSED_COMPUTED_CONSTANT_COUNT 3
#define NOT_EXPOSED_ALGEBRAIC_COUNT 9

/* The constant variables that are not exposed are then listed. */
const double neConstants[NOT_EXPOSED_CONSTANT_COUNT] = {0.0,0.3,36.0}; /* E_R, g_L and g_K */

/* I've added the exposed state for consistency in the variable declarations */

const VariableInfoWithType VOI_INFO = {"time", "millisecond", "environment", VARIABLE_OF_INTEGRATION, EXPOSED};

const VariableInfoWithType STATE_INFO[] = {
  {"V", "millivolt", "membrane", STATE, EXPOSED},
  {"h", "dimensionless", "sodium_channel_h_gate", STATE, EXPOSED},
  {"m", "dimensionless", "sodium_channel_m_gate", STATE, EXPOSED},
  {"n", "dimensionless", "potassium_channel_n_gate", STATE, EXPOSED}
};

const VariableInfoWithType CONSTANT_INFO[] = {
  {"Cm", "microF_per_cm2", "membrane", CONSTANT, EXPOSED,},
  {"E_R", "millivolt", "membrane", CONSTANT, NOT_EXPOSED,},
  {"g_L", "milliS_per_cm2", "leakage_current", CONSTANT, NOT_EXPOSED,},
  {"g_Na", "milliS_per_cm2", "sodium_channel", CONSTANT, EXPOSED,},
  {"g_K", "milliS_per_cm2", "potassium_channel", CONSTANT, NOT_EXPOSED}
};

const VariableInfoWithType COMPUTED_CONSTANT_INFO[] = {
  {"E_L", "millivolt", "leakage_current", COMPUTED_CONSTANT, NOT_EXPOSED,},
  {"E_Na", "millivolt", "sodium_channel", COMPUTED_CONSTANT, NOT_EXPOSED,},
  {"E_K", "millivolt", "potassium_channel", COMPUTED_CONSTANT, NOT_EXPOSED}
};

const VariableInfoWithType ALGEBRAIC_INFO[] = {
  {"i_Stim", "microA_per_cm2", "membrane", ALGEBRAIC, NOT_EXPOSED,}, /* I would have this a constant */
  {"i_L", "microA_per_cm2", "leakage_current", ALGEBRAIC, NOT_EXPOSED,},
  {"i_K", "microA_per_cm2", "potassium_channel", ALGEBRAIC, NOT_EXPOSED,},
  {"i_Na", "microA_per_cm2", "sodium_channel", ALGEBRAIC, EXPOSED,},
  {"alpha_m", "per_millisecond", "sodium_channel_m_gate", ALGEBRAIC, NOT_EXPOSED,},
  {"beta_m", "per_millisecond", "sodium_channel_m_gate", ALGEBRAIC, NOT_EXPOSED,},
  {"alpha_h", "per_millisecond", "sodium_channel_h_gate", ALGEBRAIC, NOT_EXPOSED,},
  {"beta_h", "per_millisecond", "sodium_channel_h_gate", ALGEBRAIC, NOT_EXPOSED,},
  {"alpha_n", "per_millisecond", "potassium_channel_n_gate", ALGEBRAIC, NOT_EXPOSED,},
  {"beta_n", "per_millisecond", "potassium_channel_n_gate", ALGEBRAIC, NOT_EXPOSED}
};


/* These can be kept although I, personally, do not required them as I would handle memory allocation externally. */

double * createStatesArray()
{
    double *res = (double *) malloc(STATE_COUNT*sizeof(double));

    for (size_t i = 0; i < STATE_COUNT; ++i) {
        res[i] = NAN;
    }

    return res;
}

double * createConstantsArray()
{
    double *res = (double *) malloc(CONSTANT_COUNT*sizeof(double));

    for (size_t i = 0; i < CONSTANT_COUNT; ++i) {
        res[i] = NAN;
    }

    return res;
}

double * createComputedConstantssArray()
{
    double *res = (double *) malloc(COMPUTED_CONSTANT_COUNT*sizeof(double));

    for (size_t i = 0; i < COMPUTED_CONSTANT_COUNT; ++i) {
        res[i] = NAN;
    }

    return res;
}

double * createAlgebraicArray()
{
    double *res = (double *) malloc(ALGEBRAIC_COUNT*sizeof(double));

    for (size_t i = 0; i < ALGEBRAIC_COUNT; ++i) {
        res[i] = NAN;
    }

    return res;
}

/* I would handle memory allocation and frees and so I wouldn't call this routine */

void deleteArray(double *array)
{
    free(array);
}

/* This routine is fine. As I would manage memory and initialisation I would try and call libCellML
   routines to work out the initial value and set it myself. I just wouldn't call this routine.
   Note the initialise routine would only init the exposed variables. */

void initialiseVariables(double *states, double *rates, double *constants)
{
    states[0] = 0.0;
    states[1] = 0.6;
    states[2] = 0.05;
    states[3] = 0.325;
    constants[0] = 1.0;
    constants[1] = 120.0;
}

/* Because there are no exposed computed constants this routine would do nothing */

void computeComputedConstants(double *computedConstants)
{
}

/* There can still be a compute rates routine. I would prefer a combined compute rates and variables routine below.
   Note 1: I have just left this routine as before for comparison with the new compute rates and variables below.
   If it was required for the separated compute routines then this would just compute the exposed variables and rates.
   Note 2: it would require the neComputedConstants to be passed in and possibly a neAlgebraic to be passed out. This
   would break the routine interface for backwards compatibility. */

void computeRates(double voi, double *states, double *rates, double *constants, double *computedConstants, double *algebraic)
{
    algebraic[0] = ((voi >= 10.0) && (voi <= 10.5))?-20.0:0.0;
    algebraic[1] = constants[2]*(states[0]-computedConstants[0]);
    algebraic[2] = constants[4]*pow(states[3], 4.0)*(states[0]-computedConstants[2]);
    algebraic[3] = constants[3]*pow(states[2], 3.0)*states[1]*(states[0]-computedConstants[1]);
    rates[0] = -(-algebraic[0]+algebraic[3]+algebraic[2]+algebraic[1])/constants[0];
    algebraic[4] = 0.1*(states[0]+25.0)/(exp((states[0]+25.0)/10.0)-1.0);
    algebraic[5] = 4.0*exp(states[0]/18.0);
    rates[2] = algebraic[4]*(1.0-states[2])-algebraic[5]*states[2];
    algebraic[6] = 0.07*exp(states[0]/20.0);
    algebraic[7] = 1.0/(exp((states[0]+30.0)/10.0)+1.0);
    rates[1] = algebraic[5]*(1.0-states[1])-algebraic[7]*states[1];
    algebraic[8] = 0.01*(states[0]+10.0)/(exp((states[0]+10.0)/10.0)-1.0);
    algebraic[9] = 0.125*exp(states[0]/80.0);
    rates[3] = algebraic[8]*(1.0-states[3])-algebraic[9]*states[3];
}

/* This would now just compute the exposed variables now. */

void computeVariables(double voi, double *states, double *rates, double *constants, double *computedConstants, double *algebraic)
{
  double neComputedConstants[1];

  neComputedConstants[0] = neConstants[0]-10.613;
  algebraic[0] = constants[1]*pow(states[2], 3.0)*states[1]*(states[0]-neComputedConstants[0]);
}

/* I would just require a compute rates and variables routine as we are just passing in exposed variables now.
   Note 1: I could live with calling separate variable an rates routines if I had to but it would require passing in the
   neComputedConstants array which would break backwards comptability etc.
   Note 2: The interface to this routine is just the same as before if we didn't turn the exposure off on any variable. */

void computeRatesAndVariables(double voi, double *states, double *rates, double *constants, double *computedConstants, double *algebraic)
{
     /* Any computed constants that are not exposed are declared locally to the routine. This is so that they are on
     the stack and not on the heap so that the computeRatesAndVariables routine can be called as part of a parallel
     loop via something like OpenMP without problems. The variables are not exposed and so they do not need to come
     out of this routine */  
    double neComputedConstants[NOT_EXPOSED_COMPUTED_CONSTANT_COUNT];
    /* Any alebraic variables that are required for the calculation of other variables that are also not exposed should be
     listed */
    double neAlgebraic[NOT_EXPOSED_ALGEBRAIC_COUNT];

    /* compute constants in this routine rather than call two routines? */
    neComputedConstants[0] = neConstants[0]-10.613;
    neComputedConstants[1] = neConstants[0]-115.0;
    neComputedConstants[2] = neConstants[0]+12.0;
    /* now compute algebraic and rates */
    neAlgebraic[0] = ((voi >= 10.0) && (voi <= 10.5))?-20.0:0.0;
    neAlgebraic[1] = neConstants[1]*(states[0]-neComputedConstants[0]);
    neAlgebraic[2] = neConstants[2]*pow(states[3], 4.0)*(states[0]-neComputedConstants[2]);
    /* The order can be whatever is required for computation. Just the positions in the arrays have been renumbered. */
    algebraic[0] = constants[1]*pow(states[2], 3.0)*states[1]*(states[0]-neComputedConstants[1]);    
    rates[0] = -(-neAlgebraic[0]+algebraic[0]+neAlgebraic[2]+neAlgebraic[1])/constants[0];
    neAlgebraic[3] = 0.1*(states[0]+25.0)/(exp((states[0]+25.0)/10.0)-1.0);
    neAlgebraic[4] = 4.0*exp(states[0]/18.0);
    rates[2] = neAlgebraic[3]*(1.0-states[2])-neAlgebraic[4]*states[2];
    /* This next equation is not actually used anywhere else.
       It could thus be removed so as to not waste FLOPs as it is not exposed. */
    neAlgebraic[5] = 0.07*exp(states[0]/20.0);    
    neAlgebraic[6] = 1.0/(exp((states[0]+30.0)/10.0)+1.0);
    rates[1] = neAlgebraic[4]*(1.0-states[1])-neAlgebraic[6]*states[1];
    neAlgebraic[7] = 0.01*(states[0]+10.0)/(exp((states[0]+10.0)/10.0)-1.0);
    neAlgebraic[8] = 0.125*exp(states[0]/80.0);
    rates[3] = neAlgebraic[7]*(1.0-states[3])-neAlgebraic[8]*states[3];
}