#include "DivisionAnalysisUtils.hpp"

#include <Triangulate.hpp>
#include <queue>
#include <MeshBuilder.hpp>

//#ifdef THRUST_BACKEND_CUDA
//#include <DivisionAnalysisCudaExport.hpp>
//#endif

using namespace std;

namespace mgx
{


  // write a 2D array to a (csv) file
  bool writeCsvFile(QString filename, QString header, vector<vector<double> >& data)
  {

    QFile file(filename);
    if(!file.open(QIODevice::WriteOnly))
    {
      //setErrorMessage(QString("File '%1' cannot be opened for writing").arg(filename));
      return false;
    }
    QTextStream out(&file);

    out << header << endl;

    forall(const vector<double> &p, data){
      int s = p.size();
      int c = 0;
      while(s>1){
        out << p[c] << ",";
        s--;
        c++;
      }
      out << p[c] << endl;
    }

  return true;
  }

   // calculates the ratio (=bigger divided by smaller value)
  double calcRatio(double value1, double value2)
  {
    if(value1 > value2)
      if(value2 == 0)
        return -1;
      else
        return value1/value2;
    else
      if(value1 == 0)
        return -1;
      else
        return value2/value1;
  }

  // returns an approximation of the 3D cell size (the average of the distance of every vertex to the cell center)
  double avgCellSizeMGX3D(const cell &c)
  {
    double sizeCell = 0;
    forall(const vertex& v, c->S){
      sizeCell += norm(v->pos - c->pos);
    }
    if(c->S.size()>0) return sizeCell/c->S.size();
    else return 0;
  }

  // returns the area of a polygon (=division plane)
  double divisionPlaneArea(Point3d cellCenter, const std::vector<Point3d>& v)
  {
    double a = 0;
    size_t prev = v.size()-1;
    Point3d cpos = cellCenter;
    Point3d dir = normalized((v[0] - cpos) ^ (v[1] - cpos));
    for(size_t cur = 0 ; cur < v.size() ; ++cur)
    {
      a += fabs(((v[prev] - cpos) ^ (v[cur] - cpos))*dir);
      prev = cur;
    }

    return fabs(a)/2;
  }

  // returns the summed area of multiple polygons (=division planes)
  double divisionPlaneAreaMulti(Point3d cellCenter, const std::vector<std::vector<Point3d> >& vMulti)
  {

    double result = 0.;
    forall(std::vector<Point3d> v, vMulti){
      result += divisionPlaneArea(cellCenter, v);
    }

    return result;
  }

  // returns whether a cell with the specified label exists, also provides the cell and the max Label of all cells
  bool findCell(cellGraph& C, cell& searchedCell, int label, int& maxLabel)
  {
    maxLabel = -1;
    bool exists = false;

    // find the cell with the specified label
    forall(const cell& c, C){
      if(c->label == label){
        searchedCell = c;
        exists = true;
      }
      if(maxLabel < c->label){
        maxLabel = c->label;
      }
    }

    return exists;

  }

 // add the neighbors of c2 to the cell graph of c1
 void calcCellGraphJoinedCells(cellGraph &C, const cell& c1, const cell& c2)
 {

    std::set<cell> neighborsC1;

    forall(const cell &c, C.neighbors(c1)){
      neighborsC1.insert(c);
    }

    forall(const cell &c, C.neighbors(c2)){
      if(neighborsC1.find(c) != neighborsC1.end()){
      } else {
        C.insertEdge(c,c1);
        C.insertEdge(c1,c);
      }
    }

 }
   // copies the cell graph
   void copyCellGraph(CellTissue &orig, CellTissue &copy){

    forall(const cell& c, copy.C){
      copy.C.erase(c);
    }

    forall(const cell& c, orig.C){
      copy.C.insert(c);
      forall(const cell& c2, orig.C.neighbors(c)){
        copy.C.insertEdge(c, c2);
      }
    }

    copy.setMeshType("MGX3D");
  }

 // count the occurances of vertices, value 1: vtx only in one cell, value 2: vtx shared between two cells
 std::map<vertex, int> countVtxs(const cell& c1, const cell& c2)
 {
   std::map<vertex, int> vtxIntMap;

    // count occurences
    forall(const vertex& v, c1->S){
      vtxIntMap[v]++;
    }
    forall(const vertex& v, c2->S){
      vtxIntMap[v]++;
    }

    return vtxIntMap;

 }


 std::vector<Point3d> getSharedVtxs(const cell& c1, const cell& c2)
 {
    std::map<vertex, int> vtxIntMap = countVtxs(c1,c2);

    std::vector<Point3d> sharedPoints;


    forall(const VtxInt& p, vtxIntMap){
      if(p.second == 2){
        sharedPoints.push_back(p.first->pos);
      }
    }

    return sharedPoints;
 }


 std::set<vertex> getVtxsOfSharedTris(const cell& c1, const cell& c2)
 {
  std::set<vertex> vtxs;
//  std::set<triVtx> sharedTris;
//  std::set<triVtx>

  forall(const vertex& v, c1->S){
    forall(const vertex& n, c1->S.neighbors(v)){
      vertex m = c1->S.nextTo(v,n);
      if(!c1->S.uniqueTri(v,n,m)) continue;
      if(!c2->S.uniqueTri(v,n,m)) continue;

      vtxs.insert(v);
      vtxs.insert(n);
      vtxs.insert(m);
    }
  }


  return vtxs;
 }

 std::set<Triangle> getSharedTris(const cell& c1, const cell& c2)
 {
  std::set<Triangle> tris;

  forall(const vertex& v, c1->S){
    forall(const vertex& n, c1->S.neighbors(v)){
      vertex m = c1->S.nextTo(v,n);
      if(!c1->S.uniqueTri(v,n,m)) continue;
      if(!c2->S.uniqueTri(v,n,m)) continue;

      tris.insert(Triangle(v,n,m));
      v->selected = true;
      n->selected = true;
      m->selected = true;
    }
  }


  return tris;
 }

 double calcTriArea(std::set<Triangle>& tris)
{
  double area = 0;

  forall(Triangle t, tris){
    area += triangleArea(t.v[0]->pos, t.v[1]->pos, t.v[2]->pos);
  }


  return area;
}


 void addTriDuplicateVertex(MeshBuilder& mb, const vertex& v, const vertex& n, const vertex& m, std::map<vertex, vertex>& oldNewMap)
 {
  if(oldNewMap.find(v) == oldNewMap.end()){
    if(oldNewMap.find(n) == oldNewMap.end()){
      if(oldNewMap.find(m) == oldNewMap.end()){
        mb.addTri(v,n,m, true);
      } else {
        mb.addTri(v,n,oldNewMap[m],true);
      }
    } else {
      if(oldNewMap.find(m) == oldNewMap.end()){
        mb.addTri(v,oldNewMap[n],m,true);
      } else {
        mb.addTri(v,oldNewMap[n],oldNewMap[m],true);
      }
    }
  } else {
    if(oldNewMap.find(n) == oldNewMap.end()){
      if(oldNewMap.find(m) == oldNewMap.end()){
        mb.addTri(oldNewMap[v],n,m,true);
      } else {
        mb.addTri(oldNewMap[v],n,oldNewMap[m],true);
      }
    } else {
      if(oldNewMap.find(m) == oldNewMap.end()){
        mb.addTri(oldNewMap[v],oldNewMap[n],m,true);
      } else {
        mb.addTri(oldNewMap[v],oldNewMap[n],oldNewMap[m],true);
      }
    }
  }
 }



    // shuffles the points of vtx list
   void shufflePoints(std::vector<Point3i> &triList, std::vector<vertex> &vtxList)
    {

    std::vector<vertex> vtxListNew;

    int nrVtxs = vtxList.size();

    std::vector<int> permVec, invPerm(nrVtxs);

    for(int i=0; i<nrVtxs; ++i){
      permVec.push_back(i);
    }

    // using built-in random generator:
    std::random_shuffle ( permVec.begin(), permVec.end() );

    for(int i = 0; i<nrVtxs; i++){
      // convert i-th point to permVec[i]-th one
      vtxListNew.push_back(vtxList[permVec[i]]);
      invPerm[permVec[i]] = i;
    }

    forall(Point3i &p, triList){
      p.x() = invPerm[p.x()];
      p.y() = invPerm[p.y()];
      p.z() = invPerm[p.z()];

    }

    vtxList = vtxListNew;

    }


   void deleteDuplicatesBorder(std::set<vertex>& vtxsBorder, std::vector<Point3i>& triList, std::vector<vertex>& vtxList, double disThreshold)
    {

    std::vector<Point3i> triListNew;
    std::vector<vertex> vtxListNew;

    std::map<int, int> triMap, joinedVtxs;
    //std::set<int> triSet;

    int vtxNr = vtxList.size();
    //double disThreshold = 0.01;

    // go though all vtxs
    for(int i=0; i<vtxNr ; i++){

      if ( triMap.find(i) == triMap.end() ) { // not found yet, check the point

        triMap[i] = i; // add to local tri map
        //triSet.insert(i);
        joinedVtxs[i]=1; // nr of joined vtxs into vtx nr i
        Point3d currentPos = vtxList[i]->pos; // vtx pos to compare with other vtx
        Point3d newPos = currentPos; // new pos of joined vtx
        for(int j=0; j<vtxNr; j++){
          if(i!=j and
             norm(currentPos - vtxList[j]->pos) <= disThreshold
             and vtxsBorder.find(vtxList[j])!=vtxsBorder.end()
             and vtxsBorder.find(vtxList[i])!=vtxsBorder.end()
             and triMap.find(j)==triMap.end()){
              triMap[j] = i;
              //triSet.insert(j);
              joinedVtxs[i]++;
              newPos += vtxList[j]->pos;
              //cout << "replace " << j << " with " << i << endl;
          }


        }
        vtxList[i]->pos = newPos / joinedVtxs[i];
      } else {
      // found
      // the point will already be replaced, do nothing
      }

    }

    typedef std::pair<int,int> IntInt;
    std::map<int, int> triMapNew;

    std::map<int, int> deleteMap;
    int shift = 0;
    int index = 0;
    for(int i=0; i<vtxNr ; i++){
      if(triMap[i] == i){ //
        vtxListNew.push_back(vtxList[i]);
        triMapNew[i] = index;
        //cout << " idx " << i << " / " << i-index;
        index++;

      }
    }

    std::map<int, int> deleteTri;

    forall(Point3i &p, triList){

      p.x() = triMapNew[triMap[p.x()]];
      p.y() = triMapNew[triMap[p.y()]];
      p.z() = triMapNew[triMap[p.z()]];

    }

    vtxList = vtxListNew;

    }

          // this will get rid of duplicated vertices
          // NOTE: this function works on a threshold base which can lead to errors
  void assembleCellFromWalls(vvGraph& S, std::vector<vertex> vtxs, std::vector<Point3i> tris,
    std::set<vertex> &borderVtxs, double mergingThreshold, int labelCounter, bool orientCheck)
  {

        if(vtxs.size() == 0 or tris.size() == 0) return;


        std::vector<Point3i> tempTris = tris;
        std::vector<vertex> tempVtxs = vtxs;

          //cout << "blubb " << tempTris.size() << "/" << tempVtxs.size() << endl;

          deleteDuplicatesBorder(borderVtxs, tempTris, tempVtxs, mergingThreshold);
          if(orientCheck) triangleOrientationCheck(tempTris);
          //triangleOrientationCheck(tempTris);
  //cout << "blubb " << tempTris.size() << "/" << tempVtxs.size() << endl;
          bool checkAgain = true;
          int attempts = 0;
          //if(testCells) checkAgain = true;
          double mergingThresholdTemp = mergingThreshold;
          while(checkAgain){
            std::set<vertex> vtxsBorder;
            attempts++;
            checkAgain = false;
            vvGraph Temp;
            Mesh::meshFromTriangles(Temp, tempVtxs, tempTris);
            markMargins(Temp, Temp, false);
            forall(const vertex &v, Temp){
              if(v->label == -1){
                checkAgain = true;
                vtxsBorder.insert(v);
              }
            }
            if(mergingThresholdTemp > 0.05 or attempts > 10){ // abort, try a last fix
              checkAgain = false;
              tempTris = tris;
              tempVtxs = vtxs;
              if(orientCheck) triangleOrientationCheck(tempTris);
              shufflePoints(tempTris, tempVtxs);
              //deleteDuplicates(tempTris, tempVtxs, mergingThreshold);
              deleteDuplicatesBorder(borderVtxs, tempTris, tempVtxs, mergingThreshold);
              cout << "last fix " << vtxs[0]->label << endl;//labelCounter << endl;
            }
            if(checkAgain){ // check again, increase merging threshold
              //shufflePoints(tempTris, tempVtxs);
              mergingThresholdTemp*=2;
              deleteDuplicatesBorder(vtxsBorder, tempTris, tempVtxs, mergingThresholdTemp);
              if(attempts == 1) cout << "checked again cell " << vtxs[0]->label << endl;//labelCounter << endl;
            }
          }
          if(!Mesh::meshFromTriangles(S, tempVtxs, tempTris)) {
            cout << "Error: There might be cells with errors. Check border vertices and consider rerunning the process. Problem with label " << labelCounter << endl;
        }


  }


 void calcJoinedCell(const cell& c1, const cell& c2, vvGraph& result, bool simple)
 {

    MeshBuilder mb;

    std::map<vertex, int> vtxIntMap = countVtxs(c1,c2);
    std::map<vertex, vertex> oldNewMap;

    // get all triangles that are not shared
    std::map<Point3d, int> pMap;
    std::vector<vertex> vtxs;
    std::vector<Point3i> tris;

    std::set<vertex> sharedVtxsAll, sharedVtxsTris;

    sharedVtxsTris = getVtxsOfSharedTris(c1,c2);

    typedef std::pair<vertex, int> VtxIntP;

    forall(VtxIntP p, vtxIntMap){
      if(p.second == 1) continue;
      if(p.second == 2){
        if(sharedVtxsTris.find(p.first) != sharedVtxsTris.end()) continue;
        // create new vertex for c2
        vertex n;
        n->pos = p.first->pos;
        oldNewMap[p.first] = n;
      }
    }

    forall(const vertex& v, c1->S){
      forall(const vertex& n, c1->S.neighbors(v)){
        vertex m = c1->S.nextTo(v,n);
        if(c1->S.uniqueTri(v,n,m) and (vtxIntMap[v]==1 or vtxIntMap[n]==1 or vtxIntMap[m]==1)){
          if(simple) mb.addTri(v,n,m);//, vtxsP, tris); //t.addTriWithLabel(vtxPMap[v], vtxPMap[n], vtxPMap[m], vtxsP, tris, Point3i(-1, -1, -1));
          else mb.addTri(v,n,m,true);//,vtxs2,tris);
        }
      }
    }
std::cout << "size " << oldNewMap.size() << std::endl;
    forall(const vertex& v, c2->S){
      forall(const vertex& n, c2->S.neighbors(v)){
        vertex m = c2->S.nextTo(v,n);
        if(c2->S.uniqueTri(v,n,m) and (vtxIntMap[v]==1 or vtxIntMap[n]==1 or vtxIntMap[m]==1)){
          if(simple) mb.addTri(v,n,m);//, vtxsP, tris); //t.addTriWithLabel(vtxPMap[v], vtxPMap[n], vtxPMap[m], vtxsP, tris, Point3i(-1, -1, -1));
          else {

          if(oldNewMap.find(v) == oldNewMap.end() and oldNewMap.find(n) == oldNewMap.end() and oldNewMap.find(m) == oldNewMap.end())
            mb.addTri(v,n,m,true);
          else
            addTriDuplicateVertex(mb, v, n, m, oldNewMap);
          }
        }
      }
    }

    if(!simple){ // advanced method: use assembleCellFromWalls to merge border vertices
      std::cout << "assemble" << std::endl;
      assembleCellFromWalls(result, mb.vtxVec, mb.triVec, sharedVtxsTris, 1E-20, c1->label, true);
      //forall(const vertex& v, result) v->label = 0;
      std::cout << "assembled " << oldNewMap.size() << std::endl;
    } else {
    // re-assemble new cell
    //if(!Mesh::meshFromTriangles(result, vtxs, tris)) {
      if(!Mesh::meshFromTriangles(result, mb.vtxVec, mb.triVec)) {
        cout << "problem joining cells " << c1->label << " and " << c2->label << endl;
        //throw(QString("Error, cannot add all the triangles"));
      }
    }



 }


void addCenterDisplacement(int centerVarSteps, double centerVarStepSize, std::vector<Point3d>& normals, std::vector<P3dDouble>& samplingDirections)
  {

    std::vector<P3dDouble> temp;

    forall(const Point3d &p, normals){
      for(int i=-centerVarSteps; i<=centerVarSteps; i++){
        temp.push_back(std::make_pair(p, (double)i*centerVarStepSize));
        //std::cout << "inside displ " << i << "/" << centerVarStepSize << "/" << (double)i*centerVarStepSize << std::endl;
      }
    }

    samplingDirections = temp;

  }


void addCenterDisplacementTestCut(int centerVarSteps, double centerVarStepSize, Point3d cellCenter, Point3d perpendicularPlane,
  std::set<Point3d> perpPlaneNormals2,
  std::unordered_map<Point3d, Point3d>& wallCenter, std::unordered_map<Point3d, std::vector<Point3d> > wallPoints,
  std::vector<P3dDouble>& samplingDirections)
  {

    //std::cout << "TestCut " << perpPlaneNormals2.size() << std::endl;

    std::vector<P3dDouble> temp;

    forall(Point3d p, perpPlaneNormals2){

      Point3d normal = p ^ perpendicularPlane;

      double sP = 0.;
      Point3d intersectP(0.,0.,0.);
      planeLineIntersect(cellCenter, normal, wallCenter[p], wallCenter[p]-normal, sP, intersectP);

      double dis = norm(wallCenter[p] - intersectP);
      if((wallCenter[p]-intersectP) * p < 0) dis *= -1;

      //std::cout << "ins " << p << "/" << wallCenter[p] << "//" << intersectP << std::endl;


      for(int i=-centerVarSteps; i<=centerVarSteps; i++){
        if(planePolygonIntersect(normal, wallCenter[p]+normal*(double)i*centerVarStepSize, wallPoints[p]))
          temp.push_back(std::make_pair(normal, (double)i*centerVarStepSize + dis));
      }
    }

    samplingDirections = temp;
    //std::cout << "TestCutend " << samplingDirections.size() << std::endl;
  }

  // with restrictions
  void generateSamplingDirections(double resolution, int centerVarSteps, double centerVarStepSize,
    Point3d perpendicularPlane, std::set<Point3d> perpPlaneNormals2, bool enforceDivThroughNeighbor,
    std::vector<P3dDouble>& samplingDirections)
  {

    std::vector<P3dDouble> temp;
    std::map<Point3d, bool > tempMap;

    samplingDirections = temp;

    std::vector<Point3d> normals;
    std::set<Point3d> orth;

    std::map<Point3d,Point3d> perpPlaneDivPlaneMap;

    std::cout << "gen " << perpendicularPlane << "/" << perpPlaneNormals2.size() << std::endl;

    if(perpendicularPlane == Point3d(0,0,0) and perpPlaneNormals2.empty()){
      // invalid input
      return;
    } else if(perpendicularPlane == Point3d(0,0,0)){
      // only perpPlaneNormals2
      forall(Point3d p, perpPlaneNormals2){
        std::cout << "perps " << p << std::endl;
        Point3d divP = orthogonalVector(p);
        orth.insert(divP);
        std::cout << "perp2s " << divP << std::endl;
        perpPlaneDivPlaneMap[p] = divP;
      }
    } else if(perpPlaneNormals2.empty()){
      // only perpendicularPlane
      Point3d divP = orthogonalVector(perpendicularPlane);
      orth.insert(divP);
      perpPlaneDivPlaneMap[divP] = perpendicularPlane;
    } else {
      // both
      if(enforceDivThroughNeighbor){
        // do that outside
        forall(Point3d p, normals){
          samplingDirections.push_back(std::make_pair(p,0.));
        }
      } else {
        forall(Point3d p, perpPlaneNormals2){
          normals.push_back(p ^ perpendicularPlane);
        }
        addCenterDisplacement(centerVarSteps, centerVarStepSize, normals, samplingDirections);
      }
      return;
    }

    // define possible division plane normals
    Point3d xDir (1,0,0);
    Point3d yDir (0,1,0);
    Point3d zDir (0,0,1);
    forall(Point3d pOrth, orth){
      std::cout << "orth! " << pOrth << "/" << perpPlaneDivPlaneMap[pOrth] << std::endl;
      for(double rx=0; rx<180; rx+=resolution){
        double rotAngleX=rx*(M_PI/180.0);
        Matrix3d rotX = Matrix3d::rotation(perpPlaneDivPlaneMap[pOrth], rotAngleX);

        Point3d p = rotX*pOrth;
        if(std::abs(p.x()) < 0.0001) p.x() = 0;
        if(std::abs(p.y()) < 0.0001) p.y() = 0;
        if(std::abs(p.z()) < 0.0001) p.z() = 0;
        if(norm(p) > 0.001) tempMap[p] = true;
      }
    }

    forall(const P3dboolPair &p, tempMap){
      int s = normals.size();
      if(s==0){
        normals.push_back(p.first);
      } else {
        if(norm(normals[s-1]-p.first)>0.001){
          normals.push_back(p.first);
        }
      }
    }
    std::cout << "normals: " << normals.size() << std::endl;
    //if(enforceDivThroughNeighbor and !perpPlaneNormals2.empty()){
      // TODO
    //} else {
      addCenterDisplacement(centerVarSteps, centerVarStepSize, normals, samplingDirections);
    //}

  }

  std::vector<Point3d> generatePointsOnSphere(int nrPoints)
  {
    std::vector<Point3d> res;

    nrPoints *= 2;

    double phi = M_PI * (3. - std::sqrt(5.));  // golden angle in radians

    for(int i = 0; i < nrPoints; i++){
      double y = 1. - ((double)(i) / (double)(nrPoints - 1)) * 2.;  // y goes from 1 to -1
      double radius = std::sqrt(1. - y * y);  // radius at y
      double theta = phi * i;  // golden angle increment

      double x = std::cos(theta) * radius;
      double z = std::sin(theta) * radius;

      if(y >= 0)
        res.push_back(Point3d(x,y,z));
    }

    return res;
  }

  // new code, no restrictions on planes
  void generateSamplingDirectionsNew(int nrPlanes, int centerVarSteps, double centerVarStepSize, std::vector<P3dDouble>& samplingDirections)
  {
    samplingDirections.clear();

    int nrPlanesPerPoint = nrPlanes / (2*centerVarSteps+1);

    std::vector<Point3d> dirs = generatePointsOnSphere(nrPlanesPerPoint);

    addCenterDisplacement(centerVarSteps, centerVarStepSize, dirs, samplingDirections);
  }

  // with center displacement but without further restrictions
  void generateSamplingDirections(double resolution, int centerVarSteps, double centerVarStepSize, std::vector<P3dDouble>& samplingDirections)
  {
    std::vector<P3dDouble> temp;
    std::map<Point3d, bool > tempMap;
    samplingDirections = temp;

    std::vector<Point3d> normals;

    // define possible division plane normals
    Point3d xDir (1,0,0);
    Point3d yDir (0,1,0);
    Point3d zDir (0,0,1);
    for(int rx=0; rx<180; rx+=resolution){
      double rotAngleX=rx*(M_PI/180.0);
      Matrix3d rotX = Matrix3d::rotation(xDir, rotAngleX);
      for(int ry=0; ry<180; ry+=resolution){
        double rotAngleY=ry*(M_PI/180.0);
        Matrix3d rotY = Matrix3d::rotation(yDir, rotAngleY);
        for(int rz=0; rz<180; rz+=resolution){
          double rotAngleZ=rz*(M_PI/180.0);
          Matrix3d rotZ = Matrix3d::rotation(zDir, rotAngleZ);
          Point3d p = rotX*rotY*rotZ*xDir;
          if(abs(p.x()) < 0.0001) p.x() = 0;
          if(abs(p.y()) < 0.0001) p.y() = 0;
          if(abs(p.z()) < 0.0001) p.z() = 0;
          tempMap[p] = true;
          //bool add = true;
          //forall(Point3d d, samplingDirections){
          //  if(norm(d-p)<0.001) add=false;
          //}
          //if(add) samplingDirections.push_back(p);
        }
      }
    }

    forall(const P3dboolPair &p, tempMap){
      int s = normals.size();
      if(s==0){
        normals.push_back(p.first);
      } else {
        if(norm(normals[s-1]-p.first)>0.001){
          normals.push_back(p.first);
        }
      }
    }

    addCenterDisplacement(centerVarSteps, centerVarStepSize, normals, samplingDirections);

  }


  double calcVolume(std::vector<Triangle> tris)
  {
  double vol = 0;
  forall(Triangle &t, tris){
    vol += signedTetraVolume(t.v[0]->pos, t.v[1]->pos, t.v[2]->pos);
  }
  return vol;

  }

  bool planeSegmentIntersect(const Point3d &p, const Point3d &nrml, const Point3d &u1, const Point3d &u2, Point3d &intersectP)
  {
    double sP = 0.;
    planeLineIntersect(p, nrml, u1, u2, sP, intersectP);

    if(sP>0 and sP<1){
      return true;
    } else {
      return false;
    }
  }

  void findNeighbTris(cellGraph& C, const cell& c, std::map<cell, std::vector<Triangle> >& neighbTris)
  {

    forall(const vertex& v, c->S){
      forall(const vertex& n, c->S.neighbors(v)){
        vertex m = c->S.nextTo(v,n); // prevto ?
        if(c->S.uniqueTri(v,n,m)){
          forall(const cell& c2, C.neighbors(c)){
            if(c2->S.uniqueTri(v,n,m)){
              Triangle t(v,n,m);
              // t[0] = v;
              // t[1] = n;
              // t[2] = m;
              neighbTris[c2].push_back(t);
            }
          }
        }
      }
    }
  }

  void cellTriNeighbors(cellGraph& C, const cell& c, std::map<Triangle, int>& trisNeighbMap)
  {

    //std::map<tri, int> triNeighbMap;

    // first get all tris of c
    forall(const vertex& v, c->S){
      forall(const vertex& n, c->S.neighbors(v)){
        vertex m = c->S.nextTo(v,n);
        if(c->S.uniqueTri(v,n,m)){
          Triangle t(v,n,m);
          // t.v[0] = v;
          // t[1] = n;
          // t[2] = m;
          trisNeighbMap[t] = -1;
        }
      }
    }
    // now go through neighbors and overwrite label
    forall(const cell& c2, C.neighbors(c)){

      forall(const vertex& v, c2->S){
        forall(const vertex& n, c2->S.neighbors(v)){
          vertex m = c2->S.prevTo(v,n);
          if(c->S.uniqueTri(v,n,m)){
            Triangle t(v,n,m);
            // tri t;
            // t[0] = v;
            // t[1] = n;
            // t[2] = m;
            trisNeighbMap[t] = c2->label;
          }
        }
      }

    }



  }

  // estimates the volume of a 3D body with 1 open side, open Edges need to be oriented right
  double estimateVolumeOpenBody(std::vector<Triangle> tris, set<vertex>& exposedVtxs)
  {

  double vol = 0;
  Point3d avgOpenPos (0,0,0);
  forall(const vertex& v, exposedVtxs){
    avgOpenPos += v->pos;
  }
  avgOpenPos/=exposedVtxs.size();

  forall(Triangle &t, tris){
    Point3d p1 = t.v[0]->pos;
    Point3d p2 = t.v[1]->pos;
    Point3d p3 = t.v[2]->pos;
    if(exposedVtxs.find(t.v[0]) != exposedVtxs.end()) p1 = avgOpenPos;
    if(exposedVtxs.find(t.v[1]) != exposedVtxs.end()) p2 = avgOpenPos;
    if(exposedVtxs.find(t.v[2]) != exposedVtxs.end()) p3 = avgOpenPos;
    vol += signedTetraVolume(p1, p2, p3);
  }

  return vol;
  }


  void estimateVolumeDaughterCells(std::vector<Triangle>& trisLeft, std::vector<Triangle>& trisRight, std::vector<Triangle>& trisSplit)
  {


  std::vector<std::pair<vertex,vertex> > openEdgesLeft, openEdgesRight;

  set<vertex> exposedLeft, exposedRight;

  int labelLeft = trisLeft[0].v[0]->label;
  //int labelRight = trisRight[0][0]->label;

  forall(const Triangle &t, trisSplit){
    std::pair<vertex,vertex> vv;

    if(t.v[0]->label == labelLeft){
      exposedLeft.insert(t.v[0]);
    } else {
      exposedRight.insert(t.v[0]);
    }
    if(t.v[1]->label == labelLeft){
      exposedLeft.insert(t.v[1]);
    } else {
      exposedRight.insert(t.v[1]);
    }
    if(t.v[2]->label == labelLeft){
      exposedLeft.insert(t.v[2]);
    } else {
      exposedRight.insert(t.v[2]);
    }
/*
    if(t[0]->label == labelLeft and t[1]->label == labelLeft){
      vv = make_pair(t[0],t[1]);
      openEdgesLeft.push_back(vv);
    } else if(t[1]->label == labelLeft and t[2]->label == labelLeft){
      vv = make_pair(t[1],t[2]);
      openEdgesLeft.push_back(vv);
    } else if(t[2]->label == labelLeft and t[0]->label == labelLeft){
      vv = make_pair(t[2],t[0]);
      openEdgesLeft.push_back(vv);

    } else if(t[0]->label == labelRight and t[1]->label == labelRight){
      vv = make_pair(t[0],t[1]);
      openEdgesRight.push_back(vv);
    } else if(t[1]->label == labelRight and t[2]->label == labelRight){
      vv = make_pair(t[1],t[2]);
      openEdgesRight.push_back(vv);
    } else if(t[2]->label == labelRight and t[0]->label == labelRight){
      vv = make_pair(t[2],t[0]);
      openEdgesRight.push_back(vv);
    }*/

  }

  double volLeft = estimateVolumeOpenBody(trisLeft, exposedLeft);
  double volRight = estimateVolumeOpenBody(trisRight, exposedRight);

  cout << "vol " << volLeft << "/" << volRight << endl;

  }



  // function to move vertices that have identical positions a tiny bit apart
  void separateIdenticalVtxs(const cell &c){

    std::map<Point3d, std::vector<vertex> > pointVerticesMap;
    typedef std::pair<Point3d, std::vector<vertex> > P3dVtxVecP;
    forall(const vertex& v, c->S){
      pointVerticesMap[v->pos].push_back(v);
    }

    forall(const P3dVtxVecP& p, pointVerticesMap){
      if(p.second.size() > 1){
        for(size_t i=0; i< p.second.size(); i++){
          Point3d nrml;
          calcNormal(c->S, p.second[i], nrml);
          cout << "found!  " << p.second[i]->pos << "/" << nrml << endl;
          p.second[i]->pos -= nrml/1000.;
        }
      }
    }

  }

  // test if all points of the polygon are above or below the plane, return true if not (plane cuts through the polygon points)
  bool planePolygonIntersect(Point3d planeNrml, Point3d planePos, std::vector<Point3d> polygonVtxs)
  {

    int count = 0;
    bool above = true;
    double eps = 1E-4;

    forall(Point3d p, polygonVtxs){
      Point3d vecPlanePoint = p - planePos;
      double scalar = vecPlanePoint * planeNrml;
      if(count == 0){ // set above
        if(scalar > eps){
          above = true;
          count++;
        } else if(scalar < -eps){
          above = false;
          count++;
        }
      } else { // test above
        if(scalar > eps and !above){
          return true;
        } else if(scalar < -eps and above){
          return true;
        }
      }
    }
    return false;

  }


  // subroutine for randwalkbetw: return a value from the C matrix (first row and col are zeros, rest identical to Cs)
  double valueC(std::vector<std::vector<double> >& Cs, int i, int j)
  {
    if(i == 0 or j == 0)
      return 0.;
    else
      return Cs[i-1][j-1];
  }

  void calcRandWalkBetweenness(std::map<IntInt, double>& neighborhoodWeights, std::map<int, int>& labelPosMap, std::vector<int>& orderedLabels,
                               std::map<int, double>& cellBetweenessMap, bool useGPU)
  {

    //std::cout << "rand here " << neighborhoodWeights.size() << "/" << labelPosMap.size() << "/" << orderedLabels.size() << std::endl;

    uint numNodes = orderedLabels.size();

    bool useEqualWeights = false;

    std::vector<std::vector<double> > L, Cs;

    std::vector<std::pair<int, int> > edges;

    //std::cout << "create L" << std::endl;

    std::map<int, std::set<int> > labelNeighborMap;

    forall(auto p, neighborhoodWeights){
      if(p.second < 1E-5) continue;
      labelNeighborMap[p.first.first].insert(p.first.second);
      labelNeighborMap[p.first.second].insert(p.first.first);
    }

    //td::cout << "fill mat " << std::endl;

    // create Laplacian
    for(uint i = 0;i<numNodes-1; i++){
      std::vector<double> newRow(numNodes-1, 0.), newRow2(numNodes-1, 0.);
      L.push_back(newRow);
      Cs.push_back(newRow2);
    }


    for(uint i = 0;i<numNodes; i++){
      //std::cout << "for " << i << "/" << numNodes << "/" << orderedLabels[i] << "/" << labelNeighborMap[orderedLabels[i]].size() << std::endl;
      int c = orderedLabels[i];
      cellBetweenessMap[c] = 0;


      //std::cout << "for2 " << std::endl;
      double nCount = 0;
      forall(int n, labelNeighborMap[c]){
         if(c < n) edges.push_back(make_pair(c,n)); // only take each edge once
         IntInt edge = make_pair(c, n);
         //std::cout << "edge " << c << "/" << n << "/"<< labelPosMap[c] << "/" << labelPosMap[n] << std::endl;
         /*if(!useEqualWeights)*/ nCount+=neighborhoodWeights[edge];
         //else nCount++;

         if(labelPosMap[c] != 0 and labelPosMap[n] != 0){
          //L[labelPosMap[c]-1][labelPosMap[n]-1] = -1.;
          /*if(!useEqualWeights)*/ L[labelPosMap[c]-1][labelPosMap[n]-1] = -1.*neighborhoodWeights[edge];
        }
      }
      if(i>0) L[i-1][i-1] = nCount;
    }

    std::cout << "mat " << L.size() << "/" << Cs.size() << std::endl;

    // generate C matrix (invert laplacian)
    if(useGPU){
//#ifdef THRUST_BACKEND_CUDA
//      matrixInverseGPU(L,Cs);
//#else
      matrixInverseGSL(L,Cs);
//#endif
    } else {
      matrixInverseGSL(L,Cs);
    }

    std::cout << "mat solved " << L.size() << "/" << Cs.size() << std::endl;

    #pragma omp parallel for schedule(guided)
    for(std::vector<IntInt>::size_type i=0; i<edges.size(); i++){ // forall edges
      IntInt p = edges[i];
      // calc Fe
      std::vector<std::pair<double, int> > rows (numNodes);
      for(uint j = 0;j<numNodes; j++){
        double row = valueC(Cs, labelPosMap[p.first], j) - valueC(Cs, labelPosMap[p.second], j);
        /*if(!useEqualWeights)*/ row*=neighborhoodWeights[make_pair(p.first, p.second)];
        rows[j] = make_pair(row, j);
      }
      // sort Fe
      std::sort(rows.begin(), rows.end());

      for(uint j=0;j<numNodes; j++){
        double fevi = rows[j].first;
        double posevi = numNodes-1-j;
        double inOrder = rows[j].second;
        #pragma omp critical
          cellBetweenessMap[p.first]+=(inOrder-posevi)*fevi;
        #pragma omp critical
          cellBetweenessMap[p.second]+=(numNodes-inOrder-1-posevi)*fevi;
      }
    }
    double nb = (numNodes - 1) * (numNodes - 2) / 2.;
    for(uint i=0;i<numNodes; i++){
      cellBetweenessMap[orderedLabels[i]] = (cellBetweenessMap[orderedLabels[i]]-i)/nb;
    }
  }


  void calcRandWalkBetweenness(cellGraph &C, std::map<IntInt, double>& weights, std::map<cell, int>& cellPosMap, std::vector<cell>& orderedCells,
                               std::map<cell, double>& cellBetweenessMap, bool useGPU)
  {

    uint numNodes = orderedCells.size();

    bool useEqualWeights = weights.size()==0 ? true : false;

    std::vector<std::vector<double> > L, Cs;

    std::vector<std::pair<cell, cell> > edges;
    typedef std::pair<cell, cell> CellCell;

    //std::cout << "create L" << std::endl;

    // create Laplacian
    for(uint i = 0;i<numNodes; i++){
      cell c = orderedCells[i];
      cellBetweenessMap[c] = 0;
      std::vector<double> newRow(numNodes-1, 0.), newRow2(numNodes-1, 0.);

      if(i>0){
        L.push_back(newRow);
        Cs.push_back(newRow2);
      }

      double nCount = 0;
      forall(const cell& n, C.neighbors(c)){
        if(c->label < n->label) edges.push_back(make_pair(c,n)); // only take each edge once
        IntInt edge = make_pair(c->label, n->label);
        if(!useEqualWeights) nCount+=weights[edge];
        else nCount++;
        if(cellPosMap[c] != 0 and cellPosMap[n] != 0){
          L[cellPosMap[c]-1][cellPosMap[n]-1] = -1.;
          if(!useEqualWeights) L[cellPosMap[c]-1][cellPosMap[n]-1] = -1.*weights[edge];
        }
      }
      if(i>0) L[i-1][i-1] = nCount;
    }

    // generate C matrix (invert laplacian)
    if(useGPU){
//#ifdef THRUST_BACKEND_CUDA
//      matrixInverseGPU(L,Cs);
//#else
      matrixInverseGSL(L,Cs);
//#endif
    } else {
      matrixInverseGSL(L,Cs);
    }

    #pragma omp parallel for schedule(guided)
    for(std::vector<CellCell>::size_type i=0; i<edges.size(); i++){ // forall edges
      CellCell p = edges[i];
      // calc Fe
      std::vector<std::pair<double, int> > rows (numNodes);
      for(uint j = 0;j<numNodes; j++){
        double row = valueC(Cs, cellPosMap[p.first], j) - valueC(Cs, cellPosMap[p.second], j);
        if(!useEqualWeights) row*=weights[make_pair(p.first->label, p.second->label)];
        rows[j] = make_pair(row, j);
      }
      // sort Fe
      std::sort(rows.begin(), rows.end());

      for(uint j=0;j<numNodes; j++){
        double fevi = rows[j].first;
        double posevi = numNodes-1-j;
        double inOrder = rows[j].second;
        #pragma omp critical
          cellBetweenessMap[p.first]+=(inOrder-posevi)*fevi;
        #pragma omp critical
          cellBetweenessMap[p.second]+=(numNodes-inOrder-1-posevi)*fevi;
      }
    }
    double nb = (numNodes - 1) * (numNodes - 2) / 2.;
    for(uint i=0;i<numNodes; i++){
      cellBetweenessMap[orderedCells[i]] = (cellBetweenessMap[orderedCells[i]]-i)/nb;
    }
  }

  void calcRandWalkBetweenness(cellGraph &C, std::map<IntInt, double>& weights, std::map<cell, double>& cellBetweenessMap, bool useGPU, bool doRevCuthill)
  {
    std::map<cell, int> cellPosMap;
    std::vector<cell> orderedCells;

    if(doRevCuthill){
      orderedCells = cuthillMckeeOrdering(C, cellPosMap);
      reverseVector(orderedCells);
    }

    int cc = 0;
    forall(const cell& c, C){
      cellPosMap[c] = cc;
      if(doRevCuthill){
        cellPosMap[c] = orderedCells.size() - cellPosMap[c] - 1;
      }
      cc++;
      orderedCells.push_back(c);
    }
    calcRandWalkBetweenness(C, weights, cellPosMap, orderedCells, cellBetweenessMap, useGPU);
  }

  void calcRandWalkBetweenness(std::map<IntInt, double>& neighborhoodGraph, std::map<int, double>& labelBetweenessMap, bool useGPU)
  {
    std::map<int, int> labelPosMap;
    std::vector<int> orderedCells;
    std::set<int> allLabels;

    forall(auto p, neighborhoodGraph){
      if(p.second < 1E-5) continue;
      allLabels.insert(p.first.first);
      allLabels.insert(p.first.second);
    }

    forall(int l, allLabels){
      labelPosMap[l] = orderedCells.size();
      orderedCells.push_back(l);
    }

    //std::cout << "rand now" << std::endl;

    calcRandWalkBetweenness(neighborhoodGraph, labelPosMap, orderedCells, labelBetweenessMap, useGPU);
  }

  std::vector<Triangle> getTrisOfCells(Mesh *mesh, int selectedLabel, bool parentMode)
  {
    vvGraph& S = mesh->graph();
    std::vector<Triangle> cellTris;
    forall(const vertex& v, S){
      int cellLabel = v->label;
      if(parentMode) cellLabel = mesh->parents()[cellLabel];
      if(cellLabel != selectedLabel) continue;
      forall(const vertex& n, S.neighbors(v)){
        vertex m = S.nextTo(v,n);
        if(!S.uniqueTri(v,n,m)) continue;
        int cellLabel = getLabel(v,n,m);
        if(parentMode) cellLabel = mesh->parents()[cellLabel];
        if(cellLabel != selectedLabel) continue;
        Triangle tr1(v,n,m);
        tr1.v[0] = v;
        tr1.v[1] = n;
        tr1.v[2] = m;
        cellTris.push_back(tr1);
      }
    }
    return cellTris;
  }

  std::vector<vertex> getVtxsOfCells(Mesh *m, int selectedLabel, bool parentMode)
  {
    std::vector<vertex> cellVtxs;
    forall(const vertex& v, m->graph()){
      int cellLabel = v->label;
      if(parentMode) cellLabel = m->parents()[v->label];
      if(cellLabel != selectedLabel) continue;
      cellVtxs.push_back(v);
    }
    return cellVtxs;
  }


  std::vector<vertex> getOrRemoveCommonWall(Mesh *m, std::vector<vertex>& cellVtxs, bool removeCommon, bool includeEdgeOfShared)
  {


    std::map<int,double> cellWallArea;
    std::map<int,double> outsideWallArea;

    std::map<Point3i, std::set<Triangle> > triVtxSetMap; // map from tri idx to vtxs of tri

    std::map<int, std::set<vertex> > posVtxMap;
    std::map<vertex, IntInt> vtxWallMap;

    // output maps
    std::map<vertex, int> vtxNewLabelMap;
    std::map<vertex, int> vtxCellCounterMap;
    std::map<IntInt, double> sharedWallArea;

    Point3dIntMap vMap;
    Point3iIntSetMap triCell; // This is set of cells for each triangle

    vvGraph& S = m->graph();
    double tolerance = 0.0001;
    forall(vertex v, cellVtxs){
      forall(const vertex& n, S.neighbors(v)) {
        const vertex& m = S.nextTo(v,n);
        if(!S.uniqueTri(v,n,m)) // Only once per triangle in the mesh
          continue;
        // If the labels are not the same, this means we have more than one label on the same cell
        if(v->label != n->label or v->label != m->label)
          std::cout << "Error, the mesh has volumetric cells with more than one label." << std::endl;
        int cell = v->label;
        //cellWallArea[cell]+=triangleArea(v->pos,n->pos, m->pos);
        // Add the cell to the triangle map, points with similar positions will get the same index
        triCell[triIndex(Point3i(vIndex(v->pos, vMap, tolerance), vIndex(n->pos, vMap, tolerance),
                               vIndex(m->pos, vMap, tolerance)))].insert(cell);

        Triangle tr1(v,n,m);
        tr1.v[0] = v;
        tr1.v[1] = n;
        tr1.v[2] = m;

        triVtxSetMap[triIndex(Point3i(vIndex(v->pos, vMap, tolerance), vIndex(n->pos, vMap, tolerance),
                               vIndex(m->pos, vMap, tolerance)))].insert(tr1);

        posVtxMap[vIndex(v->pos, vMap, tolerance)].insert(v);
        posVtxMap[vIndex(n->pos, vMap, tolerance)].insert(n);
        posVtxMap[vIndex(m->pos, vMap, tolerance)].insert(m);
      }
    }

    std::vector<vertex> vtxsSharedWall, vtxsNonSharedWall;

    forall(auto p, posVtxMap){
      if(p.second.size() == 2){
        vertex v = *p.second.begin();
        vertex v2 = *(++p.second.begin());
        vtxsSharedWall.push_back(v);
        vtxsSharedWall.push_back(v2);
      }
    }

    std::cout << "sss " << vtxsSharedWall.size() << "/" << posVtxMap.size() << std::endl;

    std::set<vertex> sharedWallSet, cellVtxsSet;
    forall(vertex v, vtxsSharedWall){
      sharedWallSet.insert(v);
    }
    forall(vertex v, cellVtxs){
      cellVtxsSet.insert(v);
    }

    if(!includeEdgeOfShared){
      std::vector<vertex> vtxsSharedWallNoEdge;

      forall(vertex v, vtxsSharedWall){
        bool nonWallNeighbor = false;
        forall(vertex n, S.neighbors(v)){
          if(sharedWallSet.find(n) == sharedWallSet.end()){
            nonWallNeighbor = true;
          }
        }
        if(!nonWallNeighbor){
          vtxsSharedWallNoEdge.push_back(v);
        }

      }
      vtxsSharedWall = vtxsSharedWallNoEdge;
    }

    sharedWallSet.clear();
    forall(vertex v, vtxsSharedWall){
      sharedWallSet.insert(v);
    }

    forall(vertex v, cellVtxs){
      if(sharedWallSet.find(v) == sharedWallSet.end()){
        vtxsNonSharedWall.push_back(v);
      }
    }


    if(removeCommon){
      return vtxsNonSharedWall;
    } else {
      return vtxsSharedWall;
    }

  }

    std::vector<Triangle> getOrRemoveCommonWall(Mesh *m, std::vector<Triangle>& cellTris, bool removeCommon, bool mergeDuplicated)
    {


      std::map<int,double> cellWallArea;
      std::map<int,double> outsideWallArea;

      std::map<Point3i, std::set<Triangle> > triVtxSetMap; // map from tri idx to vtxs of tri

      std::map<int, std::set<vertex> > posVtxMap;
      std::map<vertex, IntInt> vtxWallMap;

      // output maps
      std::map<vertex, int> vtxNewLabelMap;
      std::map<vertex, int> vtxCellCounterMap;
      std::map<IntInt, double> sharedWallArea;

      Point3dIntMap vMap;
      Point3iIntSetMap triCell; // This is set of cells for each triangle

      vvGraph& S = m->graph();
      double tolerance = 0.0001;
      forall(Triangle t, cellTris){
        //forall(const vertex& n, S.neighbors(v)) {
          //const vertex& m = S.nextTo(v,n);
          //if(!S.uniqueTri(v,n,m)) // Only once per triangle in the mesh
          //  continue;
          // If the labels are not the same, this means we have more than one label on the same cell
          vertex v = t.v[0];
          vertex n = t.v[1];
          vertex m = t.v[2];

          if(v->label != n->label or v->label != m->label)
            std::cout << "Error, the mesh has volumetric cells with more than one label." << std::endl;
          int cell = v->label;
          //cellWallArea[cell]+=triangleArea(v->pos,n->pos, m->pos);
          // Add the cell to the triangle map, points with similar positions will get the same index
          triCell[triIndex(Point3i(vIndex(v->pos, vMap, tolerance), vIndex(n->pos, vMap, tolerance),
                                 vIndex(m->pos, vMap, tolerance)))].insert(cell);

          triVtxSetMap[triIndex(Point3i(vIndex(v->pos, vMap, tolerance), vIndex(n->pos, vMap, tolerance),
                                 vIndex(m->pos, vMap, tolerance)))].insert(t);

          posVtxMap[vIndex(v->pos, vMap, tolerance)].insert(v);
          posVtxMap[vIndex(n->pos, vMap, tolerance)].insert(n);
          posVtxMap[vIndex(m->pos, vMap, tolerance)].insert(m);
        //}
      }

      std::vector<Triangle> trisSharedWall, trisNonSharedWall;

      std::set<vertex> duplicatedVtxs;

      std::map<vertex,vertex> duplicatedMapping;


            forall(auto p, posVtxMap){
              if(p.second.size() == 2){

                vertex v = *p.second.begin();
                vertex v2 = *(++p.second.begin());
                duplicatedVtxs.insert(v);
                duplicatedVtxs.insert(v2);
                duplicatedMapping[v] = v2;
              }
            }


      forall(Triangle t, cellTris){
        vertex v = t.v[0];
        if(duplicatedMapping.find(v)!=duplicatedMapping.end()) v = duplicatedMapping[v];
        vertex n = t.v[1];
        if(duplicatedMapping.find(n)!=duplicatedMapping.end()) n = duplicatedMapping[n];
        vertex m = t.v[2];
        if(duplicatedMapping.find(m)!=duplicatedMapping.end()) m = duplicatedMapping[m];

        if(duplicatedVtxs.find(v) != duplicatedVtxs.end() and duplicatedVtxs.find(n) != duplicatedVtxs.end() and duplicatedVtxs.find(m) != duplicatedVtxs.end()){
          trisSharedWall.push_back(t);
        } else {
          trisNonSharedWall.push_back(t);
        }


      }



      std::cout << "sss2 " << cellTris.size() << "/" << trisSharedWall.size() << "/" << trisNonSharedWall.size() << std::endl;

      if(removeCommon){
        return trisNonSharedWall;
      } else {
        return trisSharedWall;
      }

    }

}
