/* ----------------------------------------------------------------------------
mixture.cc

Generates a random sample or plot points from the Marron-Wand test suite.
Self documented, execute mixture with no arguments for documentation.

Reference: Marron, J. S., and M. P. Wand (1992) "Exact Mean Integrated 
Squared Error," The Annals of Statistics 20, 712--736.

See Table 1 for parameters and Figure 1 for plots.

Copyright (C) 1993, 1994, 1997, 2001, 2002.

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;

const int  max_name = 256;
const int  max_mix  = 25;

char       *progname;

void       usage_msg(void);
void       process_file(istream& cfile);
REAL       normal(REAL z, int isw);

INTEGER    sample_size;         // Set with flag -n
char       density[max_name];   // Set with flag -d
int        plot = 0;            // The flag -p sets to true
int        eformat = 0;         // The flag -e sets to true
INT_32BIT  seed = 349857;       // Can be changed with flag -s
char*      sf = 0;              // Set with flag -sf, defaults to none.

void usage_msg()
{
  cerr << '\n'
    <<"usage: " << progname << " -n value -d name [options] \n"
    <<"flags: -n value number of observations or number of plot points\n"
    <<"                the flag -m is a synonym\n"
    <<"       -d name  density to sample or plot, choices: gaussian,\n"
    <<"                skewed_unimodal, strongly_skewed, kurtotic_unimodal,\n" 
    <<"                outlier, bimodal, separated_bimodal, skewed_bimodal,\n" 
    <<"                trimodal, claw, double_claw, asymmetric_claw,\n"
    <<"                asymmetric_double_claw, smooth_comb, discrete_comb\n"
    <<"options: -s  value supply an integer seed rather than accept default\n"
    <<"         -sf name  read seed from and write next seed to file name\n"
    <<"         -p        generate plot points\n"
    <<"         -e        output observations in e format with full precision\n"
    <<"examples: "<<progname<<" -n 200 -d claw > out.dat\n"
    <<"             The output file contains y values, one per line.\n"
    <<"          "<<progname<<" -n 200 -d claw -p > out.dat\n"
    <<"             The output file contains y x pairs, one pair per line,\n"
    <<"             y first.\n"
    <<"reference: Marron, J. S., and M. P. Wand (1992) Exact Mean Integrated,\n" 
    <<"           Squared Error, The Annals of Statistics 20, 712--736.\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     ch, flag;

  sample_size = 0;
  density[0] = '\0';

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

  progname = argp[0];

  // 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"; 
            goto badflag;
          }
          argx++;
          break;

        case 's' :
	  if (memcmp(&(argp[argx][flagx]),"sf",2)==0) { 
	    flagx += 1;
	    sf = argp[argx+1];
	    argp[argx+1] = argp[argx];
	    argx++;
	    break; 
          }
	  else { 
            seed = (INTEGER) atoi(argp[argx+1]); 
            argp[argx+1] = argp[argx];
            argx++;
            break;
          }

        case 'd' : 
          strncpy(density,argp[argx+1],max_name); 
          density[max_name-1] = '\0';
          argp[argx+1] = argp[argx];
          argx++;
          break;

        badflag :
        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) || (density[0] == '\0') ) usage_msg();

      if (plot) 
        cerr<<progname<< ": "
          <<sample_size<<" plot points from density "<<density<<"\n";
      else
        cerr<<progname<< ": "
          <<sample_size<<" observations from density "<<density<<"\n";

      INTEGER mix;
      REAL    mean[max_mix];
      REAL    std_err[max_mix];
      REAL    prob[max_mix];

      INTEGER i,k;
      INTEGER j=0;   //Initialization stops a compiler warning;
      REAL tmp;

      if (strcmp(density,"gaussian")==0) {
        mix=1; 
        prob[0]=1.0;        mean[0]=0.0;         std_err[0]=1.0;
      }
      else if (strcmp(density,"skewed_unimodal")==0) {
        mix=3;
        prob[0]=1.0/5.0;    mean[0]=0.0;         std_err[0]=1.0;
        prob[1]=1.0/5.0;    mean[1]=1.0/2.0;     std_err[1]=2.0/3.0;
        prob[2]=3.0/5.0;    mean[2]=13.0/12.0;   std_err[2]=5.0/9.0;
      }
      else if (strcmp(density,"strongly_skewed")==0) {
        mix=8; 
        for (i=0, tmp=1.0; i<mix; i++) {
          prob[i]=1.0/8.0;  
          mean[i]=3.0*(tmp-1.0);         
          std_err[i]=tmp;
          tmp *= 2.0/3.0;
        }
      }
      else if (strcmp(density,"kurtotic_unimodal")==0) {
        mix=2;
        prob[0]=2.0/3.0;    mean[0]=0.0;         std_err[0]=1.0;
        prob[1]=1.0/3.0;    mean[1]=0.0;         std_err[1]=1.0/10.0;
      }
      else if (strcmp(density,"outlier")==0) {
        mix=2;
        prob[0]=1.0/10.0;   mean[0]=0.0;         std_err[0]=1.0;
        prob[1]=9.0/10.0;   mean[1]=0.0;         std_err[1]=1.0/10.0;
      }
      else if (strcmp(density,"bimodal")==0) {
        mix=2;
        prob[0]=1.0/2.0;    mean[0]=-1.0;        std_err[0]=2.0/3.0;
        prob[1]=1.0/2.0;    mean[1]=1.0;         std_err[1]=2.0/3.0;
      }
      else if (strcmp(density,"separated_bimodal")==0) {
        mix=2;
        prob[0]=1.0/2.0;    mean[0]=-3.0/2.0;    std_err[0]=1.0/2.0;
        prob[1]=1.0/2.0;    mean[1]=3.0/2.0;     std_err[1]=1.0/2.0;
      }
      else if (strcmp(density,"skewed_bimodal")==0) {
        mix=2;
        prob[0]=3.0/4.0;    mean[0]=0.0;         std_err[0]=1.0;
        prob[1]=1.0/4.0;    mean[1]=3.0/2.0;     std_err[1]=1.0/3.0;
      }
      else if (strcmp(density,"trimodal")==0) {
        mix=3;
        prob[0]=9.0/20.0;   mean[0]=-6.0/5.0;    std_err[0]=3.0/5.0;
        prob[1]=9.0/20.0;   mean[1]=6.0/5.0;     std_err[1]=3.0/5.0;
        prob[2]=2.0/20.0;   mean[2]=0.0;         std_err[2]=1.0/4.0;
      }
      else if (strcmp(density,"claw")==0) {
        mix=6; 
        prob[0]=1.0/2.0;    mean[0]=0.0;         std_err[0]=1.0;
        for (i=1; i<mix; i++) {
          tmp=i-1;
          prob[i]=1.0/10.0;  
          mean[i]=tmp/2.0 - 1.0;         
          std_err[i]=1.0/10.0;
        }
      }
      else if (strcmp(density,"double_claw")==0) {
        mix=9; 
        prob[0]=49.0/100.0; mean[0]=-1.0;        std_err[0]=2.0/3.0;
        prob[1]=49.0/100.0; mean[1]=1.0;         std_err[1]=2.0/3.0;
        for (i=2; i<mix; i++) {
          tmp=i-2;
          prob[i]=1.0/350.0;  
          mean[i]=(tmp - 3.0)/2.0;         
          std_err[i]=1.0/100.0;
        }
      }
      else if (strcmp(density,"asymmetric_claw")==0) {
        mix=6; 
        prob[0]=1.0/2.0;    mean[0]=0.0;         std_err[0]=1.0;
        for (i=1; i<mix; i++) {
          tmp=i-3;
          prob[i]=pow(2.0,1.0-tmp)/31.0;
          mean[i]=tmp+1.0/2.0;
          std_err[i]=pow(2.0,-tmp)/10.0;
        }
      }
      else if (strcmp(density,"asymmetric_double_claw")==0) {
        mix=8; 
        for (i=0; i<=1; i++) {
          tmp=i;
          prob[i]=46.0/100.0;
          mean[i]=2.0*tmp - 1.0;
          std_err[i]=2.0/3.0;
        }
        for (i=2; i<=4; i++) {
          tmp=i-1;
          prob[i]=1.0/300.0;
          mean[i]=-tmp/2.0;
          std_err[i]=1.0/100.0;
        }
        for (i=5; i<mix; i++) {
          tmp=i-4;
          prob[i]=7.0/300.0;
          mean[i]=tmp/2.0;
          std_err[i]=7.0/100.0;
        }
      }
      else if (strcmp(density,"smooth_comb")==0) {
        mix=6; 
        for (i=0; i<mix; i++) {
          tmp=i;
          prob[i]=pow(2.0,5.0-tmp)/63.0;
          mean[i]=(65.0-96.0*pow(1.0/2.0,tmp))/21.0;
          std_err[i]=(32.0/63.0)/pow(2.0,tmp);
        }
      }
      else if (strcmp(density,"discrete_comb")==0) {
        mix=6; 
        for (i=0; i<=2; i++) {
          tmp=i;
          prob[i]=2.0/7.0;
          mean[i]=(12.0*tmp-15.0)/7.0;
          std_err[i]=2.0/7.0;
        }
        for (i=3; i<mix; i++) {
          tmp=i+5;
          prob[i]=1.0/21.0;
          mean[i]=2.0*tmp/7.0;
          std_err[i]=1.0/21.0;
        }
      }
      else {
        cerr<<'\n'<<progname<< ": "
          << density << " is not one of the implemented choices\n\n"
          <<"choices: gaussian, skewed_unimodal, strongly_skewed,\n" 
          <<"         kurtotic_unimodal, outlier, bimodal,\n" 
          <<"         separated_bimodal,skewed_bimodal, trimodal,\n"
          <<"         claw, double_claw, asymmetric_claw,\n"
          <<"         asymmetric_double_claw, smooth_comb, discrete_comb\n";

        exit(1);
      }


      if (plot) {
        REAL low = -3.0;
        REAL high = 3.0;
        REAL inc = (high-low)/(sample_size-1.0);
        REAL x, y, z;
        char str[80];
	
        for (j=0; j<sample_size; j++) {
          x = low + inc*j;
          y = 0.0;
          for (i=0; i<mix; i++) {
            z = (x - mean[i])/std_err[i];
            y += prob[i]*normal(z,1)/std_err[i];
          }
	  if (eformat) { 
	    sprintf(str," %27.17e %27.17e\n",y,x);
	    cout << str;
          }
	  else {
            cout << y << "  " << x << '\n';
	  }
        }
      }
      else {
        if (sf) {
	  ifstream seedfile(sf);
	  if (!seedfile) {
	    cerr << "Error, " << progname << ", Cannot open " << sf << "\n";
	    exit(1); 
          }
	  seedfile >> seed;
	}

        INT_32BIT * ix;
        ix = &seed;

        REAL low, high;
        REAL x, y, z;
        char str[80];

        for (k=0; k<sample_size; k++) {
          x = ran(ix);
          for (low=0.0, high=0.0, i=0; i<mix; i++) { 
            high += prob[i];
            if ( (low<=x) && (x<=high) ) { j=i; break; }
            low += prob[i];
          }

          z = unsk(ix);
          y = mean[j] + std_err[j]*z;

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

	if (sf) {
	  ofstream seedfile(sf); 
	  if (!seedfile) { 
	    cerr << "Error, " << progname << ", Cannot open " << sf << "\n"; 
	    exit(1);
	  }
	  seedfile << *ix << '\n';
        }

      }
      exit(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.

  exit(0);

}

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

REAL normal(REAL z, int isw)
// Returns the value of the normal kernel at z if isw=1 and its 
// convolution with itself if isw=2.
{
  REAL const two = 2.0;
  REAL const root_two_pi = 2.5066282746310005024;
  REAL const root_two    = 1.4142135623730950488;
  REAL       h;
  REAL       tmp;

  switch (isw) {

    case 1 :
     tmp = -z*z/two;
     h = (REAL) exp( (double) tmp );
     h /= root_two_pi;
     return h;

    case 2 :
     tmp = -z*z/two;
     tmp /= two;
     h = (REAL) exp( (double) tmp );
     h /= root_two_pi;
     h /= root_two;
     return h;

    default :
      cerr << "\nError: bad data to REAL normal(REAL z, int isw)\n";
      exit(1);
   }

   // The next two statemants can't be reached but stop a compiler warning.
   cerr << "\nError: bad data to REAL normal(REAL z, int isw)\n";
   return 1;  
}
