/*
 Implementation of the 1D-model for dynamic instability
 Ludovic Brun and Francois Nedelec, Copyright EMBL 2009
 Version 1.0 on March 19 2009.
 -
 Please, refer to:
 A theory of microtubule catastrophe and their regulation
 Brun L, Rupp B, Ward J, Nedelec F
 PNAS 106 (50) 21173-21178; Dec 2009.
 
 ---
 
 To compile:
 > g++ cap_model.cc -o sim
 This should produce an executable 'sim'
 
 To run:
 > ./sim
 Simulate a dynamic filament until a catastrophe occurs.
 
 > ./sim seed=3
 Simulate one dynamic filament, seeding the random number generator with value 3.
 
 > ./sim  g=4 dt=0.005
 Simulate with specified parameters for the model
 Here, the growth rate g was set to 4 /s and the time-step was set to 0.005 sec.
 The parameters that can be changed are: g, h, N and seed, dt, cnt (see below)
 
 > ./sim  seed=3 cnt=1000 > file
 Simulate 1000 filaments, and outputs the catastrophe times in 'file'.
 You may then use your favorite program to plot the resulting histogram.
 
 */

#include <cstdio>
#include <cstdlib>
#include <vector>
#include <iostream>
#include <sstream>
#include <iomanip>
#include <ctime>


//    PARAMETERS OF THE MODEL:
double g   = 5;             // addition rate of dimers (growth)
double h   = 0.058;         // hydrolysis rate
int    N   = 2;             // coupling parameter

//    VARIABLES OF THE MODEL:
std::vector<int> fiber;     // state vector

/* Note: the state of each unit is 0 or 1
 fiber[0] = minus-end
 ...
 fiber[fiber.size()-1] = ultimate unit
*/

//    VARIABLES OF THE SIMULATION:
double         dt  = 0.005; // time step
int            cnt = 1;     // number of filaments to be simulated consecutively
unsigned long seed = 0;     // seed for random number generator
int           iter = 0;     // counter of iteration


// returns a new random number in [0, 1]
double rng()
{
    const double scale = 1.0 / static_cast<double>(RAND_MAX);
    return static_cast<double>( random() ) * scale; 
}

// initialize with N units at state '1'
void init() 
{
    iter = 0;
    fiber.clear();
    for (int i=0; i<N; ++i)
         fiber.push_back(1);
}


// examine the state to see if it correspond to a catastrophe
bool catastrophe()
{
    // catastrophe is prevented if any '1' is found in the last N units:
    for (int i = fiber.size()-N; i < fiber.size(); ++i) {
        if ( fiber[i] == 1 )
            return false;
    }
    return true;
}


// perform one simulation step, and return true if the state has changed
int step()
{
    ++iter;
    int res = 0;
    
    for (int i=fiber.size()-N; i<fiber.size(); ++i)
    {
        if ( fiber[i] == 1 ) {
            // check for hydrolysis:
            if ( rng() < dt*h )
            {
                res = 2;
                fiber[i] = 0;
            }
        }
    }
    
    // check for growth:
    if ( rng() < dt*g )
    {
        res = 1;
    	// add one unit in state '1'
        fiber.push_back(1);
        // force hydrolysis at distance 'N':
        /* Note: this is not necessary to model catastrophes,
        since units beyond depth N are not considered */
        fiber[fiber.size()-1-N] = 0;
    }
    
    return res;
}

// print the state vector to stream 'os'
void print(std::ostream & os)
{
    os << std::setw(8) << iter*dt << "s ";
    int start = 0;
    if ( fiber.size() > 73 ) {
        start = fiber.size() - 70;
        os<<"...";
    }
    for ( int i=start; i<fiber.size(); ++i )
    {
        if ( fiber[i] )
        	os<<'X';
        else
        	os << (( i % 10 ) ? '-' : '|');
    }
    os<<std::endl;    
}

// reads 'arg' checking for parameter 'par'
template<typename T>
int parse(const char arg[], const char par[], T * ptr)
{
    if ( arg == strstr(arg, par) ) {
        std::istringstream ss(arg+strlen(par));
        ss >> *ptr;
        return 1;
    }
    return 0;
}

// read string 'arg' to change various parameter values
int parse(const char arg[])
{
    if ( parse(arg, "g=",    &g) )    return 0;
    if ( parse(arg, "h=",    &h) )    return 0;
    if ( parse(arg, "N=",    &N) )    return 0;
    if ( parse(arg, "dt=",   &dt) )   return 0;
    if ( parse(arg, "seed=", &seed) ) return 0;
    if ( parse(arg, "cnt=",  &cnt) )  return 0;
    printf("Error: argument '%s' not understood\n", arg);
    return 1;
}


int main(int argc, char* argv[])
{
    // read the command line to set parameters:
    for (int a=1; a<argc; ++a)
        if ( parse(argv[a]) )
            return 1;
    
    // initialize random number generator:
    if ( seed )
    	srandom(seed);
    else
    	srandom(time(0));
    
    // initialize output stream:
    std::cout.precision(3);
    std::cout<<std::fixed;
    
    if ( cnt == 1 )
    {
    	// simulate just one fiber
        init();
        print(std::cout);
        while ( !catastrophe() )
        {
        	// perform one step, and print the state if it has changed:
            if ( step() ) 
            	print(std::cout);
        }
        std::cout<<"Catastrophe!"<<std::endl;
    }
    else
    {
        for ( int c = 0; c < cnt; ++c )
        {
            init();
            while ( !catastrophe() )
                step();
            std::cout << iter*dt << std::endl;
        }
    }
    return 0;
}