/* Copyright (c) 2015  Gerald Knizia
 * 
 * This file is part of the IboView program (see: http://www.iboview.org)
 * 
 * IboView is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, version 3.
 * 
 * IboView is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with bfint (LICENSE). If not, see http://www.gnu.org/licenses/
 * 
 * Please see IboView documentation in README.txt for:
 * -- A list of included external software and their licenses. The included
 *    external software's copyright is not touched by this agreement.
 * -- Notes on re-distribution and contributions to/further development of
 *    the IboView software
 */

#include <string>
#include <iostream>
// #include <boost/format.hpp>
// using boost::format;
#include <vector>
#include <sstream>

#include "format.h"
#include "Ir.h"
#include "CtDftFunc.h"
#include "CtIo.h"

#include "xc_meta.h"
typedef void (*FXcFn1)(double &E, double *dE, double const *D, double const Factor);


namespace ct {

struct FXcFunctionalEntry
{
   char const
      *pName;
   FXcFn1
      pFn1; // function to evaluate value + 1st derivative of xc functional.
   double
      Factor;
   FXcFunctionalEntry(char const *pName_, FXcFn1 pFn1_, double Factor_) : pName(pName_), pFn1(pFn1_), Factor(Factor_) {}
};

struct FXcFunctionalImpl
{
   FXcFunctionalImpl(std::string const &Name);
   ~FXcFunctionalImpl();

   void Eval(double *pOut, size_t nOutSt, double const *pIn, size_t nInSt, size_t nPts, size_t nDiff);

   bool
      NeedTau,
      NeedSigma;
   bool
      m_SetupFailed;

   inline uint nInputLength();

   std::vector<FXcFunctionalEntry>
      Components;
   std::string Desc() const;
};


uint FXcFunctionalImpl::nInputLength()
{
   assert(!NeedTau || NeedSigma); // tau only supported if sigma is also given.
   if (NeedTau) return 2 + 3 + 2; // rhoc, rhoo, sigmacc, sigmaco, sigmaoo, tauc, tauo
   if (NeedSigma) return 2 + 3; // rhoc, rhoo, sigmacc, sigmaco, sigmaoo
   return 2; // rhoc, rhoo
}


FXcFunctionalImpl::FXcFunctionalImpl(std::string const &Name_)
   : m_SetupFailed(false)
{
   NeedSigma = false;
   NeedTau = false;

   if ( Name_ == "LDA" ) {
      Components.push_back(FXcFunctionalEntry("DIRACX", rx::xc_diracx, 1.0));
      Components.push_back(FXcFunctionalEntry("PW92C", rx::xc_pw92c, 1.0));
   } else if ( Name_ == "PBE" ) {
      Components.push_back(FXcFunctionalEntry("PBEX", rx::xc_pbex, 1.0));
      Components.push_back(FXcFunctionalEntry("PBEC", rx::xc_pbec, 1.0));
      NeedSigma = true;
   } else if (Name_ == "None" || Name_.empty()) {
   } else {
      m_SetupFailed = true;
   }
}

std::string FXcFunctionalImpl::Desc() const
{
   std::stringstream str;
   str << " Density functional:         ";
   if (m_SetupFailed) {
      str << "ERROR/XC-SETUP-FAILED";
   } else {
      for ( uint iComp = 0; iComp < Components.size(); ++ iComp ) {
         FXcFunctionalEntry const &e = Components[iComp];
         str << fmt::format("  {} ={:6.3f}", e.pName, e.Factor);
      }
   }
   return str.str();
}



void FXcFunctionalImpl::Eval(double *pOut, size_t nOutSt, double const *pIn, size_t nInSt, size_t nPts, size_t nDiff)
{
   uint
      nInp = nInputLength(),
      nOut = 1 + nInp; // Energy + dEnergy/dInputs.
   assert(nOutSt >= nOut && nInSt >= nInp);
   assert(nDiff == 1);
   for (size_t iPt = 0; iPt < nPts; ++ iPt) {
      double const
         *pInPt = &pIn[nInSt * iPt];
      double
         *pOutPt = &pOut[nOutSt * iPt];
      for (size_t iOut = 0; iOut < nOut; ++ iOut)
         pOutPt[iOut] = 0;
      for (size_t iComp = 0; iComp < Components.size(); ++ iComp) {
         FXcFunctionalEntry &e = Components[iComp];
         e.pFn1(pOutPt[0], &pOutPt[1], pInPt, e.Factor);
         // ^- this adds to the previous results.
      }
   }
}

FXcFunctionalImpl::~FXcFunctionalImpl()
{
}

FXcFunctional::FXcFunctional(std::string const &Name)
   : p(new FXcFunctionalImpl(Name)),
     m_Name(Name)
{
}

FXcFunctional::~FXcFunctional()
{
   delete p;
}

void FXcFunctional::Eval(double *pOut, size_t nOutSt, double const *pIn, size_t nInSt, size_t nPts, size_t nDiff)
{
   return p->Eval(pOut, nOutSt, pIn, nInSt, nPts, nDiff);
}

bool FXcFunctional::NeedSigma() const
{
   return p->NeedSigma;
}

bool FXcFunctional::NeedTau() const
{
   return p->NeedTau;
}

uint FXcFunctional::nOrder() const
{
   return NeedSigma()? 1 : 0;
}

std::string FXcFunctional::Desc() const
{
   return p->Desc();
}

} // namespace ct.



// #define XC_FUNC_TEST
#ifdef XC_FUNC_TEST
#include <boost/format.hpp>
using boost::format;

static void PrintArray1(char const *pTitle, double *pValues, uint nValues, int nStep = -1, char const *pFmt = "%11.5f")
{
   using ct::xout;
   if ( nStep == -1 ) {
      // find a reasonable step size such that we get coverage of the entire
      // array but don't print out more than ~10-15 numbers.
      nStep = (12+nValues)/13;
   }
   xout << format("[%5i] %-20s") % nValues % pTitle;
   for ( uint i = 0; i < nValues; i += nStep )
      xout << format(pFmt) % pValues[i];
   xout << "\n";
}


int main_xc_func_test(int argc, char *argv[])
{
   using namespace ct;
   FXcFunctionalPtr
//       pRef = new FXcFunctional("DIRAC");
//       pRef = new FXcFunctional("PW92C");
//       pRef = new FXcFunctional("PBEX");
      pRef = new FXcFunctional("PBEC");

   double In[5] = {0.123, -0.05, 0.07, 0.01, 0.02}; // rhoc, rhoo, sigmacc, sigmaco, sigmaoo
//    double In[5] = {0.123, -0.05, 0., 0., 0.}; // rhoc, rhoo, sigmacc, sigmaco, sigmaoo
//    double In[5] = {0.123, -0.05, 1e-8, 0., 1e-8}; // rhoc, rhoo, sigmacc, sigmaco, sigmaoo
//    double In[5] = {0.123, 0., 0.05, 0., 0.}; // rhoc, rhoo, sigmacc, sigmaco, sigmaoo
//    double In[5] = {0.01, 0., 0.01, 0., 0.}; // rhoc, rhoo, sigmacc, sigmaco, sigmaoo
   double OutRef[6] = {0}; // energy and potentials.
   double OutTst[6] = {0}; // energy and potentials.

   pRef->Eval(&OutRef[0], 6, &In[0], 5, 1, 1);

   typedef void (*FXcFn1)(double &E, double *dE, double const *D, double const Factor);
   FXcFn1 xcfn1 = rx::xc_pbec;
   xcfn1(OutTst[0], &OutTst[1], &In[0], 1.0);

   uint nVar = pRef->NeedSigma()? 5 : 2;
   PrintArray1("Inputs", &In[0], nVar, 1, "%16.8f");
   PrintArray1("Output (ref)", &OutRef[0], nVar+1, 1, "%16.8f");
   PrintArray1("Output (tst)", &OutTst[0], nVar+1, 1, "%16.8f");
   return 0;
}

#endif // XC_FUNC_TEST
