/*

 *

 *  $Id$

 *

 *  This file is part of the Virtual Leaf.

 *

 *  VirtualLeaf is free software: you can redistribute it and/or modify

 *  it under the terms of the GNU General Public License as published by

 *  the Free Software Foundation, either version 3 of the License, or

 *  (at your option) any later version.

 *

 *  VirtualLeaf is distributed in the hope that it will be useful,

 *  but WITHOUT ANY WARRANTY; without even the implied warranty of

 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

 *  GNU General Public License for more details.

 *

 *  You should have received a copy of the GNU General Public License

 *  along with the Virtual Leaf.  If not, see <http://www.gnu.org/licenses/>.

 *

 *  Copyright 2010 Roeland Merks.

 *

 */

// Cell derives from Vector, where Vector is simply used as a Vertex

#ifndef _MESH_H_

#define _MESH_H_

#include <vector>

#include <algorithm>

#include <queue>

#include <iterator>

#include <functional>

#ifdef QTGRAPHICS

#include <QGraphicsScene>

#endif

#include "cell.h"

#include "node.h"

#include "simplugin.h"

#include <QVector>

#include <QPair>

#include <QDebug>

#include <QTextStream>

using namespace std;

// new queue which rejects duplicate elements

template<class T, class C = deque<T> > class unique_queue : public queue<T,C> {

 public:

 typedef typename C::value_type value_type;

 // reimplements push: reject element if it exists already

 void push(const value_type &x) {

   if (find (queue<T,C>::c.begin(),queue<T,C>::c.end(),x)==queue<T,C>::c.end()) {

     queue<T,C>::c.push_back(x);

   }

 }

 void clear(void) {

   queue<T,C>::c.clear();

 }

};

template<class P> P& deref_ptr ( P *obj) { return *obj; }

class Mesh {

  friend class Cell;

  friend class Node;

  friend class FigureEditor;

 public: 

  Mesh(void) {

    // Make sure the reserved value is large enough if a cell is added

    // in "Divide" when the capacity is insufficient, "cells" might be

    // relocated including the current Cell (i.e. the value of *this)

    // calling "Mesh::IncreaseCapacityIfNecessary" (from another

    // object than Cell, e.g. Mesh) before entering Divide will solve

    // this issue (solved now).

    cells.reserve(2);

    nodes.reserve(500);

    time = 0.;

    plugin = 0;

    boundary_polygon=0;

  };

  ~Mesh(void) {

    if (boundary_polygon) {

      delete boundary_polygon;

      boundary_polygon=0;

    }

  };

  void Clean(void);

  Cell &EllipticCell(double xc, double yc, double ra, double rb, int nnodes=10, double rotation=0);

  Cell &CircularCell(double xc, double yc, double r, int nnodes=10) {

    return EllipticCell(xc, yc, r, r, nnodes, 0);

  }

  Cell &LeafPrimordium(int n, double pet_length);

  Cell &LeafPrimordium2(int n);

  Cell *RectangularCell(const Vector ll, const Vector ur, double rotation = 0);

  void CellFiles(const Vector ll, const Vector ur);

  inline Cell &getCell(int i) {

    if ((unsigned)i<cells.size())

      return *cells[i];

    else {

#ifdef QDEBUG

      qDebug() << i << endl;

      qDebug() << "size is " << cells.size() << endl;

#endif

      abort();

    }

  }

  inline Node &getNode(int i) {

    return *nodes[i];    

  }

  //double Diffusion(void);

  inline int size(void) {

    return cells.size();

  }

  inline int nnodes(void) {

    return nodes.size();

  }

  template<class Op> void LoopCells(Op f) {

    for (vector <Cell *>::iterator i=cells.begin();

	 i!=cells.end();

	 i++) {

      f(**i);

    }

  }

  template<class Op> void LoopWalls(Op f) {

    for (list <Wall *>::iterator i=walls.begin();

	 i!=walls.end();

	 i++) {

      f(**i);

    }

  }

  // if the amount of cells might increase, during looping, use this template

  template<class Op> void LoopCurrentCells(Op f) {

    vector<Cell *> current_cells = cells;

    for (vector <Cell *>::iterator i=current_cells.begin();

	 i!=current_cells.end();

	 i++) {

      f(**i);

    }

  }

  template<class Op> void LoopNodes(Op f) {

    for (vector<Node *>::iterator i=nodes.begin();

	 i!=nodes.end();

	 i++) {

      f(**i); 

    }

  }

  template<class Op> void RandomlyLoopNodes(Op f) {

    MyUrand r(shuffled_nodes.size());

    random_shuffle(shuffled_nodes.begin(),shuffled_nodes.end(),r);

    for (vector<Node *>::const_iterator i=shuffled_nodes.begin();

	 i!=shuffled_nodes.end();

	 i++) {

    }

  }

  template<class Op> void RandomlyLoopCells(Op f) {

    MyUrand r(shuffled_cells.size());

    random_shuffle(shuffled_cells.begin(),shuffled_cells.end(),r);

    for (vector<Cell *>::const_iterator i=shuffled_cells.begin();

	 i!=shuffled_cells.end();

	 i++) {

    }

  }

  template<class Op1, class Op2> void LoopCells(Op1 f, Op2 &g) {

    for (vector<Cell *>::iterator i=cells.begin();

	 i!=cells.end();

	 i++) {

      f(**i,g); 

    }

  }

  template<class Op1, class Op2, class Op3> void LoopCells(Op1 f, Op2 &g, Op3 &h) {

    for (vector<Cell *>::iterator i=cells.begin();

	 i!=cells.end();

	 i++) {

      f(**i,g,h); 

    }

  }

  void DoCellHouseKeeping(void) {

    vector<Cell *> current_cells = cells;

    for (vector<Cell *>::iterator i = current_cells.begin();

	 i != current_cells.end();

	 i ++) {

      plugin->CellHouseKeeping(*i);

      // Call functions of Cell that cannot be called from CellBase, including Division

      if ((*i)->flag_for_divide) {

	if ((*i)->division_axis) {

	  (*i)->DivideOverAxis(*(*i)->division_axis);

	  delete (*i)->division_axis;

	  (*i)->division_axis = 0;

	} else {

	  (*i)->Divide();

	}

	(*i)->flag_for_divide=false;

      }

    }

  }

  // Apply "f" to cell i

  // i.e. this is an adapter which allows you to call a function

  // operating on Cell on its numeric index index

  template<class Op> void cell_index_adapter(Op f,int i) {

    f(cells[i]);

  }

  double DisplaceNodes(void);

  void BoundingBox(Vector &LowerLeft, Vector &UpperRight);

  int NEqs(void) {     int nwalls = walls.size();

    int ncells =cells.size();

    int nchems = Cell::NChem();

    // two eqs per chemical for each walls, and one eq per chemical for each cell

    // This is for generality. For a specific model you may optimize

    // this by removing superfluous (empty) equations.

    int neqs = 2 * nwalls * nchems + ncells * nchems;

    return neqs;

  }

  void IncreaseCellCapacityIfNecessary(void) {

    return;

    // cerr << "Entering Mesh::IncreaseCellCapacityIfNecessary \n";

    // make sure we always have enough space 

    // to have each cell divide at least once

    //

    // Note that we must do this, because Cell::Divide pushes a new Cell

    // onto Mesh::cells. As a result, Mesh::cells might be relocated 

    // if we are _within_ a Cell object: i.e. pointer "this" will be changed!!

    // 

    // An alternative solution could be to make "Mesh::cells" a list,

    // but this won't work because we need random access for 

    // the Monte Carlo algorithm.

    if (2*cells.size()>cells.capacity()) {

      cerr << "Increasing capacity to "  << 2*cells.capacity() << endl;

      cerr << "Current capacity is " << cells.capacity() << endl;

      cells.reserve(cells.capacity()*2);

    }

  }

  void ReserveMoreCells(int n) {

    if (nodes.size()+n>nodes.capacity()) {

      nodes.reserve(size()+n);

    }

  }

  double Area(void);

  double MeanArea(void) {

    double sum=0.;

    for (vector<Cell *>::const_iterator i=cells.begin();

	 i!=cells.end();

	 i++) {

      sum+=(*i)->Area();

    }

    return sum/(double)NCells();

  }

  void SetBaseArea(void);

  int NCells(void) const {

    return cells.size();

  }

  inline int NNodes(void) const {

    return nodes.size();

  }

  void PrintQueue(ostream &os) {

    while (!node_insertion_queue.empty()) {

      os << node_insertion_queue.front() << endl;

      node_insertion_queue.pop();

    }

  }

  void InsertNodes(void) {

    // insert the nodes in the insertion queue

    while (!node_insertion_queue.empty()) {

      //cerr << node_insertion_queue.front() << endl;

      InsertNode(node_insertion_queue.front());

      node_insertion_queue.pop();

    }

  }

  void Clear(); 

  void ReactDiffuse( double delta_t = 1 );

  double SumChemical(int ch);

  void SetChemical(int ch, double value) {

    for (vector<Cell *>::iterator c=cells.begin();

	 c!=cells.end();

	 c++) {

      (*c)->chem[ch]=value;

    }

  }

  // used for interacing with ODE-solvers (e.g. NRCRungeKutta)

  void setValues(double x, double *y);

  double *getValues(int *neqs);

  void Derivatives(double *derivs);

#ifdef QTGRAPHICS

  inline void DrawBoundary(QGraphicsScene *c) {

    boundary_polygon->Draw(c);

  }

  void DrawNodes(QGraphicsScene *c) const;

#endif

  double max_chem;

  void XMLSave(const char *docname, xmlNode *settings=0) const;

  void XMLRead(const char *docname, xmlNode **settings=0, bool geometry = true, bool pars = true, bool simtime = true);

  void XMLReadPars(const xmlNode * root_node);

  void XMLReadGeometry(const xmlNode *root_node);

  void XMLReadSimtime(const xmlNode *root_node);

  void XMLReadNodes(xmlNode *cur);

  void XMLReadCells(xmlNode *cur);

  void XMLParseTree(const xmlNode * root_node);

  void XMLReadWalls(xmlNode *cur, vector<Wall *> *tmp_cells);

  void XMLReadWallsToCells(xmlNode *root, vector<Wall *> *tmp_walls);

  void XMLReadNodeSets(xmlNode *root);

  void XMLReadNodeSetsToNodes(xmlNode *root);

  void PerturbChem(int chemnum, double range);

  void CleanUpCellNodeLists(void);

  void CleanUpWalls(void);

  void CutAwayBelowLine( Vector startpoint, Vector endpoint );

  void CutAwaySAM(void);

  void RepairBoundaryPolygon(void);

  void Rotate(double angle, Vector center);

  void PrintWallList( void );

  void TestIllegalWalls(void);

  Vector FirstConcMoment(int chem);

  inline Vector Centroid(void) {

    return boundary_polygon->Centroid();

  }

  inline Vector Offset(void) {

    return boundary_polygon->Offset();

  }

  inline double Factor(void) {

    return boundary_polygon->Factor();

  }

  void DeleteLooseWalls(void);

  void FitLeafToCanvas(double width, double height);

  void AddNodeSet(NodeSet *node_set) {

    node_sets.push_back(node_set);

  }

  void CleanChemicals(const vector<double> &clean_chem);

  void CleanTransporters(const vector<double> &clean_transporters);

  void RandomizeChemicals(const vector<double> &max_chem, const vector<double> &max_transporters);