// Gmsh - Copyright (C) 1997-2013 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <gmsh@geuz.org>.
//
// Contributed by Matti Pellikka <matti.pellikka@tut.fi>.

#ifndef _CELLCOMPLEX_H_
#define _CELLCOMPLEX_H_

#include <map>
#include <string.h>
#include <set>
#include <algorithm>
#include <queue>
#include <string>
#include "Cell.h"
#include "MElement.h"
#include "GModel.h"

class Cell;
class BdInfo;

class CellComplex
{
 private:


  std::pair<Cell*, double> _smallestCell;
  std::pair<Cell*, double> _biggestCell;

  GModel* _model;

  // sorted containers of unique cells in this cell complex
  // one for each dimension
  std::set<Cell*, Less_Cell> _cells[4];

  // original cells of this cell complex
  std::set<Cell*, Less_Cell> _ocells[4];

  // original cells removed during reductions
  std::vector<Cell*> _removedcells;

  // cell complex dimension
  int _dim;

  // cell are simplexes
  bool _simplicial;

  // save the original unreduced complex for another reduction run
  bool _saveorig;

  // has a relative subdomain
  bool _relative;

  int _deleteCount;
  int _createCount;

  // is the cell complex at reduced state
  bool _reduced;

  int _numRelativeCells[4];
  int _numSubdomainCells[4];


  // for constructor
  bool _insertCells(std::vector<MElement*>& elements, int domain);
  bool _removeCells(std::vector<MElement*>& elements, int domain);

  bool _immunizeCells(std::vector<MElement*>& elements);

  Cell* _omitCell(Cell* cell, bool dual);

  // enqueue cells in queue if they are not there already
  void enqueueCells(std::map<Cell*, short int, Less_Cell>& cells,
		    std::queue<Cell*>& Q, std::set<Cell*, Less_Cell>& Qset);

  // insert/remove a cell from this cell complex
  void removeCell(Cell* cell, bool other=true, bool del=false);
  void insertCell(Cell* cell);

  // queued coreduction
  int coreduction(Cell* startCell, int omit,
		  std::vector<Cell*>& omittedCells);

  static double _patience;

 public:
  CellComplex(GModel* model,
	      std::vector<MElement*>& domainElements,
	      std::vector<MElement*>& subdomainElements,
              std::vector<MElement*>& nondomainElements,
              std::vector<MElement*>& nonsubdomainElements,
              std::vector<MElement*>& immuneElements,
              bool saveOriginalComplex=true);
  ~CellComplex();


  GModel* getModel() const { return _model; }
  int getDim() const { return _dim; }
  bool simplicial() const { return _simplicial; }
  bool relative() const { return _relative; }

  // get the number of certain dimensional cells
  // if dim = -1 return the number of all cells
  int getSize(int dim, bool orig=false);

  // get domain of a cell
  // cell in domain relative to subdomain  -> domain = 0
  // cell in domain                        -> domain = 1
  // cell in subdomain                     -> domain = 2
  int getDomain(Cell* cell, std::string& str);

  // get dim-dimensional cells
  // domain = 0: cells in domain relative to subdomain
  // domain = 1: cells in domain
  // domain = 2: cells in subdomain
  void getCells(std::set<Cell*, Less_Cell>& cells, int dim, int domain=0);
  int getNumCells(int dim, int domain=0);
  Cell* getACell(int dim, int domain=0);
  //std::set<Cell*, Less_Cell> getOrigCells(int dim){ return _ocells[dim]; }

  // iterator for the cells of same dimension
  typedef std::set<Cell*, Less_Cell>::iterator citer;

  // iterators to the first and last cells of certain dimension
  citer firstCell(int dim, bool orig=false) {
    return orig ? _ocells[dim].begin() : _cells[dim].begin(); }
  citer lastCell(int dim, bool orig=false) {
    return orig ? _ocells[dim].end() : _cells[dim].end(); }

  // true if cell complex has given cell
  bool hasCell(Cell* cell, bool orig=false);

  // check whether two cells both belong to subdomain or if neither one does
  bool inSameDomain(Cell* c1, Cell* c2) const {
    return (c1->getDomain() ==  c2->getDomain()); }

  // remove cells in subdomain from this cell complex
  void removeSubdomain();

  // remove dim-dimensional cells from this cell complex
  void removeCells(int dim);

  // (co)reduction of this cell complex
  // removes (co)reduction pairs of cell of dimension dim and dim-1
  int reduction(int dim, int omit, std::vector<Cell*>& omittedCells);
  int coreduction(int dim, int omit, std::vector<Cell*>& omittedCells);

  // Cell combining for reduction and coreduction
  int combine(int dim);
  int cocombine(int dim);

  // check whether all boundary cells of a cell has this cell
  // as coboundary cell and vice versa
  // also check whether all (co)boundary cells of a cell
  // belong to this cell complex
  bool coherent();

  // full (co)reduction of this cell complex (all dimensions)
  // (combine = 1 -> with combining)
  // (omit = true -> with highest dimensional cell omitting?)
  // (homseq = true -> homology sequence splitting possible after reduction)
  // (heuristic =  0 -> no heuristic, let mesh indexing determine
  //               1 -> omit 0-cell in biggest element
  //              -1 -> omit 0-cell in smallest element)
  int reduceComplex(int combine=1, bool omit=true, bool homseq=false);
  int coreduceComplex(int combine=1, bool omit=true, int heuristic=0);

  // reduce cell complex for Betti number computation
  void bettiReduceComplex();

  bool isReduced() const { return _reduced; }

  int eulerCharacteristic() {
    return getSize(0) - getSize(1) + getSize(2) - getSize(3); }
  void printEuler() {
    printf("Euler characteristic: %d. \n", eulerCharacteristic()); }

  // restore the cell complex to its original state before (co)reduction
  bool restoreComplex();

  // print the vertices of cells of certain dimension
  void printComplex(int dim);


  // experimental
  int saveComplex(std::string filename);
  int loadComplex(std::string filename);
};

#endif
