//
// This file is part of MorphoGraphX - http://www.MorphoGraphX.org
// Copyright (C) 2012-2016 Richard S. Smith and collaborators.
//
// If you use MorphoGraphX in your work, please cite:
//   http://dx.doi.org/10.7554/eLife.05864
//
// MorphoGraphX is free software, and is licensed under under the terms of the 
// GNU General (GPL) Public License version 2.0, http://www.gnu.org/licenses.
// 

#include <DynamXProcessMorphogens.hpp>
#include <Random.hpp>

using namespace std; 

namespace mgx
{
	// Update concentrations in divided cells.  
	bool MeinhardtAI::Subdivide::updateCellData(cell c, cell cl, cell cr)
  {
    MeinhardtAI::CellData &C = c->*cellAttr;
    MeinhardtAI::CellData &Cl = cl->*cellAttr;
    MeinhardtAI::CellData &Cr = cr->*cellAttr;
    Cl.A = Cr.A = C.A;
    Cl.H = Cr.H = C.H;

    // Clear the attribute from the old cell
    cellAttr->erase(c);

    return true;
  }

  bool MeinhardtAI::Solver::setValues(const cell& c, const std::vector<double> &values)
  {
    MeinhardtAI::CellData &C = c->*cellAttr;
    C.A = values[0];
    C.H = values[1];

    return true;
  }

  std::vector<double> MeinhardtAI::Solver::getValues(const cell& c)
  {

  MeinhardtAI::CellData &C = c->*cellAttr;
  std::vector<double> values;
  values.push_back(C.A);
  values.push_back(C.H);

  return values;

  }

  std::vector<double> MeinhardtAI::Solver::getDerivatives(const cell& c)
  {

  MeinhardtAI::CellData &C = c->*cellAttr;

  std::vector<double> derivs;
  derivs.push_back(C.dA);
  derivs.push_back(C.dH);

  return derivs;

  }

  void MeinhardtAI::Solver::setValues(const cell& c, std::vector<double> values)
  {

  MeinhardtAI::CellData &C = c->*cellAttr;

  C.A = values[0];
  C.H = values[1];

  }

  void MeinhardtAI::Solver::updateDerivatives(const cell& c)
  {
    m->calcDerivatives(c);
  }
	
	// Put the signal in the graph (and the heat) and update the bounds
  void MeinhardtAI::updateSignal()
  {
    forall(const vertex &v, T->S) {
      if(v->cId) {
        cell c = T->getCell(v);
        MeinhardtAI::CellData &C = c->*cellAttr;
        v->signal = C.A;
      } else 
        v->signal = 0;
    }
		// Set signal bounds
    mesh->signalBounds() = Point2f(-aViewMax *.2, aViewMax);

		// Update mesh signal
    mesh->updateTriangles();
  }

	// Initialize model here, called from GUI thread by
	bool MeinhardtAI::initialize(QStringList &parms, QWidget* parent)
  {
		// Get the mesh, tissue and attribute map
	  mesh = currentMesh();
    T = &mesh->tissue();

    // Check the mesh type
    if(T->meshType() != "MGX2D")
      throw(QString("Invalid mesh type, must be a 2D cellular mesh"));

		// Process the parameters
    processParms(parms);

    // Set the time for statistics
    lastStat = time;

		return true;
  }

  // Process model parameters
  void MeinhardtAI::processParms(const QStringList &parms)
  {
    cellAttr = &mesh->attributes().attrMap<cell, MeinhardtAI::CellData>("MeinhardtAICellData");

    drawSteps = parms[pDrawSteps].toInt();
    statStep = parms[pStatStep].toDouble();
    dtRD = parms[pDt].toDouble();
    
		aProd = parms[pAProd].toDouble();
    aaProd = parms[pAAProd].toDouble();
    aDecay = parms[pADecay].toDouble();
    aDiff = parms[pADiff].toDouble();
    aNoise = parms[pANoise].toDouble();
    aMax = parms[pAMax].toDouble();
    aViewMax = parms[pAViewMax].toDouble();
  
    hProd = parms[pHProd].toDouble();
    haProd = parms[pHAProd].toDouble();
    hDecay = parms[pHDecay].toDouble();
    hDiff = parms[pHDiff].toDouble();
    hMax = parms[pHMax].toDouble();

    solvingMethod = stringToBool(parms[pSolvingMethod]);

    //std::cout << "done " << std::endl;
  }

  bool MeinhardtAI::rewind(QStringList &parms, QWidget *parent)
  {
    // Re-process the parms
    processParms(parms);

    // Clear attribute map for RD data
		cellAttr->clear();

    // Set simulation time to 0
    time = 0;

    // Write signal to graph 
    updateSignal();

    mesh->signalBounds() = Point2f(-aViewMax *.2, aViewMax);

    return true;
  }

  void MeinhardtAI::calcDerivatives(const cell& c)
  {
     MeinhardtAI::CellData &C = c->*cellAttr;
     // Activator
     // if cell production of activator in on
     if(C.ProdA) 
       C.dA = aaProd * C.A * C.A/(1.0 + C.H) + aProd;
     else
       C.dA = 0;
     C.dA -= aDecay * C.A;
  
     // Add Some Noise
     if(aNoise > 0.0) {
       double amt = ran(aNoise) - .5 * aNoise;
       if(amt + C.A > 0.0)
         C.dA += amt;
     }
  
     // Inhibitor
     C.dH = haProd * C.A * C.A + hProd - hDecay * C.H;
  
     // Diffusion
     forall(const cell &n, T->C.neighbors(c)) {
       MeinhardtAI::CellData &N = n->*cellAttr;

       // Activator diffusion
       C.dA += aDiff * (N.A - C.A) * T->C.edge(c,n)->length/c->area;
  
       // Inhibitor diffusion
       C.dH += hDiff * (N.H - C.H) * T->C.edge(c,n)->length/c->area;
     }

  }

  // Run reaction-diffusion step
  bool MeinhardtAI::step(const QStringList &parms)
	{

    //std::cout << "bla " << std::endl;
    // Print stats at time 0
    if(statStep > 0 and time == 0)
      checkMaxQtys(T);
      
    for(int i = 0; i < drawSteps; i++) {
      // Substance simulation
      forall(const cell &c, T->C) {
        calcDerivatives(c);
      }
      SolvingMethod s = Euler;

      //  std::cout << "bla " << solvingMethod.toLocal8Bit().constData() << std::endl;

      if(!solvingMethod) s = BackwardEuler;
      Solver solver(cellAttr, mesh, T->C, s, 2, dtRD);
      solver.m = this;

      solver.solve();
      // Advance time
      time += dtRD;
  
      // Clip to max, check < 0
      checkMaxQtys(T);
    }
    // Write signal to graph
    updateSignal();

		return true; 
	}

  // Clip max qtys and report
  void MeinhardtAI::checkMaxQtys(CellTissue *T) 
  {
    double maxA = 0.0, maxH = 0.0;
 
    forall(const cell &c, T->C) {
      MeinhardtAI::CellData &C = c->*cellAttr;

      // Can't go below 0, also check for nan
      if(C.A >= 0.0) ; else {
        Information::out << "Error Activator = " << C.A << endl;
        C.A = 0.0;
      }
      if(C.H >= 0.0) ; else {
        Information::out << "Error Inhibitor = " << C.H << endl;
        C.H = 0.0;
      }
  
      // Can't go over Max
      if(C.A > aMax)
        C.A = aMax;
      if(C.H > hMax)
        C.H = hMax;
  
      // Find Max value to print out
      if(C.A > maxA)
        maxA = C.A;
      if(C.H > maxH)
        maxH = C.H;
    }
  
    if(statStep > 0 and (time == 0 or time - lastStat >= statStep)) {
      lastStat += statStep;
      Information::out << "Time:" << time << " Max Act:" << maxA << " Max Inh: " << maxH 
                       << " Cells: " << T->C.size() << endl;
    }
  }
}
