#include <Core/Nonlinear_Equations/NONLINEAR_FUNCTION.h>
#include <Core/Krylov_Solvers/KRYLOV_SYSTEM_BASE.h>
#include <Core/Krylov_Solvers/KRYLOV_VECTOR_BASE.h>
#include <Core/Log/LOG.h>
#include <Core/Random_Numbers/RANDOM_NUMBERS.h>
#include <Core/Read_Write/OCTAVE_OUTPUT.h>
#include <cfloat>
namespace PhysBAM{
//#####################################################################
// Destructor
//#####################################################################
template<class T> NONLINEAR_FUNCTION<T(KRYLOV_VECTOR_BASE<T>&)>::
~NONLINEAR_FUNCTION()
{
}
//#####################################################################
// operator()
//#####################################################################
template<class T> T NONLINEAR_FUNCTION<T(KRYLOV_VECTOR_BASE<T>&)>::
operator()(const KRYLOV_VECTOR_BASE<T>& x) const
{
    T E=0;
    Compute(x,0,0,&E);
    return E;
}
//#####################################################################
// Function Test
//#####################################################################
template<class T> void NONLINEAR_FUNCTION<T(KRYLOV_VECTOR_BASE<T>&)>::
Test(const KRYLOV_VECTOR_BASE<T>& x,KRYLOV_SYSTEM_BASE<T>& h,
    const char* name_f,const char* name_df,bool project_dx,
    std::function<void(KRYLOV_VECTOR_BASE<T>& x)> make_feasible) const
{
    T eps=(T)1e-6;
    RANDOM_NUMBERS<T> random;
    KRYLOV_VECTOR_BASE<T> *dx=x.Clone_Default();
    KRYLOV_VECTOR_BASE<T> *g0=x.Clone_Default();
    KRYLOV_VECTOR_BASE<T> *g1=x.Clone_Default();
    KRYLOV_VECTOR_BASE<T> *a=x.Clone_Default();
    KRYLOV_VECTOR_BASE<T> *b=x.Clone_Default();

    for(int i=0,n=dx->Raw_Size();i<n;i++)
        dx->Raw_Get(i)=random.Get_Uniform_Number(-eps,eps);
    if(project_dx) h.Project(*dx);

    b->Copy(1,x,*dx);
    if(make_feasible)
    {
        PHYSBAM_ASSERT(!project_dx);
        make_feasible(*b);
        dx->Copy(-1,x,*b);
    }

    T e0=0,e1=0;
    Compute(x,&h,g0,&e0);
    h.Multiply(*dx,*a);
    h.Project(*a);

    // new x
    Compute(*b,&h,g1,&e1);
    h.Multiply(*dx,*b);
    h.Project(*b);
    T test0a=(h.Inner_Product(*g1,*dx)+h.Inner_Product(*g0,*dx))/(2*eps);
    T test0b=(e1-e0)/eps;
    T test0=abs(test0b-test0a);

    LOG::printf("diff test %s:  diff %.16P  exact %.16P   abs-error %P   rel-error %P\n",
        name_f,test0b,test0a,test0,test0/max(max(abs(test0b),abs(test0a)),(T)FLT_MIN));

    a->Copy(1,*b,*a);
    a->Copy((T).5,*a);
    T test1a=sqrt(h.Inner_Product(*a,*a))/eps;

    b->Copy(-1,*g0,*g1);
    T test1b=sqrt(h.Inner_Product(*b,*b))/eps;

    a->Copy(-1,*b,*a);
    T test1=sqrt(h.Inner_Product(*a,*a))/eps;

    LOG::printf("diff test %s:  diff %.16P  exact %.16P   abs-error %P   rel-error %P\n",
        name_df,test1b,test1a,test1,test1/max(max(test1b,test1a),(T)FLT_MIN));

    delete dx;
    delete g0;
    delete g1;
    delete a;
    delete b;
}
//#####################################################################
// Function Make_Feasible
//#####################################################################
template<class T> void NONLINEAR_FUNCTION<T(KRYLOV_VECTOR_BASE<T>&)>::
Make_Feasible(KRYLOV_VECTOR_BASE<T>& x) const
{
}
template class NONLINEAR_FUNCTION<double(KRYLOV_VECTOR_BASE<double>&)>;
template class NONLINEAR_FUNCTION<float(KRYLOV_VECTOR_BASE<float>&)>;
}
