//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
// G4SimplexDownhill inline methods implementation
//
// Author: Tatsumi Koi (SLAC/SCCS), 2007
// --------------------------------------------------------------------------

#include <cfloat>
#include <iostream>
#include <numeric>

template <class T>
void G4SimplexDownhill<T>::init()
{
  alpha = 2.0;  // refrection coefficient:  0 < alpha
  beta  = 0.5;  // contraction coefficient:   0 < beta < 1
  gamma = 2.0;  // expantion coefficient:  1 < gamma

  maximum_no_trial = 10000;
  max_se           = FLT_MIN;
  // max_ratio = FLT_EPSILON/1;
  max_ratio = DBL_EPSILON / 1;
  minimized = false;
}

/*

void G4SimplexDownhill<class T>::
SetFunction( G4int n , G4double( *afunc )( std::vector < G4double > ) )
{
   numberOfVariable = n;
   theFunction = afunc;
   minimized = false;
}

*/

template <class T>
G4double G4SimplexDownhill<T>::GetMinimum()
{
  initialize();

  // First Tryal;

  // G4cout << "Begin First Trials" << G4endl;
  doDownhill();
  // G4cout << "End First Trials" << G4endl;

  auto it_minh = std::min_element(currentHeights.cbegin(), currentHeights.cend());
  G4int imin = 0;
  G4int i    = 0;
  for(auto it = currentHeights.cbegin(); it != currentHeights.cend(); ++it)
  {
    if(it == it_minh)
    {
      imin = i;
    }
    ++i;
  }
  minimumPoint = currentSimplex[imin];

  // Second Trial

  // std::vector< G4double > minimumPoint = currentSimplex[ 0 ];
  initialize();

  currentSimplex[numberOfVariable] = minimumPoint;

  // G4cout << "Begin Second Trials" << G4endl;
  doDownhill();
  // G4cout << "End Second Trials" << G4endl;

  G4double sum =
    std::accumulate(currentHeights.begin(), currentHeights.end(), 0.0);
  G4double average = sum / (numberOfVariable + 1);
  G4double minimum = average;

  minimized = true;

  return minimum;
}

template <class T>
void G4SimplexDownhill<T>::initialize()
{
  currentSimplex.resize(numberOfVariable + 1);
  currentHeights.resize(numberOfVariable + 1);

  for(G4int i = 0; i < numberOfVariable; ++i)
  {
    std::vector<G4double> avec(numberOfVariable, 0.0);
    avec[i]           = 1.0;
    currentSimplex[i] = std::move(avec);
  }

  // std::vector< G4double > avec ( numberOfVariable , 0.0 );
  std::vector<G4double> avec(numberOfVariable, 1);
  currentSimplex[numberOfVariable] = std::move(avec);
}

template <class T>
void G4SimplexDownhill<T>::calHeights()
{
  for(G4int i = 0; i <= numberOfVariable; ++i)
  {
    currentHeights[i] = getValue(currentSimplex[i]);
  }
}

template <class T>
std::vector<G4double> G4SimplexDownhill<T>::calCentroid(G4int ih)
{
  std::vector<G4double> centroid(numberOfVariable, 0.0);

  G4int i = 0;
  for(const auto & it : currentSimplex)
  {
    if(i != ih)
    {
      for(G4int j = 0; j < numberOfVariable; ++j)
      {
        centroid[j] += it[j] / numberOfVariable;
      }
    }
    ++i;
  }

  return centroid;
}

template <class T>
std::vector<G4double> G4SimplexDownhill<T>::getReflectionPoint(
  std::vector<G4double> p, std::vector<G4double> centroid)
{
  // G4cout << "Reflection" << G4endl;

  std::vector<G4double> reflectionP(numberOfVariable, 0.0);

  for(G4int i = 0; i < numberOfVariable; ++i)
  {
    reflectionP[i] = (1 + alpha) * centroid[i] - alpha * p[i];
  }

  return reflectionP;
}

template <class T>
std::vector<G4double> G4SimplexDownhill<T>::getExpansionPoint(
  std::vector<G4double> p, std::vector<G4double> centroid)
{
  // G4cout << "Expantion" << G4endl;

  std::vector<G4double> expansionP(numberOfVariable, 0.0);

  for(G4int i = 0; i < numberOfVariable; ++i)
  {
    expansionP[i] = (1 - gamma) * centroid[i] + gamma * p[i];
  }

  return expansionP;
}

template <class T>
std::vector<G4double> G4SimplexDownhill<T>::getContractionPoint(
  std::vector<G4double> p, std::vector<G4double> centroid)
{
  std::vector<G4double> contractionP(numberOfVariable, 0.0);

  for(G4int i = 0; i < numberOfVariable; ++i)
  {
    contractionP[i] = (1 - beta) * centroid[i] + beta * p[i];
  }

  return contractionP;
}

template <class T>
G4bool G4SimplexDownhill<T>::isItGoodEnough()
{
  G4double sum =
    std::accumulate(currentHeights.begin(), currentHeights.end(), 0.0);
  G4double average = sum / (numberOfVariable + 1);

  G4double delta = 0.0;
  for(G4int i = 0; i <= numberOfVariable; ++i)
  {
    delta += std::abs(currentHeights[i] - average);
  }

  G4bool result = false;
  if (average > 0.0)
  {
    result = ((delta / (numberOfVariable + 1) / average) < max_ratio);
  }
  return result;
}

template <class T>
void G4SimplexDownhill<T>::doDownhill()
{
  G4int nth_trial = 0;

  while(nth_trial < maximum_no_trial)
  {
    calHeights();

    if(isItGoodEnough())
    {
      break;
    }

    auto it_maxh = std::max_element(currentHeights.cbegin(), currentHeights.cend());
    auto it_minh = std::min_element(currentHeights.cbegin(), currentHeights.cend());

    G4double h_H = *it_maxh;
    G4double h_L = *it_minh;

    G4int ih      = 0;
    G4int il      = 0;
    G4double h_H2 = 0.0;
    G4int i       = 0;
    for(auto it = currentHeights.cbegin(); it != currentHeights.cend(); ++it)
    {
      if(it == it_maxh)
      {
        ih = i;
      }
      else
      {
        h_H2 = std::max(h_H2, *it);
      }

      if(it == it_minh)
      {
        il = i;
      }
      ++i;
    }

    std::vector<G4double> centroidPoint = calCentroid(ih);

    // REFLECTION
    std::vector<G4double> reflectionPoint =
      getReflectionPoint(currentSimplex[ih], centroidPoint);

    G4double h = getValue(reflectionPoint);

    if(h <= h_L)
    {
      // EXPANSION
      std::vector<G4double> expansionPoint =
        getExpansionPoint(reflectionPoint, std::move(centroidPoint));
      G4double hh = getValue(expansionPoint);

      if(hh <= h_L)
      {
        // Replace
        currentSimplex[ih] = std::move(expansionPoint);
        // G4cout << "A" << G4endl;
      }
      else
      {
        // Replace
        currentSimplex[ih] = std::move(reflectionPoint);
        // G4cout << "B1" << G4endl;
      }
    }
    else
    {
      if(h <= h_H2)
      {
        // Replace
        currentSimplex[ih] = std::move(reflectionPoint);
        // G4cout << "B2" << G4endl;
      }
      else
      {
        if(h <= h_H)
        {
          // Replace
          currentSimplex[ih] = std::move(reflectionPoint);
          // G4cout << "BC" << G4endl;
        }
        // CONTRACTION
        std::vector<G4double> contractionPoint =
          getContractionPoint(currentSimplex[ih], std::move(centroidPoint));
        G4double hh = getValue(contractionPoint);
        if(hh <= h_H)
        {
          // Replace
          currentSimplex[ih] = std::move(contractionPoint);
          // G4cout << "C" << G4endl;
        }
        else
        {
          // Replace
          for(G4int j = 0; j <= numberOfVariable; ++j)
          {
            std::vector<G4double> vec(numberOfVariable, 0.0);
            for(G4int k = 0; k < numberOfVariable; ++k)
            {
              vec[k] = (currentSimplex[j][k] + currentSimplex[il][k]) / 2.0;
            }
            currentSimplex[j] = std::move(vec);
          }
        }
      }
    }

    ++nth_trial;
  }
}

template <class T>
std::vector<G4double> G4SimplexDownhill<T>::GetMinimumPoint()
{
  if(!minimized)
  {
    GetMinimum();
  }

  auto it_minh = std::min_element(currentHeights.cbegin(), currentHeights.cend());

  G4int imin = 0;
  G4int i    = 0;
  for(auto it = currentHeights.cbegin(); it != currentHeights.cend(); ++it)
  {
    if(it == it_minh)
    {
      imin = i;
    }
    ++i;
  }
  minimumPoint = currentSimplex[imin];

  return minimumPoint;
}