/* ----------------------------------------------------------------------------
orln.cc

Generates a random sample or plot points from the second order statistic in a 
sample of size r1 of a log normal random variable with location parameter r2 
and scale parameter r3.

Copyright (C) 1993, 1994, 1997, 1998, 2002, 2006.

A. Ronald Gallant
Post Office Box 659
Chapel Hill NC 27514-0659 
USA   

Permission to use, copy, modify, and distribute this software and its docu-
mentation for any purpose and without fee is hereby granted, provided that 
the above copyright notice appears in all copies and that both that copyright 
notice and this permission notice appear in supporting documentation.  This 
software is provided "as is" without any expressed or implied warranty.

-----------------------------------------------------------------------------*/

#include "libscl.h"
using namespace scl;

using namespace std;

extern REAL orden(REAL z, INTEGER n, REAL & dfwz);

char       *progname;

void       usage_msg(void);
void       process_file(istream& cfile);
void       badflag(char * text);

void badflag(char * text) {
  cerr << "Error, "<< progname <<", unknown flag or value missing: "
       << text << ".\n";
  exit(1);
}

INTEGER    sample_size;         // Set with flag -n
int        plot = 0;            // The flag -p sets to true
int        eformat = 0;         // The flag -e sets to true
char*      parm = 0;            // Set with flag -w, defaults to none.
REAL       xmin = exp(-3.0);    // Set with flag -xmin, defaults otherwise.
REAL       xmax = exp( 3.0);    // Set with flag -xmax, defaults otherwise.
INT_32BIT  seed = 349857;       // Can be changed with flag -s

void usage_msg()
{
  cerr << '\n'
    <<"usage: " << progname << " -n value -r1 value -r2 value -r3 value [options]\n\n"
    <<"flags: -n  value  number of observations or number of plot points\n"
    <<"                  the flag -m is a synonym\n"
    <<"       -r1 value  number bidders, must be larger than 3\n" 
    <<"       -r2 value  location parameter mu of lognormal valuation distribution,\n" 
    <<"       -r3 value  scale parameter sigma of lognormal valuation distribution,\n\n" 
    <<"options: -s  value     supply an integer seed rather than accept default\n"
    <<"         -p            generate plot points\n"
    <<"         -e            output observations in e format with full precision\n"
    <<"         -w  filename  write r1, r2, r3 to filename\n"
    <<"         -xmin  value  minimum plot point\n"
    <<"         -xmax  value  maximum plot point\n\n"
    <<"examples:  "<<progname<<" -n 200 -r1 4 -r2 1.0 -r3 3.0  > out.dat\n"
    <<"           The output file contains y values, one per line.\n"
    <<"           "<<progname<<" -n 200 -r1 10 -r2 1.0 -r3 3.0 -p > out.dat\n"
    <<"           The output file contains y x pairs, one pair per line, y first.\n\n"
    <<"reference: Gallant, A. R., and G. Tauchen (1999), The Relative Efficiency\n"
    <<"           of Method of Moments Estimators, J. of Econometrics.\n";
  exit(1);
}

// Template for a program that reads flags and files from the command line and 
// processes those files sequentially.  Multiple flags, grouping of flags, and 
// conditional read of next file after a flag are supported.  

int main(int argc, char *argp[], char *envp[])
{
  int      argx, flagx, filex;
  char     flag;

  sample_size = 0;

  flag = envp[0][0];  // This supresses a compiler warning.

  progname = argp[0];

  INTEGER r1=0; 
  REAL r2=REAL_MAX, r3=REAL_MIN;

  // begin: loop over flags and input files; from = file1, file2, ... 
  // If there are no files, loop once.

  argx = 0;  filex = 1;  

  while (argx < argc) {

    if (argx==0) argx++;

    // begin: read flags, if any.

    while (argx < argc && (argp[argx])[0]=='-') {

      // Read flags, if any.  Note that flags can occur singly as -a -b 
      // or they can be grouped as -ab.  Just insert code within the
      // appropriate case and don't worry about it.

      for (flagx=1; (flag=(argp[argx])[flagx])!=0; flagx++) {

        // If a flag can have more than one letter, deal with it within the 
        // appropriate case by stripping off letters and incrementing flagx.

        // If a flag indicates that argp[argx+1] is a file that should be read 
        // then, within the appropriate case, open arg[argx+1] or otherwise 
        // deal with the filename arg[argx+1], assign argp[argx] to 
        // argp[argx+1], and increment argx.  

        // If a flag indicates that argp[argx+1] is a string containing the  
        // value of a parameter, etc., then, within the appropriate case, deal 
        // with the string arg[argx+1], assign argp[argx] to argp[argx+1], and 
        // increment argx.  

        switch (flag) {

        case 'p' :
          plot = 1;
          break;

        case 'e' :
          eformat = 1;
          break;

        case 'm' :
        case 'n' :
          sample_size = (INTEGER) atoi(argp[argx+1]); 
          argp[argx+1] = argp[argx];
          if (sample_size <= 0 ) {
            cerr << "sample_size must be positive\n"; 
            badflag(&(argp[argx][flagx]));
          }
          argx++;
          break;

        case 'r' :
          if (memcmp(&(argp[argx][flagx]),"r1",2)==0) {
            r1 = (INTEGER) atoi(argp[argx+1]);
          }
          else if (memcmp(&(argp[argx][flagx]),"r2",2)==0) {
            r2 = (REAL) atof(argp[argx+1]);
          }
          else if (memcmp(&(argp[argx][flagx]),"r3",2)==0) {
            r3 = (REAL) atof(argp[argx+1]);
          }
          else {
            badflag(&(argp[argx][flagx]));
          }
          flagx++;
          argp[argx+1] = argp[argx];
          argx++;
          break;
          
        case 's' :
          seed = (INTEGER) atoi(argp[argx+1]); 
          argp[argx+1] = argp[argx];
          argx++;
          break;

        case 'w':
          parm = argp[argx+1];
          argp[argx+1] = argp[argx];
          argx++;
          break;

        case 'x' : 
          if (memcmp(&(argp[argx][flagx]),"xmin",4)==0) {
            flagx += 3;
            xmin = atof(argp[argx+1]); 
            argp[argx+1] = argp[argx];
          }
          else if (memcmp(&(argp[argx][flagx]),"xmax",4)==0) {
            flagx += 3;
            xmax = atof(argp[argx+1]); 
            argp[argx+1] = argp[argx];
          }
          else {
            badflag(&(argp[argx][flagx]));
          }
          argx++;
          break;

        default :
          cerr << "Flag was ";
          while ((flag=(argp[argx])[flagx++])!=0) 
            cerr << flag;
          cerr << "\nError: unknown flag or other error.\n";
          exit(1);
          break;

        } // end: switch

      } // end: for 

    argx++;

    } // end: while

    // end: read flags, if any.


    if (filex == 1 && argx >= argc) {

      // Put code to process cin here.  This would include the case where no  
      // files are to be read; i.e. flag processing is all that is required.

      // Example of cin processing.
      // process_file(cin);

      if ( (sample_size == 0)||(r1==0)||(r2==REAL_MAX)||(r3==REAL_MIN) ) usage_msg();
      if (r1 <= 3)  usage_msg();

      char str[80];

      if (parm) {
        ofstream parmfile(parm);
        if (!parmfile) {cerr << "Error, Cannot open " << parm << '\n'; exit(1);}
        if (eformat) {
          sprintf(str,"%30d\n",r1); parmfile << str;
          sprintf(str,"%30.18e\n",r2); parmfile << str;
          sprintf(str,"%30.18e\n",r3); parmfile << str;
        }
        else {
          cout << r1 <<'\n'<< r2 <<'\n'<< r3 << '\n';
        }
        parmfile.close();
      }

      if (plot) 
        cerr << progname << ": " << sample_size << " plot points, "
             << " r1 = " << r1 << ", r2 = " << r2 << ", r3 = " << r3 << '\n';
      else
        cerr << progname << ": " <<sample_size<<" observations, "
             << " r1 = " << r1 << ", r2 = " << r2 << ", r3 = " << r3 << '\n';

      const REAL zero = 0.0;
      const REAL one = 1.0;
      const REAL root_two = 1.414213562373095;
      const REAL root_pi  = 1.772453850905516;


      REAL    mean;
      REAL    std_err;

      mean = r2;
      std_err = r3;

      if (plot) {
        INTEGER j;
        REAL inc = (xmax-xmin)/(sample_size-1.0);
        REAL x, y, z, dummy;
        for (j=0; j<sample_size; j++) {
          x = xmin + inc*j;
          z = (log(x)-mean)/std_err;
          y = orden(z,r1,dummy)/(std_err*x);
          cout << y << "  " << x << '\n';
        }
      }
      else {
        INT_32BIT * ix;
        ix = &seed;

        REAL y;
        REAL * Z = new REAL[r1];
        if ( Z == 0) {cerr << "Error, orden, operator new failed"; return 1;}
        REAL * ztop = Z + r1;
        REAL * z;
        INTEGER j;

        for (j=0; j<sample_size; j++) {
          z = Z;

          while(z<ztop) {
            *z = unsk(ix);
            z++;
          }

          heapsort(Z,r1);

          y = exp(mean + std_err*Z[r1-2]);

          if (eformat) {
            sprintf(str,"%30.18e\n",y);
            cout << str;
          }
          else {
            cout << y << '\n';
          }
        }
      }
      return 0;
    }

    usage_msg();

    // begin: open input file, if any.

    ifstream infile(argp[argx]);
    if (!infile) {
      cerr << "Error, " << argp[0] << ", Cannot open " << argp[argx] << "\n";
      exit(1);
    } 
    // istream& from = infile;

    // end: open input file, if any.

    // Put code to process a file here.

    // Example of file processing
    // process_file(from);

    infile.close();
    filex++;
    argx++;
  } 

  // end: loop over flags and input files.

  return 0;

}

void process_file(istream& from)
{ char ch; while (from.get(ch)) cout << ch; }



/*
Function      orden - compute the density and derivative of the second highest
		      order statistic from a sample of size n from the standard 
		      normal distribution.

Syntax        REAL orden(REAL z, INTEGER n, REAL & dfwz);

Description   The return value is the value of the density at z and dfwz
	      is the derivative of the density with respect to z on return.

Return value  The value f = f(z|n).

Functions     Library: exp
called        Libcpp: pnorm
*/


REAL gaussian(REAL z, REAL & dfwz);
REAL gauscdf(REAL z);
REAL gauscdfc(REAL z);


REAL orden(REAL z, INTEGER n, REAL & dfwz)
{

// Returns the value of the density of the second largest order statistic
// in a sample of size n from the standard normal and computes its
// derivative with respect to z.

  const REAL one = 1.0;
  const REAL two = 2.0;

  REAL f;
  REAL fn, fn1, fn2;
  REAL p, pc, p2, p3;
  REAL g, g2, dgwz;

  fn  = n;
  fn1 = fn*(fn - one);
  fn2 = fn1*(fn - two);

  p = gauscdf(z);
  pc = gauscdfc(z);
  p2 = pow((double)p, (double)(n-2));
  p3 = pow((double)p, (double)(n-3));

  g  = gaussian(z, dgwz);
  g2 = pow((double)g, 2.0);

  f = fn1 * p2 * pc * g;

  dfwz  = fn2 * p3 * pc * g2;
  dfwz -= fn1 * p2 * g2;
  dfwz += fn1 * p2 * pc * dgwz;

  return f;
}


REAL gaussian(REAL z, REAL & dfwz)
{
// Returns the value of the normal density at z and computes its
// derivative with respect to z.

  const REAL two = 2.0;
  const REAL root_two_pi = 2.5066282746310005024;

  REAL f;
  REAL tmp;

  tmp = -z*z/two;

  f = (REAL) exp( (double) tmp );

  f /= root_two_pi;

  dfwz = -(z*f);

  return f;

}



REAL gauscdf(REAL z)
{
// Returns the standard normal distribution at z.

  REAL p;

  p = pnorm(z);

  return p;
}

REAL gauscdfc(REAL z)
{
// Returns the complementary standard normal distribution at z.

  REAL p;

  p = 1.0 - pnorm(z);

  return p;
}


/*
REAL gauscdf(REAL z)
{
// Returns the standard normal distribution at z.

  const REAL root_two = 1.414213562373095;
  const REAL two = 2.0;
  const REAL half = 0.5;

  REAL p;

  double zz;

  zz = (double) z/root_two;

  if (zz > 0.0) {
    p =(REAL) erf(zz); 
    p /= two;
    p += half;
  } 
  else {
    zz = fabs(zz);
    p = (REAL) erfc(zz);
    p /= two;
  }

  return p;
}

REAL gauscdfc(REAL z)
{
// Returns the complementary standard normal distribution at z.

  const REAL root_two = 1.414213562373095;
  const REAL two = 2.0;
  const REAL half = 0.5;

  REAL p;

  double zz;

  zz = (double) z/root_two;

  if (zz < 0.0) {
    zz = fabs(zz);
    p =(REAL) erf(zz); 
    p /= two;
    p += half;
  } 
  else {
    p = (REAL) erfc(zz);
    p /= two;
  }

  return p;
}
*/


