/*
 *
 *  $Id$
 *
 *  This file is part of the Virtual Leaf.
 *
 *  The Virtual Leaf 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.
 *
 *  The Virtual Leaf 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>

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;
	};
	~Mesh(void) {
		delete boundary_polygon;
	};
	
	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++) {
			f(*shuffled_nodes[*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++) {
			f(*shuffled_cells[*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);
	inline double getTime(void) const { return time; }
	string getTimeHours(void) const; 
	inline void setTime(double t) { time = t; }
	double CalcProtCellsWalls(int ch) const;  
	void SettoInitVals(void);
	QVector<qreal> VertexAngles(void);
	QVector< QPair<qreal,int> > VertexAnglesValues(void);
	void SetSimPlugin(SimPluginInterface *new_plugin) {
		plugin=new_plugin;
	}
	QString ModelID(void) { return plugin?plugin->ModelID():QString("undefined"); }
	void StandardInit(void);	

	Node* findNextBoundaryNode(Node*);

private:
	
	// Data members
	vector<Cell *> cells;
	vector<Node *> nodes;
	list<Wall *> walls; // we need to erase elements from this container frequently, hence a list.
public:
	vector<NodeSet *> node_sets;
private:
	vector<Node *> shuffled_nodes;
	vector<Cell *> shuffled_cells;
	unique_queue<Edge> node_insertion_queue;
	BoundaryPolygon *boundary_polygon;
	double time;
	SimPluginInterface *plugin;
	
	// Private member functions
	void AddNodeToCell(Cell *c, Node *n, Node *nb1 , Node *nb2);
	void AddNodeToCellAtIndex(Cell *c, Node *n, Node *nb1 , Node *nb2, list<Node *>::iterator ins_pos);
	void InsertNode(Edge &e);
	inline Node *AddNode(Node *n) {
		nodes.push_back(n);
		shuffled_nodes.push_back(n);
		n->m=this;
		return n;
	}
	
	inline Cell *AddCell(Cell *c) {
		cells.push_back(c);
		shuffled_cells.push_back(c);
		//cerr << "Shuffled cell indices:  ";
		/*copy(shuffled_cells.begin(),shuffled_cells.end(),ostream_iterator<int>(cerr," "));
		 cerr << endl;*/
		c->m=this;
		return c;
	}
	
	void CircumCircle(double x1,double y1,double x2,double y2,double x3,double y3,
					  double *xc,double *yc,double *r);
};
#endif