#include "amma/Num/NumModelData.hh"
#include "amma/Num/NumModelClosest.hh"
#include "amma/Num/NumModelPolynomial.hh"
#include "amma/Num/NumModelQuadratic.hh"
#include "amma/Num/NumModelBspline.hh"
#include "amma/Num/NumModelHybridBspline.hh"
#include "amma/Vector2d.hh"
#include <fstream.h>

main ()
{
  // Setup Data
  VectorC min = Vector2dC (0,0);
  VectorC max = Vector2dC (2,2);
  VectorC mid = Vector2dC (1,1);
  IntSArray1dC numPoints (2); numPoints.Fill (3);
  NumVVDataSetC train;
  train.Append (min,max);
  train.Append (max,max);
  train.Append (mid,min);
  train.Append (Vector2dC(0,2), Vector2dC(1,-1));
  train.Append (Vector2dC(2,0), Vector2dC(-1,1));
  train.Append (Vector2dC(1,0), Vector2dC(0.5,-0.5));
  train.Append (Vector2dC(-1,0),Vector2dC(-0.5,0.5));
  train.Append (Vector2dC(0,1), Vector2dC(0.5,0.5));
  train.Append (Vector2dC(0,-1),Vector2dC(-0.5,-0.5));
  // Setup Models
  SArray1dC<NumModelC> models (6);
  models[0] = NumModelDataC ();
  models[1] = NumModelClosestC ();
  models[2] = NumModelPolynomialC (2);
  models[3] = NumModelQuadraticC ();
  models[4] = NumModelBsplineC (numPoints,min,max);
  models[5] = NumModelHybridBsplineC (numPoints,min,max);
  // Fit and save
  {
    ofstream out ("testNumModel.data");
    FOR_SARRAY1 (models,i1) {
      models[i1].Fit (train);
      models[i1].Save (out);
    }
  }
  // Load as function and get info
  ifstream in ("testNumModel.data");
  for (IndexT i2 = 0; i2 < 6; i2++) {
    NumFuncC function (in);
    cout << "Function Name: " << function.GetName() << "\n";
    cout << "Size X and Y: " << function.SizeX() << " " << function.SizeY() << "\n";
    cout << function.GetInfo() << "\n";
    VectorC testX = Vector2dC (1.5,0.5);
    cout << "Output at X=" << testX << " is Y=" << function(testX) << "\n";
    cout << "Jacobian at X is " << function.Jacobian(testX) << "\n";
  }
}