// the following program (test2.cpp) includes the test functions
// f and neighbor, plus the random-variate generation routines,
// for Nelson's initial test problem updated to the agreed standards
// from 5/21/98


#include "test.h"

//****************************************************************************

double f(long theta)
{
// this function returns a single random value from the
// second test function of Nelson at point theta.
// it calls CalcNorm, a function to generate N(0,1)
// random variates via the inverse cdf method, by
// passing a U(0,1) generated by rand()
//
       double response, dtheta1, dtheta2;
       long m;

// convert the index theta into the vector (dtheta1, dtheta2)

       m = floor(theta/21);
       dtheta1 = -5 + m;
       dtheta2 = -10 + theta - m*21;

// evaluate the simulation

       response = pow(dtheta1, 3.0) + 12.0*dtheta1*dtheta2*dtheta2 +
                  pow(dtheta2, 3.0) + 36.0*dtheta1 - 4.0*dtheta2 +
                  4*CalcNorm(rand(2));
       return response;

}// end function f


long neighbor(long theta)
{
// this routine generates a candidate solution from the neighborhood
// of theta for the feasible region THETA1 = [-5, -4, ..., 4, 5]
// and THETA2 = [-10, -9, ..., 9, 10]
// if a bad value is entered, the function returns -9999

       double newtheta1, newtheta2, dtheta1, dtheta2;
       long m;

// convert the index theta into the vector (dtheta1, dtheta2)

       m = floor(theta/21);
       dtheta1 = -5 + m;
       dtheta2 = -10 + theta - m*21;


// input error check
	 if (dtheta1>5 || dtheta1 <-5)
	 {return -9999;}
	 if (dtheta2>10 || dtheta2 <-10)
	 {return -9999;}

// generate neighbor via rejection

      newtheta1  = floor(11.0*rand(6) - 5.0);
      while (newtheta1 == dtheta1){
               newtheta1  = floor(11.0*rand(6) - 5.0);
      }

      newtheta2  = floor(21.0*rand(6) - 10.0);
      while (newtheta2 == dtheta2){
               newtheta2  = floor(11.0*rand(6) - 5.0);
      }


      return 21*(newtheta1 + 5) + newtheta2 + 10;

}//end function neighbor


long initTheta(void)
{
// this function returns a starting solution for test function 2
// the solutions is randomly sampled from the feasible region
// THETA1 = [-5, -4, ..., 4, 5] and THETA2 = [-10, -9, ..., 9, 10] 


      double newtheta1, newtheta2;

      newtheta1  = floor(11.0*rand(6) - 5.0);
      newtheta2  = floor(21.0*rand(6) - 10.0);

      return 21*(newtheta1 + 5) + newtheta2 + 10;
}



/////////////////////////////////////////////////////////////////
double CalcNorm(double U)
{	
// Approximate inverse cdf for standard normal
// Based on vnorm in Bratley, Fox, and Schrage, p.339

	double P		= 1 - U;

	const double PLIM	= 1.0e-018; 

	const double P0		= -0.322232431088;
	const double P1		= -1.0;
	const double P2		= -0.342242088547;
	const double P3		= -0.0204231210245;
	const double P4		= -0.0000453642210148;

	const double Q0		= 0.099348462606;
	const double Q1		= 0.588581570495;
	const double Q2		= 0.531103462366;
	const double Q3		= 0.10353775285;
	const double Q4		= 0.0038560700634;

if (P > .5) {P = U;} //end if

double VTemp = 8;  //Initialize

if (P >= PLIM) //make sure P isn't tiny
{
double Y = sqrt(-log(P*P));
VTemp = Y + 
       ((((Y*P4 +P3)*Y +P2)*Y +P1)*Y +P0)/((((Y*Q4 +Q3)*Y +Q2)*Y +Q1)*Y +Q0);
}//end if

else {return 1;} //if it is tiny, return 1

if (1.0 - U < .5) {VTemp = -VTemp;}

double VNorm = VTemp;
return VNorm;
}//end CalcNorm


//////////////////////////////////////////////////////////////////////////
//Random-Number Generator form Law and Kelton, p454
//Prime modulus multiplicative linear congruential generator
//Z[i]=(630360013 *Z[i-1]) (mod(pow2,31)-1), based on Marse and
//Roberts portable Fortran random-number generator UINRAN.  Multiple
//(100) streams are supported, with seeds spaced 100,000 apart.
//throughout, input argument "stream" must be an int giving the 
//desired stream number.  The header file rand.h must be included
// in the calling program (#include "rand.h") before using these 
// functions

//Usage: (Three functions)

//1. To obtain the next U(0,1) random number from stream "stream,"
//		execute
//			u = rand(stream);
//		where rand is a float function.  The float variable u will
//		contain  the next random number.
//
//2. To set the seed for stream "stream" to a desired value zset,
//		execute
//			randst(zset, stream);
//		where randst is a void function and szet must be a long set to
//		the desired seed, a number between 1 and 2147483646 (inclusive).
//		Default seeds for all 100 streams are given in the code.
//
//3. To get the current (most recently used)integer in the sequence
//	 being generated for stream "stream" into the long variable zget,
//		execute
//		  zget = randgt(stream)
//		where randgt is a long function



/* Generate the next random number. */
double rand(int stream)
{
    long zi, lowprd, hi31;

    zi     = zrng[stream];
    lowprd = (zi & 65535) * MULT1;
    hi31   = (zi >> 16) * MULT1 + (lowprd >> 16);
    zi     = ((lowprd & 65535) - MODLUS) +
             ((hi31 & 32767) << 16) + (hi31 >> 15);
    if (zi < 0) zi += MODLUS;
    lowprd = (zi & 65535) * MULT2;
    hi31   = (zi >> 16) * MULT2 + (lowprd >> 16);
    zi     = ((lowprd & 65535) - MODLUS) +
             ((hi31 & 32767) << 16) + (hi31 >> 15);
    if (zi < 0) zi += MODLUS;
    zrng[stream] = zi;
    return (( (double) (zi >> 7 | 1) + 1 )/ 16777216.0);
}

////////////////////////////////////////////////////////////////////////////////
//Set the current zrng for stream "stream" to zset
void randst (long zset, int stream)
{
	zrng[stream] = zset;
}
/////////////////////////////////////////////////////////////////////////////////


/////////////////////////////////////////////////////////////////////////////////
//Return the current zrng for stream "stream".

long randgt (int stream)
{
	return zrng[stream];
}
//////////////////////////////////////////////////////////////////////////////////