head	1.5;access;symbols;locks	rirwin:1.5; strict;comment	@# @;1.5date	2005.06.22.18.15.00;	author rirwin;	state Exp;branches;next	1.4;1.4date	2005.05.23.20.39.51;	author rirwin;	state Exp;branches;next	1.3;1.3date	2005.03.09.03.43.18;	author patil;	state Exp;branches;next	1.2;1.2date	2005.01.20.02.41.42;	author patil;	state Exp;branches;next	1.1;1.1date	2004.12.28.00.04.32;	author patil;	state Exp;branches;next	;desc@@1.5log@Add javadoc comment to class declaration.@text@/** file: BiNormal.java** Last edited: Ryan Irwin**//*// These imports are not needed - Phil T. 6-23-03// import java.awt.*;// import java.applet.*;// import javax.swing.*;// import java.lang.*;// import java.io.*;*/import java.util.*;/** * Handles random generation of Gaussian uses */public class BiNormal {        // declare global variables    //    int snormi = 0;    double snorm = 0.0;    double snormu = 0.0;    double snorms = 0.0;    double ustar = 0.0;    double snormaa = 0.0;    double snormw = 0.0;    double snormy = 0.0;    double snormtt = 0.0;    int mltmod = 0;    int mltmoda0 = 0;    int mltmoda1 = 0;    int mltmodk = 0;    int mltmodp = 0;    int mltmodq = 0;    int mltmodqh = 0;    int mltmodrh = 0;    int mltmodh = 32768;        int Xm1 = 0;    int Xm2 = 0;    int Xa1 = 0;    int Xa2 = 0;    int Xa1w = 0;    int Xa2w = 0;    int Xa1vw = 0;    int Xa2vw = 0;    int Xcg1[] = new int[32];    int Xcg2[] = new int[32];    int Xig1[] = new int[32];    int Xig2[] = new int[32];    int Xlg1[] = new int[32];    int Xlg2[] = new int[32];    int Xqanti[] = new int[32];    Integer info = new Integer(0);    Integer qinit = new Integer(0);    Integer curntg = new Integer(1);    Integer qstate = new Integer(0);    Random random = new Random();     double[] a = {	0.0,3.917609E-2,7.841241E-2,0.11777,0.1573107,0.1970991,0.2372021,	0.2776904,0.3186394,0.36013,0.4022501,0.4450965,0.4887764,0.5334097,	0.5791322,0.626099,0.6744898,0.7245144,0.7764218,0.8305109,0.8871466,	0.9467818,1.00999,1.077516,1.150349,1.229859,1.318011,1.417797,	1.534121,1.67594,1.862732,2.153875    };     double[] d = {	0.0,0.0,0.0,0.0,0.0,0.2636843,0.2425085,0.2255674,0.2116342,0.1999243,	0.1899108,0.1812252,0.1736014,0.1668419,0.1607967,0.1553497,0.1504094,	0.1459026,0.14177,0.1379632,0.1344418,0.1311722,0.128126,0.1252791,	0.1226109,0.1201036,0.1177417,0.1155119,0.1134023,0.1114027,0.1095039    };     double[] t = {	7.673828e-4,2.30687e-3,3.860618e-3,5.438454e-3,7.0507e-3,8.708396e-3,	1.042357e-2,1.220953e-2,1.408125e-2,1.605579e-2,1.81529e-2,2.039573e-2,	2.281177e-2,2.543407e-2,2.830296e-2,3.146822e-2,3.499233e-2,	3.895483e-2,4.345878e-2,4.864035e-2,5.468334e-2,6.184222e-2,	7.047983e-2,8.113195e-2,9.462444e-2,0.1123001,0.136498,0.1716886,	0.2276241,0.330498,0.5847031    };     double[] h = {	3.920617e-2,3.932705e-2,3.951e-2,3.975703e-2,4.007093e-2,4.045533e-2,	4.091481e-2,4.145507e-2,4.208311e-2,4.280748e-2,4.363863e-2,	4.458932e-2,4.567523e-2,4.691571e-2,4.833487e-2,4.996298e-2,	5.183859e-2,5.401138e-2,5.654656e-2,5.95313e-2,6.308489e-2,	6.737503e-2,7.264544e-2,7.926471e-2,8.781922e-2,9.930398e-2,	0.11556,0.1404344,0.1836142,0.2790016,0.7010474    };    /**     * methods sets up the parameters needed to generate the multivariate      * normal deviates form the inputs given     *     *     * @@param  meanv mean vector of multivariate normal distribution     * @@param  covm covariance matrix of multivariate      *         normal distribution     * @@param  p dimension of the normal     * @@param  parm array of parameters needed to generate      *                   multivariate normal deviates     *     */    public  void setgmn(double[] meanv, double[] covm, int p, 			      double[] parm)     {	int i = 0;	int icount = 0;	int j = 0;	int D2 = 0;	int D3 = 0;	int D4 = 0;	int D5 = 0;	// test the input	//	if (p <= 0) {	    // System.out.println("P nonpositive in SETGMN");	    return;	}	parm[0] = p;	// put p and meanv into param	//	for(i=2, D2=1, D3=(p+1-i+D2)/D2; D3>0; D3--, i+=D2) {	    parm[i-1] = meanv[i-2];	}	// Cholesky decomposition to find A s.t. trans(A)*(A) = COVM	//	spofa(covm, p, p);	if(info.intValue() != 0) {	    // System.out.println(" COVM not positive definite in SETGMN");	    return;	}	icount = p+1;	// put upper half of a, which is now the cholesky factor, into parm	// covm(1,1) = parm(p+2)	// covm(1,2) = parm(p+3)	// :	// covm(1,p) = parm(2p+1)	// covm(2,2) = parm(2p+2)  ...	//	for(i=1, D4=1, D5=(p-i+D4)/D4; D5>0; D5--, i+=D4) {	    for(j=i-1; j<p; j++) {		icount += 1;		parm[icount-1] = covm[i-1+j*p];	    }	}    }        /**     * methods generates the multivariate normal deviates using the procedure:     * 1) generate p independent standard normal deviates - ei ~ n(0,1)     * 2) using cholesky decomposition find a s.t. trans(a)*a = covm     * 3) trans(a)e + meanv ~ n(meanv,covm)     *     * @@param parm array of parameters needed to generate      *             multivariate normal deviates     * @@param x vector deviate generated     * @@param work scratch array     *     */    public  void genmn(double[] parm, double[] x, double[] work)     {	int i = 0;	int j = 0;	int p = 0;	int D1 = 0;	int D2 = 0;	int D3 = 0;	int D4 = 0;	int icount = 0;	double ae = 0.0;		p = (int)parm[0];		//  generate p independent normal deviates - work ~ n(0,1)	//	for(i=1; i<=p; i++) 	{	    work[i-1] = snorm();	}	for(i=1, D3=1, D4=(p-i+D3)/D3; D4>0; D4--, i+=D3) 	{	    // parm (p+2 : p*(p+3)/2 + 1) contains a, the cholesky	    // decomposition of the desired covariance matrix.	    // trans(a)(1,1) = parm(p+2)	    // trans(a)(2,1) = parm(p+3)	    // trans(a)(2,2) = parm(p+2+p)	    // trans(a)(3,1) = parm(p+4)	    // trans(a)(3,2) = parm(p+3+p)	    // trans(a)(3,3) = parm(p+2-1+2p)  ...	    // trans(a)*work + meanv ~ n(meanv,covm)	    //	    icount = 0;	    ae = 0.0;	    for(j=1, D1=1, D2=(i-j+D1)/D1; D2>0; D2--, j+=D1) {		icount += (j-1);		ae += (parm[i+(j-1)*p-icount+p]*work[j-1]);	    }	    x[i-1] = ae+parm[i];	}    }    /**     * generates a uniform distribution over 0 - 1     *     * @@return  random floating point number from a uniform distribution     *            over 0 - 1 using the current generator     *     */    public double ranf()     {	return random.nextDouble();    }        /**     *     * linpack.  this version dated 08/14/78     * cleve moler, university of new mexico, argonne national lab     *     */    public void snorm40() {	if(ustar <= t[snormi-1]) 	{	    snorm60();	    return;	}	snormw = (ustar-t[snormi-1])*h[snormi-1];	snorm50();    }    /**     *     * linpack.  this version dated 08/14/78     * cleve moler, university of new mexico, argonne national lab     *     */    public void snorm50()     {	snormy = snormaa+snormw;	snorm = snormy;	if(snorms == 1.0) snorm = -snormy;    }    /**     *     * linpack.  this version dated 08/14/78     * cleve moler, university of new mexico, argonne national lab     *     */    public void snorm60()     {	snormu = ranf();	snormw = snormu*(a[snormi]-snormaa);	snormtt = (0.5*snormw+snormaa)*snormw;	snorm80();    }    /**     *     * linpack.  this version dated 08/14/78     * cleve moler, university of new mexico, argonne national lab     *     */    public void snorm70()     {	snormtt = snormu;	ustar = ranf();	snorm80();    }    /**     *     * linpack.  this version dated 08/14/78     * cleve moler, university of new mexico, argonne national lab     *     */    public void snorm80()     {	if(ustar > snormtt) 	{	    snorm50();	    return;	}	snormu = ranf();	if(ustar >= snormu) 	{	    snorm70();	    return;	}	ustar = ranf();	snorm40();    }    /**     *     * linpack.  this version dated 08/14/78     * cleve moler, university of new mexico, argonne national lab     *     */    public void snorm100()     {	snormi = 6;	snormaa = a[31];	snorm120();    }    /**     *     * linpack.  this version dated 08/14/78     * cleve moler, university of new mexico, argonne national lab     *     */    public void snorm110()     {	snormaa += d[snormi-1];	snormi += 1;	snorm120();    }    /**     *     * linpack.  this version dated 08/14/78     * cleve moler, university of new mexico, argonne national lab     *     */    public void snorm120()     {	snormu += snormu;	if(snormu < 1.0) 	{	    snorm110();	    return;	}	snormu -= 1.0;	snorm140();    }        /**     *     * linpack.  this version dated 08/14/78     * cleve moler, university of new mexico, argonne national lab     *     */    public void snorm140()     {	snormw = snormu*d[snormi-1];	snormtt = (0.5*snormw+snormaa)*snormw;	snorm160();    }    /**     *     * linpack.  this version dated 08/14/78     * cleve moler, university of new mexico, argonne national lab     *     */    public void snorm150()     {	snormtt = snormu;	snorm160();    }    /**     *     * linpack.  this version dated 08/14/78     * cleve moler, university of new mexico, argonne national lab     *     */    public void snorm160()     {	ustar = ranf();	if(ustar > snormtt) 	{	    snorm50();	    return;	}	snormu = ranf();	if(ustar >= snormu) 	{	    snorm150();	    return;	}	snormu = ranf();	snorm140();    }        /**     * ahrens, j.h. and dieter, u.                                 * extensions of forsythe's method for random                  * sampling from the normal distribution.                      * math. comput., 27,124 (oct. 1973), 927 - 937.         *     * @@return  standard normal distribution     *     *     */    public double snorm()     {	snormu = ranf();	snorms = 0.0;		if(snormu > 0.5) 	{	    snorms = 1.0;	}		snormu += (snormu-snorms);	snormu = 32.0*snormu;	snormi = (int)snormu;		if(snormi == 32) 	{	    snormi = 31;	}	if(snormi == 0) 	{	    snorm100();	    return snorm;	}		// start center	//	ustar = snormu-(double)snormi;	snormaa = a[snormi-1];		// start center	//	snorm40();		return snorm;    }    /**     *     * linpack.  this version dated 08/14/78     * cleve moler, university of new mexico, argonne national lab     *     * @@param n    integer     * @@param sx   array of doubles     * @@param dx   integer     * @@param incx integer     * @@param sy   array of doubles     * @@param dy   integer     * @@param incy integer     *     * @@return double value of sdot     */    public  double sdot(int n, double[] sx, int dx, int incx, 			      double[] sy, int dy, int incy)     {	int i = 0;	int ix = 0;	int iy = 0;	int m = 0;	int mp1 = 0;	double sdot = 0.0;	double stemp = 0.0;	stemp = sdot = 0.0;	if(n <= 0) return sdot;	if(incx != 1 && incy != 1) 	{	    	    ix = iy = 1;	    if(incx < 0) ix = (-n+1) * incx + 1;	    if(incy < 0) iy = (-n+1) * incy + 1;	    	    for(i=1; i<=n; i++) 	    {		stemp += (sx[dx+ix-1] * sy[dy+iy-1]);		ix += incx;		iy += incy;	    }	    	    sdot = stemp;	    return sdot;	}		m = n % 5;		if(m == 0) 	{	    mp1 = m+1;	    for(i=mp1; i<=n; i+=5) 	    {		stemp += (sx[dx+i-1] * sy[dy+i-1] +			  sx[dx+i]   * sy[dy+i]   +			  sx[dx+i+1] * sy[dy+i+1] +			  sx[dx+i+2] * sy[dy+i+2] +			  sx[dx+i+3] * sy[dy+i+3]);	    }	    sdot = stemp;	    return sdot;	} 	else 	{	    	    for(i=0; i<m; i++) 	    {		stemp += (sx[dx+i] * sy[dy+i]);	    }	    if(n < 5) 	    {		sdot = stemp;		return sdot;	    } 	    else 	    {		mp1 = m+1;		for(i=mp1; i<=n; i+=5) 		{		    stemp += (sx[dx+i-1] * sy[dy+i-1] +			      sx[dx+i]   * sy[dy+i]   +			      sx[dx+i+1] * sy[dy+i+1] +			      sx[dx+i+2] * sy[dy+i+2] +			      sx[dx+i+3] * sy[dy+i+3]);		}		sdot = stemp;		return sdot;	    }	}    }        /**     * linpack.  this version dated 08/14/78     * cleve moler, university of new mexico, argonne national lab     *     * @@param a array of doubles     * @@param lda integer     * @@param n integer     */    public  void spofa(double[] a, int lda, int n)     {    		int j = 0;	int jm1 = 0;	int k = 0;	double t = 0.0;	double s = 0.0;	int flag = 0;    S40:	for(j=1; j<=n; j++) 	{	    	    info = new Integer(j);	    s = 0.0;	    jm1 = j-1;	    	    if(jm1 >= 1)