#

#  $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.

#

CONFIG += qt

QMAKE_CXXFLAGS += -fexceptions

QMAKE_CXXFLAGS_DEBUG += -g3

QMAKE_CXXFLAGS_DEBUG += -DQDEBUG

#REACTIONS = reactions_auxin_growth.h 

#REACTIONS = reactions_meinhardt.h

#REACTIONS = reactions_pce_growth.h

DEFINES += QTGRAPHICS

DEFINES += REACTIONS_HEADER=$${REACTIONS}

DEFINES += REACTIONS_HEADER_STRING=\"$${REACTIONS}\"

DEFINES += FLEMING

PERLDIR = ./perl

BINDIR = ../bin

DESTDIR = $$BINDIR

TARGET = VirtualLeaf

TEMPLATE = app

PARTMPL = $${TARGET}par.tmpl

MAINSRC = $${TARGET}.cpp

QT -= network sql xml

QT += qt3support

!win32 {

 GRAPHICS = qt #qwt

}

win32 {

 CONFIG += console

 LIBXML2DIR = C:\libxml2

 LIBICONVDIR = C:\libiconv

 LIBZDIR = C:\libz

 system(DEL parameter.cpp parameter.h) 

 GRAPHICS = qt 

 RC_FILE = VirtualLeaf.rc

 QMAKE_CXXFLAGS += -DLIBXML_STATIC

 QMAKE_CXXFLAGS += -I$${LIBXML2DIR}\include -I$${LIBICONVDIR}\include -I$${LIBZDIR}\include

 QMAKE_POST_LINK = "\

  C:\Bin\cp release\VirtualLeaf.exe \

  C:\Qt\4.5.3\bin\Qt3Support4.dll \

  C:\Qt\4.5.3\bin\QtGui4.dll \

  C:\Qt\4.5.3\bin\QtSql4.dll \

  C:\Qt\4.5.3\bin\QtXml4.dll \

  C:\Qt\4.5.3\bin\QtCore4.dll \

  C:\Qt\4.5.3\bin\QtNetwork4.dll \

  C:\Qt\4.5.3\bin\QtSvg4.dll \

  C:\bin\iconv.dll \

#

# $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.

#

CONFIG += release

CONFIG -= debug

CONFIG += plugin

BINDIR = ../../bin

DEFINES = QTGRAPHICS # VLEAFPLUGIN

DESTDIR = $${BINDIR}/models

TARGET = test

HEADERS = ../simplugin.h $${TARGET}plugin.h  

QMAKE_CXXFLAGS += -fexceptions -I..

QMAKE_CXXFLAGS_DEBUG += -g3

QMAKE_CXXFLAGS_DEBUG += -DQDEBUG

QT += qt3support

SOURCES = $${TARGET}plugin.cpp

TEMPLATE = lib 

unix {

 LIBS += -Llib -lvleaf

 QMAKE_CXXFLAGS += -fPIC -I/usr/include/libxml2

 QMAKE_LFLAGS += -fPIC

}

win32 {

 LIBXML2DIR = C:\libxml2

 LIBICONVDIR = C:\libiconv

 LIBZDIR = C:\libz

 LIBS += -Llib -lvleaf

 QMAKE_CXXFLAGS += -DLIBXML_STATIC

 QMAKE_CXXFLAGS += -I$${LIBXML2DIR}\include -I$${LIBICONVDIR}\include -I$${LIBZDIR}\include

}

# finis

	ItList::iterator i = new_node_locations.begin();

	list< Node *>::iterator j;

	cerr << "-------------------------------\n";

	cerr << "Location of new nodes: " << (**i)->Index() << " and ";

	++i;

	j = *i; 

	if (j==nodes.begin()) j=nodes.end(); j--;

	cerr << (*j)->Index() << endl;

	cerr << "-------------------------------\n";

	if ( **new_node_locations.begin() == *j ) {

		cerr << "Rejecting proposed division (cutting off zero area).\n";

		return false;

	}

	DivideWalls(new_node_locations, v1, v2, fix_cellwall, node_set);

	return true;

}

// Core division procedure

void Cell::DivideWalls(ItList new_node_locations, const Vector from, const Vector to, bool fix_cellwall, NodeSet *node_set) {

	if (dead) return;

	bool boundary_touched_flag=false;

	// Step 0: keep some data about the parent before dividing

	ParentInfo parent_info;

	parent_info.polarization = ReduceCellAndWalls<Vector>( PINdir );

	parent_info.polarization.Normalise();

	parent_info.PINmembrane = SumTransporters(1);

	parent_info.PINendosome = Chemical(1);

	//cerr << "Parent polarization before division: " << parent_info.polarization << endl;

	// Step 1: create a daughter cell

	Cell *daughter=m->AddCell(new Cell());

	// Step 2: Copy the basics of parent cell to daughter

	for (int i=0;i<NChem();i++) {

		daughter->chem[i]=chem[i];

	}

	//extern double auxin_account;

	//auxin_account += daughter->chem[0];

	for (int i=0;i<NChem();i++) {

		daughter->new_chem[i]=new_chem[i];

	}

	daughter->boundary=boundary;

	daughter->m=m;

	daughter->target_area=target_area/2.;

	target_area/=2;

	daughter->cellvec=cellvec;

//	daughter->BaseArea()  = base_area;

	// Division currently only works for convex cells: i.e. if the division line

	// intersects the cells at two points only.

	if (new_node_locations.size()!=2) {

		// Note: if you would add the possibility of dividing non-convex

		// cells, remember to update the code below. There are some

		// fixed-size arrays over there!

		cerr << "Warning in Cell::Division: division of non-convex cells not fully implemented" << endl;

		// Reject the daughter cell and decrement the amount of cells

		// again. We can do this here because it is the last cell added.

		// Never, ever try to fully erase a cell elsewhere, because we

		// make heavy use of cell indices in this project; if you erase a

		// Cell somewhere in the middle of Mesh::Cells the indices will

		// get totally messed up...! (e.g. the indices used in Nodes::cells)

		cerr << "new_node_locations.size() = " << new_node_locations.size() <<endl;

		cerr << "daughter->index = " << daughter->index << endl;

		cerr << "cells.size() = " << m->cells.size() << endl;

		m->cells.pop_back();

		Cell::NCells()--;

		m->shuffled_cells.pop_back();

		return;

	}

	// We can be sure we only need two positions here because divisions

	// of non-convex cells are rejected above.

	Vector new_node[2];

const char* CellBase::boundary_type_names[4] = {"None", "NoFlux", "SourceSink", "SAM"};

// These statics have moved to class "CellsStaticDatamembers"

//double CellBase::static_base_area = 0.;

//int CellBase::ncells=0;

//int CellBase::NChem()=0;

#ifndef VLEAFPLUGIN

CellsStaticDatamembers *CellBase::static_data_members = new CellsStaticDatamembers();

#else

CellsStaticDatamembers *CellBase::static_data_members = 0;

#endif

CellBase::CellBase(QObject *parent) : 

QObject(parent),

Vector()

{

	chem=new double[NChem()];

	for (int i=0;i<NChem();i++) {

		chem[i]=0.;

	}

	new_chem=new double[NChem()];

	for (int i=0;i<NChem();i++) {

		new_chem[i]=0.;

	}

	boundary=None;

	index=(NCells()++);

	area=0.;

	target_area=1;

	target_length=0; //par.target_length;

	lambda_celllength = 0; //par.lambda_celllength;

	intgrl_xx=0.; intgrl_xy=0.; intgrl_yy=0.;

	intgrl_x=0.; intgrl_y=0.;

	source = false;

	source_conc = 0.;

	source_chem = 0;

	at_boundary=false;

	fixed = false;

	pin_fixed = false;

	stiffness = 0;

	marked = false;

	dead = false;

	div_counter=0;

	cell_type = 0;

	flag_for_divide = false;

}

CellBase::CellBase(double x,double y,double z) : QObject(), Vector(x,y,z)

{

#ifndef VLEAFPLUGIN

	if (static_data_members == 0) {

		static_data_members = new CellsStaticDatamembers();

	}

#endif

	chem=new double[NChem()];

	for (int i=0;i<NChem();i++) {

		chem[i]=0.;

	}

	new_chem=new double[NChem()];

	for (int i=0;i<NChem();i++) {

		new_chem[i]=0.;

	}

	boundary=None;

	area=0.;

	target_area=1;

	target_length=0; //par.target_length;

	lambda_celllength=0; // par.lambda_celllength;

	index=(NCells()++);

	intgrl_xx=0.; intgrl_xy=0.; intgrl_yy=0.;

	intgrl_x=0.; intgrl_y=0.;

	source = false;

	fixed = false;

	at_boundary=false;

	pin_fixed = false;

	stiffness = 0;

	marked=false;

	dead  = false;

	div_counter = 0;

	cell_type = 0;

	flag_for_divide = false;

}

CellBase::CellBase(const CellBase &src) :  Vector(src), QObject()

{

	chem=new double[NChem()];

	for (int i=0;i<NChem();i++) {

		chem[i]=src.chem[i];

	}

	new_chem=new double[NChem()];

	for (int i=0;i<NChem();i++) {

		new_chem[i]=src.new_chem[i];

	}

	boundary=src.boundary;

	area=src.area;

	target_length=src.target_length;

	lambda_celllength=src.lambda_celllength;

	intgrl_xx=src.intgrl_xx; intgrl_xy=src.intgrl_xy; intgrl_yy=src.intgrl_yy;

	intgrl_x=src.intgrl_x; intgrl_y=src.intgrl_y;

	target_area=src.target_area;

	index=src.index;

	nodes=src.nodes;

	neighbors=src.neighbors;

	walls=src.walls;

	source = src.source;

	fixed = src.fixed;

	source_conc = src.source_conc;

	source_chem = src.source_chem;

	cellvec = src.cellvec;

	at_boundary=src.at_boundary;

	pin_fixed = src.pin_fixed;

	stiffness = src.stiffness;

	marked = src.marked;

	dead = src.dead;

	cell_type = src.cell_type;

	div_counter = src.div_counter;

	flag_for_divide = src.flag_for_divide;

}

CellBase CellBase::operator=(const CellBase &src) {

	Vector::operator=(src);

	//  QObject::operator=(src);

	for (int i=0;i<NChem();i++) {

		chem[i]=src.chem[i];

	}

	for (int i=0;i<NChem();i++) {

		new_chem[i]=src.chem[i];

	}

	boundary=src.boundary;

	area=src.area;

	intgrl_xx=src.intgrl_xx; intgrl_xy=src.intgrl_xy; intgrl_yy=src.intgrl_yy;

	intgrl_x=src.intgrl_x; intgrl_y=src.intgrl_y;

	target_area=src.target_area;

	target_length=src.target_length;

	lambda_celllength=src.lambda_celllength;

	index=src.index;

	nodes=src.nodes;

	neighbors=src.neighbors;

	walls=src.walls;

	source = src.source;

	fixed = src.fixed;

	source_conc = src.source_conc;

	source_chem = src.source_chem;

	cellvec = src.cellvec;

	at_boundary=src.at_boundary;

	pin_fixed = src.pin_fixed;

	stiffness = src.stiffness;

	marked = src.marked;

	dead = src.dead;

	cell_type = src.cell_type;

	div_counter = src.div_counter;

	flag_for_divide = src.flag_for_divide;

	return *this;

}

void CellBase::SetChemical(int c, double conc) {

	if (c>=NChem()) {

		stringstream error;

		error << "SetChemical: value c = " << c << " is out of range\n";

		throw error.str().c_str();

	}

	chem[c]=conc;

}

ostream &CellBase::print(ostream &os) const {

	os << "[ index = " << index << " {" << x << ", " << y << ", " << z << "}: {";

	for (int i=0;i<NChem()-1;i++) {

		os << chem[i] << ", ";

	}

	os << chem[NChem()-1] << " } ]";

	os << endl << "Nodelist = { " << endl;

	for (list<Node *>::const_iterator i =  nodes.begin(); i!=nodes.end(); i++) {

		os << (*i)->Index() << "( " << *i << ") ";

	}

	os << " } ";

	for (list<Wall *>::const_iterator i =  walls.begin(); i!=walls.end(); i++) {

		(*i)->print(os);

		os << ", ";

	} 

	os << endl;

	os << " [ area = " << area << " ]";

	os << " [ walls = ";

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

		 i!=walls.end();

		 i++) {

		os << (*i)->n1->Index() << " -> " << (*i)->n2->Index() << ", " <<  (*i)->c1->Index() << " | " << (*i)->c2->Index() << ", ";

	}

	os << " ] ";

	os << "div_counter = " << div_counter << endl;

	os << "cell_type = " << cell_type << endl;

	os << endl;

	Vector polarization;

	double PINmembrane;

	double PINendosome;

};

// We need a little trick here, to make sure the plugin and the main application will see the same static datamembers

// otherwise each have their own instantation.

// My solution is as follow. I collect all original statics in a class. The main application instantiates it and

// has a static pointer to it. After loading the plugin I set a static pointer to the same class 

class CellsStaticDatamembers {

public:

	CellsStaticDatamembers(void) {

		ncells = 0;

		nchem = 0;

		base_area = 0.;

		cerr << "Constructor of CellsStaticDatamembers\n";

	}

	~CellsStaticDatamembers() {

		cerr << "Oops! Desctructor of CellsStaticDatamembers called\n";

	}

	int ncells;

	int nchem;

	double base_area;

};

class CellBase :  public QObject, public Vector 

{

	Q_OBJECT

	friend class Mesh;

	friend class CellInfo;

	friend class Node;

	friend class WallBase;

	friend class SimPluginInterface;

 public:

	CellBase(QObject *parent=0);

	CellBase(double x,double y,double z=0); // constructor

  virtual ~CellBase() {

	  delete[] chem;

	  delete[] new_chem;

	  //cerr << "CellBase " << index << " is dying. " << endl;

  }

	CellBase(const CellBase &src); // copy constructor

	virtual bool BoundaryPolP(void) const { return false; } 

	//  CellBase(const Vector &src); // not allowed (we cannot know to which mesh 

	/// the CellBase will belong...)

    CellBase operator=(const CellBase &src); // assignment operator

    CellBase operator=(const Vector &src);

    void SetChemical(int chem, double conc);

    inline void SetNewChem(int chem, double conc) { 

      new_chem[chem] = conc;

    }

    void SetSource(int chem, double conc) {

      source=true;

      source_chem = chem;

      source_conc = conc;

	}

    void UnfixNodes(void);

    void FixNodes(void);

    void UnsetSource(void) {

      source = false;

    }

	inline bool Source(void) { return source; }

    enum boundary_type {None, Noflux, SourceSink, SAM};

    static const char * boundary_type_names[4];

    inline const char *BoundaryStr(void) { return boundary_type_names[boundary]; }

    ostream &print(ostream &os) const;

    inline double Chemical(int c) const { // returns the value of chemical c

    #ifdef _undefined_

      qDebug() << endl << "Entering cellbase::chemical(" << c << "), and nchem is: " << NChem() << "." << endl;

    #endif

      int nchem = NChem();

      #ifdef _undefined_

      if ((c<0) || (c>=nchem))

	MyWarning::warning("CellBase::Chemical says: index c is: %d, but nchem is: %d. Merely return zero", c, nchem);

      #endif

    double RecalcArea(void) {

      return area = CalcArea();

    }

    Vector Centroid(void) const;

    void SetIntegrals(void) const;

    double Length(Vector *long_axis = 0, double *width = 0) const;

    double CalcLength(Vector *long_axis = 0, double *width = 0) const;

    inline int Index(void) const {

      return index;

    }

    void SetTargetArea(double tar_ar) {

      target_area=tar_ar;

    }

    inline void SetTargetLength(double tar_l) {

      target_length=tar_l;

    }

    inline void SetLambdaLength(double lambda_length) {

      lambda_celllength = lambda_length;

    }

    inline double TargetArea(void) {

      return target_area;

    }

    inline void SetStiffness(double stiff) {

      stiffness = stiff;

    }

    inline double Stiffness(void) {

      return stiffness;

    }

    inline double EnlargeTargetArea(double da) {

      return target_area+=da;

    }

    inline double Area(void) const {

      return area;

    }

	inline void Divide(void) {

		flag_for_divide = true;

	}

    inline double Circumference(void) const {

      double sum=0.;

      for (list<Wall *>::const_iterator w=walls.begin();

	   w!=walls.end();

	   w++) {

	sum +=  (*w)->Length();

      }

      return sum;

    }

	QList<WallBase *> getWalls(void) {

		QList<WallBase *> wall_list;

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

			 i!=walls.end();

			 i++) {

			wall_list << *i;

		}

		return wall_list;

	}

  //  void XFigPrint(std::ostream &os) const;

    void Dump(ostream &os) const;

	QString printednodelist(void);

   // void OnDivide(ParentInfo &parent_info, CellBase &daughter);

    inline bool DeadP(void) { return dead; }

    inline void MarkDead(void) { dead  = true; }

	static double &BaseArea(void) { 

		return static_data_members->base_area;

	}

    void CheckForDivision(void);

    // write flux from neighboring cells into "flux"

    void Flux(double *flux, double *D); 

    inline bool FixedP(void) { return fixed; }

    inline bool Fix(void) {  FixNodes(); return (fixed=true); }

    inline bool Unfix(void) { UnfixNodes(); return (fixed=false);}

    inline void setCellVec(Vector cv) { cellvec = cv; }

	bool AtBoundaryP(void) const;

    int index;

	inline void SetChemToNewchem(void) {

		for (int c=0;c<CellBase::NChem();c++) {

			chem[c]=new_chem[c];

		}

    }

    inline void SetNewChemToChem(void) {

		for (int c=0;c<CellBase::NChem();c++) {

			new_chem[c]=chem[c];

		}

    }

	inline double NewChem(int c) const { return new_chem[c]; }

 protected:

    list<Node *> nodes;

    void ConstructNeighborList(void);

	long wall_list_index (Wall *elem) const;

    // DATA MEMBERS

    // list of nodes, in clockwise order

    // a (non-ordered) list of neighboring cells (actually I think the

    // introduction of ConstructWalls() has made these

    // lists ordered (clockwise), but I am not yet 100% sure...).

    list<CellBase *> neighbors;

    list<Wall *> walls;

    double *chem;

    double *new_chem;

    boundary_type boundary;

    mutable double area;

    double target_area;

    double target_length;

    double lambda_celllength;

    double stiffness; // stiffness like in Hogeweg (2000)

    bool fixed;

    bool pin_fixed;

    bool at_boundary; 

    bool dead; 

	bool flag_for_divide;

	int cell_type;

    //double length;

    //Vector meanflux;

    //int valence;

    // for length constraint

    mutable double intgrl_xx, intgrl_xy, intgrl_yy, intgrl_x, intgrl_y;

    bool source;

    Vector cellvec;

	// STATIC DATAMEMBERS MOVED TO CLASS

	static CellsStaticDatamembers *static_data_members;

	double source_conc;

    int source_chem;

    // PRIVATE MEMBER FUNCTIONS

    inline static void ClearNCells(void) {

      NCells()=0;

    }

    bool marked;

    int div_counter;

};

ostream &operator<<(ostream &os, const CellBase &v);

inline Vector PINdir(CellBase &here, CellBase &nb, Wall &w) {

	return w.getTransporter( &here, 1)  *  w.getInfluxVector(&here);

}

#endif

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

       w!=walls.end();

       w++) {

    for (int i=0;i<Cell::NChem();i++) {

      (*w)->setTransporters1(i,max_transporters[i] * RANDOM()); (*w)->setNewTransporters1(i, (*w)->Transporters1(i) );

      (*w)->setTransporters2(i,max_transporters[i] * RANDOM()); (*w)->setNewTransporters2(i, (*w)->Transporters1(i) );

    }

  }

}

//!\brief Calculates a vector with derivatives of all variables, which

// we can pass to an ODESolver. 

void Mesh::Derivatives(double *derivs) {

  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;

  //static double *derivs = 0; 

  // derivs is allocated by RungeKutta class.

  for (int i=0;i<neqs;i++) { derivs[i]=0.;}

  // Layout of derivatives: cells [ chem1 ... chem n]  walls [ [ w1(chem 1) ... w1(chem n) ] [ w2(chem 1) ... w2(chem n) ] ]

  int i=0;

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

       c!=cells.end();

       c++) {

    //(*cr)(*c, &(derivs[i]));

	  plugin->CellDynamics(*c, &(derivs[i]));

	  i+=nchems;

  }

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

       w!=walls.end();

       w++) {

   // (*wr)(*w, &(derivs[i]), &(derivs[i+nchems]));

	  plugin->WallDynamics(*w,  &(derivs[i]), &(derivs[i+nchems]));

    // Transport function adds to derivatives of cell chemicals