/******************************************************************************
 *
 * 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 "Tools.h"	// tools::Regression

namespace AMDiS { namespace extensions {

  /**
   * strategies:
   * 0 .. get nearest point to x and eval DOFVector at this DOF-index
   * 1 .. get n nearest points, calc weighted sum of data at these points
   * 2 .. get n nearest points, calc regression plane and eval at x
   * 3 .. get n nearest points, calc weighted regression plane and eval at x
   **/
  template<typename T>
  T Box::evalAtPoint(const DOFVector<T>& vec, const PointType& x, int strategy, int nrOfPoints)
  {
    T value = 0;
    switch (strategy) {
      case 0:
	value = evalAtPoint_simple(vec, x);
	break;
      case 1:
	value = evalAtPoint_weighted_sum(vec, x, nrOfPoints);
	break;
      case 2:
	value = evalAtPoint_regression(vec, x, nrOfPoints);
	break;
      case 3:
	value = evalAtPoint_weighted_regression(vec, x, nrOfPoints);
	break;
      default:
	ERROR("ERROR: unknown strategy [%d]!\n",strategy);
	break;
    }
    return value;
  }

  
  template<typename T>
  T Box::evalAtPoint_simple(const DOFVector<T>& vec, const PointType& x)
  {
    DataType dof;
    T value; nullify(value);
    if (getNearestData(x, dof))
      value = vec[dof.first];
    return value;
  }

  
  template<typename T>
  T Box::evalAtPoint_weighted_sum(const DOFVector<T>& vec, const PointType& x, int nrOfPoints)
  {
    T value; nullify(value);

    std::vector<DataType> dofs;
    int nPoints = (nrOfPoints > 0 ? nrOfPoints : vec.getFeSpace()->getBasisFcts()->getNumber()+1);

    bool inside = getNearestData(x, dofs, nPoints);
    if (!inside)
      return value;

    // gewichtete Summe der einzelnen Funktionswerte in der Nähe von x
    std::vector<double> weights;
    for (size_t i = 0; i < dofs.size(); i++) {
      double dist = distance(x, dofs[i].second, DOW);
      weights.push_back(1.0/std::max(1.e-8, dist));
    }

    double weight = 0.0;
    for (size_t i = 0; i < dofs.size(); i++) {
      value += vec[dofs[i].first] * weights[i];
      weight += weights[i];
    }
    value *= 1.0/weight;

    return value;
  }


  template<typename T>
  T Box::evalAtPoint_regression(const DOFVector<T>& vec, const PointType& x, int nrOfPoints)
  {
    T value; nullify(value);

    std::vector<DataType> dofs;
    int nPoints = (nrOfPoints > 0 ? nrOfPoints : vec.getFeSpace()->getBasisFcts()->getNumber()+1);

    bool inside = getNearestData(x, dofs, nPoints);
    if (!inside)
      return value;

    std::vector<WorldVector<double> > points;
    std::vector<double> values;
    for (size_t i = 0; i < dofs.size(); i++) {
      points.push_back(dofs[i].second);
      values.push_back(vec[dofs[i].first]);
    }

    tools::Regression reg(DOW);
    reg.apply(points, values);
    value = reg.value(x);
    return value;
  }


  template<typename T>
  T Box::evalAtPoint_weighted_regression(const DOFVector<T>& vec, const PointType& x, int nrOfPoints)
  {
    T value; nullify(value);

    std::vector<DataType> dofs;
    int nPoints = (nrOfPoints > 0 ? nrOfPoints : vec.getFeSpace()->getBasisFcts()->getNumber()+1);

    bool inside = getNearestData(x, dofs, nPoints);
    if (!inside)
      return value;

    std::vector<WorldVector<double> > points;
    std::vector<double> values,weights;
    for (size_t i = 0; i < dofs.size(); i++) {
      points.push_back(dofs[i].second);
      values.push_back(vec[dofs[i].first]);

      double dist = distance(x, dofs[i].second, DOW);
      weights.push_back(1.0/std::max(1.e-8, dist));
    }

    tools::Regression reg(DOW);
    reg.apply(points, values, weights);
    value = reg.value(x);

    return value;
  }

} }