//
// 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 <iostream>
#include <cstdio>
#include <cmath>

#include <DynamXProcessHejnowicz.hpp>
#include <Function.hpp>
#include <Contour.hpp>
#include <Random.hpp>
#include <Information.hpp>
#include <SystemProcessLoad.hpp>

#include <gsl/gsl_vector.h>
#include <gsl/gsl_multiroots.h>

namespace mgx
{
  bool HejnowiczSurfaceGrowth::initialize(QStringList &parms, QWidget * /*parent*/)
  {
    // If there is parent, then we are in gui mode
    //if(parent) ;

    // Set the static class variables each time
    mesh = currentMesh();
    T = &mesh->tissue();

    // Read the parms from the GUI
    processParms(parms);

    // We need these two processes, but they will be defined by the particular
    // model, so we will get them by name that is passed as a parameter.
    if(!getLastParms(surfaceParmsProc, surfaceParms))
      throw(QString("Cannot get surface parameters for process %1").arg(surfaceParmsProc));
    hejnowiczSurface.initialize(vertexData, surfaceParms);

    if(!getLastParms(tissueParmsProc, tissueParms))
      throw(QString("Cannot get tissue parameters for process %1").arg(tissueParmsProc));
    T->processParms(tissueParms);

    if(T->meshType() == "MGX2D") {
      // If we have an incoming mesh, make sure coordinates are OK
      forall(vertex v, T->S)
        hejnowiczSurface.setPoint(v, v, v->pos);
    } else if(T->S.size() > 0)
      throw(QString("Invalid mesh type, must be a 2D cellular mesh"));

    Information::out << "Initialized Hejnowicz Surface Growth" << endl;

    return true;
  } 

  bool HejnowiczSurfaceGrowth::rewind(QStringList &parms, QWidget *parent)
  {
    // Get the mesh load process
    MeshLoad *meshLoad;
    QStringList meshLoadParms;
    if(!(meshLoad = getProcessParms<MeshLoad>(this, meshLoadParms)))
      throw(QString("Unable to make MeshLoad process"));

    meshLoadParms[0] = mesh->file(); // RSS: Should be MeshLoad::FileName
    return meshLoad->run(meshLoadParms);
  } 

  // Process model parameters
  bool HejnowiczSurfaceGrowth::processParms(const QStringList &parms)
  {
    // Get the attribute maps
    vertexData = &mesh->attributes().attrMap<vertex, HejnowiczSurface::VertexData>("HejnowiczVertexData");
		// Get model parms
    dt = parms[pDt].toDouble();
    drawSteps = parms[pDrawSteps].toInt();
    cellKill = parms[pCellKill].toDouble();
    surfaceParmsProc = parms[pSurfaceParmsProc];
    tissueParmsProc = parms[pTissueParmsProc];

    return true;
  }

  // Main model step
	bool HejnowiczSurfaceGrowth::step(const QStringList &parms)
  {    
    for(int i = 0; i < drawSteps; i++) {
      // Advance time
      time() += dt;
  
      // Grow the cells, avoid the centers
      forall(vertex v, T->S) 
        if(!T->getCell(v)) 
          hejnowiczSurface.growPoint(v, dt, time());

      // Update normals, area, centers, and length of cell walls
      T->updGeometry(); 

      // Now update the center's parametric coordinates and look for cells to kill
      CellVec K;
      forall(const cell &c, T->C) {
        vertex v = T->getVtx(c);
        vertex v1 = T->S.anyIn(v);
        // look for closest point on surface to v, starting from vertex already on surface (v1)
        hejnowiczSurface.setPoint(v, v1, v->pos);
        HejnowiczSurface::VertexData &V = v->*vertexData;
        if(V.nrs.x() > cellKill)
          K.push_back(c);
      }
      forall(const cell &c, K)
        T->deleteCell(c); // How to erase from attribute maps
    }
		
		// Update mesh points, edges, surfaces
    mesh->updateAll();   
    return true;
  }

  // Make initial cell, for now no parameters
  bool HejnowiczInitialCell::run(Mesh *mesh, const QString &growthProc,
                       const Point2d &top, const Point2d &bottom, const Point2d &middle)
  { 
    // Create the HejnowiczSurfaceGrowth process
    hejnowiczGrowth = dynamic_cast<HejnowiczSurfaceGrowth *>(makeProcess(growthProc));
    if(!hejnowiczGrowth)
      throw(QString("Cannot create Hejnowicz growth processs: %1").arg(growthProc));
    if(!getLastParms(growthProc, hejnowiczGrowthParms))
      throw(QString("Cannot get parameters for Hejnowicz growth processs: %1").arg(growthProc));

    // Reset the mesh
    mesh->reset();
    mesh->setMeshType("MGX2D");
    CellTissue &T = mesh->tissue();

    // Initialize the growth process
    hejnowiczGrowth->initialize(hejnowiczGrowthParms, 0);

    // Create the initial cell
    hejnowiczGrowth->hejnowiczSurface.initialCell(T, top, bottom, middle);

    // Find center and set cell area and length of cell walls
    T.updGeometry(); 

    // Propagate cell center to vertex
    forall(const cell &c, T.C) {
      vertex v = T.getVtx(c);
      hejnowiczGrowth->hejnowiczSurface.setPoint(v, v, v->pos);
    }

    // Update the GUI
    mesh->updateAll();

    return true;
  }
}
