/* demo.cc
   This is the demo main program for the Krigifier and the Kriging approximation classes
 It creates a random function, using the parameters specified in "krig.dat",
 evaluates it at several sites, and uses these values to
 create an approximation.  The approximation is then evaluated systematically and the results
 are output to the file "results.out"

 Anthony D. Padula    adpadu@wm.edu
 5/31/00
*/

/* INSTRUCTIONS FOR USE
 To compile this, you need to have downloaded the files :

 krigify.cc, krigify.h, demo5.cc rngs.c, rngs.h, rvgs.c, rvgs.h,
 krig.cc, krig.h, gamma.cc, gamma.h, ParamEstimate.cc, ParamEstimate.h,
 PatternSearch.cc, PatternSearch.h, CompassSearch.cc, CompassSearch.h,
 approx.cc, approx.h, Makefile


 You will also need the libraries CLAPACK, BLAS, and F2C, which can be obtained
 from http://www.netlib.org/
 You should modify the macros at the top of the Makefile to suit your compiler and libraries

 Then, on your command line you run

% make demo
% ./demo

*/

/* WHAT TO DO WITH THE OUTPUT 
Once you have run this program, you can plot the approximation using gnuplot.  Supposing
your CWD contains results.out, in gnuplot you should run the commmands:

  set parametric
  splot "results.out" with lines

This will produce a nice 3-d graph of the approximation  


Alternatively, you can compile and run plotpoints by

% g++ plotpoints.cc -o plotpoints
% ./plotpoints results.out results.ps

This produces a graph in postscript format, which can then be viewed or printed.
*/


#include "krig.h"
#include "krigify.h"
#include "rngs.h"

void randPoints( Matrix<double> & X, Vector<double>& lower, Vector<double> & upper);
/* Fills X with uniform random points between the bounds.  Thus
   lower[j] <= X[i][j] <= upper[j]   for all 0<=j<p
*/



int main( void ) {

  // The class object for the approximation
  krigapprox MyApproximation; 
  
  // The class object for the random function (i.e. krigifier function)
  //randfunc MyRandomFunction(960587311, 42); // choose the seed for the RNG
  randfunc MyRandomFunction;   // or take a seed from the system clock

  int m = 100;       // The number of initial design sites used to create the approximation
  int p;                       // The dimension of the space ( Given in "krig.dat" )
  Matrix<double> X(1,1);       // The matrix to store the initial design sites
 
  Vector<double> y(m);         // Vector to store the function values at the initial design sites 

  Vector<double> lower, upper; // Vectors to store the upper and lower bounds of the space 

  ofstream outf;

  outf.open( "results.out");  // open an output file 

  MyRandomFunction.setvalues("krig.dat");  // Initialize and create the random function

  cout << "Done setting values" << endl;
   /* The data will be read in the format:
   p = long integer
   n = long integer
   beta0 = double
   beta1 = p doubling-point values separated by spaces
   beta2 = p*p doubling-point values read along the rows
      Ex. the Matlab matrix [1 2; 3 4] would be input
           1 2 3 4
   x0 = p doubling-point values separated by spaces
   alpha = double                    
   theta = double                  
   sigma2 = double
   lower = p floating-point values separated by spaces
   upper = p floating-point values separated by spaces
   */

  p = MyRandomFunction.getp();  // Fetch the dimension of the space
  X.newsize(m,p);  // Make X an n by p matrix, so that each row contains one point 
  lower = MyRandomFunction.getlower();  // Fetch the lower bounds
  upper = MyRandomFunction.getupper();  // Fetch the upper bounds

  randPoints(X, lower, upper); // Randomly generate a set of initial design sites

  // Now evaluate the random function at each of these sites
  for( int i = 0; i < m; i++)
    y[i] = MyRandomFunction.evalf(X.row(i));

  cout << "creating approximation..." << endl;

  // Use these commands to fix theta in the approximation
  //double theta[2];
  //theta[0] = MyRandomFunction.gettheta();
  //theta[1] = MyRandomFunction.getalpha();
  //MyApproximation.setvalues(X,y,1,theta );  // Create the approximation

  // Otherwise, do this to estimate theta automatically
  MyApproximation.setvalues(X,y,1);  // Create the approximation

  /* Now evaluate the approximation at each point on a 100 by 100 grid in the space, and
     output these results to the file.  Results are printed in a format:

     Point Value

     So, for example, if we evaluate the function at the point x = (0.4, 0.6), and get a 
     function value y = f(x) = 23.34, then the following line would be output:

     0.4 0.6 23.34
  */

  cout << "Done creating the approximation.  Now scanning" << endl;

  scanfunction( MyApproximation, outf, 100, lower, upper );

  outf.close();
  return 0;
}  // end main


void randPoints(Matrix<double> & X, Vector<double>& lower, Vector<double>& upper)
{
  long n,p;
  double temp;
  n = X.num_rows();
  p = X.num_cols();

  if( (p != lower.dim()) || (p != upper.dim()))
    cout << "randPoints::ERROR: dimension of bounds not equal to num_cols of Matrix" << endl;
  
  SelectStream(84);
  //PutSeed(86746386); // Use this seed
  PutSeed(-1);     // Get seed from system time clock
  for( long i = 0; i < n; i++)
    for( long j = 0; j < p; j++) {
      temp = lower[j];
      X[i][j] = temp + Random()*(upper[j] - temp);
    } // end for

  return;
} // end randPoints