text:
/*
 *  vbHmmGaussIntensity.c
 *  Model-specific core functions for VB-HMM-GAUSS.
 *
 *  Created by OKAMOTO Kenji, SAKO Yasushi and RIKEN
 *  Copyright 2011-2016
 *  Cellular Informatics Laboratory, Advance Science Institute, RIKEN, Japan.
 *  All rights reserved.
 *
 *  Ver. 1.0.0
 *  Last modified on 2016.02.08
 */

#include "vbHmmGaussIntensity.h"
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_sf_psi.h>
#include <string.h>
#include "rand.h"

#ifdef _OPENMP
#include "omp.h"
#endif

#define  MAX(a,b)  ((a)>(b)?(a):(b))
#define  MIN(a,b)  ((a)<(b)?(a):(b))

static int isGlobalAnalysis = 0;

void setFunctions_gaussInt(){
    commonFunctions funcs;
    funcs.newModelParameters    = newModelParameters_gaussInt;
    funcs.freeModelParameters   = freeModelParameters_gaussInt;
    funcs.newModelStats         = newModelStats_gaussInt;
    funcs.freeModelStats        = freeModelStats_gaussInt;
    funcs.initializeVbHmm       = initializeVbHmm_gaussInt;
    funcs.pTilde_z1             = pTilde_z1_gaussInt;
    funcs.pTilde_zn_zn1         = pTilde_zn_zn1_gaussInt;
    funcs.pTilde_xn_zn          = pTilde_xn_zn_gaussInt;
    funcs.calcStatsVars         = calcStatsVars_gaussInt;
    funcs.maximization          = maximization_gaussInt;
    funcs.varLowerBound         = varLowerBound_gaussInt;
    funcs.reorderParameters     = reorderParameters_gaussInt;
    funcs.outputResults         = outputResults_gaussInt;
    setFunctions( funcs );
}

void setGFunctions_gaussInt(){
    gCommonFunctions funcs;
    funcs.newModelParameters    = newModelParameters_gaussInt;
    funcs.freeModelParameters   = freeModelParameters_gaussInt;
    funcs.newModelStats         = newModelStats_gaussInt;
    funcs.freeModelStats        = freeModelStats_gaussInt;
    funcs.newModelStatsG        = newModelStatsG_gaussInt;
    funcs.freeModelStatsG       = freeModelStatsG_gaussInt;
    funcs.initializeVbHmmG      = initializeVbHmmG_gaussInt;
    funcs.pTilde_z1             = pTilde_z1_gaussInt;
    funcs.pTilde_zn_zn1         = pTilde_zn_zn1_gaussInt;
    funcs.pTilde_xn_zn          = pTilde_xn_zn_gaussInt;
    funcs.calcStatsVarsG        = calcStatsVarsG_gaussInt;
    funcs.maximizationG         = maximizationG_gaussInt;
    funcs.varLowerBoundG        = varLowerBoundG_gaussInt;
    funcs.reorderParametersG    = reorderParametersG_gaussInt;
    funcs.outputResultsG        = outputResultsG_gaussInt;
    setGFunctions( funcs );
    isGlobalAnalysis = 1;
}

void outputResults_gaussInt( xn, gv, iv, logFP )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
FILE *logFP;
{
    outputGaussIntResults( xn, gv, iv, logFP );
}

void outputResultsG_gaussInt( xns, gv, ivs, logFP )
xnDataBundle *xns;
globalVars *gv;
indVarBundle *ivs;
FILE *logFP;
{
    outputGaussIntResultsG( xns, gv, ivs, logFP );
}

void *newModelParameters_gaussInt( xn, sNo )
xnDataSet *xn;
int sNo;
{
    int i;
    gaussIntParameters *p = (gaussIntParameters*)malloc( sizeof(gaussIntParameters) );

    p->uPiArr = (double*)malloc( sNo * sizeof(double) );
    p->sumUPi = 0.0;
    p->uAMat = (double**)malloc( sNo * sizeof(double*) );
    for( i = 0 ; i < sNo ; i++ ){
        p->uAMat[i] = (double*)malloc( sNo * sizeof(double) );
    }
    p->sumUAArr = (double*)malloc( sNo * sizeof(double) );

    p->avgPi = (double *)malloc( sNo * sizeof(double) );
    p->avgLnPi = (double *)malloc( sNo * sizeof(double) );
    p->avgA = (double **)malloc( sNo * sizeof(double*) );
    p->avgLnA = (double **)malloc( sNo * sizeof(double*) );
    for( i = 0 ; i < sNo ; i++ ){
        p->avgA[i] = (double *)malloc( sNo * sizeof(double) );
        p->avgLnA[i] = (double *)malloc( sNo * sizeof(double) );
    }

    return p;
}

void freeModelParameters_gaussInt( p, xn, sNo )
void **p;
xnDataSet *xn;
int sNo;
{
    gaussIntParameters *gp = *p;
    int i;

    free( gp->uPiArr );
    for( i = 0 ; i < sNo ; i++ ){
        free( gp->uAMat[i] );
    }
    free( gp->uAMat );
    free( gp->sumUAArr );

    free( gp->avgPi );
    free( gp->avgLnPi );
    for( i = 0 ; i < sNo ; i++ ){
        free( gp->avgA[i] );
        free( gp->avgLnA[i] );
    }
    free( gp->avgA );
    free( gp->avgLnA );

    free( gp );
    *p = NULL;
}

void *newModelStats_gaussInt( xn, gv, iv )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
{
    if( isGlobalAnalysis == 0 ){
        int sNo = gv->sNo;
        gaussIntStats *s = (gaussIntStats*)malloc( sizeof(gaussIntStats) );

        int i;
        s->Ni = (double *)malloc( sNo * sizeof(double) );
        s->Nij = (double **)malloc( sNo * sizeof(double*) );
        for( i = 0 ; i < sNo ; i++ )
        {   s->Nij[i] = (double *)malloc( sNo * sizeof(double) );   }
        s->Nii = (double *)malloc( sNo * sizeof(double) );

        return s;

    } else {

        return NULL;

    }
}

void freeModelStats_gaussInt( s, xn, gv, iv )
void **s;
xnDataSet *xn;
globalVars *gv;
indVars *iv;
{
    if( isGlobalAnalysis == 0 ){
        int sNo = gv->sNo;
        gaussIntStats *gs = *s;
        int i;

        free( gs->Ni );
        for( i = 0 ; i < sNo ; i++ )
        {   free( gs->Nij[i] );   }
        free( gs->Nij );
        free( gs->Nii );

        free( gs );
        *s = NULL;
    }
}

void *newModelStatsG_gaussInt( xns, gv, ivs)
xnDataBundle *xns;
globalVars *gv;
indVarBundle *ivs;
{
    int sNo = gv->sNo;
    gaussIntGlobalStats *gs = (gaussIntGlobalStats*)malloc( sizeof(gaussIntGlobalStats) );

    int i;
    gs->NijR = (double **)malloc( sNo * sizeof(double*) );
    for( i = 0 ; i < sNo ; i++ )
    {   gs->NijR[i] = (double *)malloc( sNo * sizeof(double) );   }
    gs->NiiR = (double *)malloc( sNo * sizeof(double) );
    gs->NiR = (double *)malloc( sNo * sizeof(double) );
    gs->z1iR = (double *)malloc( sNo * sizeof(double) );

    return gs;
}

void freeModelStatsG_gaussInt( gs, xns, gv, ivs )
void **gs;
xnDataBundle *xns;
globalVars *gv;
indVarBundle *ivs;
{
    int sNo = gv->sNo;
    gaussIntGlobalStats *ggs = *gs;
    int i;
    free( ggs->NiR );
    for( i = 0 ; i < sNo ; i++ )
    {   free( ggs->NijR[i] );   }
    free( ggs->NijR );
    free( ggs->NiiR );
    free( ggs->z1iR );

    free( *gs );
    *gs = NULL;
}

void initializeVbHmm_gaussInt( xn, gv, iv )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
{
    gaussIntData *d = xn->data;
    size_t dLen = xn->N;
    int sNo = gv->sNo;
    gaussIntParameters *p = gv->params;

    int i, j;
    double totalX = 0.0, precX, varX = 0.0;
    for( i = 0 ; i < dLen ; i++ ){
        totalX += d->v[i];
    }
    double meanX = totalX / (double)dLen;
    for( i = 0 ; i < dLen ; i++ ){
        varX += pow(d->v[i] - meanX, 2.0);
    }
    varX /= (double)(dLen - 1);
    precX = 1.0 / varX;

    // hyper parameter for p( pi(i) )
    p->sumUPi = 0.0;
    for( i = 0 ; i < sNo ; i++ ){
        p->uPiArr[i] = 1.0;
        p->sumUPi += p->uPiArr[i];
    }

    // hyper parameter for p( A(i,j) )
    for( i = 0 ; i < sNo ; i++ ){
        p->sumUAArr[i] = 0.0;
        for( j = 0 ; j < sNo ; j++ ){
            if( j == i ){
                p->uAMat[i][j] = 5.0;
            } else {
                p->uAMat[i][j] = 1.0;
            }
            p->sumUAArr[i] += p->uAMat[i][j];
        }
    }

    // hyper parameter for p( mu(k), lm(k) )
    for( i = 0 ; i < sNo ; i++ ){
        p->uBt = 1.0;
        p->uMu = meanX * 2.0 / (double)sNo;
        p->uA  = 1.0;
        p->uB  = p->uA / precX ;
    }

    initialize_indVars_gaussInt( xn, gv, iv );

    calcStatsVars_gaussInt( xn, gv, iv );
    maximization_gaussInt( xn, gv, iv );

}

void initializeVbHmmG_gaussInt( xns, gv, ivs )
xnDataBundle *xns;
globalVars *gv;
indVarBundle *ivs;
{
    gaussIntData *d;
    size_t dLen, totalN;
    int sNo = gv->sNo, rNo = xns->R;
    gaussIntParameters *p = gv->params;

    int i, j, r;
    double totalX, precX, varX, meanX;

    totalX = 0.0;
    totalN = 0;
    for( r = 0 ; r < rNo ; r++ ){
        d = xns->xn[r]->data;
        dLen = xns->xn[r]->N;
        for( i = 0 ; i < dLen ; i++ ){
            totalX += d->v[i];
        }
        totalN += dLen;
    }
    meanX = totalX / (double)totalN;
    varX = 0.0;
    for( r = 0 ; r < rNo ; r++ ){
        d = xns->xn[r]->data;
        dLen = xns->xn[r]->N;
        for( i = 0 ; i < dLen ; i++ ){
            varX += pow(d->v[i] - meanX, 2.0);
        }
    }
    varX /= (double)(totalN - 1);
    precX = 1.0 / varX;

    p->sumUPi = 0.0;
    for( i = 0 ; i < sNo ; i++ ){
        p->uPiArr[i] = 1.0;
        p->sumUPi += p->uPiArr[i];
    }

    for( i = 0 ; i < sNo ; i++ ){
        p->sumUAArr[i] = 0.0;
        for( j = 0 ; j < sNo ; j++ ){
            if( j == i ){
                p->uAMat[i][j] = 5.0;
            } else {
                p->uAMat[i][j] = 1.0;
            }
            p->sumUAArr[i] += p->uAMat[i][j];
        }
    }

    // hyper parameter for p( mu(k), lm(k) )
    for( i = 0 ; i < sNo ; i++ ){
        p->uBt = 1.0;
        p->uMu = meanX * 2.0 / (double)sNo;
        p->uA  = 1.0;
        p->uB  = p->uA / precX ;
    }

    for( r = 0 ; r < rNo ; r++ ){
        initialize_indVars_gaussInt( xns->xn[r], gv, ivs->indVars[r] );
    }

    calcStatsVarsG_gaussInt( xns, gv, ivs );
    maximizationG_gaussInt( xns, gv, ivs );
}

void initialize_indVars_gaussInt( xn, gv, iv )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
{
    size_t dLen = xn->N;
    int sNo = gv->sNo;
    double **gmMat = iv->gmMat;

    int i;
    size_t n;
    double sumPar;
    for( n = 0 ; n < dLen ; n++ ){
        sumPar = 0.0;
        for( i = 0 ; i < sNo ; i++ ){
            gmMat[n][i] = enoise(1.0) + 1.0;
            sumPar += gmMat[n][i];
        }
        for( i = 0 ; i < sNo ; i++ ){
            gmMat[n][i] /= sumPar;
        }
    }
}

xnDataSet *newXnDataSet_gaussInt( filename )
const char *filename;
{
    xnDataSet *xn = (xnDataSet*)malloc( sizeof(xnDataSet) );
    xn->name = (char*)malloc( strlen(filename) + 2 );
    strncpy( xn->name, filename, strlen(filename)+1 );
    xn->data = (gaussIntData*)malloc( sizeof(gaussIntData) );
    gaussIntData *d = (gaussIntData*)xn->data;
    d->v = NULL;
    return xn;
}

void freeXnDataSet_gaussInt( xn )
xnDataSet **xn;
{
    gaussIntData *d = (gaussIntData*)(*xn)->data;
    free( d->v );
    free( (*xn)->data );
    free( (*xn)->name );
    free( *xn );
    *xn = NULL;
}

double pTilde_z1_gaussInt( i, params )
int i;
void *params;
{
    gaussIntParameters *p = (gaussIntParameters*)params;
    return exp( p->avgLnPi[i] );
}

double pTilde_zn_zn1_gaussInt( i, j, params )
int i, j;
void *params;
{
    gaussIntParameters *p = (gaussIntParameters*)params;
    return exp( p->avgLnA[i][j] );
}

double pTilde_xn_zn_gaussInt( xnWv, n, i, params )
xnDataSet *xnWv;
size_t n;
int i;
void *params;
{
    gaussIntParameters *p = (gaussIntParameters*)params;
    gaussIntData *xn = (gaussIntData*)xnWv->data;
    double val, di = (double)(i+1);
    val  = p->avgLnLm - log(2.0 * M_PI * di);
    val -= di/p->btMu + p->aLm * pow(xn->v[n] - di * p->mu0, 2.0) / di / p->bLm;
    return exp(val / 2.0);
}

void calcStatsVars_gaussInt( xn, gv, iv )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
{
    gaussIntData *d = (gaussIntData*)xn->data;
    gaussIntStats *s = (gaussIntStats*)iv->stats;
    size_t dLen = xn->N;
    int sNo = gv->sNo;
    double **gmMat = iv->gmMat, ***xiMat = iv->xiMat;
    double *Nii = s->Nii, **Nij = s->Nij, *Ni = s->Ni;
    size_t n;
    int i, j;

    s->N0 = (double)xn->N;
    s->N0i = 1e-10;
    s->Nx  = 1e-10;
    for( i = 0 ; i < sNo ; i++ ){
        Ni[i]   = 1e-10;
        Nii[i]  = 1e-10;
        for( j = 0 ; j < sNo ; j++ ){
            Nij[i][j] = 1e-10;
        }

        for( n = 0 ; n < dLen ; n++ ){
            Ni[i] += gmMat[n][i];
            s->Nx += gmMat[n][i] * d->v[n];
            for( j = 0 ; j < sNo ; j++ ){
                Nii[i]    += xiMat[n][i][j];
                Nij[i][j] += xiMat[n][i][j];
            }
        }
        s->N0i += (double)(i+1) * Ni[i];
    }
}

void calcStatsVarsG_gaussInt( xns, gv, ivs )
xnDataBundle *xns;
globalVars *gv;
indVarBundle *ivs;
{
    gaussIntData *d;
    gaussIntGlobalStats *gs = (gaussIntGlobalStats*)ivs->stats;
    int sNo = gv->sNo, rNo = xns->R;
    double **gmMat, ***xiMat;
    double *NiiR = gs->NiiR, **NijR = gs->NijR, *NiR = gs->NiR, *z1iR = gs->z1iR;
    size_t dLen, n;
    int i, j, r;

    gs->N0R = 0.0;
    gs->N0iR = 1e-10;
    gs->NxR  = 1e-10;
    for( i = 0 ; i < sNo ; i++ ){
        NiR[i]   = 1e-10;
        NiiR[i]  = 1e-10;
        for( j = 0 ; j < sNo ; j++ ){
            NijR[i][j] = 1e-10;
        }
        z1iR[i] = 1e-10;
    }
    for( r = 0 ; r < rNo ; r++ ){
        d = (gaussIntData*)xns->xn[r]->data;
        dLen = xns->xn[r]->N;
        gmMat = ivs->indVars[r]->gmMat;
        xiMat = ivs->indVars[r]->xiMat;
        for( i = 0 ; i < sNo ; i++ ){
            z1iR[i] += gmMat[0][i];
            for( n = 0 ; n < dLen ; n++ ){
                NiR[i]  += gmMat[n][i];
                gs->NxR += gmMat[n][i] * d->v[n];
                for( j = 0 ; j < sNo ; j++ ){
                    NiiR[i]    += xiMat[n][i][j];
                    NijR[i][j] += xiMat[n][i][j];
                }
            }
        }
        gs->N0R += (double)xns->xn[r]->N;
    }
    for( i = 0 ; i < sNo ; i++ ){
        gs->N0iR += (double)(i+1) * NiR[i];
    }
}

void maximization_gaussInt( xn, gv, iv )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
{
    gaussIntParameters *p = (gaussIntParameters*)gv->params;
    gaussIntStats *s = (gaussIntStats*)iv->stats;
    gaussIntData *d = (gaussIntData*)xn->data;
    size_t dLen = xn->N;
    int sNo = gv->sNo;
    double **gmMat = iv->gmMat;
    double *uPiArr = p->uPiArr, sumUPi = p->sumUPi;
    double **uAMat = p->uAMat, *sumUAArr = p->sumUAArr;
    double uBt = p->uBt, uMu = p->uMu, uA = p->uA, uB = p->uB;
    double *avgPi = p->avgPi, *avgLnPi = p->avgLnPi, **avgA = p->avgA, **avgLnA = p->avgLnA;
    double *Nii = s->Nii, **Nij = s->Nij, N0 = s->N0, N0i = s->N0i, Nx = s->Nx;
    int i, j, m;

    double bLm = 0.0, di;
    for( i = 0 ; i < sNo ; i++ ){
        avgPi[i] = ( uPiArr[i] + gmMat[0][i] ) / ( sumUPi + 1.0 );
        avgLnPi[i] = gsl_sf_psi( uPiArr[i] + gmMat[0][i] ) - gsl_sf_psi( sumUPi + 1.0 );

        for( j = 0 ; j < sNo ; j++ ){
            avgA[i][j] = ( uAMat[i][j] + Nij[i][j] ) / ( sumUAArr[i] + Nii[i] );
            avgLnA[i][j] = gsl_sf_psi( uAMat[i][j] + Nij[i][j] ) - gsl_sf_psi( sumUAArr[i] + Nii[i] );
        }

        di = (double)(i+1);
        for( m = 0 ; m < dLen ; m++ ){
            bLm += gmMat[m][i] * d->v[m] * d->v[m] / di;
        }
    }
    bLm /= 2.0;

    p->btMu = uBt + N0i;
    p->mu0  = (uBt * uMu + Nx) / p->btMu;
    p->aLm  = uA + (N0 / 2.0);
    bLm += uB + (uBt * uMu * uMu / 2.0);
    bLm -= p->btMu * p->mu0 * p->mu0 / 2.0;
    p->bLm = bLm;

    p->avgMu    = p->mu0;
    p->avgLm    = p->aLm / p->bLm;
    p->avgLnLm  = gsl_sf_psi( p->aLm ) - log( p->bLm );
}

void maximizationG_gaussInt( xns, gv, ivs )
xnDataBundle *xns;
globalVars *gv;
indVarBundle *ivs;
{
    gaussIntParameters *p = (gaussIntParameters*)gv->params;
    double *uPiArr = p->uPiArr, sumUPi = p->sumUPi;
    double **uAMat = p->uAMat, *sumUAArr = p->sumUAArr;
    double uBt = p->uBt, uMu = p->uMu, uA = p->uA, uB = p->uB;
    double *avgPi = p->avgPi, *avgLnPi = p->avgLnPi, **avgA = p->avgA, **avgLnA = p->avgLnA;
    gaussIntGlobalStats *gs = (gaussIntGlobalStats*)ivs->stats;
    double *NiiR = gs->NiiR, **NijR = gs->NijR, N0R = gs->N0R, N0iR = gs->N0iR, NxR = gs->NxR;
    double *z1iR = gs->z1iR, dR = (double)(xns->R), di, bLm = 0.0;
    int sNo = gv->sNo, rNo = xns->R;
    int i, j, m, r;

    for( i = 0 ; i < sNo ; i++ ){
        avgPi[i] = ( uPiArr[i] + z1iR[i] ) / ( sumUPi + dR );
        avgLnPi[i] = gsl_sf_psi( uPiArr[i] + z1iR[i] ) - gsl_sf_psi( sumUPi + dR );

        for( j = 0 ; j < sNo ; j++ ){
            avgA[i][j] = ( uAMat[i][j] + NijR[i][j] ) / ( sumUAArr[i] + NiiR[i] );
            avgLnA[i][j] = gsl_sf_psi( uAMat[i][j] + NijR[i][j] ) - gsl_sf_psi( sumUAArr[i] + NiiR[i] );
        }

        di = (double)(i+1);
        for( r = 0 ; r < rNo ; r++ ){
            size_t dLen = xns->xn[r]->N;
            gaussIntData *d = (gaussIntData*)xns->xn[r]->data;
            double **gmMat = ivs->indVars[r]->gmMat;
            for( m = 0 ; m < dLen ; m++ ){
                bLm += gmMat[m][i] * d->v[m] * d->v[m] / di;
            }
        }

    }
    bLm /= 2.0;

    p->btMu = uBt + N0iR;
    p->mu0  = (uBt * uMu + NxR) / p->btMu;
    p->aLm  = uA + N0R / 2.0;
    bLm += uB + uBt * uMu * uMu / 2.0;
    bLm -= p->btMu * p->mu0 * p->mu0 / 2.0;
    p->bLm = bLm;

    p->avgMu    = p->mu0;
    p->avgLm    = p->aLm / p->bLm;
    p->avgLnLm  = gsl_sf_psi( p->aLm ) - log( p->bLm );
}

double varLowerBound_gaussInt( xn, gv, iv )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
{
    gaussIntParameters *p = (gaussIntParameters*)gv->params;
    size_t dLen = xn->N;
    int sNo = gv->sNo;
    double **gmMat = iv->gmMat, *cn = iv->cn;
    double *uPiArr = p->uPiArr, sumUPi = p->sumUPi;
    double **uAMat = p->uAMat, *sumUAArr = p->sumUAArr;
    double uBt = p->uBt, uMu = p->uMu, uA = p->uA, uB = p->uB;
    double *avgLnPi = p->avgLnPi, **avgLnA = p->avgLnA;
    double avgLm = p->avgLm, avgLnLm = p->avgLnLm;
    double mu0 = p->mu0, btMu = p->btMu, aLm = p->aLm, bLm = p->bLm;
    gaussIntStats *s = (gaussIntStats*)iv->stats;
    double *Nii = s->Nii, **Nij = s->Nij;
    size_t n;
    int i, j;

    double lnpPi = gsl_sf_lngamma(sumUPi);
    double lnpA = 0.0;
    double lnpMuLm = log(uBt) / 2.0;
    lnpMuLm -= uBt * (1.0 / btMu + aLm * pow(mu0 - uMu, 2.0) / bLm ) / 2.0;
    lnpMuLm += - gsl_sf_lngamma(uA) + uA * log(uB);
    lnpMuLm += (uA - 0.5) * avgLnLm - uB * avgLm;

    double lnqPi = gsl_sf_lngamma(sumUPi + 1.0);
    double lnqA = 0.0;
    double lnqMuLm = (log(btMu) - 1.0) / 2.0 - gsl_sf_lngamma(aLm) + aLm * log(bLm);
    lnqMuLm += (aLm - 0.5) * avgLnLm - aLm;

    for( i = 0 ; i < sNo ; i++ ){
        lnpPi += (uPiArr[i]-1.0) * avgLnPi[i] - gsl_sf_lngamma(uPiArr[i]);

        lnqPi += (uPiArr[i]+gmMat[0][i]-1.0) * (gsl_sf_psi(uPiArr[i]+gmMat[0][i]) - gsl_sf_psi(sumUPi+1.0));
        lnqPi -= gsl_sf_lngamma(uPiArr[i] + gmMat[0][i]);

        lnpA += gsl_sf_lngamma(sumUAArr[i]);
        lnqA += gsl_sf_lngamma(sumUAArr[i] + Nii[i]);
        for( j = 0 ; j < sNo ; j++ ){
            lnpA += (uAMat[i][j]-1.0)*avgLnA[i][j] - gsl_sf_lngamma(uAMat[i][j]);

            lnqA += (uAMat[i][j] + Nij[i][j] - 1.0) * (gsl_sf_psi(uAMat[i][j]+Nij[i][j]) - gsl_sf_psi(sumUAArr[i]+Nii[i]));
            lnqA -= gsl_sf_lngamma( uAMat[i][j] + Nij[i][j] );
        }
    }

    double lnpX = 0.0;
    for( n = 0 ; n < dLen ; n++ ){
        lnpX += log( cn[n] );
    }

    double val;
    val  = lnpPi + lnpA + lnpMuLm;
    val -= lnqPi + lnqA + lnqMuLm;
    val += lnpX;
    val += log(gsl_sf_fact(sNo));

    return val;
}

double varLowerBoundG_gaussInt( xns, gv, ivs )
xnDataBundle *xns;
globalVars *gv;
indVarBundle *ivs;
{
    gaussIntParameters *p = (gaussIntParameters*)gv->params;
    double *uPiArr = p->uPiArr, sumUPi = p->sumUPi;
    double **uAMat = p->uAMat, *sumUAArr = p->sumUAArr;
    double uBt = p->uBt, uMu = p->uMu, uA = p->uA, uB = p->uB;
    double *avgLnPi = p->avgLnPi, **avgLnA = p->avgLnA;
    double avgLm = p->avgLm, avgLnLm = p->avgLnLm;
    double mu0 = p->mu0, btMu = p->btMu, aLm = p->aLm, bLm = p->bLm;
    gaussIntGlobalStats *gs = (gaussIntGlobalStats*)ivs->stats;
    double *NiiR = gs->NiiR, **NijR = gs->NijR, *z1iR = gs->z1iR, dR = (double)xns->R;
    size_t n;
    int sNo = gv->sNo, rNo = xns->R;
    int i, j, r;

    double lnpPi = gsl_sf_lngamma(sumUPi);
    double lnpA = 0.0;
    double lnpMuLm = log(uBt) / 2.0;
    lnpMuLm -= uBt * (1.0 / btMu + aLm * pow(mu0 - uMu, 2.0) / bLm ) / 2.0;
    lnpMuLm += - gsl_sf_lngamma(uA) + uA * log(uB);
    lnpMuLm += (uA - 0.5) * avgLnLm - uB * avgLm;

    double lnqPi = gsl_sf_lngamma(sumUPi + dR);
    double lnqA = 0.0;
    double lnqMuLm = (log(btMu) - 1.0) / 2.0 - gsl_sf_lngamma(aLm) + aLm * log(bLm);
    lnqMuLm += (aLm - 0.5) * avgLnLm - aLm;

    for( i = 0 ; i < sNo ; i++ ){
        lnpPi += (uPiArr[i]-1.0) * avgLnPi[i] - gsl_sf_lngamma(uPiArr[i]);

        lnqPi += (uPiArr[i]+z1iR[i]-1.0) * (gsl_sf_psi(uPiArr[i]+z1iR[i]) - gsl_sf_psi(sumUPi+dR));
        lnqPi -= gsl_sf_lngamma(uPiArr[i] + z1iR[i]);

        lnpA += gsl_sf_lngamma(sumUAArr[i]);
        lnqA += gsl_sf_lngamma(sumUAArr[i] + NiiR[i]);
        for( j = 0 ; j < sNo ; j++ ){
            lnpA += (uAMat[i][j]-1.0)*avgLnA[i][j] - gsl_sf_lngamma(uAMat[i][j]);

            lnqA += (uAMat[i][j] + NijR[i][j] - 1.0) * (gsl_sf_psi(uAMat[i][j]+NijR[i][j]) - gsl_sf_psi(sumUAArr[i]+NiiR[i]));
            lnqA -= gsl_sf_lngamma( uAMat[i][j] + NijR[i][j] );
        }
    }

    double lnpX = 0.0;
    for( r = 0 ; r < rNo ; r++ ){
        size_t dLen = xns->xn[r]->N;
        for( n = 0 ; n < dLen ; n++ ){
            lnpX += log( ivs->indVars[r]->cn[n] );
        }
    }

    double val;
    val  = lnpPi + lnpA + lnpMuLm;
    val -= lnqPi + lnqA + lnqMuLm;
    val += lnpX;
    val += log(gsl_sf_fact(sNo));

    return val;
}

void reorderParameters_gaussInt( xn, gv, iv )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
{
}

void reorderParametersG_gaussInt( xns, gv, ivs )
xnDataBundle *xns;
globalVars *gv;
indVarBundle *ivs;
{
}

void outputGaussIntResults( xn, gv, iv, logFP )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
FILE *logFP;
{
    gaussIntParameters *p = (gaussIntParameters*)gv->params;
    int sNo = gv->sNo;

    int i, j;
    fprintf(logFP, "  results: K = %d \n", sNo);

    fprintf(logFP, "   mean: ( %g ) \n", p->avgMu);

    fprintf(logFP, "   lambda: ( %g ) \n", p->avgLm);

    fprintf(logFP, "   pi: ( %g", p->avgPi[0]);
    for( i = 1 ; i < sNo ; i++ ){
        fprintf(logFP, ", %g", p->avgPi[i]);
    }
    fprintf(logFP, " ) \n");

    fprintf(logFP, "   A_matrix: [");
    for( i = 0 ; i < sNo ; i++ ){
            fprintf(logFP, " ( %g", p->avgA[i][0]);
            for( j = 1 ; j < sNo ; j++ )
            {   fprintf(logFP, ", %g", p->avgA[i][j]);   }
            fprintf(logFP, ")");
    }
    fprintf(logFP, " ] \n\n");

    char fn[256];
    FILE *fp;
    size_t n;

    sprintf( fn, "%s.param%03d", xn->name, sNo );
    if( (fp = fopen( fn, "w")) != NULL ){
        fprintf(fp, "mu, lambda, pi");
        for( i = 0 ; i < sNo ; i++ )
        {   fprintf(fp, ", A%dx", i);   }
        fprintf(fp, "\n");

        for( i = 0 ; i < sNo ; i++ ){
            fprintf(fp, "%g, %g, %g", (double)(i+1)*p->avgMu, (double)(i+1)*p->avgLm, p->avgPi[i]);
            for( j = 0 ; j < sNo ; j++ )
            {   fprintf(fp, ", %g", p->avgA[j][i]);   }
            fprintf(fp, "\n");
        }
        fclose(fp);
    }

    sprintf( fn, "%s.Lq%03d", xn->name, sNo );
    if( (fp = fopen( fn, "w")) != NULL ){
        for( n = 0 ; n < gv->iteration ; n++ ){
            fprintf( fp, "%24.20e\n", gv->LqArr[n] );
        }
        fclose(fp);
    }

    sprintf( fn, "%s.maxS%03d", xn->name, sNo );
    if( (fp = fopen( fn, "w")) != NULL ){
        for( n = 0 ; n < xn->N ; n++ ){
            fprintf( fp, "%d\n", iv->stateTraj[n] );
        }
        fclose(fp);
    }

}

void outputGaussIntResultsG( xns, gv, ivs, logFP )
xnDataBundle *xns;
globalVars *gv;
indVarBundle *ivs;
FILE *logFP;
{
    int sNo = gv->sNo, rNo = xns->R;
    gaussIntParameters *p = (gaussIntParameters*)gv->params;
    int i, j, r;

    fprintf(logFP, "  results: K = %d \n", sNo);

    fprintf(logFP, "   mean: ( %g ) \n", p->avgMu);

    fprintf(logFP, "   lambda: ( %g ) \n", p->avgLm);

    fprintf(logFP, "   pi: ( %g", p->avgPi[0]);
    for( i = 1 ; i < sNo ; i++ ){
        fprintf(logFP, ", %g", p->avgPi[i]);
    }
    fprintf(logFP, " ) \n");

    fprintf(logFP, "   A_matrix: [");
    for( i = 0 ; i < sNo ; i++ ){
        fprintf(logFP, " ( %g", p->avgA[i][0]);
        for( j = 1 ; j < sNo ; j++ )
        {   fprintf(logFP, ", %g", p->avgA[i][j]);   }
        fprintf(logFP, ")");
    }
    fprintf(logFP, " ] \n\n");

    char fn[256];
    FILE *fp;
    size_t n;

    sprintf( fn, "%s.param%03d", xns->xn[0]->name, sNo );
    if( (fp = fopen( fn, "w")) != NULL ){
        fprintf(fp, "mu, lambda, pi");
        for( i = 0 ; i < sNo ; i++ )
        {   fprintf(fp, ", A%dx", i);   }
        fprintf(fp, "\n");

        for( i = 0 ; i < sNo ; i++ ){
            fprintf(fp, "%g, %g, %g", (double)(i+1)*p->avgMu, (double)(i+1)*p->avgLm, p->avgPi[i]);
            for( j = 0 ; j < sNo ; j++ )
            {   fprintf(fp, ", %g", p->avgA[j][i]);   }
            fprintf(fp, "\n");
        }
        fclose(fp);
    }

    sprintf( fn, "%s.Lq%03d", xns->xn[0]->name, sNo );
    if( (fp = fopen( fn, "w")) != NULL ){
        for( n = 0 ; n < gv->iteration ; n++ ){
            fprintf( fp, "%24.20e\n", gv->LqArr[n] );
        }
        fclose(fp);
    }

    sprintf( fn, "%s.maxS%03d", xns->xn[0]->name, sNo );
    int flag = 0;
    if( (fp = fopen( fn, "w")) != NULL ){
        n = 0;
        do{
            flag = 1;
            for( r = 0 ; r < rNo ; r++ ){
                xnDataSet *xn = xns->xn[r];
                indVars *iv = ivs->indVars[r];

                if( r > 0 ){
                    fprintf( fp, "," );
                }
                if( n < xn->N ){
                    fprintf( fp, "%d", iv->stateTraj[n] );
                }
                flag &= (n >= (xn->N - 1));
            }
            fprintf( fp, "\n" );
            n++;
        }while( !flag );
        fclose(fp);
    }
}

//

max_stars_repo_name: okamoto-kenji/varBayes-HMM

max_issues_repo_name: okamoto-kenji/varBayes-HMM

max_forks_repo_name: okamoto-kenji/varBayes-HMM


text:
/* ieee-utils/fp-gnuc99.c
 *
 * Copyright (C) 2003, 2004, 2007 Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#define _GNU_SOURCE 1

#include <math.h>
#include <stdio.h>
#include <fenv.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_errno.h>

int
gsl_ieee_set_mode (int precision, int rounding, int exception_mask)
{
  int mode;

  switch (precision)
    {
    case GSL_IEEE_SINGLE_PRECISION:
      GSL_ERROR ("single precision rounding is not supported by <fenv.h>",
                 GSL_EUNSUP) ;
      break ;
    case GSL_IEEE_DOUBLE_PRECISION:
      GSL_ERROR ("double precision rounding is not supported by <fenv.h>",
                 GSL_EUNSUP) ;
      break ;
    case GSL_IEEE_EXTENDED_PRECISION:
      GSL_ERROR ("extended precision rounding is not supported by <fenv.h>",
                 GSL_EUNSUP) ;
      break ;
    }

  switch (rounding)
    {
    case GSL_IEEE_ROUND_TO_NEAREST:
#ifdef FE_TONEAREST
      fesetround (FE_TONEAREST) ;
#else
      GSL_ERROR ("round-to-nearest is not supported by <fenv.h>", GSL_EUNSUP) ;
#endif
      break ;
    case GSL_IEEE_ROUND_DOWN:
#ifdef FE_DOWNWARD
      fesetround (FE_DOWNWARD) ;
#else
      GSL_ERROR ("round-down is not supported by <fenv.h>", GSL_EUNSUP) ;
#endif
      break ;
    case GSL_IEEE_ROUND_UP:
#ifdef FE_UPWARD
      fesetround (FE_UPWARD) ;
#else
      GSL_ERROR ("round-up is not supported by <fenv.h>", GSL_EUNSUP) ;
#endif
      break ;
    case GSL_IEEE_ROUND_TO_ZERO:
#ifdef FE_TOWARDZERO
      fesetround (FE_TOWARDZERO) ;
#else
      GSL_ERROR ("round-toward-zero is not supported by <fenv.h>", GSL_EUNSUP) ;
#endif
      break ;
    default:
#ifdef FE_TONEAREST
      fesetround (FE_TONEAREST) ;
#else
      GSL_ERROR ("default round-to-nearest mode is not supported by <fenv.h>", GSL_EUNSUP) ;
#endif
    }

  /* Turn on all the exceptions apart from 'inexact' */

  mode = 0;

#ifdef FE_INVALID
  mode |= FE_INVALID;
#endif

#ifdef FE_DIVBYZERO
  mode |= FE_DIVBYZERO;
#endif

#ifdef FE_OVERFLOW
  mode |= FE_OVERFLOW ;
#endif

#ifdef FE_UNDERFLOW
  mode |= FE_UNDERFLOW ;
#endif

  if (exception_mask & GSL_IEEE_MASK_INVALID)
    {
#ifdef FE_INVALID
    mode &= ~ FE_INVALID ;
#else
    GSL_ERROR ("invalid operation exception not supported by <fenv.h>",
               GSL_EUNSUP);
#endif
    }

  if (exception_mask & GSL_IEEE_MASK_DENORMALIZED)
    {
      /* do nothing */
    }
  else
    {
      GSL_ERROR ("denormalized operand exception not supported by <fenv.h>. "
                 "Use 'mask-denormalized' to work around this.", GSL_EUNSUP) ;
    }

  if (exception_mask & GSL_IEEE_MASK_DIVISION_BY_ZERO)
    {
#ifdef FE_DIVBYZERO
      mode &= ~ FE_DIVBYZERO ;
#else
      GSL_ERROR ("division by zero exception not supported by <fenv.h>",
                 GSL_EUNSUP);
#endif
    }

  if (exception_mask & GSL_IEEE_MASK_OVERFLOW)
    {
#ifdef FE_OVERFLOW
      mode &= ~ FE_OVERFLOW ;
#else
      GSL_ERROR ("overflow exception not supported by <fenv.h>", GSL_EUNSUP);
#endif
    }

  if (exception_mask & GSL_IEEE_MASK_UNDERFLOW)
    {
#ifdef FE_UNDERFLOW
      mode &=  ~ FE_UNDERFLOW ;
#else
      GSL_ERROR ("underflow exception not supported by <fenv.h>", GSL_EUNSUP);
#endif
    }

  if (exception_mask & GSL_IEEE_TRAP_INEXACT)
    {
#ifdef FE_INEXACT
      mode |= FE_INEXACT ;
#else
      GSL_ERROR ("inexact exception not supported by <fenv.h>", GSL_EUNSUP);
#endif
    }
  else
    {
#ifdef FE_INEXACT
      mode &= ~ FE_INEXACT ;
#else
      /* do nothing */
#endif
    }

#if HAVE_DECL_FEENABLEEXCEPT
  feenableexcept (mode) ;
#elif HAVE_DECL_FESETTRAPENABLE
  fesettrapenable (mode);
#else
  GSL_ERROR ("unknown exception trap method", GSL_EUNSUP)
#endif

  return GSL_SUCCESS ;
}

ext: c

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
/* Modified from SUPOLAR, Copyright (c) Colorado School of Mines, 2011.*/
/* All rights reserved.                                                */

/* SUPOLAR_PS: $Revision: 1.9 $ ; $Date: Sat, 23 Jul 2016 00:10:26 -0700 */

#include <omp.h>
#include <stdio.h>
#include <gsl/gsl_statistics_float.h>
#include "su.h"
#include "segy.h"
#include "header.h"

/*********************** self documentation *****************************/
char *sdoc[] = {
"                                                                       ",
" SUPOLAR_PS - POLarization analysis of three-component data            ",
"            - to provide P and S arrivals                              ",
"            - Modification of SUPOLAR                                  ",
"                                                                       ",
" supolar_PS <stdin [optional parameters]                               ",
"                                                                       ",
" Required parameters:                                                  ",
"    -ntr (number of traces) needs to be specified in the SU file       ",
"                                                                       ",
" Optional parameters:                                                  ",
"    dt=(from header)  time sampling interval in seconds                ",
"    wl=0.1            correlation window length in seconds             ",
"    win=boxcar        correlation window shape, choose \"boxcar\",     ",
"                      \"hanning\", \"bartlett\", or \"welsh\"          ",
"    file=polar        base of output file name(s)                      ",
"    rl=1              1 = rectilinearity evaluating 2 eigenvalues,     ",
"                      2, 3 = rectilinearity evaluating 3 eigenvalues   ",
"    rlq=1.0           contrast parameter for rectilinearity            ",
"    dir=1             1 = 3 components of direction of polarization    ",
"                          (the only three-component output file)       ",
"    tau=0             1 = global polarization parameter                ",
"    ellip=0           1 = principal, subprincipal, and transverse      ",
"                          ellipticities e21, e31, and e32              ",
"    pln=0             1 = planarity measure                            ",
"    f1=0              1 = flatness or oblateness coefficient           ",
"    l1=0              1 = linearity coefficient                        ",
"    amp=0             1 = amplitude parameters: instantaneous,         ",
"                          quadratic, and eigenresultant ir, qr, and er ",
"    theta=0           1, 2, 3 = incidence angle of principal axis      ",
"    phi=0             1, 2, 3 = horizontal azimuth of principal axis   ",
"    angle=rad         unit of angles theta and phi, choose \"rad\",    ",
"                      \"deg\", or \"gon\"                              ",
"    all=0             1, 2, 3 = set all output flags to that value     ",
"    verbose=0         1 = echo additional information                  ",
"    kwl=5*(1/dt)      kurtosis window length, samples                  ",
"                                                                       ",
" Notes:                                                                ",
"    Three adjacent traces are considered as one three-component        ",
"    dataset.                                                           ",
"    Correct calculation of angles theta and phi requires the first of  ",
"    these traces to be the vertical component, followed by the two     ",
"    horizontal components (e.g. Z, N, E, or Z, inline, crossline).     ",
"    Significant signal energy on Z is necessary to resolve the 180 deg ",
"    ambiguity of phi (options phi=2,3 only).                           ",
"                                                                       ",
"    Each calculated polarization attribute is written into its own     ",
"    SU file. These files get the same base name (set with \"file=\")   ",
"    and the parameter flag as an extension (e.g. polar.rl).            ",
"                                                                       ",
"    In case of a tapered correlation window, the window length wl may  ",
"    have to be increased compared to the boxcar case, because of their ",
"    smaller effective widths (Bartlett, Hanning: 1/2, Welsh: 1/3).     ",
"                                                                       ",
" Range of values:                                                      ",
"    parameter     option  interval                                     ",
"    rl            1, 2    0.0 ... 1.0   (1.0: linear polarization)     ",
"    rl            3      -1.0 ... 1.0                                  ",
"    tau, l1       1       0.0 ... 1.0   (1.0: linear polarization)     ",
"    pln, f1       1       0.0 ... 1.0   (1.0: planar polarization)     ",
"    e21, e31, e32 1       0.0 ... 1.0   (0.0: linear polarization)     ",
"    theta         1      -pi/2... pi/2  rad                            ",
"    theta         2, 3    0.0 ... pi/2  rad                            ",
"    phi           1      -pi/2... pi/2  rad                            ",
"    phi           2      -pi  ... pi    rad   (see notes above)        ",
"    phi           3       0.0 ... 2 pi  rad   (see notes above)        ",
"                                                                       ",
"                                                                       ",
NULL};

/*
 * Author: Nils Maercklin,
 *         GeoForschungsZentrum (GFZ) Potsdam, Germany, 1998-2001.
 *         E-mail: nils@gfz-potsdam.de
 *
 *
 * References:
 *    Jurkevics, A., 1988: Polarization analysis of three-component
 *       array data. Bulletin of the Seismological Society of America,
 *       vol. 78, no. 5.
 *    Kanasewich, E. R., 1981: Time Sequence Analysis in Geophysics.
 *       The University of Alberta Press.
 *    Kanasewich, E. R., 1990: Seismic Noise Attenuation.
 *       Handbook of Geophysical Exploration, Pergamon Press, Oxford.
 *    Meyer, J. H. 1988: First Comparative Results of Integral and
 *       Instantaneous Polarization Attributes for Multicomponent Seismic
 *       Data. Institut Francais du Petrole.
 *    Press, W. H., Teukolsky, S. A., Vetterling, W. T., and Flannery, B. P.
 *       1996: Numerical Recipes in C - The Art of Scientific Computing.
 *       Cambridge University Press, Cambridge.
 *    Samson, J. C., 1973: Description of the Polarisation States of Vector
 *       Processes: Application to ULF Electromagnetic Fields.
 *       Geophysical Journal vol. 34, p. 403-419.
 *    Sheriff, R. E., 1991: Encyclopedic Dictionary of Exploration
 *       Geophysics. 3rd ed., Society of Exploration Geophysicists, Tulsa.
 *
 * Trace header fields accessed: ns, dt
 * Trace header fields modified: none
 */
/**************** end self doc *******************************************/

/* prototypes of functions used internally */
float calc_ellip(float *d, int d1, int d2);
float calc_er(float *d);
float calc_f1(float *d);
float calc_l1(float *d);
float calc_phi(float **v, int opt);
float calc_plan(float *d);
float calc_rl(float d[], float q, int opt);
float calc_tau(float *d);
float calc_theta(float **v, int opt);
float covar(float *data1, float *data2, int istart, int iwl, float *w);
float calc_pfilt(float rl, float theta);
float calc_sfilt(float rl, float theta);
float kurtosiswindow(float *data, float *data_kwl, int it, int kwl, int nt);
void ampparams(float **indata, float *data_ir, float *data_qr, int nt, int iwl);
void calc_dir(float **data3c_dir, float **v, int it);
void calc_window(float *w, int iwl, int iwin);
void fputdata(FILE *fileptr, FILE *headerptr, float outdata[], int nt);
void fputdata3c(FILE *fileptr, FILE *headerptr, float **outdata3c, int nt);

/* window shape identifiers */
#define WBOXCAR 0
#define WBARTLETT 1
#define WHANNING 2
#define WWELSH 3

segy tr;                /* SEG-Y record (trace and header) */

int main(int argc, char **argv)
{
    /* output file pointers */
    FILE *rlfp=NULL, *taufp=NULL, *e21fp=NULL;
    FILE *e31fp=NULL, *e32fp=NULL, *plnfp=NULL;
    FILE *f1fp=NULL, *l1fp=NULL;
    FILE *thetafp=NULL, *phifp=NULL, *dirfp=NULL;
    FILE *erfp=NULL, *irfp=NULL, *qrfp=NULL;
    FILE *headerfp=NULL, *pfiltfp=NULL, *sfiltfp=NULL;
    FILE *nfiltfp=NULL, *efiltfp=NULL;
 //   FILE *pkur=NULL, *skur=NULL;
        /* temporary file for trace headers */
			 /* (one 3C station only) */

    char *file=NULL;         /* base of output file name(s) */
    char *fname=NULL;        /* complete output file name */
    char *angle=NULL;        /* unit used for angles theta and phi */
    char *win=NULL;          /* shape of used time window */
    float fangle=0.0;   /* unit conversion factor applied to angles theta and phi */
    int iwin=0;           /* time window shape identifier */

    /* flags (see selfdoc) */
    int rl,theta,phi,tau,ellip,pln,f1,l1,dir,amp,verbose,all;

    int i,j,icomp;      /* indices for components (in loops) */
    int it;             /* index for time sample in main loop */
    int iwl;            /* correlation window length in samples */
    int nstat;          /* number of 3-component datasets */
    int nt;             /* number of time samples in one trace */
//    int kwl;            /* kurtosis window length in seconds */

    float **data3c;     /* three-component data ([1..3][0..nt-1]) */
    float **a;          /* covariance matrix (a[1..3][1..3]) */
    float **v;          /* eigenvectors of covariance matrix (v[1..3][1..3]) */
    float *d;           /* the corresponding eigenvalues (d[1..3]) */
    float *w;           /* time window weights for correlation window */
    float dt;           /* sampling interval in seconds */
    float rlq;          /* contrast factor of rectilinearity */
    float wl;           /* correlation window length in seconds */

    float *data_e21=NULL;    /* main ellipticity */
    float *data_e31=NULL;    /* second ellipticity */
    float *data_e32=NULL;    /* transverse ellipticity */
    float *data_er=NULL;     /* eigenresultant */
    float *data_f1=NULL;     /* flatness coefficient */
    float *data_ir=NULL;     /* instantaneous resultant */
    float *data_l1=NULL;     /* linearity coefficient */
    float *data_phi=NULL;    /* horizontal azimuth phi */
    float *data_pln=NULL;    /* planarity */
    float *data_qr=NULL;     /* quadratic resultant */
    float *data_rl=NULL;     /* rectilinearity factor */
    float *data_tau=NULL;    /* polarization parameter tau */
    float *data_theta=NULL;  /* inclination angle theta */
    float *data_pfilt=NULL;  /* P (vertical) polarization filter */
    float *data_sfilt=NULL;  /* S (horizontal) polarization filter */
    float *data_zfilt=NULL;  /* Z Filtered Trace */
    float *data_nfilt=NULL;  /* N Filtered Trace */
    float *data_efilt=NULL;  /* E Filtered Trace */
//    float *data_kwl=NULL;   /* Data for Kurtosis Window */
//    float *data_pkur=NULL;   /* Kurtosis detector for P */
//    float *data_skur=NULL;   /* Kurtosis detector for S */
    float **data3c_dir=NULL; /* 3 components of direction of polarization ([1..3][0..nt-1]) */

    /* initialize */
    initargs(argc, argv);
    requestdoc(1);

    /* get info from first trace */
    if(!gettr(&tr)) err("can't get first trace");
    nt = tr.ns;

    /* get parameters ... */
    if (!getparstring("file", &file)) file="polar";
    if (!getparstring("angle", &angle)) angle="rad";
    if (!getparstring("win", &win)) win="boxcar";
    if (!getparfloat("wl", &wl)) wl = 0.1;
    if (!getparfloat("dt", &dt)) dt = ((double) tr.dt)/1000000.0;
    if (!getparfloat("rlq", &rlq)) rlq = 1.0;
    if (!getparint("verbose", &verbose)) verbose = 0;
//    if (!getparint("kwl", &kwl)) kwl = 5 * ((int) 1/dt);

    /* ... and output flags */
    if (!getparint("all", &all))       all = 0;
    if (!getparint("rl", &rl))          rl = (all) ? all : 1;
    if (!getparint("dir", &dir))       dir = (all) ? 1 : 1;
    if (!getparint("theta", &theta)) theta = (all) ? all : 0;
    if (!getparint("phi", &phi))       phi = (all) ? all : 0;
    if (!getparint("tau", &tau))       tau = (all) ? 1 : 0;
    if (!getparint("ellip", &ellip)) ellip = (all) ? 1 : 0;
    if (!getparint("pln", &pln))       pln = (all) ? 1 : 0;
    if (!getparint("f1", &f1))          f1 = (all) ? 1 : 0;
    if (!getparint("l1", &l1))          l1 = (all) ? 1 : 0;
    if (!getparint("amp", &amp))       amp = (all) ? 1 : 0;

    checkpars();

    /* get time window shape */
    if      (STREQ(win, "boxcar"))   iwin=WBOXCAR;
    else if (STREQ(win, "bartlett")) iwin=WBARTLETT;
    else if (STREQ(win, "hanning"))  iwin=WHANNING;
    else if (STREQ(win, "welsh"))    iwin=WWELSH;
    else err("unknown win=%s", win);

    /* get unit conversion factor for angles */
    if      (STREQ(angle, "rad")) fangle=1.0;
    else if (STREQ(angle, "deg")) fangle=180.0/PI;
    else if (STREQ(angle, "gon")) fangle=200.0/PI;
    else err("unknown angle=%s", angle);

    /* convert seconds to samples */
    if (!dt) {
        dt = 0.004;
        warn("dt not set, assuming dt=0.004");
    }

    iwl = NINT(wl/dt);

    /* data validation */
    if (iwl<1) err("wl=%g must be positive", wl);
    if (iwl>nt) err("wl=%g too long for trace", wl);
    if (!strlen(file)) err("file= not set and default overridden");

    /* echo some information */
    if (verbose && (theta || phi)) warn("computing angles in %s", angle);
    if (verbose) warn("%s window length = %d samples\n", win, iwl);

    if (rl && theta) warn("computing filtered phase");

    /* open temporary file for trace headers */
    headerfp = etmpfile();

    /* set filenames and open files */
    fname = malloc( strlen(file)+7 );
    sprintf(fname, "%s.rl", file);    if (rl) rlfp = efopen(fname, "w");
    sprintf(fname, "%s.theta", file); if (theta) thetafp = efopen(fname, "w");
    sprintf(fname, "%s.phi", file);   if (phi) phifp = efopen(fname, "w");
    sprintf(fname, "%s.tau", file);   if (tau) taufp = efopen(fname, "w");
    sprintf(fname, "%s.e21", file);   if (ellip) e21fp = efopen(fname, "w");
    sprintf(fname, "%s.e31", file);   if (ellip) e31fp = efopen(fname, "w");
    sprintf(fname, "%s.e32", file);   if (ellip) e32fp = efopen(fname, "w");
    sprintf(fname, "%s.pln", file);   if (pln) plnfp = efopen(fname, "w");
    sprintf(fname, "%s.f1", file);    if (f1) f1fp = efopen(fname, "w");
    sprintf(fname, "%s.l1", file);    if (l1) l1fp = efopen(fname, "w");
    sprintf(fname, "%s.dir", file);   if (dir) dirfp = efopen(fname, "w");
    sprintf(fname, "%s.er", file);    if (amp) erfp = efopen(fname, "w");
    sprintf(fname, "%s.ir", file);    if (amp) irfp = efopen(fname, "w");
    sprintf(fname, "%s.qr", file);    if (amp) qrfp = efopen(fname, "w");
    sprintf(fname, "%s.pfilt", file); if (rl && theta) pfiltfp = efopen(fname, "w");
    sprintf(fname, "%s.sfilt", file); if (rl && theta) sfiltfp = efopen(fname, "w");
    sprintf(fname, "%s.nfilt", file); if (rl && theta) nfiltfp = efopen(fname, "w");
    sprintf(fname, "%s.efilt", file); if (rl && theta) efiltfp = efopen(fname, "w");
  //  sprintf(fname, "%s.pkur", file); if (rl && theta) pkur = efopen(fname, "w");
  //  sprintf(fname, "%s.skur", file); if (rl && theta) skur = efopen(fname, "w");
    free(fname);

    /* allocate space for input data and analysis matrices */
    /* index ranges used here: data3c[1..3][0..nt-1],      */
    /* a[1..3][1..3], v[1..3][1..3], d[1..3]               */
    data3c = ealloc2float(nt,3); data3c-=1;
    a = ealloc2float(3,3); a[0]-=1; a-=1;
    v = ealloc2float(3,3); v[0]-=1; v-=1;
    d = ealloc1float(3); d-=1;

    /* calculate time window weights */
    w = ealloc1float(iwl);
    memset((void *) w, 0, iwl*FSIZE);
    calc_window(w, iwl, iwin);

    /* allocate and zero out space for output data */
    if (rl) {
        data_rl = ealloc1float(nt);
        memset((void *) data_rl, 0, nt*FSIZE);
    }
    if (theta) {
        data_theta = ealloc1float(nt);
        memset((void *) data_theta, 0, nt*FSIZE);
    }
    if (phi) {
        data_phi = ealloc1float(nt);
        memset((void *) data_phi, 0, nt*FSIZE);
    }
    if (tau) {
        data_tau = ealloc1float(nt);
        memset((void *) data_tau, 0, nt*FSIZE);
    }
    if (ellip) {
        data_e21 = ealloc1float(nt);
        data_e31 = ealloc1float(nt);
        data_e32 = ealloc1float(nt);
        memset((void *) data_e21, 0, nt*FSIZE);
        memset((void *) data_e31, 0, nt*FSIZE);
        memset((void *) data_e32, 0, nt*FSIZE);
    }
    if (pln) {
        data_pln = ealloc1float(nt);
        memset((void *) data_pln, 0, nt*FSIZE);
    }
    if (f1) {
        data_f1 = ealloc1float(nt);
        memset((void *) data_f1, 0, nt*FSIZE);
    }
    if (l1) {
        data_l1 = ealloc1float(nt);
        memset((void *) data_l1, 0, nt*FSIZE);
    }
    if (amp) {
        data_er = ealloc1float(nt);
        data_ir = ealloc1float(nt);
        data_qr = ealloc1float(nt);
        memset((void *) data_er, 0, nt*FSIZE);
        memset((void *) data_ir, 0, nt*FSIZE);
        memset((void *) data_qr, 0, nt*FSIZE);
    }
    if (dir) {
        data3c_dir = ealloc2float(nt,3); data3c_dir-=1;
        for (i=1;i<=3;i++) memset((void *) data3c_dir[i], 0, nt*FSIZE);
    }
    if (rl && theta) {
        data_pfilt = ealloc1float(nt);
        memset((void *) data_pfilt, 0, nt*FSIZE);
        data_sfilt = ealloc1float(nt);
        memset((void *) data_pfilt, 0, nt*FSIZE);
/*        data_3Cfilt = ealloc2float(nt,3);
        for (i=1;i<=3;i++) memset((void *) data_3Cfilt[i], 0, nt*FSIZE); */
        data_zfilt = ealloc1float(nt);
        data_nfilt = ealloc1float(nt);
        data_efilt = ealloc1float(nt);
        memset((void *) data_zfilt, 0, nt*FSIZE);
        memset((void *) data_nfilt, 0, nt*FSIZE);
        memset((void *) data_efilt, 0, nt*FSIZE);

 //       data_pkur = ealloc1float(nt);
//        data_skur = ealloc1float(nt);
//        memset((void *) data_pkur, 0, nt*FSIZE);
//        memset((void *) data_skur, 0, nt*FSIZE);
        /* Allocate data for kurtosis window arrays */
//        data_kwl = ealloc1float(iwl);
//        memset((void *) data_kwl, 0, kwl*FSIZE);
    }

/* ************************ BEGIN CALCULATION ******************************* */

    /* loop over traces */
    icomp=0;
    nstat=0;
    // Need to convert this do while loop into a for loop so as to be easier
    // to parallelize
    warn("Trace Start Time: %d %d %d %d %d", tr.year, tr.day, tr.hour, tr.minute, tr.sec);
    do {
        /* store trace header in temporary file and read data */
        efwrite(&tr, HDRBYTES, 1, headerfp);
        icomp++;
        memcpy((void *)data3c[icomp], (const void *) tr.data, nt*FSIZE);

        /* process 3-component dataset */
        if (icomp==3) {
            erewind(headerfp);
            icomp = 0;
            nstat++;
            if (verbose)
                fprintf(stderr,"%s: analyzing station %d \r",argv[0], nstat);

            /* start loop over samples */

            for (it=iwl/2;it<nt-iwl/2;it++) {
                //warn("Sample %d", it);
                /* covariance matrix */

                for (i=1;i<=3;i++)
                {
                    for (j=i;j<=3;j++)
                    {
                        a[i][j]=a[j][i]=covar(data3c[i], data3c[j], it-iwl/2, iwl, w);
                    }
                }

                /* compute eigenvalues and vectors */
                eig_jacobi(a,d,v,3);
                sort_eigenvalues(d,v,3);

                /* polarization parameters */
                if (rl) data_rl[it]=calc_rl(d,rlq,rl);
                if (theta) data_theta[it]=calc_theta(v, theta) * fangle;
                if (phi) data_phi[it]=calc_phi(v, phi) * fangle;
                if (tau) data_tau[it]=calc_tau(d);
                if (ellip) {
                    data_e21[it]=calc_ellip(d,2,1);
                    data_e31[it]=calc_ellip(d,3,1);
                    data_e32[it]=calc_ellip(d,3,2);
                }
                if (pln) data_pln[it]=calc_plan(d);
                if (f1) data_f1[it]=calc_f1(d);
                if (l1) data_l1[it]=calc_l1(d);
                if (amp) data_er[it]=calc_er(d);
                if (dir) calc_dir(data3c_dir,v,it);
                if (rl && theta)
                {
                    data_zfilt[it] = data3c[1][it] * calc_pfilt(rl, theta);
                    data_nfilt[it] = data3c[2][it] * calc_sfilt(rl, theta);
                    data_efilt[it] = data3c[3][it] * calc_sfilt(rl, theta);
                    data_pfilt[it] = data_zfilt[it];
                    data_sfilt[it] = (data_nfilt[it] + data_efilt[it]) / 2;
//                    data_pkur[it] = kurtosiswindow(data_pfilt,data_kwl,it - kwl/2,kwl,nt);
//                    data_skur[it] = kurtosiswindow(data_sfilt,data_kwl,it - kwl/2,kwl,nt);
               }

            } /* end loop over samples */

            /* compute amplitude parameters */
            if (amp) ampparams(data3c, data_ir, data_qr, nt, iwl);

/* *************************** END CALCULATION ****************************** */
/* ***************************** BEGIN WRITE ******************************** */

            /* write polarization attributes to files */
            if (rl) fputdata(rlfp, headerfp, data_rl, nt);
            if (theta) fputdata(thetafp, headerfp, data_theta, nt);
            if (phi) fputdata(phifp, headerfp, data_phi, nt);
            if (tau) fputdata(taufp, headerfp, data_tau, nt);
            if (ellip) {
                fputdata(e21fp, headerfp, data_e21, nt);
                fputdata(e31fp, headerfp, data_e31, nt);
                fputdata(e32fp, headerfp, data_e32, nt);
            }
            if (pln) fputdata(plnfp, headerfp, data_pln, nt);
            if (f1) fputdata(f1fp, headerfp, data_f1, nt);
            if (l1) fputdata(l1fp, headerfp, data_l1, nt);
            if (amp) {
                fputdata(erfp, headerfp, data_er, nt);
                fputdata(irfp, headerfp, data_ir, nt);
                fputdata(qrfp, headerfp, data_qr, nt);
            }
            if (dir) fputdata3c(dirfp, headerfp, data3c_dir, nt);
            if (rl && theta)
            {
                fputdata(pfiltfp, headerfp, data_pfilt, nt);
                fputdata(sfiltfp, headerfp, data_sfilt, nt);
                fputdata(nfiltfp, headerfp, data_nfilt, nt);
                fputdata(efiltfp, headerfp, data_efilt, nt);
               // fputdata(pkur, headerfp, data_pkur, nt);
               // fputdata(skur, headerfp, data_skur, nt);
            }

/* ****************************** END WRITE ********************************* */

        } /* end of processing three-component dataset */

    } while (gettr(&tr)); /* end loop over traces */

    if (verbose) {
        fprintf(stderr,"\n");
        if (icomp) warn("last %d trace(s) skipped", icomp);
    }

    /* close files */
    efclose(headerfp);
    if (rl) efclose(rlfp);
    if (theta) efclose(thetafp);
    if (phi) efclose(phifp);
    if (tau) efclose(taufp);
    if (ellip) {
        efclose(e21fp);
        efclose(e31fp);
        efclose(e32fp);
    }
    if (pln) efclose(plnfp);
    if (f1) efclose(f1fp);
    if (l1) efclose(l1fp);
    if (amp) {
        efclose(erfp);
        efclose(irfp);
        efclose(qrfp);
    }
    if (dir) efclose(dirfp);
    if (rl && theta) {
        efclose(pfiltfp);
        efclose(sfiltfp);
//        efclose(pkur);
//        efclose(skur);
        }
    return(CWP_Exit());
}

/**********************************************************************/
/* Functions used internally                                          */
/**********************************************************************/

/* calculate time window weights for a smooth time window,        */
/* after Sheriff, 1991, p. 335; tapered windows are suggested by  */
/* Jurkevics, 1988                                                */

void calc_window(float *w, int iwl, int iwin)
{
    int i,j;      /* loop indices within window */
    float m;      /* half of time window (iwl/2) */

    m = (float) iwl/2;

    switch (iwin) {
        case WBOXCAR :
            for (i=-iwl/2,j=0;i<iwl/2;i++,j++) {
                w[j] = 1.0;
            }
            if (j<iwl) w[iwl-1] = w[0];
            break;
        case WBARTLETT :
            for (i=-iwl/2,j=0;i<iwl/2;i++,j++) {
                w[j] = 1.0 - (fabs((float) i) / m);

            }
            if (j<iwl) w[iwl-1] = w[0];
            break;
        case WHANNING :
            for (i=-iwl/2,j=0;i<iwl/2;i++,j++) {
                w[j] = 0.5 + 0.5*cos(PI*fabs((float) i) / m);
            }
            if (j<iwl) w[iwl-1] = w[0];
            break;
        case WWELSH :
            for (i=-iwl/2,j=0;i<iwl/2;i++,j++) {
                w[j] = fabs( fabs((float) i) / m - 1.0);
                w[j] *= w[j];
            }
            if (j<iwl) w[iwl-1] = w[0];
            break;
    }
}

/* covariance of two time sequences *data1 and *data2,   */
/* calculated in a window of length iwl,                 */
/* weighted with factors *w and starting at index istart */

float covar(float *data1, float *data2, int istart, int iwl, float *w)
{
    register int i, j;
    float cov=0.0;
    float mean1=0.0;
    float mean2=0.0;

    for (i=istart,j=0;i<(istart+iwl);i++,j++) {
        mean1 += data1[i] * w[j];
        mean2 += data2[i] * w[j];
    }
    mean1 = mean1 / (float) iwl;
    mean2 = mean2 / (float) iwl;

    for (i=istart,j=0;i<(istart+iwl);i++,j++) {
	/* modified by Mega G. Baitoff 6 May 2015 */
        /* cov += (data1[i]-mean1) * (data2[i]-mean2) * w[j]; */
        cov += (data1[i] * w[j] - mean1) * (data2[i] * w[j] - mean2);
    }
    cov = cov / (float) iwl;
    return cov;
}

float kurtosiswindow(float *data, float *data_kwl, int it, int kwl, int nt)
    {
        int i, j;
        for (i=it, j=0; i<(it+kwl); i++, j++)
        {
            data_kwl[j] = data[i];
        }
        return gsl_stats_float_kurtosis(data_kwl,1,kwl);
    }

/**********************************************************************/
/* Functions defining polarization attributes and angles              */
/**********************************************************************/

/* rectilinearity factor rl (3 different definitions) */

float calc_rl(float *d, float rlq, int opt)
{
    float rl;

    if (d[1]!=0.0) {
        switch (opt) {
            default:
                /* case 1: */
                /* rl definition after Kanasewich, 1981 */
                rl = 1 - pow(fabs(d[2]/d[1]), rlq);
                break;
            case 2:
                /* rl definition after Jurkevics, 1988 */
                /* case used by Z. Ross et. al., 2014 */
                rl = 1 - pow( 0.5*(fabs(d[2]/d[1]) + fabs(d[3]/d[1]) ), rlq);
                break;
            case 3:
                /* rl definition after Meyer, 1988 */
                rl = 1 - pow(fabs(d[2]/d[1]) + fabs(d[3]/d[1]), rlq);
                break;
        }
        return rl;
    }
    else {
        return 0.0;
    }
}

/* vertical (incidence) angle theta (different definitions) */
/* assumes a right handed coordinate system starting        */
/* with the vertical component Z (Z=1, N,R=2, E,T=3)        */

float calc_theta(float **v, int opt)
{
    float theta, horiz;

    switch (opt) {
        default:
            /* case 1: */
            /* definition after Kanasewich, 1981 */
            /* interval -pi/2 <= theta <= pi/2 */
            if (v[1][1]) {
                horiz = sqrt( v[2][1]*v[2][1] + v[3][1]*v[3][1] );
                theta = atan( horiz / v[1][1] );
            }
            else {
                theta = 0.0;
            }
            break;
        case 2:
        case 3:
            /* definition after Jurkevics, 1988 */
            /* interval 0.0 <= theta <= pi/2 */
            /* case used by Z. Ross et. al., 2014 */
            theta = acos( fabs(v[1][1]) );
            break;
    }
    return theta;
}

/* horizontal azimuth phi (different definitions)    */
/* assumes a right handed coordinate system starting */
/* with the vertical component Z (Z=1, N,R=2, E,T=3) */

/* Note: The SIGN function is introduced to resolve  */
/* the 180 deg ambiguity by taking the positive      */
/* vertical component of v[1][1] (Jurkevics, 1988).  */
#define VSIGN ( (v[1][1]<0.0) ? -1.0 : 1.0 )

float calc_phi(float **v, int opt)
{
    float phi;

    switch (opt) {
        default:
            /* case 1: */
            /* definition after Kanasewich, 1981 */
            /* interval -pi/2 <= phi <= pi/2 */
            if (v[1][2]) {
                phi = atan( v[3][1] / v[2][1] );
            }
            else {
                phi = (v[3][1]>0.0) ? 0.5*PI : -0.5*PI;
            }
            break;
        case 2:
        case 3:
            /* definitions after Jurkevics, 1988 */
            /* interval -pi <= phi <= pi */
            if (v[2][1]) {
                phi = atan2( v[3][1]*VSIGN, v[2][1]*VSIGN);
            }
            else {
                phi = (v[3][1]>0.0) ? 0.5*PI*VSIGN : -0.5*PI*VSIGN;
            }

            /* interval 0.0 <= phi <= 2*pi */
            if (phi<0.0 && opt==3) phi += 2.0*PI;
            break;
    }
    return phi;
}
#undef VSIGN

/* global polarization parameter tau (Samson, 1973) */

float calc_tau(float *d)
{
    float x1, x2, x3, x4, tau;

    if (d[1]) {
        x1 = pow(( 1 - fabs(d[2]/d[1])), 2.);
        x2 = pow(( 1 - fabs(d[3]/d[1])), 2.);
        x3 = pow(( fabs(d[2]/d[1]) - fabs(d[3]/d[1]) ), 2.);
        x4 = pow(( 1 + fabs(d[2]/d[1]) + fabs(d[3]/d[1])), 2.);
        tau = sqrt( (x1+x2+x3) / (2*x4) );
        return tau;
    }
    else
        return 0.0;
}

/* ellipticities e_ik */

float calc_ellip(float *d, int d1, int d2)
{
    float ellip;

    if (d[d2]) {
        ellip = sqrt( fabs(d[d1]/d[d2]) );
        return ellip;
    }
    else
        return 0.0;
}

/* planarity after Jurkevics, 1988 */

float calc_plan(float *d)
{
    float pln;

    if (d[1]+d[2]) {
        pln = 1.0 - 2.0*d[3] / (d[1] + d[2]);
        return pln;
    }
    else
        return 0.0;
}

/* flatness coefficient f1 after Benhama et. al, 1988 */

float calc_f1(float *d)
{
    float f1,x1,x2;

    x1 = 3.0 * calc_ellip(d,3,1);
    x2 = 1.0 + calc_ellip(d,2,1) + calc_ellip(d,3,1);
    f1 = 1.0 - x1 / x2;
    return f1;
}

/* linearity coefficient l1 */

float calc_l1(float *d)
{
    float l1,x1,x2;

    x1 = 3. * ( calc_ellip(d,2,1) + calc_ellip(d,3,1) );
    x2 = 2. * ( 1 + calc_ellip(d,2,1) + calc_ellip(d,3,1));
    l1 = 1. - x1 / x2;
    return l1;
}

/* direction cosines or directivity functions (3 components) */

void calc_dir(float **data3c_dir, float **v, int it)
{
    int i;
    for (i=1;i<=3;i++) {
        data3c_dir[i][it]=v[i][1];
    }
}

/* amplitude parameters */

/* eigenresultant */

float calc_er(float *d)
{
    float er;
    er = sqrt(fabs(d[1]));
    return er;
}

/* instantaneous and quadratic resultant (Meyer, 1988) */

void ampparams(float **indata, float *data_ir, float *data_qr, int nt, int iwl)
{
    int i,it;
    float sqrsum;

    for (it=0;it<nt;it++) {
        /* instantaneous resultant */
        data_ir[it] = sqrt(indata[1][it]*indata[1][it]+ \
            indata[2][it]*indata[2][it]+indata[3][it]*indata[3][it]);

        /* quadratic resultant */
        if ( (it >= iwl/2) && (it < nt-iwl/2) ) {
            sqrsum=0.0;
            for (i=it-iwl/2;i<(it+iwl);i++) {
                sqrsum += indata[1][i]*indata[1][i]+ \
                    indata[2][i]*indata[2][i]+indata[3][i]*indata[3][i];
            }
            data_qr[it] = sqrsum / (float) iwl;
        }
    }
}

/* P, S detection parameters */

/* P-Wave energy filter */
float calc_pfilt(float rl, float theta)
{
    float pfilt;
    pfilt = rl*cos(theta);
    return pfilt;
}

/* S-Wave energy filter */
float calc_sfilt(float rl, float theta)
{
    float sfilt;
    sfilt = rl * (1 - cos(theta));
    return sfilt;
}

/**********************************************************************/
/* Functions for data output                                          */
/**********************************************************************/

/* write one-component data into file */

void fputdata(FILE *fileptr, FILE *headerptr, float *outdata, int nt)
{
    efread(&tr, 1, HDRBYTES, headerptr);
    erewind(headerptr);

    memcpy((void *)tr.data, (const void *) outdata, nt*FSIZE);

    fputtr(fileptr, &tr);
}

/* write three-component data into file */

void fputdata3c(FILE *fileptr, FILE *headerptr, float **outdata3c, int nt)
{
    int i;

    for(i=1;i<=3;i++) {
        efread(&tr, 1, HDRBYTES, headerptr);

        memcpy((void *)tr.data, (const void *) outdata3c[i], nt*FSIZE);

        fputtr(fileptr, &tr);
    }
    erewind(headerptr);
}

/* END OF FILE */


text:
#include <gsl/gsl_math.h>
#include <gsl/gsl_cblas.h>
#include "cblas.h"
#include "error_cblas_l3.h"

void
cblas_ssymm (const enum CBLAS_ORDER Order, const enum CBLAS_SIDE Side,
             const enum CBLAS_UPLO Uplo, const int M, const int N,
             const float alpha, const float *A, const int lda, const float *B,
             const int ldb, const float beta, float *C, const int ldc)
{
#define BASE float
#include "source_symm_r.h"
#undef BASE
}

max_stars_repo_licenses: Net-SNMP Xnet

max_issues_repo_licenses: Net-SNMP Xnet

max_forks_repo_licenses: Net-SNMP Xnet


text:
static char help[] = "Compute ln 2 in serial with PETSc, using random\n"
"permutation of sum order.  Shows that floating-point arithmetic is not\n"
"associative.\n\n";

#include <petsc.h>
#include <time.h>

int main(int argc, char **args) {
  PetscErrorCode  ierr;
  PetscMPIInt     size;
  PetscInt        i, j, n=10;
  PetscReal       v, tmp, *a, sum;
  PetscRandom     r;

  PetscInitialize(&argc,&args,NULL,help);

  ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size); CHKERRQ(ierr);
  if (size > 1) {
      SETERRQ(PETSC_COMM_SELF,1,"lntwo only works in serial\n");
  }

  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","options for lntwo",""); CHKERRQ(ierr);
  ierr = PetscOptionsInt("-n","number of terms in sum",
                          "lntwo.c",n,&n,NULL); CHKERRQ(ierr);
  ierr = PetscOptionsEnd(); CHKERRQ(ierr);

  ierr = PetscRandomCreate(PETSC_COMM_WORLD,&r); CHKERRQ(ierr);
  ierr = PetscRandomSetType(r,PETSCRAND48); CHKERRQ(ierr);
  ierr = PetscRandomSetSeed(r,(int)time(NULL)); CHKERRQ(ierr);
  ierr = PetscRandomSeed(r); CHKERRQ(ierr);

  // fill array with the terms   (-1)^i / (i+1)  for i=0 .. n
  ierr = PetscMalloc1(n,&a); CHKERRQ(ierr);
  for (i=0; i<n; i++) {
      a[i] = pow(-1.0,(double)i) / (double)(i+1);
  }

  // shuffle the terms
  for (i=n-1; i>0; i--) {
      ierr = PetscRandomGetValueReal(r,&v); CHKERRQ(ierr);
      j = (int)floor(i*v);
      tmp = a[i];
      a[i] = a[j];
      a[j] = tmp;
  }

  // sum the terms
  sum = 0.0;
  for (i=0; i<n; i++) {
      sum += a[i];
  }

  // print the result
  ierr = PetscPrintf(PETSC_COMM_WORLD,"ln 2 is approximately %18.16f\n",sum); CHKERRQ(ierr);

  PetscRandomDestroy(&r);
  PetscFree(a);
  return PetscFinalize();
}


text:
/* interpolation/interp2d.c
 *
 * Copyright 2012 David Zaslavsky
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_interp.h>
#include <gsl/gsl_interp2d.h>

/**
 * Triggers a GSL error if the argument is not equal to GSL_SUCCESS.
 * If the argument is GSL_SUCCESS, this does nothing.
 */
#define DISCARD_STATUS(s) if ((s) != GSL_SUCCESS) { GSL_ERROR_VAL("interpolation error", (s),  GSL_NAN); }

#define IDX2D(i, j, w) ((j) * ((w)->xsize) + (i))

gsl_interp2d *
gsl_interp2d_alloc(const gsl_interp2d_type * T, const size_t xsize,
                   const size_t ysize)
{
  gsl_interp2d * interp;

  if (xsize < T->min_size || ysize < T->min_size)
    {
      GSL_ERROR_NULL ("insufficient number of points for interpolation type",
                      GSL_EINVAL);
    }

  interp = (gsl_interp2d *) calloc(1, sizeof(gsl_interp2d));
  if (interp == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate space for gsl_interp2d struct",
                      GSL_ENOMEM);
    }

  interp->type = T;
  interp->xsize = xsize;
  interp->ysize = ysize;

  if (interp->type->alloc == NULL)
    {
      interp->state = NULL;
      return interp;
    }

  interp->state = interp->type->alloc(xsize, ysize);
  if (interp->state == NULL)
    {
      free(interp);
      GSL_ERROR_NULL ("failed to allocate space for gsl_interp2d state",
                      GSL_ENOMEM);
    }

  return interp;
} /* gsl_interp2d_alloc() */

void
gsl_interp2d_free (gsl_interp2d * interp)
{
  RETURN_IF_NULL(interp);

  if (interp->type->free)
    interp->type->free(interp->state);

  free(interp);
} /* gsl_interp2d_free() */

int
gsl_interp2d_init (gsl_interp2d * interp, const double xarr[], const double yarr[],
                   const double zarr[], const size_t xsize, const size_t ysize)
{
  size_t i;

  if (xsize != interp->xsize || ysize != interp->ysize)
    {
      GSL_ERROR("data must match size of interpolation object", GSL_EINVAL);
    }

  for (i = 1; i < xsize; i++)
    {
      if (xarr[i-1] >= xarr[i])
        {
          GSL_ERROR("x values must be strictly increasing", GSL_EINVAL);
        }
    }

  for (i = 1; i < ysize; i++)
    {
      if (yarr[i-1] >= yarr[i])
        {
          GSL_ERROR("y values must be strictly increasing", GSL_EINVAL);
        }
    }

  interp->xmin = xarr[0];
  interp->xmax = xarr[xsize - 1];
  interp->ymin = yarr[0];
  interp->ymax = yarr[ysize - 1];

  {
    int status = interp->type->init(interp->state, xarr, yarr, zarr,
                                    xsize, ysize);
    return status;
  }
} /* gsl_interp2d_init() */

/*
 * A wrapper function that checks boundary conditions, calls an evaluator
 * which implements the actual calculation of the function value or
 * derivative etc., and checks the return status.
 */
static int
interp2d_eval(int (*evaluator)(const void *, const double xa[], const double ya[],
                               const double za[], size_t xsize, size_t ysize,
                               double x, double y, gsl_interp_accel *,
                               gsl_interp_accel *, double * z),
              const gsl_interp2d * interp, const double xarr[],
              const double yarr[], const double zarr[],
              const double x, const double y,
              gsl_interp_accel * xa, gsl_interp_accel * ya,
              double * result)
{
  if (x < interp->xmin || x > interp->xmax)
    {
      GSL_ERROR ("interpolation x value out of range", GSL_EDOM);
    }
  else if (y < interp->ymin || y > interp->ymax)
    {
      GSL_ERROR ("interpolation y value out of range", GSL_EDOM);
    }

  return evaluator(interp->state, xarr, yarr, zarr,
                   interp->xsize, interp->ysize,
                   x, y, xa, ya, result);
}

/*
 * Another wrapper function that serves as a drop-in replacement for
 * interp2d_eval but does not check the bounds. This can be used
 * for extrapolation.
 */
static int
interp2d_eval_extrap(int (*evaluator)(const void *, const double xa[],
                                      const double ya[], const double za[],
                                      size_t xsize, size_t ysize,
                                      double x, double y,
                                      gsl_interp_accel *,
                                      gsl_interp_accel *, double * z),
                     const gsl_interp2d * interp, const double xarr[],
                     const double yarr[], const double zarr[],
                     const double x, const double y,
                     gsl_interp_accel * xa, gsl_interp_accel * ya,
                     double * result)
{
  return evaluator(interp->state, xarr, yarr, zarr,
                   interp->xsize, interp->ysize, x, y, xa, ya, result);
}

double
gsl_interp2d_eval (const gsl_interp2d * interp, const double xarr[],
                   const double yarr[], const double zarr[],
                   const double x, const double y,
                   gsl_interp_accel * xa, gsl_interp_accel * ya)
{
  double z;
  int status = gsl_interp2d_eval_e(interp, xarr, yarr, zarr, x, y, xa, ya, &z);
  DISCARD_STATUS(status)
  return z;
} /* gsl_interp2d_eval() */

double
gsl_interp2d_eval_extrap (const gsl_interp2d * interp,
                          const double xarr[],
                          const double yarr[],
                          const double zarr[],
                          const double x,
                          const double y,
                          gsl_interp_accel * xa,
                          gsl_interp_accel * ya)
{
  double z;
  int status =
    interp2d_eval_extrap(interp->type->eval, interp,
                         xarr, yarr, zarr, x, y, xa, ya, &z);
  DISCARD_STATUS(status)
  return z;
}

int
gsl_interp2d_eval_e (const gsl_interp2d * interp, const double xarr[],
                     const double yarr[], const double zarr[],
                     const double x, const double y,
                     gsl_interp_accel * xa, gsl_interp_accel * ya, double * z)
{
  return interp2d_eval(interp->type->eval, interp,
                       xarr, yarr, zarr, x, y, xa, ya, z);
} /* gsl_interp2d_eval_e() */

int
gsl_interp2d_eval_e_extrap (const gsl_interp2d * interp,
                            const double xarr[], const double yarr[],
                            const double zarr[], const double x,
                            const double y, gsl_interp_accel * xa,
                            gsl_interp_accel * ya, double * z)
{
  return interp2d_eval_extrap(interp->type->eval, interp,
                              xarr, yarr, zarr, x, y, xa, ya, z);
}

double
gsl_interp2d_eval_deriv_x (const gsl_interp2d * interp, const double xarr[],
                           const double yarr[], const double zarr[],
                           const double x, const double y,
                           gsl_interp_accel * xa, gsl_interp_accel * ya)
{
  double z;
  int status = gsl_interp2d_eval_deriv_x_e(interp, xarr, yarr, zarr, x, y, xa, ya, &z);
  DISCARD_STATUS(status)
  return z;
}

int
gsl_interp2d_eval_deriv_x_e (const gsl_interp2d * interp, const double xarr[],
                             const double yarr[], const double zarr[],
                             const double x, const double y,
                             gsl_interp_accel * xa, gsl_interp_accel * ya, double * z)
{
  return interp2d_eval(interp->type->eval_deriv_x, interp,
                       xarr, yarr, zarr, x, y, xa, ya, z);
}

double
gsl_interp2d_eval_deriv_y (const gsl_interp2d * interp, const double xarr[],
                           const double yarr[], const double zarr[],
                           const double x, const double y,
                           gsl_interp_accel * xa, gsl_interp_accel * ya)
{
  double z;
  int status = gsl_interp2d_eval_deriv_y_e(interp, xarr, yarr, zarr, x, y, xa, ya, &z);
  DISCARD_STATUS(status)
  return z;
}

int
gsl_interp2d_eval_deriv_y_e (const gsl_interp2d * interp, const double xarr[],
                             const double yarr[], const double zarr[],
                             const double x, const double y,
                             gsl_interp_accel * xa, gsl_interp_accel * ya, double * z)
{
  return interp2d_eval(interp->type->eval_deriv_y, interp,
                       xarr, yarr, zarr, x, y, xa, ya, z);
}

double
gsl_interp2d_eval_deriv_xx (const gsl_interp2d * interp, const double xarr[],
                            const double yarr[], const double zarr[],
                            const double x, const double y,
                            gsl_interp_accel * xa, gsl_interp_accel * ya)
{
  double z;
  int status = gsl_interp2d_eval_deriv_xx_e(interp, xarr, yarr, zarr, x, y, xa, ya, &z);
  DISCARD_STATUS(status)
  return z;
}

int
gsl_interp2d_eval_deriv_xx_e (const gsl_interp2d * interp, const double xarr[],
                              const double yarr[], const double zarr[],
                              const double x, const double y,
                              gsl_interp_accel * xa, gsl_interp_accel * ya, double * z)
{
  return interp2d_eval(interp->type->eval_deriv_xx, interp,
                       xarr, yarr, zarr, x, y, xa, ya, z);
}

double
gsl_interp2d_eval_deriv_yy (const gsl_interp2d * interp, const double xarr[],
                            const double yarr[], const double zarr[],
                            const double x, const double y,
                            gsl_interp_accel * xa, gsl_interp_accel * ya)
{
  double z;
  int status = gsl_interp2d_eval_deriv_yy_e(interp, xarr, yarr, zarr, x, y, xa, ya, &z);
  DISCARD_STATUS(status)
  return z;
}

int
gsl_interp2d_eval_deriv_yy_e (const gsl_interp2d * interp, const double xarr[],
                              const double yarr[], const double zarr[],
                              const double x, const double y,
                              gsl_interp_accel * xa, gsl_interp_accel * ya, double * z)
{
  return interp2d_eval(interp->type->eval_deriv_yy, interp,
                       xarr, yarr, zarr, x, y, xa, ya, z);
}

double
gsl_interp2d_eval_deriv_xy (const gsl_interp2d * interp, const double xarr[],
                            const double yarr[], const double zarr[],
                            const double x, const double y,
                            gsl_interp_accel * xa, gsl_interp_accel * ya)
{
  double z;
  int status = gsl_interp2d_eval_deriv_xy_e(interp, xarr, yarr, zarr, x, y, xa, ya, &z);
  DISCARD_STATUS(status)
  return z;
}

int
gsl_interp2d_eval_deriv_xy_e (const gsl_interp2d * interp, const double xarr[],
                              const double yarr[], const double zarr[],
                              const double x, const double y,
                              gsl_interp_accel * xa, gsl_interp_accel * ya, double * z)
{
  return interp2d_eval(interp->type->eval_deriv_xy, interp,
                       xarr, yarr, zarr, x, y, xa, ya, z);
}

size_t
gsl_interp2d_type_min_size(const gsl_interp2d_type * T)
{
  return T->min_size;
}

size_t
gsl_interp2d_min_size(const gsl_interp2d * interp)
{
  return interp->type->min_size;
}

const char *
gsl_interp2d_name(const gsl_interp2d * interp)
{
  return interp->type->name;
}

size_t
gsl_interp2d_idx(const gsl_interp2d * interp,
                 const size_t i, const size_t j)
{
  if (i >= interp->xsize)
    {
      GSL_ERROR_VAL ("x index out of range", GSL_ERANGE, 0);
    }
  else if (j >= interp->ysize)
    {
      GSL_ERROR_VAL ("y index out of range", GSL_ERANGE, 0);
    }
  else
    {
      return IDX2D(i, j, interp);
    }
} /* gsl_interp2d_idx() */

int
gsl_interp2d_set(const gsl_interp2d * interp, double zarr[],
                 const size_t i, const size_t j, const double z)
{
  if (i >= interp->xsize)
    {
      GSL_ERROR ("x index out of range", GSL_ERANGE);
    }
  else if (j >= interp->ysize)
    {
      GSL_ERROR ("y index out of range", GSL_ERANGE);
    }
  else
    {
      zarr[IDX2D(i, j, interp)] = z;
      return GSL_SUCCESS;
    }
} /* gsl_interp2d_set() */

double
gsl_interp2d_get(const gsl_interp2d * interp, const double zarr[],
                 const size_t i, const size_t j)
{
  if (i >= interp->xsize)
    {
      GSL_ERROR_VAL ("x index out of range", GSL_ERANGE, 0);
    }
  else if (j >= interp->ysize)
    {
      GSL_ERROR_VAL ("y index out of range", GSL_ERANGE, 0);
    }
  else
    {
      return zarr[IDX2D(i, j, interp)];
    }
} /* gsl_interp2d_get() */

#undef IDX2D


text:
/**
 *
 * @generated c Tue Jan  7 11:45:26 2014
 *
 **/
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <cblas.h>
#include <lapacke.h>
#include <plasma.h>
#include <core_blas.h>
#include "auxiliary.h"

/*-------------------------------------------------------------------
 * Check the orthogonality of Q
 */

int c_check_orthogonality(int M, int N, int LDQ, PLASMA_Complex32_t *Q)
{
    float alpha, beta;
    float normQ;
    int info_ortho;
    int i;
    int minMN = min(M, N);
    float eps;
    float *work = (float *)malloc(minMN*sizeof(float));

    eps = LAPACKE_slamch_work('e');
    alpha = 1.0;
    beta  = -1.0;

    /* Build the idendity matrix USE DLASET?*/
    PLASMA_Complex32_t *Id = (PLASMA_Complex32_t *) malloc(minMN*minMN*sizeof(PLASMA_Complex32_t));
    memset((void*)Id, 0, minMN*minMN*sizeof(PLASMA_Complex32_t));
    for (i = 0; i < minMN; i++)
        Id[i*minMN+i] = (PLASMA_Complex32_t)1.0;

    /* Perform Id - Q'Q */
    if (M >= N)
        cblas_cherk(CblasColMajor, CblasUpper, CblasConjTrans, N, M, alpha, Q, LDQ, beta, Id, N);
    else
        cblas_cherk(CblasColMajor, CblasUpper, CblasNoTrans, M, N, alpha, Q, LDQ, beta, Id, M);

    normQ = LAPACKE_clansy_work(LAPACK_COL_MAJOR, 'i', 'u', minMN, Id, minMN, work);

    printf("============\n");
    printf("Checking the orthogonality of Q \n");
    printf("||Id-Q'*Q||_oo / (N*eps) = %e \n",normQ/(minMN*eps));

    if ( isnan(normQ / (minMN * eps)) || (normQ / (minMN * eps) > 10.0) ) {
        printf("-- Orthogonality is suspicious ! \n");
        info_ortho=1;
    }
    else {
        printf("-- Orthogonality is CORRECT ! \n");
        info_ortho=0;
    }

    free(work); free(Id);

    return info_ortho;
}

/*------------------------------------------------------------
 *  Check the factorization QR
 */

int c_check_QRfactorization(int M, int N, PLASMA_Complex32_t *A1, PLASMA_Complex32_t *A2, int LDA, PLASMA_Complex32_t *Q)
{
    float Anorm, Rnorm;
    PLASMA_Complex32_t alpha, beta;
    int info_factorization;
    int i,j;
    float eps;

    eps = LAPACKE_slamch_work('e');

    PLASMA_Complex32_t *Ql       = (PLASMA_Complex32_t *)malloc(M*N*sizeof(PLASMA_Complex32_t));
    PLASMA_Complex32_t *Residual = (PLASMA_Complex32_t *)malloc(M*N*sizeof(PLASMA_Complex32_t));
    float *work              = (float *)malloc(max(M,N)*sizeof(float));

    alpha=1.0;
    beta=0.0;

    if (M >= N) {
        /* Extract the R */
        PLASMA_Complex32_t *R = (PLASMA_Complex32_t *)malloc(N*N*sizeof(PLASMA_Complex32_t));
        memset((void*)R, 0, N*N*sizeof(PLASMA_Complex32_t));
        LAPACKE_clacpy_work(LAPACK_COL_MAJOR,'u', M, N, A2, LDA, R, N);

        /* Perform Ql=Q*R */
        memset((void*)Ql, 0, M*N*sizeof(PLASMA_Complex32_t));
        cblas_cgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, M, N, N, CBLAS_SADDR(alpha), Q, LDA, R, N, CBLAS_SADDR(beta), Ql, M);
        free(R);
    }
    else {
        /* Extract the L */
        PLASMA_Complex32_t *L = (PLASMA_Complex32_t *)malloc(M*M*sizeof(PLASMA_Complex32_t));
        memset((void*)L, 0, M*M*sizeof(PLASMA_Complex32_t));
        LAPACKE_clacpy_work(LAPACK_COL_MAJOR,'l', M, N, A2, LDA, L, M);

    /* Perform Ql=LQ */
        memset((void*)Ql, 0, M*N*sizeof(PLASMA_Complex32_t));
        cblas_cgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, M, N, M, CBLAS_SADDR(alpha), L, M, Q, LDA, CBLAS_SADDR(beta), Ql, M);
        free(L);
    }

    /* Compute the Residual */
    for (i = 0; i < M; i++)
        for (j = 0 ; j < N; j++)
            Residual[j*M+i] = A1[j*LDA+i]-Ql[j*M+i];

    Rnorm = LAPACKE_clange_work(LAPACK_COL_MAJOR, 'i', M, N, Residual, M, work);
    Anorm = LAPACKE_clange_work(LAPACK_COL_MAJOR, 'i', M, N, A2, LDA, work);

    if (M >= N) {
        printf("============\n");
        printf("Checking the QR Factorization \n");
        printf("-- ||A-QR||_oo/(||A||_oo.N.eps) = %e \n",Rnorm/(Anorm*N*eps));
    }
    else {
        printf("============\n");
        printf("Checking the LQ Factorization \n");
        printf("-- ||A-LQ||_oo/(||A||_oo.N.eps) = %e \n",Rnorm/(Anorm*N*eps));
    }

    if (isnan(Rnorm / (Anorm * N *eps)) || (Rnorm / (Anorm * N * eps) > 10.0) ) {
        printf("-- Factorization is suspicious ! \n");
        info_factorization = 1;
    }
    else {
        printf("-- Factorization is CORRECT ! \n");
        info_factorization = 0;
    }

    free(work); free(Ql); free(Residual);

    return info_factorization;
}

/*------------------------------------------------------------------------
 *  Check the factorization of the matrix A2
 */

int c_check_LLTfactorization(int N, PLASMA_Complex32_t *A1, PLASMA_Complex32_t *A2, int LDA, int uplo)
{
    float Anorm, Rnorm;
    PLASMA_Complex32_t alpha;
    int info_factorization;
    int i,j;
    float eps;

    eps = LAPACKE_slamch_work('e');

    PLASMA_Complex32_t *Residual = (PLASMA_Complex32_t *)malloc(N*N*sizeof(PLASMA_Complex32_t));
    PLASMA_Complex32_t *L1       = (PLASMA_Complex32_t *)malloc(N*N*sizeof(PLASMA_Complex32_t));
    PLASMA_Complex32_t *L2       = (PLASMA_Complex32_t *)malloc(N*N*sizeof(PLASMA_Complex32_t));
    float *work              = (float *)malloc(N*sizeof(float));

    memset((void*)L1, 0, N*N*sizeof(PLASMA_Complex32_t));
    memset((void*)L2, 0, N*N*sizeof(PLASMA_Complex32_t));

    alpha= 1.0;

    LAPACKE_clacpy_work(LAPACK_COL_MAJOR,' ', N, N, A1, LDA, Residual, N);

    /* Dealing with L'L or U'U  */
    if (uplo == PlasmaUpper){
        LAPACKE_clacpy_work(LAPACK_COL_MAJOR,'u', N, N, A2, LDA, L1, N);
        LAPACKE_clacpy_work(LAPACK_COL_MAJOR,'u', N, N, A2, LDA, L2, N);
        cblas_ctrmm(CblasColMajor, CblasLeft, CblasUpper, CblasConjTrans, CblasNonUnit, N, N, CBLAS_SADDR(alpha), L1, N, L2, N);
    }
    else{
        LAPACKE_clacpy_work(LAPACK_COL_MAJOR,'l', N, N, A2, LDA, L1, N);
        LAPACKE_clacpy_work(LAPACK_COL_MAJOR,'l', N, N, A2, LDA, L2, N);
        cblas_ctrmm(CblasColMajor, CblasRight, CblasLower, CblasConjTrans, CblasNonUnit, N, N, CBLAS_SADDR(alpha), L1, N, L2, N);
    }

    /* Compute the Residual || A -L'L|| */
    for (i = 0; i < N; i++)
        for (j = 0; j < N; j++)
           Residual[j*N+i] = L2[j*N+i] - Residual[j*N+i];

    Rnorm = LAPACKE_clange_work(LAPACK_COL_MAJOR, 'i', N, N, Residual, N, work);
    Anorm = LAPACKE_clange_work(LAPACK_COL_MAJOR, 'i', N, N, A1, LDA, work);

    printf("============\n");
    printf("Checking the Cholesky Factorization \n");
    printf("-- ||L'L-A||_oo/(||A||_oo.N.eps) = %e \n",Rnorm/(Anorm*N*eps));

    if ( isnan(Rnorm/(Anorm*N*eps)) || (Rnorm/(Anorm*N*eps) > 10.0) ){
        printf("-- Factorization is suspicious ! \n");
        info_factorization = 1;
    }
    else{
        printf("-- Factorization is CORRECT ! \n");
        info_factorization = 0;
    }

    free(Residual); free(L1); free(L2); free(work);

    return info_factorization;
}

/*--------------------------------------------------------------
 * Check the gemm
 */
float c_check_gemm(PLASMA_enum transA, PLASMA_enum transB, int M, int N, int K,
                   PLASMA_Complex32_t alpha, PLASMA_Complex32_t *A, int LDA,
                   PLASMA_Complex32_t *B, int LDB,
                   PLASMA_Complex32_t beta, PLASMA_Complex32_t *Cplasma,
                   PLASMA_Complex32_t *Cref, int LDC,
                   float *Cinitnorm, float *Cplasmanorm, float *Clapacknorm )
{
    PLASMA_Complex32_t beta_const = -1.0;
    float Rnorm;
    float *work = (float *)malloc(max(K,max(M, N))* sizeof(float));

    *Cinitnorm   = LAPACKE_clange_work(LAPACK_COL_MAJOR, 'i', M, N, Cref,    LDC, work);
    *Cplasmanorm = LAPACKE_clange_work(LAPACK_COL_MAJOR, 'i', M, N, Cplasma, LDC, work);

    cblas_cgemm(CblasColMajor, (CBLAS_TRANSPOSE)transA, (CBLAS_TRANSPOSE)transB, M, N, K,
                CBLAS_SADDR(alpha), A, LDA, B, LDB, CBLAS_SADDR(beta), Cref, LDC);

    *Clapacknorm = LAPACKE_clange_work(LAPACK_COL_MAJOR, 'i', M, N, Cref, LDC, work);

    cblas_caxpy(LDC * N, CBLAS_SADDR(beta_const), Cplasma, 1, Cref, 1);

    Rnorm = LAPACKE_clange_work(LAPACK_COL_MAJOR, 'i', M, N, Cref, LDC, work);

    free(work);

    return Rnorm;
}

/*--------------------------------------------------------------
 * Check the trsm
 */
float c_check_trsm(PLASMA_enum side, PLASMA_enum uplo, PLASMA_enum trans, PLASMA_enum diag,
                   int M, int NRHS, PLASMA_Complex32_t alpha,
                   PLASMA_Complex32_t *A, int LDA,
                   PLASMA_Complex32_t *Bplasma, PLASMA_Complex32_t *Bref, int LDB,
                   float *Binitnorm, float *Bplasmanorm, float *Blapacknorm )
{
    PLASMA_Complex32_t beta_const = -1.0;
    float Rnorm;
    float *work = (float *)malloc(max(M, NRHS)* sizeof(float));
    /*float eps = LAPACKE_slamch_work('e');*/

    *Binitnorm   = LAPACKE_clange_work(LAPACK_COL_MAJOR, 'i', M, NRHS, Bref,    LDB, work);
    *Bplasmanorm = LAPACKE_clange_work(LAPACK_COL_MAJOR, 'm', M, NRHS, Bplasma, LDB, work);

    cblas_ctrsm(CblasColMajor, (CBLAS_SIDE)side, (CBLAS_UPLO)uplo,
                (CBLAS_TRANSPOSE)trans, (CBLAS_DIAG)diag, M, NRHS,
                CBLAS_SADDR(alpha), A, LDA, Bref, LDB);

    *Blapacknorm = LAPACKE_clange_work(LAPACK_COL_MAJOR, 'm', M, NRHS, Bref, LDB, work);

    cblas_caxpy(LDB * NRHS, CBLAS_SADDR(beta_const), Bplasma, 1, Bref, 1);

    Rnorm = LAPACKE_clange_work(LAPACK_COL_MAJOR, 'm', M, NRHS, Bref, LDB, work);
    Rnorm = Rnorm / *Blapacknorm;
    /* max(M,NRHS) * eps);*/

    free(work);

    return Rnorm;
}

/*--------------------------------------------------------------
 * Check the solution
 */

float c_check_solution(int M, int N, int NRHS, PLASMA_Complex32_t *A, int LDA,
                      PLASMA_Complex32_t *B,  PLASMA_Complex32_t *X, int LDB,
                      float *anorm, float *bnorm, float *xnorm )
{
/*     int info_solution; */
    float Rnorm = -1.00;
    PLASMA_Complex32_t zone  =  1.0;
    PLASMA_Complex32_t mzone = -1.0;
    float *work = (float *)malloc(max(M, N)* sizeof(float));

    *anorm = LAPACKE_clange_work(LAPACK_COL_MAJOR, 'i', M, N,    A, LDA, work);
    *xnorm = LAPACKE_clange_work(LAPACK_COL_MAJOR, 'i', M, NRHS, X, LDB, work);
    *bnorm = LAPACKE_clange_work(LAPACK_COL_MAJOR, 'i', N, NRHS, B, LDB, work);

    cblas_cgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, M, NRHS, N, CBLAS_SADDR(zone), A, LDA, X, LDB, CBLAS_SADDR(mzone), B, LDB);

    Rnorm = LAPACKE_clange_work(LAPACK_COL_MAJOR, 'i', N, NRHS, B, LDB, work);

    free(work);

    return Rnorm;
}


text:
/* specfunc/laguerre.c
 *
 * Copyright (C) 2007 Brian Gough
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

/* Author:  G. Jungman */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_exp.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_sf_laguerre.h>

#include "error.h"

/*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/

/* based on the large 2b-4a asymptotic for 1F1
 * [Abramowitz+Stegun, 13.5.21]
 * L^a_n(x) = (a+1)_n / n! 1F1(-n,a+1,x)
 *
 * The second term (ser_term2) is from Slater,"The Confluent
 * Hypergeometric Function" p.73.  I think there may be an error in
 * the first term of the expression given there, comparing with AS
 * 13.5.21 (cf sin(a\pi+\Theta) vs sin(a\pi) + sin(\Theta)) - but the
 * second term appears correct.
 *
 */
static
int
laguerre_large_n(const int n, const double alpha, const double x,
                 gsl_sf_result * result)
{
  const double a = -n;
  const double b = alpha + 1.0;
  const double eta    = 2.0*b - 4.0*a;
  const double cos2th = x/eta;
  const double sin2th = 1.0 - cos2th;
  const double eps = asin(sqrt(cos2th));  /* theta = pi/2 - eps */
  const double pre_h  = 0.25*M_PI*M_PI*eta*eta*cos2th*sin2th;
  gsl_sf_result lg_b;
  gsl_sf_result lnfact;
  int stat_lg = gsl_sf_lngamma_e(b+n, &lg_b);
  int stat_lf = gsl_sf_lnfact_e(n, &lnfact);
  double pre_term1 = 0.5*(1.0-b)*log(0.25*x*eta);
  double pre_term2 = 0.25*log(pre_h);
  double lnpre_val = lg_b.val - lnfact.val + 0.5*x + pre_term1 - pre_term2;
  double lnpre_err = lg_b.err + lnfact.err + GSL_DBL_EPSILON * (fabs(pre_term1)+fabs(pre_term2));

  double phi1 = 0.25*eta*(2*eps + sin(2.0*eps));
  double ser_term1 = -sin(phi1);

  double A1 = (1.0/12.0)*(5.0/(4.0*sin2th)+(3.0*b*b-6.0*b+2.0)*sin2th - 1.0);
  double ser_term2 = -A1 * cos(phi1)/(0.25*eta*sin(2.0*eps));

  double ser_val = ser_term1 + ser_term2;
  double ser_err = ser_term2*ser_term2 + GSL_DBL_EPSILON * (fabs(ser_term1) + fabs(ser_term2));
  int stat_e = gsl_sf_exp_mult_err_e(lnpre_val, lnpre_err, ser_val, ser_err, result);
  result->err += 2.0 * GSL_SQRT_DBL_EPSILON * fabs(result->val);
  return GSL_ERROR_SELECT_3(stat_e, stat_lf, stat_lg);
}

/* Evaluate polynomial based on confluent hypergeometric representation.
 *
 * L^a_n(x) = (a+1)_n / n! 1F1(-n,a+1,x)
 *
 * assumes n > 0 and a != negative integer greater than -n
 */
static
int
laguerre_n_cp(const int n, const double a, const double x, gsl_sf_result * result)
{
  gsl_sf_result lnfact;
  gsl_sf_result lg1;
  gsl_sf_result lg2;
  double s1, s2;
  int stat_f = gsl_sf_lnfact_e(n, &lnfact);
  int stat_g1 = gsl_sf_lngamma_sgn_e(a+1.0+n, &lg1, &s1);
  int stat_g2 = gsl_sf_lngamma_sgn_e(a+1.0, &lg2, &s2);
  double poly_1F1_val = 1.0;
  double poly_1F1_err = 0.0;
  int stat_e;
  int k;

  double lnpre_val = (lg1.val - lg2.val) - lnfact.val;
  double lnpre_err = lg1.err + lg2.err + lnfact.err + 2.0 * GSL_DBL_EPSILON * fabs(lnpre_val);

  for(k=n-1; k>=0; k--) {
    double t = (-n+k)/(a+1.0+k) * (x/(k+1));
    double r = t + 1.0/poly_1F1_val;
    if(r > 0.9*GSL_DBL_MAX/poly_1F1_val) {
      /* internal error only, don't call the error handler */
      INTERNAL_OVERFLOW_ERROR(result);
    }
    else {
      /* Collect the Horner terms. */
      poly_1F1_val  = 1.0 + t * poly_1F1_val;
      poly_1F1_err += GSL_DBL_EPSILON + fabs(t) * poly_1F1_err;
    }
  }

  stat_e = gsl_sf_exp_mult_err_e(lnpre_val, lnpre_err,
                                    poly_1F1_val, poly_1F1_err,
                                    result);

  return GSL_ERROR_SELECT_4(stat_e, stat_f, stat_g1, stat_g2);
}

/* Evaluate the polynomial based on the confluent hypergeometric
 * function in a safe way, with no restriction on the arguments.
 *
 * assumes x != 0
 */
static
int
laguerre_n_poly_safe(const int n, const double a, const double x, gsl_sf_result * result)
{
  const double b = a + 1.0;
  const double mx = -x;
  const double tc_sgn = (x < 0.0 ? 1.0 : (GSL_IS_ODD(n) ? -1.0 : 1.0));
  gsl_sf_result tc;
  int stat_tc = gsl_sf_taylorcoeff_e(n, fabs(x), &tc);

  if(stat_tc == GSL_SUCCESS) {
    double term = tc.val * tc_sgn;
    double sum_val = term;
    double sum_err = tc.err;
    int k;
    for(k=n-1; k>=0; k--) {
      term *= ((b+k)/(n-k))*(k+1.0)/mx;
      sum_val += term;
      sum_err += 4.0 * GSL_DBL_EPSILON * fabs(term);
    }
    result->val = sum_val;
    result->err = sum_err + 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(stat_tc == GSL_EOVRFLW) {
    result->val = 0.0; /* FIXME: should be Inf */
    result->err = 0.0;
    return stat_tc;
  }
  else {
    result->val = 0.0;
    result->err = 0.0;
    return stat_tc;
  }
}

/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*/

int
gsl_sf_laguerre_1_e(const double a, const double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  {
    result->val = 1.0 + a - x;
    result->err = 2.0 * GSL_DBL_EPSILON * (1.0 + fabs(a) + fabs(x));
    return GSL_SUCCESS;
  }
}

int
gsl_sf_laguerre_2_e(const double a, const double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(a == -2.0) {
    result->val = 0.5*x*x;
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    double c0 = 0.5 * (2.0+a)*(1.0+a);
    double c1 = -(2.0+a);
    double c2 = -0.5/(2.0+a);
    result->val  = c0 + c1*x*(1.0 + c2*x);
    result->err  = 2.0 * GSL_DBL_EPSILON * (fabs(c0) + 2.0 * fabs(c1*x) * (1.0 + 2.0 * fabs(c2*x)));
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
}

int
gsl_sf_laguerre_3_e(const double a, const double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(a == -2.0) {
    double x2_6  = x*x/6.0;
    result->val  = x2_6 * (3.0 - x);
    result->err  = x2_6 * (3.0 + fabs(x)) * 2.0 * GSL_DBL_EPSILON;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(a == -3.0) {
    result->val = -x*x/6.0;
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    double c0 = (3.0+a)*(2.0+a)*(1.0+a) / 6.0;
    double c1 = -c0 * 3.0 / (1.0+a);
    double c2 = -1.0/(2.0+a);
    double c3 = -1.0/(3.0*(3.0+a));
    result->val  = c0 + c1*x*(1.0 + c2*x*(1.0 + c3*x));
    result->err  = 1.0 + 2.0 * fabs(c3*x);
    result->err  = 1.0 + 2.0 * fabs(c2*x) * result->err;
    result->err  = 2.0 * GSL_DBL_EPSILON * (fabs(c0) + 2.0 * fabs(c1*x) * result->err);
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
}

int
gsl_sf_laguerre_n_e(const int n, const double a, const double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(n < 0) {
    DOMAIN_ERROR(result);
  }
  else if(n == 0) {
    result->val = 1.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else if(n == 1) {
    result->val = 1.0 + a - x;
    result->err = 2.0 * GSL_DBL_EPSILON * (1.0 + fabs(a) + fabs(x));
    return GSL_SUCCESS;
  }
  else if(x == 0.0) {
    double product = a + 1.0;
    int k;
    for(k=2; k<=n; k++) {
      product *= (a + k)/k;
    }
    result->val = product;
    result->err = 2.0 * (n + 1.0) * GSL_DBL_EPSILON * fabs(product) + GSL_DBL_EPSILON;
    return GSL_SUCCESS;
  }
  else if(x < 0.0 && a > -1.0) {
    /* In this case all the terms in the polynomial
     * are of the same sign. Note that this also
     * catches overflows correctly.
     */
    return laguerre_n_cp(n, a, x, result);
  }
  else if(n < 5 || (x > 0.0 && a < -n-1)) {
    /* Either the polynomial will not lose too much accuracy
     * or all the terms are negative. In any case,
     * the error estimate here is good. We try both
     * explicit summation methods, as they have different
     * characteristics. One may underflow/overflow while the
     * other does not.
     */
    if(laguerre_n_cp(n, a, x, result) == GSL_SUCCESS)
      return GSL_SUCCESS;
    else
      return laguerre_n_poly_safe(n, a, x, result);
  }
  else if(n > 1.0e+07 && x > 0.0 && a > -1.0 && x < 2.0*(a+1.0)+4.0*n) {
    return laguerre_large_n(n, a, x, result);
  }
  else if(a >= 0.0 || (x > 0.0 && a < -n-1)) {
    gsl_sf_result lg2;
    int stat_lg2 = gsl_sf_laguerre_2_e(a, x, &lg2);
    double Lkm1 = 1.0 + a - x;
    double Lk   = lg2.val;
    double Lkp1;
    int k;

    for(k=2; k<n; k++) {
      Lkp1 = (-(k+a)*Lkm1 + (2.0*k+a+1.0-x)*Lk)/(k+1.0);
      Lkm1 = Lk;
      Lk   = Lkp1;
    }
    result->val = Lk;
    result->err = (fabs(lg2.err/lg2.val) + GSL_DBL_EPSILON) * n * fabs(Lk);
    return stat_lg2;
  }
  else {
    /* Despair... or magic? */
    return laguerre_n_poly_safe(n, a, x, result);
  }
}

/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/

#include "eval.h"

double gsl_sf_laguerre_1(double a, double x)
{
  EVAL_RESULT(gsl_sf_laguerre_1_e(a, x, &result));
}

double gsl_sf_laguerre_2(double a, double x)
{
  EVAL_RESULT(gsl_sf_laguerre_2_e(a, x, &result));
}

double gsl_sf_laguerre_3(double a, double x)
{
  EVAL_RESULT(gsl_sf_laguerre_3_e(a, x, &result));
}

double gsl_sf_laguerre_n(int n, double a, double x)
{
  EVAL_RESULT(gsl_sf_laguerre_n_e(n, a, x, &result));
}

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
#include <gsl/gsl_math.h>
#include "gsl_cblas.h"
#include "cblas.h"

void
cblas_zhpmv (const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
	     const int N, const void *alpha, const void *Ap, const void *X,
	     const int incX, const void *beta, void *Y, const int incY)
{
#define BASE double
#include "source_hpmv.h"
#undef BASE
}