// example of using multiroot finder 
// based on GSL algorithm 
//
// The MultiRootFinder is based on GSL and it requires the MathMore library 
// installed 
//
// Usage: 
//  >.x exampleMultiRoot.C()
// or 
// >.x exampleMultiRoot(algoname,printlevel)
//
// where algoname is for an algorithm not using the derivatives:
//  hybridS (default) , hybrid, dnewton, broyden
// 
#include "RConfigure.h"
#ifdef R__HAS_MATHMORE
#include "Math/MultiRootFinder.h"
#endif
#include "Math/WrappedMultiTF1.h"
//#include "TF2.h"
#include "TError.h"
// example of using multi root finder based on GSL 
// need to use an algorithm not requiring the derivative
//like hybrids (default), hybrid, dnewton, broyden
double MyFunc1 (double *x, double *p1 ) {
  return ( p1[0] + x[0] + p1[2]*x[1] + p1[3]*x[2] + p1[4]*x[3] + p1[5]*x[4] + p1[6]*x[5] - p1[1] * x[6] * x[7] );
}
double MyFunc2 (double *x, double *p2 ) {
  return ( p2[0] + p2[2]*x[0] + p2[3]*x[1] + p2[4]*x[2] + x[3] + p2[5]*x[4] + p2[6]*x[5] + p2[1] * x[6] * x[7] );
}
double MyFunc3 (double *x, double *p3 ) {
  return ( p3[0] + p3[2]*x[0] + x[1] + p3[3]*x[2] + p3[4]*x[3] + p3[5]*x[4] + p3[6]*x[5] - p3[1] * x[6] * sqrt(1 - x[7]*x[7] ) );
}
double MyFunc4 (double *x, double *p4 ) {
  return ( p4[0] + p4[2]*x[0] + p4[3]*x[1] + p4[4]*x[2] + p4[5]*x[3] + x[4] + p4[6]*x[5] + p4[1] * x[6] * sqrt(1 - x[7]*x[7] ) );
}
double MyFunc5 (double *x, double *p5 ) {
  return ( p5[0] + p5[2]*x[0] + p5[3]*x[1] + x[2] + p5[4]*x[3] + p5[5]*x[4] + p5[6]*x[5] - p5[1] * sqrt(1 - x[6]*x[6]) + p5[7]*x[7] );
}
double MyFunc6 (double *x, double *p6 ) {
  return ( p6[0] + p6[2]*x[0] + p6[3]*x[1] + p6[4]*x[2] + p6[5]*x[3] + p6[6]*x[4] + x[5] + p6[1] * sqrt(1 - x[6]*x[6]) + p6[7]*x[7] );
}
double MyFunc7 (double *x, double *p7 ) {
  return ( sqrt(x[0]*x[0] + x[1]*x[1] + x[2]*x[2]) + p7[1]*x[3] + p7[2]*x[4] + p7[3]*x[5] + p7[4]*x[6] + p7[5]*x[7] - p7[0] );
}
double MyFunc8 (double *x, double *p8 ) {
  return ( p8[1]*x[0] + p8[2]*x[1] + p8[3]*x[2] + sqrt(x[3]*x[3] + x[4]*x[5] + x[5]*x[5]) + p8[4]*x[6] + p8[5]*x[7] - p8[0] );
}
void exampleMultiRoot(const char * algo = 0, int printlevel = 1) {
#ifndef R__HAS_MATHMORE
  Error("exampleMultiRoot","libMathMore is not available - cannot run this tutorial");
#else 
  ROOT::Math::MultiRootFinder r(algo);
  //defining the function
  TF1 * MyFunc1 = new TF1("MyFunc1", MyFunc1, 0, 0, 7);
  TF1 * MyFunc2 = new TF1("MyFunc2", MyFunc2, 0, 0, 7);
  TF1 * MyFunc3 = new TF1("MyFunc3", MyFunc3, 0, 0, 7);
  TF1 * MyFunc4 = new TF1("MyFunc4", MyFunc4, 0, 0, 7);
  TF1 * MyFunc5 = new TF1("MyFunc5", MyFunc5, 0, 0, 8);
  TF1 * MyFunc6 = new TF1("MyFunc6", MyFunc6, 0, 0, 8);
  TF1 * MyFunc7 = new TF1("MyFunc7", MyFunc7, 0, 0, 6);
  TF1 * MyFunc8 = new TF1("MyFunc8", MyFunc8, 0, 0, 6);


  double mom = 0.267;
  MyFunc1->SetParameters( -0.34,  mom, 0., 0., 0., 0., 0.);
  MyFunc2->SetParameters( -0.28,  mom, 0., 0., 0., 0., 0.);
  MyFunc3->SetParameters(  0.39,  mom, 0., 0., 0., 0., 0.);
  MyFunc4->SetParameters(  0.037, mom, 0., 0., 0., 0., 0.);
  MyFunc5->SetParameters( -0.74,  mom, 0., 0., 0., 0., 0., 0.);
  MyFunc6->SetParameters( -1.02,  mom, 0., 0., 0., 0., 0., 0.);
  MyFunc7->SetParameters(  1.11, 0., 0., 0., 0., 0. );
  MyFunc8->SetParameters(  0.97, 0., 0., 0., 0., 0. );

  // wrap the functions
  ROOT::Math::WrappedMultiTF1 g1(*MyFunc1, 8);
  ROOT::Math::WrappedMultiTF1 g2(*MyFunc2, 8);
  ROOT::Math::WrappedMultiTF1 g3(*MyFunc3, 8);
  ROOT::Math::WrappedMultiTF1 g4(*MyFunc4, 8);
  ROOT::Math::WrappedMultiTF1 g5(*MyFunc5, 8);
  ROOT::Math::WrappedMultiTF1 g6(*MyFunc6, 8);
  ROOT::Math::WrappedMultiTF1 g7(*MyFunc7, 8);
  ROOT::Math::WrappedMultiTF1 g8(*MyFunc8, 8);

  r.AddFunction(g1);
  r.AddFunction(g2);
  r.AddFunction(g3);
  r.AddFunction(g4);
  r.AddFunction(g5);
  r.AddFunction(g6);
  r.AddFunction(g7);
  r.AddFunction(g8);

  r.SetPrintLevel(printlevel);
  // starting point - This starting point is the solution and when used here, program converges.
  double x0[8]={ 0.31, -0.62, 0.86, 0.31, 0.2, 0.9, -0.89, 0.12};
  /*
 root [0] .x exampleMultiRoot.C("kDNewton",1)
 GSLMultiRootFinder::Solve:dnewton max iterations 100 and  tolera1e-06
 Info in <ROOT::Math::GSLMultiRootFinder::Solve>: The iteration converged
 GSL Algorithm used is :  dnewton
 Number of iterations  =  4
 Root values     = x[0] =     0.249173   x[1] =    -0.557852   x[2] =     0.926722   x[3] =     0.370827   x[4] =     0.130852   x[5] =     0.833278   x[6] =    -0.714794   x[7] =     0.475908
 Function values = f[0] = -4.58058e-10   f[1] =  4.58043e-10   f[2] = -5.12292e-10   f[3] =  5.12289e-10   f[4] =   2.1217e-10   f[5] = -2.12151e-10   f[6] =  8.80505e-11   f[7] =  1.00293e-10
  */
  // With this starting value, program does not converge.
  //    double x0[8]={ 0.01, -0.02, 0.06, 0.01, 0., 0.01, -0.09, 0.02};
  /*
 root [0] .x exampleMultiRoot.C("kDNewton",1)
 GSLMultiRootFinder::Solve:dnewton max iterations 100 and  tolerance 1e-06
 Info in <ROOT::Math::GSLMultiRootFinder::Solve>: exceeded max iterations, reached tolerance is not sufficient; absTol = 1e-06
  */
  r.Solve(x0);
#endif
}