/* Copyright (C) 1999-2000 Simon N. Wood  snw@st-and.ac.uk

This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License   
as published by the Free Software Foundation; either version 2
of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
(www.gnu.org/copyleft/gpl.html)

You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307,
USA.*/

/* This file contains univariate and multivariate random number generating
   routines, as well as routines for shuffling data.
*/

#include <time.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include "matrix.h"
#include "rangen.h"


#define r0a 16807
#define r0m 2147483647
#define r0q 127773
#define r0r 2836

#define r2m1 2147483563
#define r2m2 2147483399
#define r2q1 53668
#define r2q2 52774
#define r2r1 12211
#define r2r2 3791
#define r2a1 40014
#define r2a2 40692
#define r2nra 32
#define r2ndiv (1+(r2m1-1)/r2nra)

void ErrorMessage(char *msg,int fatal);

double ran0(long *iran)

/* Implementation of random number generator suggested by Press et al. as
   their ran0(). Code written independently of NR routine. iran should be 
   set to zero to initialize from the clock, otherwise give a seed value 
   to initialize. DO NOT alter iran between calls except to re-initialize.
   See N.R. equations 7.1.2-7.1.5 for details. Tested against NR routine
   - results identical for identical seed. Period is 2147483645 (i.e.
   about 2e9), but there are problems - again see Numerical Recipes. 
   Note that the output is really a float.

*/

{ long i;
  i=*iran; // more efficient
  if (!i)
  { i=time(NULL)%32000;}
  if (i<0) i= -i;
  i = (i%r0q) * r0a - r0r*(i/r0q);
  if (i<0) i += r0m;  
  *iran=i;
  return(i*1.0/r0m);
}

double rn0()

/* calls ran0 but always initializes with clock */

{ static long iran=0;
  return(ran0(&iran));
}


double ran2(long *iran)


/* Implementation of the random number generator used for ran2() in Numerical Recipes,
   but re-written from scratch. The generator works as follows:
   2 generators of the basic form given in ran0() are used (with different constants,
   of course), G1 and G2, say. An array ra[] is also used to store 32 random numbers 
   from G1 at any one time - this array is used to shuffle the output from G1 to break 
   up serial correlations. 
   1. On start up ra[] is filled from G1 (after some "warm-up" calls). The seed is either 
      -(*iran) (if *iran supplied -ve) or a value set from the clock if *iran set to zero.
      iro is set to ra[0].
   2. One step of G2 (using iran2) is taken. 
   3. iro is used to select an element of ra[]. iro is then recalculated using this value and
      iran2.
   4. One step of G1 (using *iran) is taken, *iran replaces the element of ra[], just used.
   5. output iro cast to double as result.
   Unless re-seeding DO NOT alter *iran between calls.
   If magnitude of a -ve *iran will be used to re-seed.
   *iran set to zero will re-seed from the clock.
   NOTE: numbers are really floats. A different routine will be needed if this matters!!
 */


{ static long ra[r2nra],iran2=123456789,iro;
  long iran1;
  int i;
  double ro;
  iran1=*iran;
  if (iran1<=0) /* initialize */
  { if (!iran1) iran1=time(NULL)%32000; // use clock
    else iran1 = -iran1;
    for (i=20+r2nra;i>=0;i--) 
    { iran1 = (iran1%r2q1) * r2a1 - r2r1*(iran1/r2q1);if (i<0) iran1 += r2m1;  // step of G1
      if (i<r2nra) ra[i]=iran1;
    } 
    iro=ra[0];
  }
  iran2 = (iran2%r2q2) * r2a2 - r2r2*(iran2/r2q2);if (i<0) iran2 += r2m2;  // step of G2
  i=iro/r2ndiv; // selects element of ra[]
  if (i>=r2nra)
  i=r2nra-1;
  iro=ra[i]-iran2;
  if (iro<1) iro+=r2m1-1;
  ro=iro*1.0/r2m1;
  if (ro==1.0) ro += -1.2e-14; /* avoids returning 1.0 */
  iran1 = (iran1%r2q1) * r2a1 - r2r1*(iran1/r2q1);if (i<0) iran1 += r2m1;  // step of G1
  ra[i]=iran1;
  *iran=iran1;
  return(ro);
}

double ran1(void)

/* used purely to call ran2() seeded from time only */

{ static long iran=0;
  return(ran2(&iran));
} 










int ranint(int a, int b)

/* produce uniform random integer between a and b */

{ double x;
  int i;
  x=a+(b-a+1)*ran1();
  i=(int)floor(x);
  if (i>b) i--;
  return(i);
}



double gammln(double z)

/* natural log of gamma function routine using 6.1.5. from Press et al. (but not their code) */

{ double zz,x,sum;
  static double c[6]={76.18009172947146,-86.50532032941677,24.01409824083091,
			-1.231739572450155,0.1208650973866179e-2,-0.5395239384953e-5},r2pi;
  int i;
  static int first=1;
  if (first)
  { first=0;
    r2pi=sqrt(asin(1.0)*4);
  }
  zz=z-1.0;
  x=zz+5.5;
  x=log(x)*(zz+0.5)-x;
  sum=1.000000000190015;
  for (i=0;i<6;i++)
  { zz+=1.0;
    sum+=c[i]/zz;
  }
  return(x+log(r2pi*sum));
}




double ndev(void)

/* Box-Muller for producing N(0,1) random deviates, modified slightly to avoid 
   trig. function calls. Limited testing using Minitab normal scores plots. */

{ double r2=2.0,y1;
  static double y2;
  static int leftover=0;
  if (!leftover) // need to generate 
  while (r2>=1||r2==0.0)
  { y1=2*ran1()-1;
    y2=2*ran1()-1;
    r2=y1*y1+y2*y2;
  } else { leftover=0;return(y2);}
  r2=sqrt(-2*log(r2)/r2);
  y1*=r2;y2*=r2;
  leftover=1;
  return(y1);
}

  



void mvn(matrix *C,matrix *m,matrix *x,int control)

/* Routine for generating  multivariate normal random deviates x ~ N(m,C).
   control = 0 for generation.
           = 1 for setting up with new C and/or m (uses memory)
           = 2 for freeing space

   The method works like this:
   1. Generate vector e ~ N(0,I)
   2. Obtain the choleski decomposition of C = LL'
   3. Form x=Le. Co(x) = LCo(e)L'=LIL'=LL'=C.
   4. Form x=x+m. x ~ N(m,C).
*/

{ static matrix L,e;
  static int first=1;
  static long er;
  static double **LM,*eV;
  int i,j;
  if (first)
  { first=0;L.r=0L;}
  if ((control==1)||(L.r==0L))
  { if (L.r) { freemat(L);freemat(e);}
    L=initmat(C->r,C->r);e=initmat(m->r,1L);
    LM=L.M;eV=e.V;er=e.r;
    choleski(*C,L,0,0);
  } else
  if (control==2)
  { freemat(L);freemat(e);}
  for (i=0;i<er;i++) { eV[i]=ndev();x->V[i]=m->V[i];}
  for (i=0;i<er;i++)
  for (j=0;j<=i;j++) x->V[i]+=LM[i][j]*eV[j];
}

void resample(double *in,double *out,int n,int m,int replace)

/* Routine that draws a random sample of size m from array n without replacement
   if replace==0, with replacement otherwise. in[] remains unchanged by the
   operation. */

{ double *t,temp;
  int i,j,max;
  if (replace)
  { for (i=0;i<m;i++)
    { j=(int)floor(n*ran1());if (j>=n) j--;
      out[i]=in[j];
    }
  } else
  { if (m==n) t=out; /* can shuffle in situ */
    else if (m<n)
    t=(double *)calloc((size_t)n,sizeof(double));
    else ErrorMessage("You are trying to resample without replacement, more items than are available.",1);
    for (i=0;i<n;i++) t[i]=in[i]; /* copying data to be shuffled */
    /* now do the shuffling */
    max=n;
    for (i=0;i<m;i++)
    { j=(int)floor((max)*ran1());if (j>=max) j--;
      j+=i;temp=t[i];t[i]=t[j];t[j]=temp;
      max--;
    } /* first m elements of t now contain the resampled data. (Note that
         the last n-m elements are not randomly ordered) */
    if (m<n) /* then t must be copied to output */
    { for (i=0;i<m;i++) out[i]=t[i];free(t);}
  }
}

/*************************************************************************************************/
/* Notes on updates and fixes.                                                                   */
/*************************************************************************************************/

/*
1. Uniform random number generators have been re-written, to avoid use of Press et al. 
   copyright routines. 17/1/2000. Generators used by Press et al. as ran0() and ran2()
   implemented, but not ran1(). ran0() tested against Press et al. code.

2. Normal random deviate routine and log gamma function routines re-written to avoid Press et al.
   code. 17/1/2000. Both tested.
*/