/******************************************************************************
 *
 * Extension of AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
 * Authors: Simon Praetorius et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * See also license.opensource.txt in the distribution.
 * 
 ******************************************************************************/


#include "BackgroundMesh.h"
#include "VectorOperations.h"	// getMin, sort, CompairPair
#include "Initfile.h"

#include <boost/math/special_functions/round.hpp>

namespace experimental {
  
  std::map<const FiniteElemSpace*, std::pair<int, Box*> > Box::boxMap;


  Box::Box(int DOW_, std::vector<int> N_)
  : DOW(DOW_), N(N_), boxFilled(false)
  {
    min_corner.set(-1.0);
    max_corner.set(1.0);
    init();
  }


  Box::Box(int DOW_, PointType min_corner_, PointType max_corner_, std::vector<int> N_)
  : DOW(DOW_), N(N_),
    boxFilled(false), 
    min_corner(min_corner_),
    max_corner(max_corner_)
  {
    init();
  }


  void Box::fillBox(const FiniteElemSpace* feSpace)
  {
    DOFVector<WorldVector<double> > coords(feSpace, "coords");
    feSpace->getMesh()->dofCompress();
    feSpace->getMesh()->getDofIndexCoords(coords);

    for (DegreeOfFreedom i = 0; i < coords.getUsedSize(); i++) {
      PointType location(coords[i]);
      addData(location, std::make_pair(i, location));
    }
    boxFilled = true;
  }


  void Box::clearBox()
  {
    boxData.clear();
  }


  void Box::clearBoxData()
  {
    for (size_t i = 0; i < boxData.size(); i++) {
      boxData[i].clear();
    }
  }
  

  int Box::getMaxBoxSize()
  {
    size_t maxSize = 0;
    for (size_t i = 0; i < boxData.size(); i++) {
      maxSize = std::max(maxSize, boxData[i].size());
    }
    return maxSize;
  }


  bool Box::inBox(const PointType& x)
  {
    for (int i = 0; i < DOW; i++)
      if (min_corner[i] > x[i] || max_corner[i] < x[i])
	return false;
    return true;
  }


  /**
   *
   * Ny=2
   * ^----------------
   * | 4 | 5 | 6 | 7 |
   * ----+---+---+----
   * | 0 | 1 | 2 | 3 |
   * -----------------> Nx=4
   **/
  int Box::getBox(const PointType& x)
  {
    std::vector<int> idx(DOW);
    for (int i = 0; i < DOW; i++)
      idx[i] = static_cast<int>((x[i]-min_corner[i])*N[i]/(max_corner[i]-min_corner[i]));
//       idx[i] = static_cast<int>(floor((x[i]-min_corner[i])*N[i]/(max_corner[i]-min_corner[i])));
    return idx2nr(idx);
  }


  bool Box::getNearestData(const PointType& x, DataType& data)
  {
    if (!inBox(x))
      return false;

    int nr = getBox(x);
    unsigned int n = boxData[nr].size();
    // Box, in der x liegt, ist leer
    if (n == 0) {
      std::vector<DataType> datas;
      getNearestData(x, datas, 1);
      data = datas[0];
      return true;
    }
    // suche Punkt in aktueller Box, der am nächsten zu x liegt
    std::vector<double> distances;
    std::vector<DataType>::iterator dataIter;
    for (dataIter = boxData[nr].begin(); dataIter != boxData[nr].end(); dataIter++) {
      distances.push_back(distance2(x, (*dataIter).second, DOW));
    }
    double minDist = 1.e15;
    size_t minIdx = 0;
    vector_operations::getMin(distances, minDist, minIdx);
    data = boxData[nr][minIdx];
    return true;
  }


  bool Box::getNearestData(const PointType& x, std::vector<DataType>& data, int nData)
  {
    if (!inBox(x))
      return false;
    int center_box = getBox(x);

    std::set<int> nrs; nrs.insert(center_box);
    double boxBoundaryDist = getBoxBoundaryDist(center_box, x);
    std::vector<double> distances;
    std::set<int>::iterator nrIter;
    for (int level = 0; static_cast<int>(data.size()) < nData && nrs.size() > 0; level++) {
      int oldDataSize = data.size();
      for (nrIter = nrs.begin(); nrIter != nrs.end(); nrIter++) {
	data.insert(data.end(), boxData[*nrIter].begin(), boxData[*nrIter].end());
      }
      int addon = data.size()-oldDataSize;
      if (addon > 0) {
	// 1.) bestimme Abstand zu x
	int n = distances.size();
	TEST_EXIT(n + addon ==static_cast<int>(data.size()) && oldDataSize == n)("hier ist ein Index falsch: distances.size = %d, addon = %d, data.size = %d, oldDataSize = %d!\n",n,addon,data.size(),oldDataSize);
	for (int i = 0; i < addon; i++) {
	  distances.push_back(distance2(x, data[n+i].second, DOW));
	}
	// 2.) sortiere Paare aus dist und data
	comparePair<double, DataType, 0> comp;
	vector_operations::sort(distances, data, comp);
      }

      // 3.) wenn noch nicht genügend Punkte, dann umgebende Boxen mit durchsuchen
      if (static_cast<int>(data.size()) < nData
	|| (level == 0 && distances[distances.size()-1] > boxBoundaryDist))
      {
	getSurroundingBoxes(center_box, level+1, nrs);
      }
    }
    data.resize(nData);
    return true;
  }


  Box* Box::provideBackgroundMesh(const FiniteElemSpace* feSpace)
  {
    if (boxMap.find(feSpace) != boxMap.end()) {
      if (boxMap[feSpace].first != feSpace->getMesh()->getChangeIndex()) {
	boxMap[feSpace].second->clearBox();
	boxMap[feSpace].second->init();
	boxMap[feSpace].second->fillBox(feSpace);
	boxMap[feSpace].first = feSpace->getMesh()->getChangeIndex();
	std::cout<<"max #data/box = "<<boxMap[feSpace].second->getMaxBoxSize()<<"\n";
      }
    } else {
      int N_const = 1<<Global::getGeo(WORLD);
      PointType min_corner_; min_corner_.set(-1.0);
      PointType max_corner_; max_corner_.set(1.0);
      Parameters::get("backgroundMesh->N", N_const); // Gesamtanzahl an Boxen
      Parameters::get("backgroundMesh->min_corner", min_corner_);
      Parameters::get("backgroundMesh->max_corner", max_corner_);

      // Anzahl an Box je Dimension anhand der Größe des gesamten Box bestimmen
      double sum_lengths = 0.0;
      std::vector<double> lengths;
      for (int i = 0; i < Global::getGeo(WORLD); i++) {
	lengths.push_back(abs(max_corner_[i]-min_corner_[i]));
	sum_lengths += lengths[i];
      }
      std::vector<int> N_;
      int N_sum = 0;
      for (int i = 0; i < Global::getGeo(WORLD); i++) {
	N_.push_back(boost::math::iround(pow(static_cast<double>(N_const),lengths[i]/sum_lengths)));
	N_sum += N_[i];
      }
      
      Box* box = new Box(Global::getGeo(WORLD), min_corner_, max_corner_, N_);
      box->fillBox(feSpace);

      boxMap[feSpace] = std::make_pair(feSpace->getMesh()->getChangeIndex(), box);
      std::cout<<"max #data/box = "<<boxMap[feSpace].second->getMaxBoxSize()<<"\n";
    }

    return boxMap[feSpace].second;
  }


  void Box::delete_all()
  {
    std::map<const FiniteElemSpace*, std::pair<int, Box*> >::iterator boxIter;
    for (boxIter = boxMap.begin(); boxIter != boxMap.end(); boxIter++) {
      (*boxIter).second.second->clearBox();
      delete (*boxIter).second.second;
    }
    boxMap.clear();
  }


  void Box::init()
  {
    int maxNr = 1;
    for (int i = 0; i < DOW; i++)
      maxNr *= N[i];
    boxData.resize(maxNr);

    std::cout<<"\n background mesh\n";
    std::cout<<"=================\n";
    std::cout<<"min_corner: "<<min_corner<<"\n";
    std::cout<<"max_corner: "<<max_corner<<"\n";
    std::cout<<"N: "<<N<<"\n\n";

    PointType tolerance; tolerance.set(2.0*DBL_TOL);
    min_corner -= tolerance;
    max_corner += tolerance;
  }


/// protected methods
///_____________________________________________________________________________________________


  inline int Box::idx2nr(std::vector<int>& idx)
  {
    int nr = idx.back();
    std::vector<int>::reverse_iterator iIter, nIter;
    for (iIter = idx.rbegin()+1, nIter = N.rbegin()+1; iIter < idx.rend(); iIter++,nIter++) {
      nr = *iIter + nr*(*nIter);
    }
    return nr;
  }


  inline void Box::nr2idx(int nr, std::vector<int>& idx)
  {
    idx.resize(DOW);
    switch (DOW) {
      case 1:
	idx[0] = nr % N[0];
	break;
      case 2:
	idx[1] = nr / N[0];
	idx[0] = nr % N[0];
	break;
      case 3:
	idx[2] = nr / (N[1]*N[0]);
	nr = nr % (N[1]*N[0]);
	idx[1] = nr / N[0];
	idx[0] = nr % N[0];
	break;
      default:
	throw(std::runtime_error("DOW="+boost::lexical_cast<std::string>(DOW)+", nur die Version für DOW=1 verwendet!"));
	break;
    }
  }


  void Box::getSurroundingBoxes(int center, int radius, std::set<int> &nrs)
  {
    std::set<int> newNrs;
    std::vector<int> idx, newIdx;
    nr2idx(center, idx); // --> idx
    newIdx.resize(idx.size());
    double incr = m_pi/3.0/radius;
    switch(DOW) {
      case 1:
	if (idx[0]-radius >= 0) { newIdx[0] = idx[0]-radius; newNrs.insert(idx2nr(newIdx)); }
	if (idx[0]+radius < N[0]) { newIdx[0] = idx[0]+radius; newNrs.insert(idx2nr(newIdx)); }
	break;
      case 2:
	for (double phi = 0.0; phi < 2.0*m_pi; phi += incr) {
	  newIdx[0] = idx[0] + boost::math::iround(radius*cos(phi));
	  newIdx[1] = idx[1] + boost::math::iround(radius*sin(phi));
	  if (newIdx[0] >= 0 && newIdx[0] < N[0]
	   && newIdx[1] >= 0 && newIdx[1] < N[1])
	    newNrs.insert(idx2nr(newIdx));
	}
	break;
      case 3:
	for (double phi = -m_pi; phi < m_pi; phi += incr) {
	for (double theta = 0.0; theta < m_pi; theta += incr) {
	  newIdx[0] = idx[0] + boost::math::iround(radius*sin(theta)*cos(phi));
	  newIdx[1] = idx[1] + boost::math::iround(radius*sin(theta)*sin(phi));
	  newIdx[2] = idx[2] + boost::math::iround(radius*cos(theta));
	  if (newIdx[0] >= 0 && newIdx[0] < N[0]
	   && newIdx[1] >= 0 && newIdx[1] < N[1]
	   && newIdx[2] >= 0 && newIdx[2] < N[2])
	    newNrs.insert(idx2nr(newIdx));
	}
	}
	break;
      default:
	ERROR("unknown dimension!\n");
    }
    
    std::set<int>::iterator nrsIter;
    for (nrsIter = nrs.begin(); nrsIter != nrs.end(); nrsIter++)
      newNrs.erase(*nrsIter);
    swap(nrs, newNrs);
  }
  

  void Box::getSurroundingBoxes(int nr, std::set<int> &surrounding_nrs)
  {
    std::vector<int> idx;
    nr2idx(nr, idx); // --> idx

    for (int i = 0; i < DOW; i++) {
      std::vector<int> idx_i = idx;
      idx_i[i] = std::min(N[i]-1, idx_i[i]+1);
      surrounding_nrs.insert(idx2nr(idx_i));
      idx_i[i] = std::max(0, idx_i[i]-2);
      surrounding_nrs.insert(idx2nr(idx_i));
    }
    surrounding_nrs.erase(nr);
  }


  void Box::getSurroundingBoxes(std::set<int> &nrs, std::set<int> &surrounding_nrs)
  {
    std::set<int>::iterator nrsIter;
    for (nrsIter = nrs.begin(); nrsIter != nrs.end(); nrsIter++)
      getSurroundingBoxes(*nrsIter, surrounding_nrs);

    for (nrsIter = nrs.begin(); nrsIter != nrs.end(); nrsIter++)
      surrounding_nrs.erase(*nrsIter);
  }


  /// berechne die minimale Entfernung der Ränder der Box zu x
  double Box::getBoxBoundaryDist(int nr, const PointType& x)
  {
    std::vector<int> idx;
    nr2idx(nr, idx);
    PointType min_c, max_c;
    for (int i = 0; i < DOW; i++) {
      min_c[i] = min_corner[i] + (max_corner[i]-min_corner[i])*idx[i]/N[i];
      max_c[i] = min_corner[i] + (max_corner[i]-min_corner[i])*(idx[i]+1)/N[i];
    }
    double dist = 1.e15;
    for (int i = 0; i < DOW; i++)
      dist = std::min(dist, std::min(abs(max_c[i]-x[i]), abs(min_c[i]-x[i])));
    return dist;
  }


  void Box::addData(PointType x, DataType data)
  {
    int nr = getBox(x);
    if (nr < 0 || nr >= static_cast<int>(boxData.size()))
      throw(std::runtime_error("box-nr out of range: nr = "+boost::lexical_cast<std::string>(nr)+", x = "+boost::lexical_cast<std::string>(x)+", boxSize = "+boost::lexical_cast<std::string>(boxData.size())));
    boxData[nr].push_back(data);
  }
  
}