pvmica.tar

dsum.c100644   3665     17        3304  6571322457  10642 0ustar  enghoffcnl
/*  -- translated by f2c (version 19940927).
   You must link the resulting object file with the libraries:
	-lf2c -lm   (in that order)
*/

#include "f2c.h"

doublereal dsum_(integer *n, doublereal *dx, integer *incx)
{


    /* System generated locals */
    integer i__1, i__2;
    doublereal ret_val, d__1, d__2, d__3, d__4, d__5, d__6;

    /* Local variables */
    static integer i, m;
    static doublereal dtemp;
    static integer nincx, mp1;


/*     takes the sum of the values.   
       jack dongarra, linpack, 3/11/78.   
       modified 3/93 to return if incx .le. 0.   
       modified 12/3/93, array(1) declarations changed to array(*)   


    
   Parameter adjustments   
       Function Body */
#define DX(I) dx[(I)-1]


    ret_val = 0.;
    dtemp = 0.;
    if (*n <= 0 || *incx <= 0) {
	return ret_val;
    }
    if (*incx == 1) {
	goto L20;
    }

/*        code for increment not equal to 1 */

    nincx = *n * *incx;
    i__1 = nincx;
    i__2 = *incx;
    for (i = 1; *incx < 0 ? i >= nincx : i <= nincx; i += *incx) {
	dtemp += (d__1 = DX(i), d__1);
/* L10: */
    }
    ret_val = dtemp;
    return ret_val;

/*        code for increment equal to 1   


          clean-up loop */

L20:
    m = *n % 6;
    if (m == 0) {
	goto L40;
    }
    i__2 = m;
    for (i = 1; i <= m; ++i) {
	dtemp += (d__1 = DX(i), d__1);
/* L30: */
    }
    if (*n < 6) {
	goto L60;
    }
L40:
    mp1 = m + 1;
    i__2 = *n;
    for (i = mp1; i <= *n; i += 6) {
	dtemp = dtemp + (d__1 = DX(i), d__1) + (d__2 = DX(i + 1), 
		d__2) + (d__3 = DX(i + 2), d__3) + (d__4 = DX(i + 3), 
		d__4) + (d__5 = DX(i + 4), d__5) + (d__6 = DX(i + 5)
		, d__6);
/* L50: */
    }
L60:
    ret_val = dtemp;
    return ret_val;
} /* dsum_ */

ica.c100644   3665     17       77075  6710205660  10455 0ustar  enghoffcnl/******************************************************************************/
/* For all support and information please contact:                            */
/*                                                                            */
/*   Sigurd Enghoff                                                           */
/*   The Salk Institute, CNL                                                  */
/*   enghoff@salk.edu                                                         */
/*                                                                            */
/* Additional ICA software:                                                   */
/*   http://www.cnl.salk.edu/~enghoff/                                        */
/*                                                                            */
/******************************************************************************/

#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <math.h>
#include "ica.h"

#ifdef R250
#include "r250.h"
#endif

/* Global and externally accessible variables */
integer extended, extblocks, pdfsize, nsub;
integer verbose, block, maxsteps;

doublereal lrate, annealstep, annealdeg, nochange, momentum;


/***************************** Zero-fill a array ******************************/
/* Zero-fill array A of n elements.                                           */
/*                                                                            */
/* n: integer (input)                                                         */
/* A: double array [n] (output)                                                 */

void zero(integer n, doublereal *A) {
	doublereal dzero = 0.0;
	integer izero = 0, ione = 1;

	dcopy_(&n,&dzero,&izero,A,&ione);
}

/************************ Construct an identity matrix ************************/
/* Construct an identity matrix A of size m x m.                              */
/*                                                                            */
/* m: integer (input)                                                         */
/* A: double array [m,m] (output)                                             */

void eye(integer m, doublereal *A) {
	doublereal dzero = 0.0, done = 1.0;
	integer izero = 0, ione = 1;
	integer mxm = m*m, m1 = m+1;

	dcopy_(&mxm,&dzero,&izero,A,&ione);
	dcopy_(&m,&done,&izero,A,&m1);
}


/*********************** Initialize a permutation vector **********************/
/* Initialize a vector perm of n elements for constructing a random.          */
/*                                                                            */
/* permutation of the number in the range [0..n-1]                            */
/* perm: integer array [n] (output)                                           */
/* n:    integer (input)                                                      */

void initperm(integer *perm, integer n) {
	integer i;
	for (i=0 ; i<n ; i++) perm[i]=i;
}


/******************* Obtain and project a random permutation ******************/
/* n elements of vector perm are picked at random and used as indecies into   */
/* data. The designated elements of data are projected by trsf employing the  */
/* bias, if not NULL. k denotes the number of elements of perm left to be     */
/* extracted; perm is updated while calling routine must update k. The        */
/* resulting projections are returned in proj.                                */
/*                                                                            */
/* data: double array [m,?] (input)                                           */
/* trsf: double array [m,m] (input)                                           */
/* bias: double array [m] or NULL (input)                                     */
/* perm: integer array [k>] (input/output)                                    */
/* m:    integer (input)                                                      */
/* n:    integer (input)                                                      */
/* k:    integer (input)                                                      */
/* proj: double array [m,n] (output)                                          */

void randperm(doublereal *data, doublereal *trsf, doublereal *bias, integer *perm, integer m, integer n, integer k, doublereal *proj) {
	integer i, im, swap, inc = 1;
	doublereal alpha = 1.0, beta;
	char trans='N';
	
	if (bias) {
		for (i=0,im=0 ; i<n ; i++,im+=m) dcopy_(&m,bias,&inc,&proj[im],&inc);
		beta = 1.0;
	}
	else
		beta = 0.0;

	for (i=0,im=0 ; i<n ; i++,im+=m) {
		swap = RAND()%k;
		dgemv_(&trans,&m,&m,&alpha,trsf,&m,&data[perm[swap]*m],&inc,&beta,&proj[im],&inc);
		perm[swap] = perm[--k];
	}
}


/********* Obtain and project data points from Laplacian distribution *********/
/* n samples are picked from a Laplacian distribution centered at the center  */
/* point of a length frames window. Data points are extracted from frames in  */
/* data according to the windowed distribution while chosen uniformly over    */
/* epochs. The designated elements of data are projected by trsf employing    */
/* the bias, if not NULL. The resulting projections are returned in proj.     */
/*                                                                            */
/* data:   double array [m,?] (input)                                         */
/* trsf:   double array [m,m] (input)                                         */
/* bias:   double array [m] or NULL (input)                                   */
/* m:      integer (input)                                                    */
/* n:      integer (input)                                                    */
/* frames: integer (input)                                                    */
/* epochs: integer (input)                                                    */
/* proj:   double array [m,n] (output)                                        */

void probperm(doublereal *data, doublereal *trsf, doublereal *bias, integer m, integer n, integer frames, integer epochs, doublereal *proj) {
	integer i, im, idx, norm, prob, inc = 1;
	doublereal alpha = 1.0, beta;
	char trans='N';

	if (bias) {
		for (i=0,im=0 ; i<n ; i++,im+=m) dcopy_(&m,bias,&inc,&proj[im],&inc);
		beta = 1.0;
	}
	else
		beta = 0.0;

	norm = (1<<((frames+1)/2+1))-2*((frames%2)-1);

	for (i=0,im=0 ; i<n ; i++,im+=m) {
		prob = RAND()%norm;

		if (prob & 1) {
			idx = 0;
			prob = (prob>>1) + 1;
			while (prob > 1) {
				prob >>= 1;
				idx++;
			}
		}
		else {
			idx = frames-1;
			prob = (prob>>1) + 1;
			while (prob > 1) {
				prob >>= 1;
				idx--;
			}
		}

		idx += frames*(RAND()%epochs);
		dgemv_(&trans,&m,&m,&alpha,trsf,&m,&data[idx*m],&inc,&beta,&proj[im],&inc);
	}
}


/**************************** Compute PDF estimate ****************************/
/* If perm is not NULL, n elements of vector perm are picked at random and    */
/* used as indecies into data. If perm is NULL, the first n elements of data  */
/* are referenced. The designated n elements of data used for the Kurosis     */
/* estimate. k denotes the number of elements of perm left to be extracted;   */
/* perm is updated while calling routine must update k. The resulting         */
/* PDF estimate is returned in kk.                                            */
/*                                                                            */
/* data: double array [m,? or n] (input)                                      */
/* trsf: double array [m,m] (input)                                           */
/* perm: integer array [k>] (input/output) or NULL                            */
/* m:    integer (input)                                                      */
/* n:    integer (input)                                                      */
/* k:    integer (input)                                                      */
/* kk:   double array [m,n] (output)                                          */

void pdf(doublereal *data, doublereal *trsf, integer *perm, integer m, integer n, integer k, doublereal *kk) {
	integer i, j, im, swap, inc = 1;
	doublereal alpha = 1.0, beta = 0.0, tmp;
	char trans='N';

#ifdef KURTOSIS_ESTIMATE
    zero(2*m,&kk[m]);

    for (i=0,im=0 ; i<n ; i++,im+=m) {
        if (perm) {
            swap = RAND()%k;
            dgemv_(&trans,&m,&m,&alpha,trsf,&m,&data[perm[swap]*m],&inc,&beta,kk,&inc);
            perm[swap] = perm[--k];
        }
        else
            dgemv_(&trans,&m,&m,&alpha,trsf,&m,&data[im],&inc,&beta,kk,&inc);
        
        for (j=0 ; j<m ; j++) {
            tmp = kk[j]*kk[j];
            kk[j+m] += tmp;
            kk[j+2*m] += tmp*tmp;
        }
    }

    for (i=0 ; i<m ; i++)
        kk[i] = kk[i+2*m]*(doublereal)n / (kk[i+m]*kk[i+m]) - 3.0;

#else

	zero(3*m,&kk[m]);

	for (i=0,im=0 ; i<n ; i++,im+=m) {
		if (perm) {
			swap = RAND()%k;
			dgemv_(&trans,&m,&m,&alpha,trsf,&m,&data[perm[swap]*m],&inc,&beta,kk,&inc);
			perm[swap] = perm[--k];
		}
		else
			dgemv_(&trans,&m,&m,&alpha,trsf,&m,&data[im],&inc,&beta,kk,&inc);
		
		for (j=0 ; j<m ; j++) {
			tmp = 2.0 / (exp(kk[j])+exp(-kk[j]));
			kk[j+m] += tmp*tmp;
			kk[j+2*m] += kk[j]*kk[j];
			kk[j+3*m] += tanh(kk[j])*kk[j];
		}
	}

	for (i=0 ; i<m ; i++)
		kk[i] = (kk[i+m]*kk[i+2*m]/((doublereal)n)-kk[i+3*m])/((doublereal)n);
#endif
}


/******************* Project data using a general projection ******************/
/* Project data using trsf. Return projections in proj.                       */
/*                                                                            */
/* data: double array [m,n] (input)                                           */
/* trsf: double array [m,m] (input)                                           */
/* m:    integer (input)                                                      */
/* n:    integer (input)                                                      */
/* proj: double array [m,n] (output)                                          */

void geproj(doublereal *data, doublereal *trsf, integer m, integer n, doublereal *proj) {
	doublereal alpha = 1.0, beta = 0.0;
	char trans='N';

	dgemm_(&trans,&trans,&m,&n,&m,&alpha,trsf,&m,data,&m,&beta,proj,&m);
}


/***************** Project data using a symmetrical projection ****************/
/* Project data using the upper triangular part of trsf, assuming trsf is     */
/* symmetrical. Return projections in proj.                                   */
/*                                                                            */
/* data: double array [m,n] (input)                                           */
/* trsf: double array [m,m] (input)                                           */
/* m:    integer (input)                                                      */
/* n:    integer (input)                                                      */
/* proj: double array [m,n] (output)                                          */

void syproj(doublereal *data, doublereal *trsf, integer m, integer n, doublereal *proj) {
	doublereal alpha = 1.0, beta = 0.0;
	char uplo='U', side = 'L';

	dsymm_(&side,&uplo,&m,&n,&alpha,trsf,&m,data,&m,&beta,proj,&m);
}


/******************* Project data using a reduced projection ******************/
/* Project data using the first k rows of trsf only. Return projections       */
/* in proj.                                                                   */
/*                                                                            */
/* data: double array [k,n] (input)                                           */
/* trsf: double array [m,k] (input)                                           */
/* m:    integer (input)                                                      */
/* n:    integer (input)                                                      */
/* k:    integer (input)                                                      */
/* proj: double array [m,n] (output)                                          */

void pcaproj(doublereal *data, doublereal *eigv, integer m, integer n, integer k, doublereal *proj) {
	doublereal alpha = 1.0, beta = 0.0;
	char transn='N', transt='T';

	dgemm_(&transt,&transn,&m,&n,&k,&alpha,eigv,&k,data,&k,&beta,proj,&m);
}


/******************************* Used by varsort ******************************/

int compar(const void *x, const void *y) {
	if (((idxelm*)x)->val < ((idxelm*)y)->val) return 1;
	if (((idxelm*)x)->val > ((idxelm*)y)->val) return -1;
	return 0;
}

/****************** Sort data according to projected variance *****************/
/* Compute back-projected variances for each component based on inverse of    */
/* weights and sphere or pseudoinverse of weights, sphere, and eigv. Reorder  */
/* data and weights accordingly. Also if not NULL, reorder bias and signs.    */
/*                                                                            */
/* data:    double array [m,n] (input/output)                                 */
/* weights: double array [m,m] (input/output)                                 */
/* sphere:  double array [m,m] (input)                                        */
/* eigv:    double array [m,k] (input) or NULL                                */
/* bias:    double array [m] (input/output) or NULL                           */
/* signs:   integer array [m] (input/output) or NULL                          */
/* m:       integer (input)                                                   */
/* n:       integer (input)                                                   */
/* k:       integer (input)                                                   */

void varsort(doublereal *data, doublereal *weights, doublereal *sphere, doublereal *eigv, doublereal *bias, integer *signs, integer m, integer n, integer k) {
	char name[] = "DGETRI\0", opts[] = "\0";
	doublereal alpha = 1.0, beta = 0.0;
	integer i, j, l, jm, ik, info = 0, ispec = 1, na = -1;
	char transn='N', transt='T', uplo='U', side = 'R';
	
	integer nb = ilaenv_(&ispec,name,opts,&m,&na,&na,&na);
	integer itmp, lwork = m*nb, inc = 1;
	doublereal act, dtmp, *wcpy;

	integer    *ipiv = (integer*)malloc(m*sizeof(integer));
	doublereal *work = (doublereal*)malloc(lwork*sizeof(doublereal));
	doublereal *winv = (doublereal*)malloc(m*k*sizeof(doublereal));
	doublereal *sum  = (doublereal*)malloc(n*sizeof(doublereal));
	idxelm  *meanvar = (idxelm*)malloc(m*sizeof(idxelm));
	
	if (eigv) {
/* Compute pseudoinverse of weights*sphere*eigv */
		wcpy = (doublereal*)malloc(m*m*sizeof(doublereal));
		dsymm_(&side,&uplo,&m,&m,&alpha,sphere,&m,weights,&m,&beta,wcpy,&m);
		dgemm_(&transn,&transt,&m,&k,&m,&alpha,wcpy,&m,eigv,&k,&beta,weights,&m);

		dgetrf_(&m,&m,wcpy,&m,ipiv,&info);
		dgetri_(&m,wcpy,&m,ipiv,work,&lwork,&info);
		dgemm_(&transn,&transn,&k,&m,&m,&alpha,eigv,&k,wcpy,&m,&beta,winv,&k);
		free(wcpy);
	}
	else {
/* Compute inverse of weights*sphere */
		dsymm_(&side,&uplo,&m,&m,&alpha,sphere,&m,weights,&m,&beta,winv,&m);

		dgetrf_(&m,&m,winv,&m,ipiv,&info);
		dgetri_(&m,winv,&m,ipiv,work,&lwork,&info);
	}


/* Compute mean variances for back-projected components */
	for (i=0 ; i<m*k ; i++) winv[i] = winv[i]*winv[i];

	for (i=0,ik=0 ; i<m ; i++,ik+=k) {
		for (j=0,jm=0 ; j<n ; j++,jm+=m) {
			sum[j] = 0;
			act = data[i+jm]*data[i+jm];
			for(l=0 ; l<k ; l++) sum[j] += act*winv[l+ik];
		}
		
		meanvar[i].idx = i;
		meanvar[i].val = dsum_(&n,sum,&inc);
		if (verbose) printf("%d ",(int)(i+1));
	}
	if (verbose) printf("\n");
	
/* Sort meanvar */
	qsort(meanvar,m,sizeof(idxelm),compar);

	if (verbose) printf("Permuting the activation wave forms ...\n");

/* Perform in-place reordering of weights, data, bias, and signs */
	for (i=0 ; i<m-1 ; i++) {
		j = meanvar[i].idx;
		if (i != j) {
			dswap_(&k,&weights[i],&m,&weights[j],&m);
			dswap_(&n,&data[i],&m,&data[j],&m);

			if (bias) {
				dtmp = bias[i];
				bias[i] = bias[j];
				bias[j] = dtmp;
			}

			if (signs) {
				itmp = signs[i];
				signs[i] = signs[j];
				signs[j] = itmp;
			}
			
			for (l=i+1 ; i!=meanvar[l].idx ; l++);
			meanvar[l].idx = j;
		}
	}
	
	free(ipiv);
	free(work);
	free(winv);
	free(sum);
	free(meanvar);
}


/******************************* Remove row mean ******************************/
/* Remove row means from data.                                                */
/*                                                                            */
/* data: double array [m,n] (input/output)                                    */
/* m:    integer (input)                                                      */
/* n:    integer (input)                                                      */

void rmmean(doublereal *data, integer m, integer n) {
	doublereal mean;
	integer i, j;
	
	for (i=0 ; i<m ; i++) {
		mean = dsum_(&n,&data[i],&m) / (doublereal)n;
		for (j=i ; j<n*m ; j+=m) data[j] -= mean;
	}
}


/*************************** Compute sphering matrix **************************/
/* Compute sphering matrix based on data.                                     */
/*                                                                            */
/* data: double array [m,n] (input)                                           */
/* m:    integer (input)                                                      */
/* n:    integer (input)                                                      */
/* sphe: double array [m,m] (output)                                          */

void do_sphere(doublereal *data, integer m, integer n, doublereal *sphe) {
	char name[] = "SSYTRD\0", opts[] = "U\0";
	doublereal alpha = 1.0/(doublereal)(n-1);
	doublereal beta = 0.0;
	integer info = 0, ispec = 1, na = -1;
	
	integer nb = ilaenv_(&ispec,name,opts,&m,&na,&na,&na);
	integer i, im, lwork = (nb+2)*m, inc = 1, mxm = m*m;

	char uplo='U', transn='N', jobz='V';
	doublereal *eigv = (doublereal*)malloc(m*m*sizeof(doublereal));
	doublereal *eigd = (doublereal*)malloc(m*sizeof(doublereal));
	integer    *ipiv = (integer*)malloc(m*sizeof(integer));
	doublereal *work = (doublereal*)malloc(lwork*sizeof(doublereal));

/*
   sphere = 2*inv(sqrtm(cov(data'))) is computes as:

   [v d] = eig(data*data')
   sphere = inv(diag(0.5 * sqrt(d))*v')*v'
*/

	dsyrk_(&uplo,&transn,&m,&n,&alpha,data,&m,&beta,sphe,&m);
	dsyev_(&jobz,&uplo,&m,sphe,&m,eigd,work,&lwork,&info);

/* Store transpose of sphe in eigv */
	for (i=0,im=0 ; i<m ; i++,im+=m)
		dcopy_(&m,&sphe[im],&inc,&eigv[i],&m);
	
/* Copy eigv to sphe */
	dcopy_(&mxm,eigv,&inc,sphe,&inc);
	
/* Scale rows of eigv by corresponding eigd values */
	for (i=0 ; i<m ; i++) {
		eigd[i] = 0.5 * sqrt(eigd[i]);
		dscal_(&m,&eigd[i],&eigv[i],&m);
	}

/* Solve eigv * sphere = sphe  ~  diag(0.5*sqrt(d))*v' * sphere = v' */
	dgesv_(&m,&m,eigv,&m,ipiv,sphe,&m,&info);

	free(eigv);
	free(eigd);
	free(ipiv);
	free(work);
}


/***************************** Compute PCA matrix *****************************/
/* Compute PCA decomposition matrix based on data.                            */
/*                                                                            */
/* data: double array [m,n] (input)                                           */
/* m:    integer (input)                                                      */
/* n:    integer (input)                                                      */
/* trsf: double array [m,m] (output)                                          */

void pca(doublereal *data, integer m, integer n, doublereal *trsf) {
	char name[] = "SSYTRD\0", opts[] = "U\0";
	doublereal alpha = 1.0/(doublereal)(n-1);
	doublereal beta = 0.0;
	integer info = 0, ispec = 1, na = -1;
	
	integer nb = ilaenv_(&ispec,name,opts,&m,&na,&na,&na);
	integer lwork = (nb+2)*m;

	char uplo='U', transn='N', jobz='V';
	doublereal *eigd = (doublereal*)malloc(m*sizeof(doublereal));
	doublereal *work = (doublereal*)malloc(lwork*sizeof(doublereal));

	dsyrk_(&uplo,&transn,&m,&n,&alpha,data,&m,&beta,trsf,&m);
	dsyev_(&jobz,&uplo,&m,trsf,&m,eigd,work,&lwork,&info);
	
	free(eigd);
	free(work);
}


/*************** Orient components toward positive activations ****************/
/* Project data using trsf. Negate components and their corresponding weights */
/* to assure positive RMS. Returns projections in proj.                       */
/*                                                                            */
/* data: double array [m,n] (input)                                           */
/* trsf: double array [m,m] (input/output)                                    */
/* m:    integer (input)                                                      */
/* n:    integer (input)                                                      */
/* proj: double array [m,n] (output)                                          */

void posact(doublereal *data, doublereal *trsf, integer m, integer n, doublereal *proj) {
	char trans='N';
	doublereal alpha = 1.0;
	doublereal beta = 0.0;
	doublereal posrms, negrms;
	integer pos, neg, i, j;
	
	dgemm_(&trans,&trans,&m,&n,&m,&alpha,trsf,&m,data,&m,&beta,proj,&m);

	if (verbose) printf("Inverting negative activations: ");
	for (i=0 ; i<m ; i++) {
		posrms = 0.0; negrms = 0.0;
		pos = 0; neg = 0;
		for (j=i ; j<m*n ; j+=m)
			if (proj[j] >= 0) {
				posrms += proj[j]*proj[j];
				pos++;
			}
			else {
				negrms += proj[j]*proj[j];
				neg++;
			}

		if (negrms*(doublereal)pos > posrms*(doublereal)neg) {
			if (verbose) printf("-");
			for (j=i ; j<m*n ; j+=m) proj[j] = -proj[j];
			for (j=i ; j<m*m ; j+=m) trsf[j] = -trsf[j];
		}
		if (verbose) printf("%d ",(int)(i+1));
	}
	printf("\n");
}


/******************************* Perform infomax ******************************/
/* Perform infomax or extended-infomax on data. If any elements of weights    */
/* are none-zero use weights as starting weights else use identity matrix.    */
/* If bias not NULL employ biasing. The following externally accessible       */
/* variables are assumed initialized: extended, extblocks, pdfsize, nsub,     */
/* verbose, block, maxsteps, lrate, annealstep, annealdeg, nochange, and      */
/* momentum. If the boolean variable extended is set, signs must be defined   */
/* (i.e. not NULL)                                                            */
/*                                                                            */
/* data:    double array [chans,frames*epoch] (input)                         */
/* weights: double array [chans,chans] (input/output)                         */
/* chans:   integer (input)                                                   */
/* frames:  integer (input)                                                   */
/* epochs:  integer (input)                                                   */
/* bias:    double array [chans] (output) or NULL                             */
/* signs:   integer array [chans] (output) or NULL                            */

void runica(doublereal *data, doublereal *weights, integer chans, integer frames, integer epochs, doublereal *bias, integer *signs) {
	doublereal alpha = 1.0, beta = 0.0, gamma = -1.0, epsilon;
	doublereal change, oldchange, angledelta;
	char uplo='U', transn='N', transt='T';
	integer i, j, t = 0, inc = 1;
	
	integer datalength = frames*epochs;
	integer chxch = chans*chans;
	integer norm = (1<<((frames+1)/2+1))-2*(frames%2)-2;

	doublereal signsbias = DEFAULT_SIGNSBIAS;
	doublereal extmomentum = DEFAULT_EXTMOMENTUM;
	
	integer wts_blowup = 0;
	integer signcount = 0;
	integer pleft = 0;
	integer step = 0;
	integer blockno = 1;
	integer urextblocks;

/* Allocate arrays */
	doublereal *startweights = (doublereal*)malloc(chxch*sizeof(doublereal));
	doublereal *oldweights = (doublereal*)malloc(chxch*sizeof(doublereal));
	doublereal *tmpweights = (doublereal*)malloc(chxch*sizeof(doublereal));
	doublereal *delta = (doublereal*)malloc(chxch*sizeof(doublereal));
	doublereal *olddelta = (doublereal*)malloc(chxch*sizeof(doublereal));
	doublereal *u = (doublereal*)malloc(chans*block*sizeof(doublereal));
	doublereal *y = (doublereal*)malloc(chans*block*sizeof(doublereal));
	doublereal *yu = (doublereal*)malloc(chxch*sizeof(doublereal));
	doublereal *prevweights, *prevwtchange;
	doublereal *bsum, *kk, *old_kk;
	integer *oldsigns, *pdfperm;

#ifndef PROB_WINDOW
	integer *dataperm = (integer*)malloc(datalength*sizeof(integer));
#endif

/* Initialize weights if zero array */
	if (weights[idamax_(&chxch,weights,&inc)-1] == 0.0) eye(chans,weights);

/* Allocate and initialize arrays and variables needed for momentum */
	if (momentum > 0.0) {
		prevweights = (doublereal*)malloc(chxch*sizeof(doublereal));
		prevwtchange = (doublereal*)malloc(chxch*sizeof(doublereal));
		dcopy_(&chxch,weights,&inc,prevweights,&inc);
		zero(chxch,prevwtchange);
	}
	else
		prevweights = prevwtchange = NULL;

/* Allocate and initialize arrays and variables needed for biasing */
	if (bias) {
		bsum = (doublereal*)malloc(chans*sizeof(doublereal));
		zero(chans,bias);
	}
	else
		bsum = NULL;

/* Allocate and initialize arrays and variables needed for extended-ICA */
	if (extended) {
		oldsigns = (integer*)malloc(chans*sizeof(integer));
		for (i=0 ; i<chans ; i++) oldsigns[i] = -1;
		for (i=0 ; i<nsub ; i++) signs[i] = 1;
		for (i=nsub ; i<chans ; i++) signs[i] = 0;
		
		if (pdfsize < datalength)
			pdfperm = (integer*)malloc(datalength*sizeof(integer));
		else {
			pdfsize = datalength;
			pdfperm = NULL;
		}
		
#ifdef KURTOSIS_ESTIMATE
		kk = (doublereal*)malloc(3*chans*sizeof(doublereal));
#else
		kk = (doublereal*)malloc(4*chans*sizeof(doublereal));
#endif
		old_kk = (doublereal*)malloc(chans*sizeof(doublereal));
		zero(chans,old_kk);

		if (extblocks<0 && verbose) {
			printf("Fixed extended-ICA sign assignments:  ");
			for (i=0 ; i<chans ; i++) printf("%d ",(int)(signs[i]));
			printf("\n");
		}
	}
	else {
		oldsigns = pdfperm = NULL;
		old_kk = kk = NULL;
	}
	urextblocks = extblocks;

/*************************** Initialize ICA training **************************/

	dcopy_(&chxch,weights,&inc,startweights,&inc);
	dcopy_(&chxch,weights,&inc,oldweights,&inc);
	
	if (verbose) {
		printf("Beginning ICA training ...");
		if (extended) printf(" first training step may be slow ...\n");
		else printf("\n");
	}

#ifdef FIX_SEED	
	SRAND(1);
#else
	SRAND((int)time(NULL));
#endif

	while (step < maxsteps) {

#ifndef PROB_WINDOW
		initperm(dataperm,datalength);
#endif

/***************************** ICA training block *****************************/
		for (t=0 ; t<datalength-block && !wts_blowup ; t+=block) {

#ifdef PROB_WINDOW
			probperm(data,weights,bias,chans,block,frames,epochs,u);
#else
			randperm(data,weights,bias,dataperm,chans,block,datalength-t,u);
#endif

			if (!extended) {
/************************* Logistic ICA weight update *************************/
				for (i=0 ; i<chans*block ; i++)
					y[i] = 1.0 - 2.0 / (1.0+exp(-u[i]));

/*Bias sum for logistic ICA */
				if (bias)
					for (i=0 ; i<chans ; i++)
						bsum[i] = dsum_(&block,&y[i],&chans);

/* Compute: (1-2*y) * u' */
				dgemm_(&transn,&transt,&chans,&chans,&block,&alpha,y,&chans,u,&chans,&beta,yu,&chans);
			}
			else {
/************************* Extended-ICA weight update *************************/
				for (i=0 ; i<chans*block ; i++)
					y[i] = tanh(u[i]);

/* Bias sum for extended-ICA */
				if (bias)
					for (i=0 ; i<chans ; i++)
						bsum[i] = -2*dsum_(&block,&y[i],&chans);

/* Apply sign matrix */
				for (i=0 ; i<chans ; i++)
					if (signs[i])
						for (j=i ; j<block*chans ; j+=chans)
							y[j] = -y[j];

/* Compute: u * u' */
				dsyrk_(&uplo,&transn,&chans,&block,&alpha,u,&chans,&beta,yu,&chans);

				j = chxch - 2;
				for (i=1 ; i<chans ; i++) {
					dcopy_(&i,&yu[j],&chans,&yu[j-chans+1],&inc);
					j -= chans+1;
				}

/* Compute: -y * u' -u*u' */
				dgemm_(&transn,&transt,&chans,&chans,&block,&gamma,y,&chans,u,&chans,&gamma,yu,&chans);
			}
			
/* Add block identity matix */
			for (i=0 ; i<chxch ; i+=(chans+1))
				yu[i] += (doublereal)block;

/* Apply weight change */
			dcopy_(&chxch,weights,&inc,tmpweights,&inc);
			dgemm_(&transn,&transn,&chans,&chans,&chans,&lrate,yu,&chans,tmpweights,&chans,&alpha,weights,&chans);

/* Apply bias change */
			if (bias) daxpy_(&chans,&lrate,bsum,&inc,bias,&inc);
			
/******************************** Add momentum ********************************/
			if (momentum > 0.0) {
				daxpy_(&chxch,&momentum,prevwtchange,&inc,weights,&inc);
				for (i=0 ; i<chxch ; i++)
					prevwtchange[i] = weights[i] - prevweights[i];

				dcopy_(&chxch,weights,&inc,prevweights,&inc);
			}
			
			if (abs(weights[idamax_(&chxch,weights,&inc)-1]) > MAX_WEIGHT)
				wts_blowup = 1;

			if (extended && !wts_blowup && extblocks>0 && blockno%extblocks==0) {
				if (pdfperm && pleft < pdfsize) {
					initperm(pdfperm,datalength);
					pleft = datalength;
				}
				
				pdf(data,weights,pdfperm,chans,pdfsize,pleft,kk);
				pleft -= pdfsize;
				
				if (extmomentum > 0.0) {
					epsilon = 1.0-extmomentum;
					dscal_(&chans,&epsilon,kk,&inc);
					daxpy_(&chans,&extmomentum,old_kk,&inc,kk,&inc);
					dcopy_(&chans,kk,&inc,old_kk,&inc);
				}
				
				for (i=0 ; i<chans ; i++)
					signs[i] = kk[i] < -signsbias;
				
				for (i=0 ; i<chans && signs[i]==oldsigns[i] ; i++);
				if (i==chans) signcount++;
				else signcount = 0;
				
				for (i=0 ; i<chans ; i++)
					oldsigns[i] = signs[i];
								
				if (signcount >= SIGNCOUNT_THRESHOLD) {
					extblocks = (integer)(extblocks * SIGNCOUNT_STEP);
					signcount = 0;
				}
			}
			blockno++;			
		}

		if (!wts_blowup) {
			step++;
			angledelta = 0.0;
			
			for (i=0 ; i<chxch ; i++)
				delta[i] = weights[i]-oldweights[i];
			
			change = ddot_(&chxch,delta,&inc,delta,&inc);
		}
		
/************************* Restart if weights blow up *************************/
		if (wts_blowup) {
			printf("\a");
			step = 0;
			change = nochange;
			wts_blowup = 0;
			blockno = 1;
			extblocks = urextblocks;
			lrate = lrate*DEFAULT_RESTART_FAC;
			dcopy_(&chxch,startweights,&inc,weights,&inc);
			zero(chxch,delta);
			zero(chxch,olddelta);

			if (bias) zero(chans,bias);
			
			if (momentum > 0.0) {
				dcopy_(&chxch,startweights,&inc,oldweights,&inc);
				dcopy_(&chxch,startweights,&inc,prevweights,&inc);
				zero(chxch,prevwtchange);
			}
						
			if (extended) {
				for (i=0 ; i<chans ; i++) oldsigns[i] = -1;
				for (i=0 ; i<nsub ; i++) signs[i] = 1;
				for (i=nsub ; i<chans ; i++) signs[i] = 0;
			}
			
			if (lrate > MIN_LRATE)
				printf("Lowering learning rate to %g and starting again.\n",lrate);
			else
				error("QUITTING - weight matrix may not be invertible!\n");
		}
		else {
/*********************** Print weight update information **********************/
			if (step > 2) {
				epsilon = ddot_(&chxch,delta,&inc,olddelta,&inc);
				angledelta = acos(epsilon/sqrt(change*oldchange));

				if (verbose) {
					if (!extended)
						printf("step %d - lrate %5f, wchange %7.6f, angledelta %4.1f deg\n",step,lrate,change,DEGCONST*angledelta);
					else {
						for (i=0,j=0 ; i<chans ; i++) j += signs[i];
						printf("step %d - lrate %5f, wchange %7.6f, angledelta %4.1f deg, %d subgauss\n",step,lrate,change,DEGCONST*angledelta,j);
					}
				}
			}
			else
				if (verbose) {
					if (!extended)
						printf("step %d - lrate %5f, wchange %7.6f\n",step,lrate,change);
					else {
						for (i=0,j=0 ; i<chans ; i++) j += signs[i];
						printf("step %d - lrate %5f, wchange %7.6f, %d subgauss\n",step,lrate,change,j);
					}
				}
		}

/**************************** Save current values *****************************/
		dcopy_(&chxch,weights,&inc,oldweights,&inc);
		
		if (DEGCONST*angledelta > annealdeg) {
			dcopy_(&chxch,delta,&inc,olddelta,&inc);
			lrate = lrate*annealstep;
			oldchange = change;
		}
		else
			if (step == 1) {
				dcopy_(&chxch,delta,&inc,olddelta,&inc);
				oldchange = change;
			}

		if (step > 2 && change < nochange)
			step = maxsteps;
		else
			if (change > DEFAULT_BLOWUP)
				lrate = lrate*DEFAULT_BLOWUP_FAC;
	}

	if (bsum) free(bsum);
	if (oldsigns) free(oldsigns);
	if (prevweights) free(prevweights);
	if (prevwtchange) free(prevwtchange);
	if (startweights) free(startweights);
	if (oldweights) free(oldweights);
	if (tmpweights) free(tmpweights);
	if (delta) free(delta);
	if (olddelta) free(olddelta);
	if (kk) free(kk);
	if (old_kk) free(old_kk);
	if (yu) free(yu);
	if (y) free(y);
	if (u) free(u);
	if (pdfperm) free(pdfperm);

#ifndef PROB_WINDOW
	if (dataperm) free(dataperm);
#endif
}

interfc.c100644   3665     17       65274  6714110025  11343 0ustar  enghoffcnl/******************************************************************************/
/* For all support and information please contact:                            */
/*                                                                            */
/*   Sigurd Enghoff                                                           */
/*   The Salk Institute, CNL                                                  */
/*   enghoff@salk.edu                                                         */
/*                                                                            */
/* Additional ICA software:                                                   */
/*   http://www.cnl.salk.edu/~enghoff/                                        */
/*                                                                            */
/******************************************************************************/

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <ctype.h>

#ifdef PVM
#include "pvmica.h"
#endif

#ifdef _WIN32
#include <io.h>
#else
#include <unistd.h>
#include "memap.h"
#endif

#include "ica.h"

#ifdef _WIN32
#define isascii __isascii
#define isatty  _isatty
#define fileno  _fileno
#define strdup  _strdup
#endif

#define TOKENLEN 160      /* Number of characters allocated for tokens. */
#define BUFSIZE 0x10      /* Help file buffer size. */
#define HELPFILE "ica.sc" /* File containing help message. */
#define COMM "!#%\0"      /* Characters to preceed comments. */

/* Globally defined variables */
integer pcaflag, sphering, posactflag, biasflag;


/**************************** Initialize variables ****************************/
/* Initialize global and external variables to their default values as        */
/* defined in ica.h                                                           */

void initdefaults () {
	pcaflag    = DEFAULT_PCAFLAG;
	sphering   = DEFAULT_SPHEREFLAG;
	posactflag = DEFAULT_POSACT;
	verbose    = DEFAULT_VERBOSE;
	biasflag   = DEFAULT_BIASFLAG;
	
	block      = 0;
	lrate      = 0.0;
	annealdeg  = DEFAULT_ANNEALDEG;
	annealstep = 0.0;
	nochange   = DEFAULT_STOP;
	momentum   = DEFAULT_MOMENTUM;
	maxsteps   = DEFAULT_MAXSTEPS;
	
	extended   = DEFAULT_EXTENDED;
	extblocks  = DEFAULT_EXTBLOCKS;
	pdfsize    = MAX_PDFSIZE;
	nsub       = DEFAULT_NSUB;
}


/***************************** Print help message *****************************/
/* Print help message retrieved from HELPFILE to standard output.             */

void help() {
	FILE *file = fopen(HELPFILE,"r");
	char *buffer;
	int items;

	if (!file) error("help file is missing");
	
	buffer = (char*)malloc((BUFSIZE+1)*sizeof(char));
	
	do {
		items = (int)fread(buffer,sizeof(char),BUFSIZE,file);
		buffer[items] = '\0';
		printf("%s",buffer);
	} while (items == BUFSIZE);
	
	fclose(file);
}


/************************ Convert string to lower case ************************/
/* Convert a null-terminated string to lower case.                            */
/*                                                                            */
/* str: char array (input/output)                                             */

void lower(char *str) {
	int i;
	for (i=0 ; str[i]!=0 ; i++) str[i] = (char)tolower(str[i]);
}


/**************************** Free linked key list ****************************/
/* Recursively dismantel a linked list of key elements.                       */
/*                                                                            */
/* keys: key pointer (input)                                                  */

void rmkeys(key *keys) {
	if (keys->prev) {
		rmkeys((key*)keys->prev);
		free(keys->prev);
	}
	free(keys->token);
}


/******************************* Parse a string *******************************/
/* Parse str for the tokens: on (1), off (0), and none (-2) and return the    */
/* corresponding values. Else, str is assumed to contain a boolean value      */
/* which state is returned.                                                   */
/*                                                                            */
/* str: char array (input)                                                    */

integer swtch(char *str) {
	lower(str);

	if (!strcmp(str,"on")) return 1;
	if (!strcmp(str,"off")) return 0;
	if (!strcmp(str,"none")) return -2;
	return atoi(str) > 0;
}


/**************************** Check file contents *****************************/
/* Checks file for binary (none-ascii) contents, i.e. bytes with values >127. */
/* Returns 1 if file contains binary data else 0.                             */
/*                                                                            */
/* file: FILE pointer (input)                                                 */

int isbin(FILE *file) {
	int ch;
	do ch = fgetc(file);
	while (isascii(ch) && ch!=EOF);

	fseek(file,0,SEEK_SET);
	return ch != EOF;
}


/***************************** Check file access ******************************/
/* Checks the access permission for the file named fname. The file is created */
/* with no contence.                                                          */
/*                                                                            */
/* fname: char array (input)                                                  */

int faccess(char *fname) {
	FILE *file = fopen(fname,"wb");

	if (!file) return 0;

	fclose(file);
	return 1;
}


/*************************** Print message and exit ***************************/
/* Prints the error message contain in str and exit with exit code -1.        */
/*                                                                            */
/* str: char array (input)                                                    */

void error(char *str) {
	fprintf(stderr,"\n\aica: %s\n\n",str);

#ifdef PVM
	fflush(stdout);
	pvm_exit();
#endif

	exit(-1);
}


/*************************** Read binary float file ***************************/
/* Read a total of size floting point values from the file specified by       */
/* fname. Subsequently convert and store the values as doublereal in mat.     */
/*                                                                            */
/* fname: char array (input)                                                  */
/* size:  int (input)                                                         */
/* mat:   doublereal array [size] (output)                                    */

void fb_matread(char *fname, int size, doublereal *mat) {
	FILE *file = fopen(fname,"rb");
	float *buffer;
	int i, items;

	if (!file) error("open failed");

#ifdef MMAP
	buffer = (float*)mapmalloc(size*sizeof(float));
#else
	buffer = (float*)malloc(size*sizeof(float));
#endif

	items = (int)fread(buffer,sizeof(float),size,file);
	if (items != size) error("invalid number of elements");
	
	for (i=0 ; i<size ; i++) mat[i] = (doublereal)buffer[i];

#ifdef MMAP
	mapfree(buffer,size*sizeof(float));
#else
	free(buffer);
#endif

	fclose(file);
}


/*************************** Write binary float file **************************/
/* Convert size doublereal elements of mat to floating point values and store */
/* the values in the binary file specified by fname.                          */
/*                                                                            */
/* fname: char array (input)                                                  */
/* size:  int (input)                                                         */
/* mat:   doublereal array [size] (input)                                     */

void fb_matwrite(char *fname, int size, doublereal *mat) {
	FILE *file = fopen(fname,"wb");
	float *buffer;
	int i, items;

	if (!file) error("open failed");

#ifdef MMAP
	buffer = (float*)mapmalloc(size*sizeof(float));
#else
	buffer = (float*)malloc(size*sizeof(float));
#endif

	for (i=0 ; i<size ; i++) buffer[i] = (float)mat[i];
	
	items = (int)fwrite(buffer,sizeof(float),size,file);
	if (items != size) error("invalid number of elements");
	
#ifdef MMAP
	mapfree(buffer,size*sizeof(float));
#else
	free(buffer);
#endif

	fclose(file);
}


/************************** Write ascii integer file **************************/
/* Write size integer values of mat as ascii numbers to the file specified by */
/* fname.                                                                     */
/*                                                                            */
/* fname: char array (input)                                                  */
/* size:  int (input)                                                         */
/* mat:   integer array [size] (input)                                        */

void ia_matwrite(char *fname, int size, integer *mat) {
	FILE *file = fopen(fname,"w");
	int i;

	if (!file) error("open failed");
	for (i=0 ; i<size ; i++) fprintf(file,"%d\n",(int)(mat[i]));
	fclose(file);
}


/*************************** Process linked key list **************************/
/* Process the linked key list pointed to by keys. The processing includes    */
/* interpreting tokens, pre-processing, ICA-decomposition, and post-          */
/* processing.                                                                */
/*                                                                            */
/* keys: key pointer (input)                                                  */

void doit(key *keys) {
	char *keyword, *value;
	
	char *data_f = NULL;
	char *act_f = NULL;
	char *weights_in_f = NULL;
	char *weights_out_f = NULL;
	char *sphere_f = NULL;
	char *bias_f = NULL;
	char *sign_f = NULL;
	
	int i, datasize, chans = 0, frames = 0, epochs = 0, datalength = 0, ncomps = 0;
	doublereal *dataA, *dataB, *weights, *sphere, *bias, *eigv;
	integer *signs;

#ifdef PVM
	int   window[NWINDOW];
	char *fnames[3];
	
	for (i=0 ; i<NWINDOW ; i++) window[i] = -1;

	fprintf(stderr,"\nICA/PVM %s\n",VER_INFO);
#else
	fprintf(stderr,"\nICA %s\n",VER_INFO);
#endif

	initdefaults();

/********************* Argument parsing and error checking ********************/
	while (keys) {
/* Extract next value from linked keys list */
		value = keys->token;
		keys = (key*)(keys->prev);
		
		if (!keys) error("even number of input arguments");
		
/* Extract next keyword from linked keys list */
		keyword = keys->token;
		keys = (key*)(keys->prev);
		lower(keyword);
		
/* Keyword: datafile */
		if (!strcmp(keyword,"datafile"))
			data_f = value;
				
/* Keyword: activationsfile */
		else if (!strcmp(keyword,"activationsfile"))
			act_f = value;
				
/* Keyword: weightsinfile */
		else if (!strcmp(keyword,"weightsinfile"))
			weights_in_f = value;
		
/* Keyword: weightsoutfile */
		else if (!strcmp(keyword,"weightsoutfile"))
			weights_out_f = value;
		
/* Keyword: spherefile */
		else if (!strcmp(keyword,"spherefile"))
			sphere_f = value;
		
/* Keyword: biasfile */
		else if (!strcmp(keyword,"biasfile"))
			bias_f = value;

/* Keyword: signfile */
		else if (!strcmp(keyword,"signfile"))
			sign_f = value;

/* Keywords: chans, chan */
		else if (!strcmp(keyword,"chans") || !strcmp(keyword,"chan")) {
			chans = atoi(value);
			if (chans < 2) error("chans value should be the number of input channels");
		}

#ifdef PVM
/* Keywords: frame, frames */
		else if (!strcmp(keyword,"frame") || !strcmp(keyword,"frames")) {
			frames = atoi(value);
			if (frames < 2) error("frames value should be the number of data points per epoch");
		}

/* Keywords: epoch, epochs */
		else if (!strcmp(keyword,"epoch") || !strcmp(keyword,"epochs")) {
			epochs = atoi(value);
			if (epochs < 2) error("epochs value should be the total number of epochs");
		}
#else
/* Keywords: frames, datalength */
		else if (!strcmp(keyword,"frames") || !strcmp(keyword,"datalength")) {
			datalength = atoi(value);
			if (datalength < 2) error("frames value should be the number of data points");
		}
#endif

/* Keyword: pca */
		else if (!strcmp(keyword,"pca")) {
			ncomps = atoi(value);
			if (ncomps != 0) {
				pcaflag = 1;
				if (ncomps < 1)
					error("pca value should be the number of principal components to retain");
			}
			else
				pcaflag = 0;
		}

/* Keyword: lrate */
		else if (!strcmp(keyword,"lrate")) {
			lrate = atof(value);
			if (lrate>MAX_LRATE || lrate<MIN_LRATE)
				error("lrate value is out of bounds");
		}

/* Keywords: block, blocksize */
		else if (!strcmp(keyword,"block") || !strcmp(keyword,"blocksize")) {
			block = atoi(value);
			if (block < 0)
				error("block size value must be positive");
		}
		
/* Keywords: stop, nochange */
		else if (!strcmp(keyword,"stop") || !strcmp(keyword,"nochange") || !strcmp(keyword,"stopping")) {
			nochange = atof(value);
			if (nochange < 0.0)
				error("stop wchange value must be positive");
		}

/* Keywords: maxsteps, steps */
		else if (!strcmp(keyword,"maxsteps") || !strcmp(keyword,"steps")) {
				maxsteps = atoi(value);
				if (maxsteps < 0)
					error("maxsteps value must be a positive integer");
		}

/* Keywords: anneal, annealstep */
		else if (!strcmp(keyword,"anneal") || !strcmp(keyword,"annealstep")) {
			annealstep = atof(value);
			if (annealstep<=0 || annealstep>1)
				error("anneal step value must be (0,1]");
		}

/* Keywords: annealdeg, degrees */
		else if (!strcmp(keyword,"annealdeg") || !strcmp(keyword,"degrees")) {
			annealdeg = atof(value);
			if (annealdeg>180 || annealdeg<0)
				error("annealdeg value is out of bounds [0,180]");
		}

/* Keyword: momentum */
		else if (!strcmp(keyword,"momentum")) {
			momentum = atof(value);
			if (momentum>1.0 || momentum<0.0)
				error("momentum value is out of bounds [0,1]");
		}

		else if (!strcmp(keyword,"sphering") || !strcmp(keyword,"sphereing") || !strcmp(keyword,"sphere")) {
			sphering = swtch(value);
			if (sphering == -1) error("sphering value must be on, off, or none");
		}

/* Keyword: bias */
		else if (!strcmp(keyword,"bias")) {
			biasflag = swtch(value);
			if (biasflag < 0) error("bias value must be on or off");
		}

/* Keywords: extended, extend */
		else if (!strcmp(keyword,"extended") || !strcmp(keyword,"extend")) {
			extblocks = atoi(value);
			extended = 1;
					
			if (extblocks == 0) extended = 0;
			else
				if (extblocks < 0) nsub = -extblocks;
		}

/* Keyword: posact */
		else if (!strcmp(keyword,"posact")) {
			posactflag = swtch(value);
			if (posactflag < 0) error("posact value must be on or off");
		}
			
/* Keyword: verbose */
		else if (!strcmp(keyword,"verbose")) {
			verbose = swtch(value);
			if (verbose < 0) error("verbose flag value must be on or off");
		}
		
#ifdef PVM
/* Keyword: framewindow */
		else if (!strcmp(keyword,"framewindow")) {
			window[FRAMEWINDOW] = atoi(value);
			if (window[FRAMEWINDOW] <= 0) error("framewindow value must be positive");
		}
/* Keyword: framestep */
		else if (!strcmp(keyword,"framestep")) {
			window[FRAMESTEP] = atoi(value);
			if (window[FRAMESTEP] <= 0) error("framestep value must be positive");
		}
/* Keyword: epochwindow */
		else if (!strcmp(keyword,"epochwindow")) {
			window[EPOCHWINDOW] = atoi(value);
			if (window[EPOCHWINDOW] <= 0) error("epochwindow value must be positive");
		}
/* Keyword: epochstep */
		else if (!strcmp(keyword,"epochstep")) {
			window[EPOCHSTEP] = atoi(value);
			if (window[EPOCHSTEP] <= 0) error("epochstep value must be positive");
		}
/* Keyword: baseline */
		else if (!strcmp(keyword,"baseline")) {
			window[BASELINE] = atoi(value);
			if (window[BASELINE] <= 0) error("Length of baseline must be positive");
		}
#endif
		else
			error("unknown flag");
	}


#ifdef PVM
	datalength = frames*epochs;

	if (window[FRAMEWINDOW] < 0) window[FRAMEWINDOW] = frames;
	if (window[FRAMESTEP]   < 0) window[FRAMESTEP]   = 1;
	if (window[EPOCHWINDOW] < 0) window[EPOCHWINDOW] = epochs;
	if (window[EPOCHSTEP]   < 0) window[EPOCHSTEP]   = 1;
	if (window[BASELINE]    < 0) window[BASELINE]    = 25;

	if (window[FRAMEWINDOW] > frames) error("window frame length must be <= frames");
	if (window[FRAMESTEP]   > frames) error("frame step size must be <= frames");
	if (window[EPOCHWINDOW] > epochs) error("window epoch length must be <= epochs");
	if (window[EPOCHSTEP]   > epochs) error("epoch step size must be <= epochs");
	if (window[BASELINE]    > frames) error("length of baseline must be <= frames");

	datasize = window[FRAMEWINDOW] * window[EPOCHWINDOW];
#else
	datasize = datalength;
#endif

	if (chans < 2) error("invalid number of channels");
	if (datasize < 3) error("invalid data length");

	if (lrate == 0.0) lrate = DEFAULT_LRATE(chans);
	if (block == 0) block = DEFAULT_BLOCK(datasize);
	if (ncomps == 0) ncomps = chans;

	if (ncomps > chans || ncomps < 1) error("invalid number of components");
	if (datasize < chans) error("data length less than data channels");
	if (block < 2) error("block size too small!");
	if (block > datasize) error("block size exceeds data length!");
	if (nsub > ncomps) error("sub-Gaussian components exceeds total number of components!");

	if (annealstep == 0.0)
		annealstep = (extended) ? DEFAULT_EXTANNEAL : DEFAULT_ANNEALSTEP;


	if (extended && extblocks>0) {
		pdfsize = MIN(pdfsize,datalength);
		if (pdfsize < MIN_PDFSIZE)
			fprintf(stderr,"ica: warning, PDF values are inexact\n");
	}


	if (data_f==NULL)
		error("input data file required");

	if (weights_out_f==NULL)
		error("output weights file required");
		
	if (sphere_f==NULL)
		error("output sphering file required");


	if (!faccess(weights_out_f))
		error("weights output file not writable");

	if (!faccess(sphere_f))
		error("sphere file not writable");


	if (act_f!=NULL && !faccess(act_f))
		error("activations file not writable");

	if (bias_f!=NULL && !faccess(bias_f))
		error("bias file not writable");

	if (sign_f!=NULL && !faccess(sign_f))
		error("sign file not writable");

/****************************** Process the data ******************************/
	if (verbose) {
#ifdef PVM
		printf("\nInput data size [%d,%d] = %d channels, %d epoch of %d frames.\n",chans,datalength,chans,epochs,frames);
#else
		printf("\nInput data size [%d,%d] = %d channels, %d frames.\n",chans,datalength,chans,datalength);
#endif
		if (pcaflag) printf("After PCA dimension reduction,\n  finding ");
		else printf("Finding ");
	
		if (!extended)
			printf("%d ICA components using logistic ICA.\n",ncomps);
		else {
			printf("%d ICA components using extended ICA.\n",ncomps);
			
			if (extblocks > 0)
				printf("PDF will be calculated initially every %d blocks using %d data points.\n",extblocks,pdfsize);
			else
				printf("PDF will not be calculated. Exactly %d sub-Gaussian components assumed.\n",nsub);
		}
		
		printf("Initial learning rate will be %g, block size %d.\n",lrate,block);
		
		if (momentum > 0.0)
			printf("Momentum will be %g.\n",momentum);
			
		printf("Learning rate will be multiplied by %g whenever angledelta >= %g deg.\n",annealstep,annealdeg);
		printf("Training will end when wchange < %g or after %d steps.\n",nochange,maxsteps);
		
		if (biasflag)
			printf("Online bias adjustment will be used.\n");
		else
			printf("Online bias adjustment will not be used.\n");
	}
	
/******************************* Allocate memory ******************************/
	if (verbose) printf("\nLoading data from %s\n",data_f);

#ifdef MMAP
	dataA = (doublereal*)mapmalloc(chans*datalength*sizeof(doublereal));
#else
	dataA = (doublereal*)malloc(chans*datalength*sizeof(doublereal));
#endif

	fb_matread(data_f,chans*datalength,dataA);

	weights = (doublereal*)malloc(ncomps*chans*sizeof(doublereal));
	if (weights_in_f!=NULL) {
		if (verbose) printf("Loading weights from %s\n",weights_in_f);
		fb_matread(weights_in_f,ncomps*ncomps,weights);
	}
	else
		zero(ncomps*chans,weights);

	sphere = (doublereal*)malloc(chans*chans*sizeof(doublereal));

	if (biasflag)
		bias = (doublereal*)malloc(ncomps*sizeof(doublereal));
	else
		bias = NULL;
	
	if (extended)
		signs = (integer*)malloc(ncomps*sizeof(integer));
	else
		signs = NULL;


/************************** Remove overall row means **************************/
	if (verbose) printf("Removing mean of each channel ...\n");
	rmmean(dataA,(integer)chans,(integer)datalength);


/**************************** Perform PCA reduction ***************************/
	if (pcaflag) {
		if (verbose) printf("Reducing the data to %d principal dimensions...\n",ncomps);

		eigv = (doublereal*)malloc(chans*chans*sizeof(doublereal));
		pca(dataA,(integer)chans,(integer)datalength,eigv);

#ifdef MMAP
		dataB = (doublereal*)mapmalloc(ncomps*datalength*sizeof(doublereal));
		pcaproj(dataA,&eigv[chans*(chans-ncomps)],(integer)ncomps,(integer)datalength,(integer)chans,dataB);
		mapfree(dataA,chans*datalength*sizeof(doublereal));
#else
		dataB = (doublereal*)malloc(ncomps*datalength*sizeof(doublereal));
		pcaproj(dataA,&eigv[chans*(chans-ncomps)],(integer)ncomps,(integer)datalength,(integer)chans,dataB);
		free(dataA);
#endif
		dataA = dataB;
	}
	else
		eigv = NULL;
	
/**************************** Apply sphering matrix ***************************/
	if (sphering == 1) {
		if (verbose) printf("Computing the sphering matrix...\n");
		do_sphere(dataA,(integer)ncomps,(integer)datalength,sphere);
		
		if (verbose) printf("Sphering the data ...\n");

#ifdef MMAP
		dataB = (doublereal*)mapmalloc(ncomps*datalength*sizeof(doublereal));
		syproj(dataA,sphere,(integer)ncomps,(integer)datalength,dataB);
		mapfree(dataA,ncomps*datalength*sizeof(doublereal));
#else
		dataB = (doublereal*)malloc(ncomps*datalength*sizeof(doublereal));
		syproj(dataA,sphere,(integer)ncomps,(integer)datalength,dataB);
		free(dataA);
#endif
		dataA = dataB;
	}
	else if (sphering == 0) {
		if (weights_in_f==NULL) {
			if (verbose) printf("Using the sphering matrix as the starting weight matrix ...\n");
			do_sphere(dataA,(integer)ncomps,(integer)datalength,weights);
		}

		if (verbose) printf("Returning the identity matrix in variable \"sphere\" ...\n");
		eye((integer)ncomps,sphere);
	}
	else if (sphering == -2) {
		if (verbose) printf("Returning the identity matrix in variable \"sphere\" ...\n");
		eye((integer)ncomps,sphere);
	}

#ifdef PVM
	fnames[0] = weights_out_f;
	fnames[1] = bias_f;
	fnames[2] = sign_f;
	pvmica(dataA,weights,sphere,eigv,(integer)chans,(integer)ncomps,(integer)frames,(integer)epochs,window,bias,signs,fnames);
#else
	runica(dataA,weights,(integer)ncomps,(integer)datalength,1,bias,signs);

/*************** Orient components toward positive activations ****************/
#ifdef MMAP
	dataB = (doublereal*)mapmalloc(ncomps*datalength*sizeof(doublereal));

	if (posactflag) posact(dataA,weights,(integer)ncomps,(integer)datalength,dataB);
	else geproj(dataA,weights,(integer)ncomps,(integer)datalength,dataB);

	mapfree(dataA,ncomps*datalength*sizeof(doublereal));
#else
	dataB = (doublereal*)malloc(ncomps*datalength*sizeof(doublereal));
		
	if (posactflag) posact(dataA,weights,(integer)ncomps,(integer)datalength,dataB);
	else geproj(dataA,weights,(integer)ncomps,(integer)datalength,dataB);

	free(dataA);
#endif
	dataA = dataB;


/******* Sort components in descending order of max projected variance ********/
	if (verbose) {
		if (pcaflag) {
			printf("Composing the eigenvector, weights, and sphere matrices\n");
			printf("  into a single rectangular weights matrix; sphere=eye(%d)\n",chans);
		}
		printf("Sorting components in descending order of mean projected variance ...\n");
	}
	
	if (eigv) {
		varsort(dataA,weights,sphere,&eigv[chans*(chans-ncomps)],bias,signs,(integer)ncomps,(integer)datalength,(integer)chans);
		eye((integer)chans,sphere);
	}
	else
		varsort(dataA,weights,sphere,NULL,bias,signs,(integer)ncomps,(integer)datalength,(integer)chans);


/**************************** Write results to disk ***************************/
	if (verbose) printf("Storing weights in %s\n",weights_out_f);
	fb_matwrite(weights_out_f,ncomps*chans,weights);
	
	if (act_f!=NULL) {
		if (verbose) printf("Storing activations in %s\n",act_f);
		fb_matwrite(act_f,ncomps*datalength,dataA);
	}
	
	if (bias_f!=NULL && bias) {
		if (verbose) printf("Storing bias vector in %s\n",bias_f);
		fb_matwrite(bias_f,ncomps,bias);
	}
	
	if (sign_f!=NULL && signs) {
		if (verbose) printf("Storing sign vector in %s\n",sign_f);
		for (i=0 ; i<ncomps ; i++) signs[i] = (signs[i]) ? (-1) : 1;
		ia_matwrite(sign_f,ncomps,signs);
	}
#endif

	if (verbose) printf("Storing sphering matrix in %s\n",sphere_f);
	fb_matwrite(sphere_f,chans*chans,sphere);

#ifdef MMAP
	if (dataA) mapfree(dataA,ncomps*datalength*sizeof(doublereal));
#else
	if (dataA) free(dataA);
#endif

	if (weights) free(weights);
	if (sphere) free(sphere);
	if (bias) free(bias);
	if (signs) free(signs);
	if (eigv) free(eigv);
}


/******************************** Script parser *******************************/
/* Creates a linked list of keywords and values based on the redirected       */
/* script file.                                                               */
/*                                                                            */
/* argc: int (input)                                                          */
/* argv: char array pointer (input)                                           */

int master(int argc, char **argv) {
	char token[TOKENLEN];
	key  *pkey, *keys = NULL;
	int  gchar;

/* Print help message and exit if arguments were supplied. */
	if (argc > 1) {
		help();
		return 0;
	}

/* Print help and error message if no file was redirected. */
	if (isatty(fileno(stdin))) {
		help();
		error("takes an ascii script file as redirected input");
	}

/* Print help and error message if a redirected file contains binary data. */
	if (isbin(stdin)) {
		help();
		error("takes an ascii script file as input, not a binary file");
	}

/* Decompose script file into a linked list of tokens. */
	while (1) {
		if (scanf("%s",token) < 0) break;

		if (strchr(COMM,(int)token[0]))
			do gchar = getchar();
			while (gchar!=EOL && gchar!=EOF);
		else {
			pkey = keys;
			keys = (key*)malloc(sizeof(key));
			keys->prev = (void*)pkey;
			keys->token = strdup(token);
		}
	}

/* Process token list - do everything */
	doit(keys);

/* Dismantle token list */
	rmkeys(keys);
	
	return 0;
}


#ifdef PVM

/************************ OS entry point (PVM version) ************************/
/* Separate master and slave processes and calls their respective subroutines.*/
/*                                                                            */
/* argc: int (input)                                                          */
/* argv: char array pointer (input)                                           */

int main(int argc, char **argv) {
	int rtnval, parent;

	parent = pvm_parent();

	if (parent < 0 && parent != PvmNoParent) error("Invalid parent process");
	rtnval = (parent == PvmNoParent) ? master(argc,argv) : slave();

	pvm_exit();
	return rtnval;
}

#else

/********************** OS entry point (non-PVM version) **********************/
/* Dummy routine. Allows PVM and non-PVM versions to both access master (old  */
/* main)                                                                      */
/*                                                                            */
/* argc: int (input)                                                          */
/* argv: char array pointer (input)                                           */

int main(int argc, char **argv) {
	return master(argc,argv);
}

#endif
memap.c100644   3665     17        2651  6570622637  10777 0ustar  enghoffcnl/******************************************************************************/
/* For all support and information please contact:                            */
/*                                                                            */
/*   Sigurd Enghoff                                                           */
/*   The Salk Institute, CNL                                                  */
/*   enghoff@salk.edu                                                         */
/*                                                                            */
/* Additional ICA software:                                                   */
/*   http://www.cnl.salk.edu/~enghoff/                                        */
/*                                                                            */
/******************************************************************************/

#include <sys/types.h>
#include <sys/mman.h>
#include <sys/stat.h>
#include <fcntl.h>
#include <stddef.h>

#ifndef ALPHA

void *mapmalloc(int size) {
	int fd;
	caddr_t base;

	fd = open("/dev/zero",O_RDWR);
	base = mmap(NULL,(size_t)(size),PROT_READ|PROT_WRITE,MAP_PRIVATE,fd,0);
	close(fd);
	return (void*)base;
}

#else

void *mapmalloc(int size) {
	caddr_t base;
	
	base = mmap(NULL,(size_t)(size),PROT_READ|PROT_WRITE,MAP_ANONYMOUS|MAP_VARIABLE,-1,0);
	return (void*)base;
}

#endif

void mapfree(void *addr, int size) {
	munmap((caddr_t)addr,size);
}

pvmica.c100644   3665     17       56521  6714114641  11173 0ustar  enghoffcnl/******************************************************************************/
/* For all support and information please contact:                            */
/*                                                                            */
/*   Sigurd Enghoff                                                           */
/*   The Salk Institute, CNL                                                  */
/*   enghoff@salk.edu                                                         */
/*                                                                            */
/* Additional ICA software:                                                   */
/*   http://www.cnl.salk.edu/~enghoff/                                        */
/*                                                                            */
/******************************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <sys/time.h>
#include <sys/resource.h>
#include "pvmica.h"


/********************* Concatenate binary float matricies *********************/
/* Convert size doublereal elements of mat to floating point values and       */
/* appends values to the binary file specified by file.                       */
/*                                                                            */
/* file: FILE pointer (input)                                                 */
/* size: int (input)                                                          */
/* mat:  doublereal array [size] (input)                                      */

void fbc_matwrite(FILE *file, int size, doublereal *mat) {
	float *buffer;
	int i, items;

	if (!file) error("open failed");

	buffer = (float*)malloc(size*sizeof(float));

	for (i=0 ; i<size ; i++) buffer[i] = (float)mat[i];
	
	items = (int)fwrite(buffer,sizeof(float),size,file);
	if (items != size) error("invalid number of elements");
	
	free(buffer);
}


/********************* Concatenate ascii integer matricies ********************/
/* Append size integer values of mat as ascii numbers to the file specified   */
/* by file.                                                                   */
/*                                                                            */
/* file: FILE pointer (input)                                                 */
/* size: int (input)                                                          */
/* mat:  integer array [size] (input)                                         */

void iac_matwrite(FILE *file, int size, integer *mat) {
	int i;

	if (!file) error("open failed");
	for (i=0 ; i<size-1 ; i++) fprintf(file,"%d ",(int)(mat[i]));
	fprintf(file,"%d\n",(int)(mat[i]));
}


/************************** Send assignment to slave **************************/
/* Send assignment structure assign to slave process tid.                     */
/*                                                                            */
/* tid:    int (input)                                                        */
/* assign: tassign pointer (input)                                            */

void send_assign(int tid, tassign *assign) {
	int     bufid, datasize, i, toint[ASSIGN_NINT];
	double *tdbl, todbl[ASSIGN_NDBL];
	float  *tflt;

	toint[ASSIGN_ID]         = (int)(assign->id);
	toint[ASSIGN_CHANS]      = (int)(assign->chans);
	toint[ASSIGN_FRAMES]     = (int)(assign->frames);
	toint[ASSIGN_EPOCHS]     = (int)(assign->epochs);
	toint[ASSIGN_BIAS]       = (int)(assign->bias);
	toint[ASSIGN_SIGNS]      = (int)(assign->signs);
	toint[ASSIGN_EXTENDED]   = (int)(assign->extended);
	toint[ASSIGN_EXTBLOCKS]  = (int)(assign->extblocks);
	toint[ASSIGN_PDFSIZE]    = (int)(assign->pdfsize);
	toint[ASSIGN_NSUB]       = (int)(assign->nsub);
	toint[ASSIGN_VERBOSE]    = (int)(assign->verbose);
	toint[ASSIGN_BLOCK]      = (int)(assign->block);
	toint[ASSIGN_MAXSTEPS]   = (int)(assign->maxsteps);
	
	todbl[ASSIGN_LRATE]      = (double)(assign->lrate);
	todbl[ASSIGN_ANNEALSTEP] = (double)(assign->annealstep);
	todbl[ASSIGN_ANNEALDEG]  = (double)(assign->annealdeg);
	todbl[ASSIGN_NOCHANGE]   = (double)(assign->nochange);
	todbl[ASSIGN_MOMENTUM]   = (double)(assign->momentum);

	bufid = pvm_initsend(PvmDataDefault);
	pvm_pkint(toint,ASSIGN_NINT,1);
	pvm_pkdouble(todbl,ASSIGN_NDBL,1);

	datasize = assign->chans * assign->frames * assign->epochs;
	if (sizeof(float) != sizeof(doublereal)) {
		tflt = (float*)malloc(datasize*sizeof(float));
		for (i=0 ; i<datasize ; i++) tflt[i] = (float)(assign->data[i]);
		pvm_pkfloat(tflt,datasize,1);
		free(tflt);
	}
	else
		pvm_pkfloat((float*)(assign->data),datasize,1);
		
	datasize = assign->chans * assign->chans;
	if (sizeof(double) != sizeof(doublereal)) {
		tdbl = (double*)malloc(datasize*sizeof(double));
		for (i=0 ; i<datasize ; i++) tdbl[i] = (double)(assign->weights[i]);
		pvm_pkdouble(tdbl,datasize,1);
		free(tdbl);
	}
	else
		pvm_pkdouble((double*)(assign->weights),datasize,1);

	pvm_send(tid,1);
	pvm_freebuf(bufid);
}


/********************************* Kill slave *********************************/
/* Send kill signal to slave process tid.                                     */
/*                                                                            */
/* tid: int (input)                                                           */

void send_kill(int tid) {
	int bufid;

	bufid = pvm_initsend(PvmDataDefault);
	pvm_send(tid,2);
	pvm_freebuf(bufid);
}


/*********************** Receive assignment from master ***********************/
/* Receive assignment structure assign from master.                           */
/*                                                                            */
/* assign: tassign pointer (output)                                           */

int receive_assign(tassign *assign) {
	int     bufid, tag, dummy, datasize, i, toint[ASSIGN_NINT];
	double *tdbl, todbl[ASSIGN_NDBL];
	float  *tflt;

	bufid = pvm_recv(-1,-1);
	pvm_bufinfo(bufid,&dummy,&tag,&dummy);
	if (tag == 2) return 0;
	
	pvm_upkint(toint,ASSIGN_NINT,1);
	pvm_upkdouble(todbl,ASSIGN_NDBL,1);

	assign->id         = (integer)toint[ASSIGN_ID];
	assign->chans      = (integer)toint[ASSIGN_CHANS];
	assign->frames     = (integer)toint[ASSIGN_FRAMES];
	assign->epochs     = (integer)toint[ASSIGN_EPOCHS];
	assign->bias       = (integer)toint[ASSIGN_BIAS];
	assign->signs      = (integer)toint[ASSIGN_SIGNS];
	assign->extended   = (integer)toint[ASSIGN_EXTENDED];
	assign->extblocks  = (integer)toint[ASSIGN_EXTBLOCKS];
	assign->pdfsize    = (integer)toint[ASSIGN_PDFSIZE];
	assign->nsub       = (integer)toint[ASSIGN_NSUB];
	assign->verbose    = (integer)toint[ASSIGN_VERBOSE];
	assign->block      = (integer)toint[ASSIGN_BLOCK];
	assign->maxsteps   = (integer)toint[ASSIGN_MAXSTEPS];
	
	assign->lrate      = (doublereal)todbl[ASSIGN_LRATE];
	assign->annealstep = (doublereal)todbl[ASSIGN_ANNEALSTEP];
	assign->annealdeg  = (doublereal)todbl[ASSIGN_ANNEALDEG];
	assign->nochange   = (doublereal)todbl[ASSIGN_NOCHANGE];
	assign->momentum   = (doublereal)todbl[ASSIGN_MOMENTUM];

	datasize = assign->chans * assign->frames * assign->epochs;
	assign->data = (doublereal*)malloc(datasize*sizeof(doublereal));
	if (sizeof(float) != sizeof(doublereal)) {
		tflt = (float*)malloc(datasize*sizeof(float));
		pvm_upkfloat(tflt,datasize,1);
		for (i=0 ; i<datasize ; i++) assign->data[i] = (doublereal)tflt[i];
		free(tflt);
	}
	else
		pvm_upkfloat((float*)(assign->data),datasize,1);

	datasize = assign->chans * assign->chans;
	assign->weights = (doublereal*)malloc(datasize*sizeof(doublereal));
	if (sizeof(double) != sizeof(doublereal)) {
		tdbl = (double*)malloc(datasize*sizeof(double));
		pvm_upkdouble(tdbl,datasize,1);
		for (i=0 ; i<datasize ; i++) assign->weights[i] = (doublereal)tdbl[i];
		free(tdbl);
	}
	else
		pvm_upkdouble((double*)(assign->weights),datasize,1);

	pvm_freebuf(bufid);
	return 1;
}


/*************************** Send results to master ***************************/
/* Send result structure result to master process tid.                        */
/*                                                                            */
/* tid:    int (input)                                                        */
/* result: tresult pointer (input)                                            */

void send_result(int tid, tresult *result) {
	int    *tint, bufid, datasize, i, toint[RESULT_NINT];
	double *tdbl, todbl[RESULT_NDBL];

	toint[RESULT_ID]    = (int)(result->id);
	toint[RESULT_CHANS] = (int)(result->chans);
	toint[RESULT_BIAS]  = (int)(result->bias != NULL);
	toint[RESULT_SIGNS] = (int)(result->signs != NULL);

	todbl[RESULT_LRATE] = (double)(result->lrate);
	
	bufid = pvm_initsend(PvmDataDefault);
	pvm_pkint(toint,RESULT_NINT,1);
	pvm_pkdouble(todbl,RESULT_NDBL,1);

	datasize = result->chans * result->chans;
	if (sizeof(double) != sizeof(doublereal)) {
		tdbl = (double*)malloc(datasize*sizeof(double));
		for (i=0 ; i<datasize ; i++) tdbl[i] = (double)(result->weights[i]);
		pvm_pkdouble(tdbl,datasize,1);
		free(tdbl);
	}
	else
		pvm_pkdouble((double*)(result->weights),datasize,1);

	if (result->bias != NULL) {
		datasize = result->chans;
		if (sizeof(double) != sizeof(doublereal)) {
			tdbl = (double*)malloc(datasize*sizeof(double));
			for (i=0 ; i<datasize ; i++) tdbl[i] = (double)(result->bias[i]);
			pvm_pkdouble(tdbl,datasize,1);
			free(tdbl);
		}
		else
			pvm_pkdouble((double*)(result->bias),datasize,1);
	}

	if (result->signs != NULL) {
		datasize = result->chans;
		if (sizeof(int) != sizeof(integer)) {
			tint = (int*)malloc(datasize*sizeof(int));
			for (i=0 ; i<datasize ; i++) tint[i] = (int)(result->signs[i]);
			pvm_pkint(tint,datasize,1);
			free(tint);
		}
		else
			pvm_pkint((int*)(result->signs),datasize,1);
	}
	
	pvm_send(tid,1);
	pvm_freebuf(bufid);
}


/************************** Receive result from slave *************************/
/* Receive result structure result from the first slave to responding.        */
/*                                                                            */
/* result: tresult pointer (output)                                           */

int receive_result(tresult *result) {
	int    *tint, bufid, datasize, bias, signs, toint[RESULT_NINT];
	int     i, buflen, tag, tid;
	double *tdbl, todbl[RESULT_NDBL];

	bufid = pvm_recv(-1,-1);
	pvm_upkint(toint,RESULT_NINT,1);
	pvm_upkdouble(todbl,RESULT_NDBL,1);

	result->id         = (integer)toint[RESULT_ID];
	result->chans      = (integer)toint[RESULT_CHANS];
	bias               = (integer)toint[RESULT_BIAS];
	signs              = (integer)toint[RESULT_SIGNS];

	result->lrate      = (doublereal)todbl[RESULT_LRATE];

	datasize = result->chans * result->chans;
	result->weights = (doublereal*)malloc(datasize*sizeof(doublereal));
	if (sizeof(double) != sizeof(doublereal)) {
		tdbl = (double*)malloc(datasize*sizeof(double));
		pvm_upkdouble(tdbl,datasize,1);
		for (i=0 ; i<datasize ; i++) result->weights[i] = (doublereal)tdbl[i];
		free(tdbl);
	}
	else
		pvm_upkdouble((double*)(result->weights),datasize,1);

	if (bias) {
		datasize = result->chans;
		result->bias = (doublereal*)malloc(datasize*sizeof(doublereal));
		if (sizeof(double) != sizeof(doublereal)) {
			tdbl = (double*)malloc(datasize*sizeof(double));
			pvm_upkdouble(tdbl,datasize,1);
			for (i=0 ; i<datasize ; i++) result->bias[i] = (doublereal)tdbl[i];
			free(tdbl);
		}
		else
			pvm_upkdouble((double*)(result->bias),datasize,1);
	}
	else
		result->bias = NULL;

	if (signs) {
		datasize = result->chans;
		result->signs = (integer*)malloc(datasize*sizeof(integer));
		if (sizeof(int) != sizeof(integer)) {
			tint = (int*)malloc(datasize*sizeof(int));
			pvm_upkint(tint,datasize,1);
			for (i=0 ; i<datasize ; i++) result->signs[i] = (integer)tint[i];
			free(tint);
		}
		else
			pvm_upkint((int*)(result->signs),datasize,1);
	}
	else
		result->signs = NULL;

	pvm_bufinfo(bufid,&buflen,&tag,&tid);
	pvm_freebuf(bufid);
	
	return tid;
}


/***************************** Main slave routine *****************************/
/* Repeatedly performs ICA on data segments. Assignments are received by the  */
/* process which decomposes transmitted data and returns the results to the   */
/* assigning master process. The sequence is repeated until a kill signal is  */
/* received.                                                                  */

int slave() {
	int         i, datasize, id, ptid;
	tassign     assign;
	tresult     result;
	doublereal *bias, *data, *weights;
	integer    *signs, chans, frames, epochs;

#ifndef CRAY
	setpriority(PRIO_PROCESS,0,SLAVE_NICE);
#endif

	while (receive_assign(&assign)) {
		id         = assign.id;
		data       = assign.data;
		weights    = assign.weights;
		chans      = assign.chans;
		frames     = assign.frames;
		epochs     = assign.epochs;

		extended   = assign.extended;
		extblocks  = assign.extblocks;
		pdfsize    = assign.pdfsize;
		nsub       = assign.nsub;
		verbose    = assign.verbose;
		block      = assign.block;
		maxsteps   = assign.maxsteps;
		lrate      = assign.lrate;
		annealstep = assign.annealstep;
		annealdeg  = assign.annealdeg;
		nochange   = assign.nochange;
		momentum   = assign.momentum;

		if (assign.bias)
			bias = (doublereal*)malloc(chans*sizeof(doublereal));
		else
			bias = NULL;

		if (assign.signs)
			signs = (integer*)malloc(chans*sizeof(integer));
		else
			signs = NULL;

		runica(data,weights,chans,frames,epochs,bias,signs);
	
		result.id      = id;
		result.weights = weights;
		result.chans   = chans;
		result.bias    = bias;
		result.signs   = signs;
		result.lrate   = lrate;

		ptid = pvm_parent();
		send_result(ptid,&result);

		if (data != NULL) free(data);
		if (weights != NULL) free(weights);
		if (bias != NULL) free(bias);
		if (signs != NULL) free(signs);
	}
	return 0;
}


/*************************** Extracts a data segment **************************/
/* Exctracts the data segment corresponding to step from matrix data. The     */
/* resulting segment is stored in windata. Chans, frames and epochs denote    */
/* the dimensions of data. The window array specifies dimensions for grid     */
/* from which to extract windata                                              */
/*                                                                            */
/* data:     doublereal array [chans,frames*epoch] (input)                    */
/* datawin:  doublereal array -                                               */
/*                   [chans,window(FRAMEWINDOW)*window(EPOCHWINDOW)] (output) */
/* chans:    int (input)                                                      */
/* frames:   int (input)                                                      */
/* epochs:   int (input)                                                      */
/* step:     int (input)                                                      */
/* window:   int array [NWINDOW] (input)                                      */

doublereal *extract(doublereal *data, doublereal *windata, int chans, int frames, int epochs, int step, int *window) {
	int i, j, nep, ofs, size;

	nep = (epochs-window[EPOCHWINDOW])/window[EPOCHSTEP] + 1;
	ofs = frames * window[EPOCHSTEP] * (step%nep);
	ofs += window[FRAMESTEP] * (step/nep);
	size = window[FRAMEWINDOW]*chans;
	
	for (i=0,j=0 ; i<window[EPOCHWINDOW] ; i++,j+=size,ofs+=frames)
		memcpy(&windata[j],&data[ofs*chans],size*sizeof(doublereal));

	return windata;
}


/*************************** Extracts data baseline ***************************/
/* Exctracts and concatenates all baseline data contained in matrix data. The */
/* resulting data block is stored in basedata. Chans, frames and epochs       */
/* denote the dimensions of data. The window array contains, amongst others,  */
/* the baseline length.                                                       */
/*                                                                            */
/* data:     doublereal array [chans,frames*epoch] (input)                    */
/* basedata: doublereal array [chans,window(BASELINE)*epoch] (output)         */
/* chans:    int (input)                                                      */
/* frames:   int (input)                                                      */
/* epochs:   int (input)                                                      */
/* window:   int array [NWINDOW] (input)                                      */

doublereal *baseline(doublereal *data, doublereal *basedata, int chans, int frames, int epochs, int *window) {
	int i, j, k, size;

	size = window[BASELINE]*chans;
	for (i=0,j=0,k=0 ; i<epochs ; i++,j+=size,k+=frames*chans)
		memcpy(&basedata[j],&data[k],size*sizeof(doublereal));

	return basedata;
}


/********************** Converts speed to CPU quantities **********************/
/* Converts speed values as represented by host info to CPU quantities in     */
/* multi-processor systems.                                                   */
/*                                                                            */
/* speed: int (input)                                                         */

int speed2proc(int speed) {
	if (speed >= 500) return speed/500;
	return 2;
}


/************************** Executes ICA in parallel **************************/
/* A serial ICA decomposition is initially performed on the baseline portions */
/* of data, all succeeding decompositions are use this result as initial      */
/* weight estimate. A number of tasks are spawned on each node according to   */
/* host info speed values. Each node is repeatedly assigned jobs until all    */
/* decompositions in the frame/epoch grid specified by window are distributed.*/
/* Weights returned by slave processes are continously sorted by projected    */
/* variance stored and to the files specified by fnames. The following        */
/* externally accessible variables are assumed initialized: extended,         */
/* extblocks, pdfsize, nsub, verbose, block, maxsteps, lrate, annealstep,     */
/* annealdeg, nochange, and momentum. If the boolean variable extended is set,*/
/* signs must be defined (i.e. not NULL)                                      */
/*                                                                            */
/* data:    double array [ncomps,frames*epoch] (input)                        */
/* weights: double array [ncomps,ncomps] (input)                              */
/* sphere:  double array [ncomps,ncomps] (input)                              */
/* eigv:    double array [ncomps,chans] (input)                               */
/* chans:   integer (input)                                                   */
/* ncomps:  integer (input)                                                   */
/* frames:  integer (input)                                                   */
/* epochs:  integer (input)                                                   */
/* window:  integer array [NWINDOW] (input)                                   */
/* bias:    double array [ncomps] (dummy) or NULL                             */
/* signs:   integer array [ncomps] (dummy) or NULL                            */
/* fnames:  char array [3] of array                                           */

void pvmica(doublereal *data, doublereal *weights, doublereal *sphere, doublereal *eigv, integer chans, integer ncomps, integer frames, integer epochs, int *window, doublereal *bias, integer *signs, char **fnames) {
	struct pvmhostinfo *hinfo;
	int                 i, j, datasize, maxep, maxfr, segs, speed, spwnd = 0, id = 0;
	int                *tids, tid, nhost, narch, nproc = 0, ntask = 0, last = 0;
	char               *name;
	FILE               *fids[3];
	integer           **srec;
	doublereal        **wrec, **brec, *basedata, *windata, *prjdata, deflr = lrate;
	tresult             result;
	tassign             assign;

/*	pvm_catchout(stdout);*/
	pvm_config(&nhost,&narch,&hinfo);
	for (i=0 ; i<nhost ; i++) nproc += speed2proc(hinfo[i].hi_speed);

	basedata = (doublereal*)malloc(ncomps*epochs*window[BASELINE]*sizeof(doublereal));
	baseline(data,basedata,ncomps,frames,epochs,window);
	runica(basedata,weights,ncomps,1,window[BASELINE]*epochs,bias,signs);
	free(basedata);

	assign.weights = weights;
	assign.chans = ncomps;
	assign.frames = window[FRAMEWINDOW];
	assign.epochs = window[EPOCHWINDOW];
	assign.bias = (int)(bias!=NULL);
	assign.signs = (int)(signs!=NULL);
	assign.extended = extended;
	assign.extblocks = extblocks;
	assign.pdfsize = pdfsize;
	assign.nsub = nsub;
	assign.verbose = verbose;
	assign.block = block;
	assign.maxsteps = maxsteps;
	assign.lrate = deflr;
	assign.annealstep = annealstep;
	assign.annealdeg = annealdeg;
	assign.nochange = nochange;
	assign.momentum = momentum;
	
	datasize = window[FRAMEWINDOW] * window[EPOCHWINDOW] * ncomps;
	windata = (doublereal*)malloc(datasize*sizeof(doublereal));
	prjdata = (doublereal*)malloc(datasize*sizeof(doublereal));
	
	maxep = (epochs-window[EPOCHWINDOW])/window[EPOCHSTEP] + 1;
	maxfr = (frames-window[FRAMEWINDOW])/window[FRAMESTEP] + 1;
	segs = maxep*maxfr;

	if (segs > nproc) {
		tids = (int*)malloc(nproc*sizeof(int));
		for (i=0 ; i<nhost ; i++) {
			name = hinfo[i].hi_name;
			speed = hinfo[i].hi_speed;
			spwnd = pvm_spawn(SPAWN_ICA,NULL,1,name,speed2proc(speed),&tids[ntask]);
			if (spwnd <= 0) error("Failed to spawn processes");
			ntask += spwnd;
		}
	}
	else {
		tids = (int*)malloc(segs*sizeof(int));
		spwnd += pvm_spawn(SPAWN_ICA,NULL,0,"",segs,tids);
		if (spwnd <= 0) error("Failed to spawn processes");
		ntask += spwnd;
	}

	for (i=0,j=0 ; i<ntask ; i++,j+=2,id++) {
		if (j >= ntask) j=1;
		assign.id   = id;
		assign.data = extract(data,windata,(int)ncomps,(int)frames,(int)epochs,(int)id,window);
		send_assign(tids[j],&assign);
	}

	wrec = (doublereal**)malloc(segs*sizeof(doublereal*));
	brec = (doublereal**)malloc(segs*sizeof(doublereal*));
	srec = (integer**)malloc(segs*sizeof(integer*));
	for (i=0 ; i<segs ; i++) {
		wrec[i] = NULL;
		brec[i] = NULL;
		srec[i] = NULL;
	}

	for (i=0 ; i<3 ; i++) fids[i] = NULL;
	if (fnames[0] != NULL) fids[0] = fopen(fnames[0],"wb");
	if (fnames[1] != NULL) fids[1] = fopen(fnames[1],"wb");
	if (fnames[2] != NULL) fids[2] = fopen(fnames[2],"wt");

	for (i=0 ; i<segs ; i++) {
		tid = receive_result(&result);
		
		printf("Received id %d\n",result.id);
		wrec[result.id] = result.weights;
		brec[result.id] = result.bias;
		srec[result.id] = result.signs;

		if (ntask < segs) {
			assign.id = id;
			assign.data = extract(data,windata,(int)ncomps,(int)frames,(int)epochs,(int)id,window);
			send_assign(tid,&assign);
			ntask++;
			id++;
		}
		else {
			send_kill(tid);
		}

		datasize = window[FRAMEWINDOW] * window[EPOCHWINDOW];
		extract(data,windata,(int)ncomps,(int)frames,(int)epochs,(int)result.id,window);
		geproj(windata,result.weights,(integer)ncomps,(integer)datasize,prjdata);
		if (eigv)
			varsort(prjdata,result.weights,sphere,&eigv[chans*(chans-ncomps)],result.bias,result.signs,(integer)ncomps,(integer)datasize,(integer)chans);
		else
			varsort(prjdata,result.weights,sphere,NULL,result.bias,result.signs,(integer)ncomps,(integer)datasize,(integer)chans);

		while (last<segs && wrec[last]!=NULL) {
			if (fids[0]!=NULL && wrec[last]!=NULL) fbc_matwrite(fids[0],chans*ncomps,wrec[last]);
			if (fids[1]!=NULL && brec[last]!=NULL) fbc_matwrite(fids[1],ncomps,brec[last]);
			if (fids[2]!=NULL && srec[last]!=NULL) iac_matwrite(fids[2],ncomps,srec[last]);
			if (wrec[last] != NULL) free(wrec[last]);
			if (brec[last] != NULL) free(brec[last]);
			if (srec[last] != NULL) free(srec[last]);
			last++;
		}
	}

	for (i=0 ; i<3 ; i++) 
		if (fids[i] != NULL) fclose(fids[i]);

	if (wrec != NULL) free(wrec);
	if (brec != NULL) free(brec);
	if (srec != NULL) free(srec);
	if (windata != NULL) free(windata);
	if (prjdata != NULL) free(prjdata);
	if (tids != NULL)    free(tids);
}
r250.c100644   3665     17       10451  6670057345  10404 0ustar  enghoffcnl/* r250.c	the r250 uniform random number algorithm

		Kirkpatrick, S., and E. Stoll, 1981; "A Very Fast
		Shift-Register Sequence Random Number Generator",
		Journal of Computational Physics, V.40

		also:

		see W.L. Maier, DDJ May 1991



*/

static char rcsid[] = "@(#)r250.c	1.2 15:50:31 11/21/94   EFC";


#include <limits.h>

#include "r250.h"

/* set the following if you trust rand(), otherwise the minimal standard
   generator is used
*/
/* #define TRUST_RAND */


#ifndef TRUST_RAND
#include "randlcg.h"
#endif

/* defines to allow for 16 or 32 bit integers */
#define BITS 31


#if WORD_BIT == 32
#ifndef BITS
#define BITS	32
#endif
#else
#ifndef BITS
#define BITS    16
#endif
#endif

#if BITS == 31
#define MSB          0x40000000L
#define ALL_BITS     0x7fffffffL
#define HALF_RANGE   0x20000000L
#define STEP         7
#endif

#if BITS == 32
#define MSB          0x80000000L
#define ALL_BITS     0xffffffffL
#define HALF_RANGE   0x40000000L
#define STEP         7
#endif

#if BITS == 16
#define MSB         0x8000
#define ALL_BITS    0xffff
#define HALF_RANGE  0x4000
#define STEP        11
#endif

static unsigned int r250_buffer[ 250 ];
static int r250_index;

#ifdef NO_PROTO
void r250_init(sd)
int seed;
#else
void r250_init(int sd)
#endif
{
	int j, k;
	unsigned int mask, msb;

#ifdef TRUST_RAND        

#if BITS == 32 || BITS == 31       
	srand48( sd );
#else
	srand( sd );
#endif	


#else
	set_seed( sd );
#endif
	
	r250_index = 0;
	for (j = 0; j < 250; j++)      /* fill r250 buffer with BITS-1 bit values */
#ifdef TRUST_RAND
#if BITS == 32 || BITS == 31
		r250_buffer[j] = (unsigned int)lrand48();
#else
		r250_buffer[j] = rand();
#endif
#else
		r250_buffer[j] = randlcg();
#endif


	for (j = 0; j < 250; j++)	/* set some MSBs to 1 */
#ifdef TRUST_RAND
		if ( rand() > HALF_RANGE )
			r250_buffer[j] |= MSB;
#else
		if ( randlcg() > HALF_RANGE )
			r250_buffer[j] |= MSB;
#endif


	msb = MSB;	        /* turn on diagonal bit */
	mask = ALL_BITS;	/* turn off the leftmost bits */

	for (j = 0; j < BITS; j++)
	{
		k = STEP * j + 3;	/* select a word to operate on */
		r250_buffer[k] &= mask; /* turn off bits left of the diagonal */
		r250_buffer[k] |= msb;	/* turn on the diagonal bit */
		mask >>= 1;
		msb  >>= 1;
	}

}

unsigned int r250()		/* returns a random unsigned integer */
{
	register int	j;
	register unsigned int new_rand;

	if ( r250_index >= 147 )
		j = r250_index - 147;	/* wrap pointer around */
	else
		j = r250_index + 103;

	new_rand = r250_buffer[ r250_index ] ^ r250_buffer[ j ];
	r250_buffer[ r250_index ] = new_rand;

	if ( r250_index >= 249 )	/* increment pointer for next time */
		r250_index = 0;
	else
		r250_index++;

	return new_rand;

}


double dr250()		/* returns a random double in range 0..1 */
{
	register int	j;
	register unsigned int new_rand;

	if ( r250_index >= 147 )
		j = r250_index - 147;	/* wrap pointer around */
	else
		j = r250_index + 103;

	new_rand = r250_buffer[ r250_index ] ^ r250_buffer[ j ];
	r250_buffer[ r250_index ] = new_rand;

	if ( r250_index >= 249 )	/* increment pointer for next time */
		r250_index = 0;
	else
		r250_index++;

	return (double)new_rand / ALL_BITS;

}

#ifdef MAIN

/* test driver	prints out either NMR_RAND values or a histogram	*/

#include <stdio.h>

#define NMR_RAND	5000
#define MAX_BINS	500

#ifdef NO_PROTO
void main(argc, argv)
int argc;
char **argv;
#else
void main(int argc, char **argv)
#endif
{
	int j,k,nmr_bins,seed;
	int bins[MAX_BINS];
	double randm, bin_inc;
	double bin_limit[MAX_BINS];

	if ( argc != 3 )
	{
		printf("Usage -- %s nmr_bins seed\n", argv[0]);
		exit(1);
	}

	nmr_bins = atoi( argv[1] );
	if ( nmr_bins > MAX_BINS )
	{
		printf("ERROR -- maximum number of bins is %d\n", MAX_BINS);
		exit(1);
	}

	seed = atoi( argv[2] );

	r250_init( seed );

	if ( nmr_bins < 1 )	/* just print out the numbers */
	{
		for (j = 0; j < NMR_RAND; j++)
			printf("%f\n", dr250() );
		exit(0);
	}
	
	bin_inc = 1.0 / nmr_bins;
	for (j = 0; j < nmr_bins; j++)	/* initialize bins to zero */
	{
		bins[j] = 0;
		bin_limit[j] = (j + 1) * bin_inc;
	}

	bin_limit[nmr_bins-1] = 1.0e7;	/* make sure all others are in last bin */

	for (j = 0; j < NMR_RAND; j++)
	{
		randm = r250() / (double)ALL_BITS;
		for (k = 0; k < nmr_bins; k++)
			if ( randm < bin_limit[k] )
			{
				bins[k]++;
				break;
			}
	}


	for (j = 0; j < nmr_bins; j++)
		printf("%d\n", bins[j]);
	
}

#endif

randlcg.c100644   3665     17        2136  6671542625  11310 0ustar  enghoffcnl/* rndlcg            Linear Congruential Method, the "minimal standard generator"
                     Park & Miller, 1988, Comm of the ACM, 31(10), pp. 1192-1201

*/

static char rcsid[] = "@(#)randlcg.c	1.1 15:48:15 11/21/94   EFC";

#include <math.h>
#include <limits.h>

#define ALL_BITS     0xffffffff

static long int quotient  = LONG_MAX / 16807L;
static long int remain    = LONG_MAX % 16807L;

static long int seed_val = 1L;

long set_seed(long int sd)
{
        return seed_val = sd;
}

long get_seed()
{
        return seed_val;
}


unsigned long int randlcg()       /* returns a random unsigned integer */
{
        if ( seed_val <= quotient )
                seed_val = (seed_val * 16807L) % LONG_MAX;
        else
        {
                long int high_part = seed_val / quotient;
                long int low_part  = seed_val % quotient;

                long int test = 16807L * low_part - remain * high_part;

                if ( test > 0 )
                        seed_val = test;
                else
                        seed_val = test + LONG_MAX;

        }

        return seed_val;
}





blas2mkl.h100644   3665     17        4127  6676315102  11405 0ustar  enghoffcnl/******************************************************************************/
/* For all support and information please contact:                            */
/*                                                                            */
/*   Sigurd Enghoff                                                           */
/*   The Salk Institute, CNL                                                  */
/*   enghoff@salk.edu                                                         */
/*                                                                            */
/* Additional ICA software:                                                   */
/*   http://www.cnl.salk.edu/~enghoff/                                        */
/*                                                                            */
/******************************************************************************/

#ifndef BLAS2MKL_INCLUDE
#define BLAS2MKL_INCLUDE

/* BLAS Level 1 */
#define drotg_  DROTG
#define drotmg_ DROTMG
#define drot_   DROT
#define drotm_  DROTM
#define dswap_  DSWAP
#define dscal_  DSCAL
#define dcopy_  DCOPY
#define daxpy_  DAXPY
#define ddot_   DDOT
#define ddotu_  DDOTU
#define ddotc_  DDOTC
#define dnrm2_  DNRM2
#define dasum_  DASUM
#define idamax_ IDAMAX

/* BLAS Level 2 */
#define dgemv_  DGEMV
#define dgbmv_  DGBMV
#define dhemv_  DHEMV
#define dhbmv_  DHBMV
#define dhpmv_  DHPMV
#define dsymv_  DSYMV
#define dsbmv_  DSBMV
#define dspmv_  DSPMV
#define dtrmv_  DTRMV
#define dtbmv_  DTBMV
#define dtpmv_  DTPMV
#define dtrsv_  DTRMV
#define dtbsv_  DTBMV
#define dtpsv_  DTPMV

#define dger_   DGER
#define dgeru_  DGERU
#define dgerc_  DGERC
#define dher_   DHER
#define dhrp_   DHRP
#define dher2_  DHER2
#define dhrp2_  DHRP2
#define dsyr_   DSYR
#define dspr_   DSPR
#define dsyr2_  DSYR2
#define dspr2_  DSPR2

/* BLAS Level 3 */
#define dgemm_  DGEMM
#define dsymm_  DSYMM
#define dhemm_  DHEMM
#define dsyrk_  DSYRK
#define dherk_  DHERK
#define dsyr2k_ DSYR2K
#define dher2k_ DHER2K
#define dtrmm_  DTRMM
#define dtrsm_  DTRSM

#endif
f2c.h120777   3665     17           0  6555473346  15626 2/home/enghoff/CLAPACK/F2CLIBS/f2c.hustar  enghoffcnlica.h100644   3665     17       13752  6701245155  10455 0ustar  enghoffcnl/******************************************************************************/
/* For all support and information please contact:                            */
/*                                                                            */
/*   Sigurd Enghoff                                                           */
/*   The Salk Institute, CNL                                                  */
/*   enghoff@salk.edu                                                         */
/*                                                                            */
/* Additional ICA software:                                                   */
/*   http://www.cnl.salk.edu/~enghoff/                                        */
/*                                                                            */
/******************************************************************************/

#ifndef __defs_h
#define __defs_h

#include "f2c.h"

#define VER_INFO               "Version 1.2  (Dec 31, 1998)"

/*#ifdef PVM
#define PROB_WINDOW
#endif*/

/* If R250 is defined random number are generated using the r250 algorithm    */
/* otherwise standard srand and rand function calls are used                  */
#define R250

/* If KURTOSIS_ESTIMATE is defined sub- and super-Gaussian distributions are  */
/* estimated using:                                                           */
/*    k = E{u^2}^2 / E{x^4} - 3                                               */
/* else:                                                                      */
/*    k = E{sech(u)^2} * E{u^2} - E{tanh(u) * u}                              */

#define KURTOSIS_ESTIMATE

#define MAX_WEIGHT             1e8
#define DEFAULT_STOP           0.000001
#define DEFAULT_ANNEALDEG      60
#define DEFAULT_ANNEALSTEP     0.90
#define DEFAULT_EXTANNEAL      0.98
#define DEFAULT_MAXSTEPS       512
#define DEFAULT_MOMENTUM       0.0

#define DEFAULT_BLOWUP         1000000000.0
#define DEFAULT_BLOWUP_FAC     0.8
#define DEFAULT_RESTART_FAC    0.9

#define MIN_LRATE              0.000001
#define MAX_LRATE              0.1
#define DEFAULT_LRATE(chans)   0.015/log((doublereal)chans)

#define DEFAULT_BLOCK(length)  (integer)sqrt((doublereal)length/3.0)

/* Extended-ICA option: */
#define DEFAULT_EXTENDED       0
#define DEFAULT_EXTBLOCKS      1
#define DEFAULT_NSUB           1

#define DEFAULT_EXTMOMENTUM    0.5
#define MAX_PDFSIZE            6000
#define MIN_PDFSIZE            2000
#define SIGNCOUNT_THRESHOLD    25

#ifdef KURTOSIS_ESTIMATE
#define DEFAULT_SIGNSBIAS      0.02
#else
#define DEFAULT_SIGNSBIAS      0.0
#endif

#define SIGNCOUNT_STEP         2

#define DEFAULT_BIASFLAG       1

#define MIN(x,y) (((x) > (y)) ? (y) : (x))
#define DEGCONST               180.0/ICA_PI
#define ICA_PI                 3.14159265358979324

#define DEFAULT_PCAFLAG        0
#define DEFAULT_SPHEREFLAG     1
#define DEFAULT_POSACT         1
#define DEFAULT_VERBOSE        1

#define NWINDOW                5
#define FRAMEWINDOW            0
#define FRAMESTEP              1
#define EPOCHWINDOW            2
#define EPOCHSTEP              3
#define BASELINE               4

#define EOL 0x0A
/* Fixing the seed causes multiple runs of a given */
/* configuration to yield identical results */
/* #define FIX_SEED */

#ifdef R250
#define SRAND(seed) r250_init(seed)
#define RAND        r250
#else
#define SRAND(seed) srand(seed)
#define RAND        rand
#endif

typedef struct {
    integer    idx;
    doublereal val;
} idxelm;

typedef struct {
    char *token;
    void *prev;
} key;

extern integer extended, extblocks, pdfsize, nsub;
extern integer verbose, block, maxsteps;

extern doublereal lrate, annealstep, annealdeg, nochange, momentum;

extern void error(char*);

extern void zero(integer, doublereal*);
extern void eye(integer, doublereal*);
extern void rmmean(doublereal*, integer, integer);
extern void syproj(doublereal*, doublereal*, integer, integer, doublereal*);
extern void geproj(doublereal*, doublereal*, integer, integer, doublereal*);
extern void pcaproj(doublereal*, doublereal*, integer, integer, integer, doublereal*);
extern void varsort(doublereal*, doublereal*, doublereal*, doublereal*, doublereal*, integer*, integer, integer, integer);
extern void do_sphere(doublereal*, integer, integer, doublereal*);
extern void pca(doublereal*, integer, integer, doublereal*);
extern void posact(doublereal*, doublereal*, integer, integer, doublereal*);
extern void runica(doublereal*, doublereal*, integer, integer, integer, doublereal*, integer*);

extern doublereal dsum_(integer*, doublereal*, integer*);
extern doublereal ddot_(integer*, doublereal*, integer*, doublereal*, integer*);
extern integer idamax_(integer*, doublereal*, integer*);
extern integer ilaenv_(integer*, char*, char*, integer*, integer*, integer*, integer*);
extern void dscal_(integer*, doublereal*, doublereal*, integer*);
extern void dcopy_(integer*, doublereal*, integer*, doublereal*, integer*);
extern void dswap_(integer*, doublereal*, integer*, doublereal*, integer*);
extern void daxpy_(integer*, doublereal*, doublereal*, integer*, doublereal*, integer*);
extern void dgemv_(char*, integer*, integer*, doublereal*, doublereal*, integer*, doublereal*, integer*, doublereal*, doublereal*, integer*);
extern void dgemm_(char*, char*, integer*, integer*, integer*, doublereal*, doublereal*, integer*, doublereal*, integer*, doublereal*, doublereal*, integer*);
extern void dsymm_(char*, char*, integer*, integer*, doublereal*, doublereal*, integer*, doublereal*, integer*, doublereal*, doublereal*, integer*);
extern void dsyrk_(char*, char*, integer*, integer*, doublereal*, doublereal*, integer*, doublereal*, doublereal*, integer*);
extern void dsyev_(char*, char*, integer*, doublereal*, integer*, doublereal*, doublereal*, integer*, integer*);
extern void dgesv_(integer*, integer*, doublereal*, integer*, integer*, doublereal*, integer*, integer*);
extern void dgetri_(integer*, doublereal*, integer*, integer*, doublereal*, integer*, integer*);
extern void dgetrf_(integer*, integer*, doublereal*, integer*, integer*, integer*);

#endif
memap.h100644   3665     17        1766  6570622657  11014 0ustar  enghoffcnl/******************************************************************************/
/* For all support and information please contact:                            */
/*                                                                            */
/*   Sigurd Enghoff                                                           */
/*   The Salk Institute, CNL                                                  */
/*   enghoff@salk.edu                                                         */
/*                                                                            */
/* Additional ICA software:                                                   */
/*   http://www.cnl.salk.edu/~enghoff/                                        */
/*                                                                            */
/******************************************************************************/

#ifndef __memap_h
#define __memap_h

extern void *mapmalloc(int size);
extern void mapfree(void *addr, int size);

#endif
pvmica.h100644   3665     17        6210  6714105771  11152 0ustar  enghoffcnl/******************************************************************************/
/* For all support and information please contact:                            */
/*                                                                            */
/*   Sigurd Enghoff                                                           */
/*   The Salk Institute, CNL                                                  */
/*   enghoff@salk.edu                                                         */
/*                                                                            */
/* Additional ICA software:                                                   */
/*   http://www.cnl.salk.edu/~enghoff/                                        */
/*                                                                            */
/******************************************************************************/

#ifndef __pvmica_h
#define __pvmica_h

#include "pvm3.h"
#include "ica.h"

/* Command to be excuted to spawn new ICA processes */
#define SPAWN_ICA "ica"

/* Nice value used by slave processes */
#define SLAVE_NICE 20

/* Number of integer elements in the assign buffer */
#define ASSIGN_NINT 13

/* Integer elements in the assign buffer */
#define ASSIGN_ID         0
#define ASSIGN_CHANS      1
#define ASSIGN_FRAMES     2
#define ASSIGN_EPOCHS     3
#define ASSIGN_BIAS       4
#define ASSIGN_SIGNS      5
#define ASSIGN_EXTENDED   6
#define ASSIGN_EXTBLOCKS  7
#define ASSIGN_PDFSIZE    8 
#define ASSIGN_NSUB       9
#define ASSIGN_VERBOSE    10 
#define ASSIGN_BLOCK      11
#define ASSIGN_MAXSTEPS   12


/* Number of double elements in the PVM assign buffer */
#define ASSIGN_NDBL       5

/* Double elements in the PVM assign buffer */
#define ASSIGN_LRATE      0
#define ASSIGN_ANNEALSTEP 1
#define ASSIGN_ANNEALDEG  2
#define ASSIGN_NOCHANGE   3
#define ASSIGN_MOMENTUM   4

/* Structure processed by the assign send and receive routines */
typedef struct {
	integer     id;

	doublereal *data;
	doublereal *weights;
	integer     chans;
	integer     frames;
	integer     epochs;
	integer     bias;      /* boolean value */
	integer     signs;     /* boolean value */

	integer     extended;
	integer     extblocks;
	integer     pdfsize;
	integer     nsub;
	integer     verbose;
	integer     block;
	integer     maxsteps;
	doublereal  lrate;
	doublereal  annealstep;
	doublereal  annealdeg;
	doublereal  nochange;
	doublereal  momentum;
} tassign;


/* Number of integer elements in the PVM result buffer */
#define RESULT_NINT  4

/* Integer elements in the PVM result buffer */
#define RESULT_ID    0
#define RESULT_CHANS 1
#define RESULT_BIAS  2
#define RESULT_SIGNS 3


/* Number of double elements in the result buffer */
#define RESULT_NDBL  1

/* Double elements in the result buffer */
#define RESULT_LRATE 0

/* Structure processed by the result send and receive routines */
typedef struct {
	integer     id;

	doublereal *weights;
	integer     chans;
	doublereal *bias;
	integer    *signs;

	doublereal  lrate;
} tresult;

extern int slave();
extern void pvmica(doublereal*, doublereal*, doublereal*, doublereal*, integer, integer, integer, integer, int*, doublereal*, integer*, char**);

#endif
r250.h100644   3665     17        1047  6670057362  10371 0ustar  enghoffcnl/* r250.h	prototypes for r250 random number generator,

		Kirkpatrick, S., and E. Stoll, 1981; "A Very Fast
		Shift-Register Sequence Random Number Generator",
		Journal of Computational Physics, V.40

		also:

		see W.L. Maier, DDJ May 1991


*/

#ifndef _R250_H_
#define _R250_H_ 1.2

#ifdef __cplusplus
extern "C" {
#endif

#ifdef NO_PROTO
void         r250_init();
unsigned int r250();
double      dr250();

#else
void         r250_init(int seed);
unsigned int r250( void );
double       dr250( void );
#endif

#ifdef __cplusplus
}
#endif

#endif
randlcg.h100644   3665     17        1104  6670057362  11305 0ustar  enghoffcnl
/* randlcg.h	prototypes for the minimal standard random number generator,

   Linear Congruential Method, the "minimal standard generator"
   Park & Miller, 1988, Comm of the ACM, 31(10), pp. 1192-1201


  rcsid: @(#)randlcg.h	1.1 15:48:09 11/21/94   EFC

*/

#ifndef _RANDLCG_H_
#define _RANDLCG_H_ 1.1

#ifdef __cplusplus
extern "C" {
#endif

#ifdef NO_PROTO
long         set_seed();
long         get_seed();
unsigned long int randlcg();

#else
long         set_seed(long);
long         get_seed(long);
unsigned long int randlcg();

#endif

#ifdef __cplusplus
}
#endif

#endif
ica.sc100644   3665     17        6606  6702233504  10607 0ustar  enghoffcnl# ica - Perform Independent Component Analysis, standalone-version
#
#   Run the ICA algorithm of Bell & Sejnowski (1996) or the extended-ICA 
#   of Lee, Girolami & Sejnowski (1998). Original Matlab code: Scott Makeig,
#   Tony Bell, et al.; C++ code: Sigurd Enghoff, CNL / Salk Institute 7/98
#
#   Usage:   % ica < my.sc
#
#   Leading # -> use default values
#   Edit a copy of this file to run an ica decomposition
#   Contacts: {enghoff,scott,terry,tony,tewon}@salk.edu

# Required variables:
 
    DataFile     berger/modeldata # Input data to decompose (floats multiplexed
                           #   by channel (i.e., chan1, chan2, ...))
    chans        31        # Number of data channels (= data rows) 
    frames       768       # Number of data points per epoch (= data columns)
    epochs       436       # Number of epochs

#	FrameWindow  20        # Number of frames per window
#	FrameStep    4         # Number of frames to step per window
#	EpochWindow  100       # Number of epochs per window
#	EpochStep    25        # Number of epochs to step per window
#	Baseline     25        # Number of data points contained in baseline

    WeightsOutFile berger/data.wts # Output ICA weight matrix (floats)
    SphereFile   berger/data.sph  # Output sphering matrix (floats)
 
# Processing options:
 
#   sphering     on        # Flag sphering of data (on/off)   {default: on}
#   bias         on        # Perform bias adjustment (on/off) {default: on}
    extended     1         # Perform "extended-ICA" using tnah() with kurtosis
                           #  estimation every N training blocks. If N < 0,
                           #  fix number of sub-Gaussian components to -N 
                           #  {default|0: off}
#   pca          0         # Decompose a principal component subspace of
                           #  the data. Retain this many PCs. {default|0: all}
# Optional input variables:
 
#  WeightsInFile input.wts # Starting ICA weight matrix (nchans,ncomps)
                           #  {default: identity or sphering matrix}
    lrate        2.0e-3    # Initial ICA learning rate (float << 1)
                           #  {default: heuristic ~5e-4}
#   blocksize    20        # ICA block size (integer << datalength) 
                           #  {default: heuristic fraction of log data length}
#   stop         1.0e-6    # Stop training when weight-change < this value
                           #  {default: heuristic ~0.000001}
    maxsteps     512       # Max. number of ICA training steps {default: 128}
#   posact       on        # Make each component activation net-positive
                           # (on/off) {default: on}
#   annealstep   0.98      # Annealing factor (range (0,1]) - controls 
                           #  the speed of convergence.
#   annealdeg    60        # Angledelta threshold for annealing {default: 60}
#   momentum     0.0       # Momentum gain (range [0,1])      {default: 0}
#   verbose      off        # Give ascii messages (on/off) {default: on}
 
# Optional outputs:
 
#  ActivationsFile data.act # Activations of each component (ncomps,points)
#  BiasFile      data.bs   # Bias weights (ncomps,1)
#  SignFile      data.sgn  # Signs designating (-1) sub- and (1) super-Gaussian 
                           #  components (ncomps,1)

# This script, "ica.sc" is a sample ica script file. Copy and modify it as
# desired. Note that the input data file(s) must be native floats.
Makefile120755   3665     17           0  6667546146  13321 2Makefile.bsdustar  enghoffcnlMakefile.alpha100644   3665     17        2157  6677062402  12255 0ustar  enghoffcnlPLAT = _alpha

CC        = cc
CFLAGS    = -fast -DALPHA -DMMAP # -DPVM
BLASLIB   = ./CLAPACK/blas$(PLAT).a
LAPACKLIB = ./CLAPACK/lapack$(PLAT).a
F2CLIB    = ./CLAPACK/F2CLIBS/libF77$(PLAT).a ./CLAPACK/F2CLIBS/libI77$(PLAT).a
PVMLIB    = -L$(PVM_ROOT)/lib/$(PVM_ARCH) -lpvm3
LIBS      = $(LAPACKLIB) $(BLASLIB) $(F2CLIB) $(PVMLIB) -lm
INCLUDE   = -I$(PVM_ROOT)/include
OBJS      = ica.o pvmica.o interfc.o memap.o r250.o randlcg.o dsum.o

ica_alpha: Makefile.alpha ica.h pvmica.h memap.h $(OBJS)
	$(CC) -o ica_alpha $(OBJS) $(CFLAGS) $(LIBS)

ica.o: Makefile.alpha ica.h memap.h ica.c
	$(CC) -c ica.c $(CFLAGS) $(INCLUDE)

pvmica.o: Makefile.alpha ica.h pvmica.h pvmica.c
	$(CC) -c pvmica.c $(CFLAGS) $(INCLUDE)

interfc.o: Makefile.alpha ica.h memap.h interfc.c
	$(CC) -c interfc.c $(CFLAGS) $(INCLUDE)

memap.o: Makefile.alpha memap.h memap.c
	$(CC) -c memap.c $(CFLAGS) $(INCLUDE)

r250.o: Makefile.alpha r250.h r250.c
	$(CC) -c r250.c $(CFLAGS) $(INCLUDE)

randlcg.o: Makefile.alpha randlcg.h randlcg.c
	$(CC) -c randlcg.c $(CFLAGS) $(INCLUDE)

dsum.o: Makefile.alpha dsum.c
	$(CC) -c dsum.c $(CFLAGS) $(INCLUDE)

clean:
	rm -f *.o
Makefile.bsd100644   3665     17        2203  6714633145  11730 0ustar  enghoffcnlPLAT = _bsd

CC        = gcc
CFLAGS    = -O3 -DMMAP -ansi # -DPVM
BLASLIB   = ./CLAPACK/lsblasppro1.1n_01.99.a
#BLASLIB   = ./CLAPACK/blas$(PLAT).a
LAPACKLIB = ./CLAPACK/lapack$(PLAT).a
F2CLIB    = ./CLAPACK/F2CLIBS/libF77$(PLAT).a ./CLAPACK/F2CLIBS/libI77$(PLAT).a
PVMLIB    = -L$(PVM_ROOT)/lib/$(PVM_ARCH) -lpvm3
LIBS      = $(LAPACKLIB) $(BLASLIB) $(F2CLIB) $(PVMLIB) -lm
INCLUDE   = -I$(PVM_ROOT)/include
OBJS      = ica.o pvmica.o interfc.o memap.o r250.o randlcg.o dsum.o

ica_bsd: Makefile.bsd ica.h pvmica.h memap.h $(OBJS)
	$(CC) -o ica_bsd $(OBJS) $(LIBS)

ica.o: Makefile.bsd ica.h pvmica.h memap.h ica.c
	$(CC) -c ica.c $(CFLAGS) $(INCLUDE)

pvmica.o: Makefile.bsd ica.h pvmica.h pvmica.c
	$(CC) -c pvmica.c $(CFLAGS) $(INCLUDE)

interfc.o: Makefile.bsd ica.h memap.h interfc.c
	$(CC) -c interfc.c $(CFLAGS) $(INCLUDE)

memap.o: Makefile.bsd memap.h memap.c
	$(CC) -c memap.c $(CFLAGS) $(INCLUDE)

r250.o: Makefile.bsd r250.h r250.c
	$(CC) -c r250.c $(CFLAGS) $(INCLUDE)

randlcg.o: Makefile.bsd randlcg.h randlcg.c
	$(CC) -c randlcg.c $(CFLAGS) $(INCLUDE)

dsum.o: Makefile.bsd dsum.c
	$(CC) -c dsum.c $(CFLAGS) $(INCLUDE)

clean:
	rm -f *.o
Makefile.sgi100644   3665     17        2104  6714160026  11733 0ustar  enghoffcnlPLAT = _sgi

CC        = cc
CFLAGS    = -O2 -DMMAP -w # -DPVM
BLASLIB   = -lblas
#BLASLIB   = ./CLAPACK/blas$(PLAT).a
LAPACKLIB = ./CLAPACK/lapack$(PLAT).a
F2CLIB    = ./CLAPACK/F2CLIBS/libF77$(PLAT).a
PVMLIB    = -L$(PVM_ROOT)/lib/$(PVM_ARCH) -lpvm3
LIBS      = $(LAPACKLIB) $(BLASLIB) $(F2CLIB) $(PVMLIB) -lm
INCLUDE   = -I$(PVM_ROOT)/include
OBJS      = ica.o interfc.o memap.o r250.o randlcg.o dsum.o

ica_sgi: Makefile.sgi ica.h pvmica.h memap.h $(OBJS)
	$(CC) -o ica_sgi $(OBJS) $(CFLAGS) $(LIBS)

ica.o: Makefile.sgi ica.h pvmica.h memap.h ica.c
	$(CC) -c ica.c $(CFLAGS) $(INCLUDE)

pvmica.o: Makefile.sgi ica.h pvmica.h pvmica.c
	$(CC) -c pvmica.c $(CFLAGS) $(INCLUDE)

interfc.o: Makefile.sgi ica.h memap.h interfc.c
	$(CC) -c interfc.c $(CFLAGS) $(INCLUDE)

memap.o: Makefile.sgi memap.h memap.c
	$(CC) -c memap.c $(CFLAGS) $(INCLUDE)

r250.o: Makefile.sgi r250.h r250.c
	$(CC) -c r250.c $(CFLAGS) $(INCLUDE)

randlcg.o: Makefile.sgi randlcg.h randlcg.c
	$(CC) -c randlcg.c $(CFLAGS) $(INCLUDE)

dsum.o: Makefile.sgi dsum.c
	$(CC) -c dsum.c $(CFLAGS) $(INCLUDE)

clean:
	rm -f *.o
Makefile.sun100644   3665     17        2120  6676312255  11766 0ustar  enghoffcnlPLAT = _sun

CC        = acc
CFLAGS    = -fast -DMMAP -DPVM
BLASLIB   = ./CLAPACK/blas$(PLAT).a
LAPACKLIB = ./CLAPACK/lapack$(PLAT).a
F2CLIB    = ./CLAPACK/F2CLIBS/libF77$(PLAT).a ./CLAPACK/F2CLIBS/libI77$(PLAT).a
PVMLIB    = -L$(PVM_ROOT)/lib/$(PVM_ARCH) -lpvm3
LIBS      = $(LAPACKLIB) $(BLASLIB) $(F2CLIB) $(PVMLIB) -lm
INCLUDE   = -I$(PVM_ROOT)/include
OBJS      = ica.o pvmica.o interfc.o memap.o r250.o randlcg.o dsum.o

ica_sun: Makefile.sun ica.h pvmica.h memap.h $(OBJS)
	$(CC) -o ica_sun $(OBJS) $(CFLAGS) $(LIBS)

ica.o: Makefile.sun ica.h memap.h ica.c
	$(CC) -c ica.c $(CFLAGS) $(INCLUDE)

pvmica.o: Makefile.sun ica.h pvmica.h pvmica.c
	$(CC) -c pvmica.c $(CFLAGS) $(INCLUDE)

interfc.o: Makefile.sun ica.h memap.h interfc.c
	$(CC) -c interfc.c $(CFLAGS) $(INCLUDE)

memap.o: Makefile.sun memap.h memap.c
	$(CC) -c memap.c $(CFLAGS) $(INCLUDE)

r250.o: Makefile.sun r250.h r250.c
	$(CC) -c r250.c $(CFLAGS) $(INCLUDE)

randlcg.o: Makefile.sun randlcg.h randlcg.c
	$(CC) -c randlcg.c $(CFLAGS) $(INCLUDE)

dsum.o: Makefile.sun dsum.c
	$(CC) -c dsum.c $(CFLAGS) $(INCLUDE)

clean:
	rm -f *.o