text:
/*
 *  vbHmm_template.c
 *  Template for original model to be analyzed by VB-HMM.
 *
 *  Created by OKAMOTO Kenji and SAKO Yasushi
 *  Copyright 2011
 *  Cellular Informatics Laboratory, Advance Science Institute, RIKEN, Japan.
 *  All rights reserved.
 *
 *  Ver. 1.0.0
 *  Last modified on 2011.04.19
 */

#include <math.h>
#include <string.h>

// Special functions is probably necessary.
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_sf_psi.h>

#include "vbHmm_template.h"

// For random number generation.
#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__(){
    commonFunctions funcs;
    funcs.newModelParameters    = newModelParameters__;
    funcs.freeModelParameters   = freeModelParameters__;
    funcs.newModelStats         = newModelStats__;
    funcs.freeModelStats        = freeModelStats__;
    funcs.initializeVbHmm       = initializeVbHmm__;
    funcs.pTilde_z1             = pTilde_z1__;
    funcs.pTilde_zn_zn1         = pTilde_zn_zn1__;
    funcs.pTilde_xn_zn          = pTilde_xn_zn__;
    funcs.calcStatsVars         = calcStatsVars__;
    funcs.maximization          = maximization__;
    funcs.varLowerBound         = varLowerBound__;
    funcs.reorderParameters     = reorderParameters__;
    funcs.outputResults         = outputResults__;
    setFunctions( funcs );
}

void setGFunctions__(){
    gCommonFunctions funcs;
    funcs.newModelParameters    = newModelParameters__;
    funcs.freeModelParameters   = freeModelParameters__;
    funcs.newModelStats         = newModelStats__;
    funcs.freeModelStats        = freeModelStats__;
    funcs.newModelStatsG        = newModelStatsG__;
    funcs.freeModelStatsG       = freeModelStatsG__;
    funcs.initializeVbHmmG      = initializeVbHmmG__;
    funcs.pTilde_z1             = pTilde_z1__;
    funcs.pTilde_zn_zn1         = pTilde_zn_zn1__;
    funcs.pTilde_xn_zn          = pTilde_xn_zn__;
    funcs.calcStatsVarsG        = calcStatsVarsG__;
    funcs.maximizationG         = maximizationG__;
    funcs.varLowerBoundG        = varLowerBoundG__;
    funcs.reorderParametersG    = reorderParametersG__;
    funcs.outputResultsG        = outputResultsG__;
    setGFunctions( funcs );
    isGlobalAnalysis = 1;
}

// Redirect output function to that for model-specific data.
void outputResults__( xn, gv, iv, logFP )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
FILE *logFP;
{
    output__Results( xn, gv, iv, logFP );
}

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

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

    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__( p, xn, sNo )
void **p;
xnDataSet *xn;
int sNo;
{
    tempParameters *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( *p );
    *p = NULL;
}

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

        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__( s, xn, gv, iv )
void **s;
xnDataSet *xn;
globalVars *gv;
indVars *iv;
{
    if( isGlobalAnalysis == 0 ){
        int sNo = gv->sNo;
        tempStats *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__( xns, gv, ivs)
xnDataBundle *xns;
globalVars *gv;
indVarBundle *ivs;
{
    int sNo = gv->sNo;
    tempGlobalStats *gs = (tempGlobalStats*)malloc( sizeof(tempGlobalStats) );

    int i;
    gs->NiR = (double *)malloc( sNo * sizeof(double) );
    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->z1iR = (double *)malloc( sNo * sizeof(double) );

    return gs;
}

void freeModelStatsG__( gs, xns, gv, ivs )
void **gs;
xnDataBundle *xns;
globalVars *gv;
indVarBundle *ivs;
{
    int sNo = gv->sNo;
    tempGlobalStats *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__( xn, gv, iv )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
{
    // an example of hyperparameter initialization

    int sNo = gv->sNo;
    tempParameters *p = gv->params;

    int i, j;
    // 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];
        }
    }

    initialize_indVars__( xn, gv, iv );

    calcStatsVars__( xn, gv, iv );
    maximization__( xn, gv, iv );
}

void initializeVbHmmG__( xns, gv, ivs )
xnDataBundle *xns;
globalVars *gv;
indVarBundle *ivs;
{
    // an example of hyperparameter initialization

    int sNo = gv->sNo, rNo = xns->R;
    tempParameters *p = gv->params;

    int i, j, r;

    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];
        }
    }

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

    calcStatsVarsG__( xns, gv, ivs );
    maximizationG__( xns, gv, ivs );
}

void initialize_indVars__( xn, gv, iv )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
{
    // an example of latent variable distribution

    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__( 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 = (tempData*)malloc( sizeof(tempData) );
    tempData *d = (tempData*)xn->data;
    return xn;
}

void freeXnDataSet__( xn )
xnDataSet **xn;
{
    tempData *d = (tempData*)(*xn)->data;
    free( (*xn)->data );
    free( (*xn)->name );
    free( *xn );
    *xn = NULL;
}

double pTilde_z1__( i, params )
int i;
void *params;
{
    // Return the value for p(z_{1i}|pi_i) for given i.

    // for common HMM
    tempParameters *p = (tempParameters*)params;
    return exp( p->avgLnPi[i] );
}

double pTilde_zn_zn1__( i, j, params )
int i, j;
void *params;
{
    // Return the value for p(z_{nj}|z_{n-1,i},theta) for given i and j.

    // for common HMM
    tempParameters *p = (tempParameters*)params;
    return exp( p->avgLnA[i][j] );
}

double pTilde_xn_zn__( xnWv, n, i, params )
xnDataSet *xnWv;
size_t n;
int i;
void *params;
{
    // Return the value for p(x_{mi}|z_{ni},phi) for given i.

    // return emission probability for the model.
    tempParameters *p = (tempParameters*)params;
    return 0.0;
}

void calcStatsVars__( xn, gv, iv )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
{
    // Execute calculation of statistics, which may be necessary for M-step.

    // calculation of stats for common HMM
    tempStats *s = (tempStats*)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;

    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];
            for( j = 0 ; j < sNo ; j++ ){
                Nii[i]    += xiMat[n][i][j];
                Nij[i][j] += xiMat[n][i][j];
            }
        }
    }
}

void calcStatsVarsG__( xns, gv, ivs )
xnDataBundle *xns;
globalVars *gv;
indVarBundle *ivs;
{
    // Execute calculation of statistics for global analysis, which may be necessary for M-step.

    // calculation of stats for global analysis of common HMM
    tempGlobalStats *gs = (tempGlobalStats*)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;

    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++ ){
        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];
                for( j = 0 ; j < sNo ; j++ ){
                    NiiR[i]  += xiMat[n][i][j];
                    NijR[i][j] += xiMat[n][i][j];
                }
            }
        }
    }
}

void maximization__( xn, gv, iv )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
{
    // Calculation of M-step, depending on the model.

    // calculation of M-step for common HMM
    tempParameters *p = (tempParameters*)gv->params;
    tempStats *s = (tempStats*)iv->stats;
    int sNo = gv->sNo;
    double **gmMat = iv->gmMat;
    double *uPiArr = p->uPiArr, sumUPi = p->sumUPi;
    double **uAMat = p->uAMat, *sumUAArr = p->sumUAArr;
    double *avgPi = p->avgPi, *avgLnPi = p->avgLnPi, **avgA = p->avgA, **avgLnA = p->avgLnA;
    double *Nii = s->Nii, **Nij = s->Nij;
    int i, j;

    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] );
        }
    }
}

void maximizationG__( xns, gv, ivs )
xnDataBundle *xns;
globalVars *gv;
indVarBundle *ivs;
{
    // Calculation of M-step for global analysis, depending on the model.

    // calculation of stats for global analysis of common HMM
    tempParameters *p = (tempParameters*)gv->params;
    double *uPiArr = p->uPiArr, sumUPi = p->sumUPi;
    double **uAMat = p->uAMat, *sumUAArr = p->sumUAArr;
    double *avgPi = p->avgPi, *avgLnPi = p->avgLnPi, **avgA = p->avgA, **avgLnA = p->avgLnA;
    tempGlobalStats *gs = (tempGlobalStats*)ivs->stats;
    double *NiiR = gs->NiiR, **NijR = gs->NijR, *z1iR = gs->z1iR, dR = (double)(xns->R);
    int sNo = gv->sNo;
    int i, j;

    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] );
        }
    }
}

double varLowerBound__( xn, gv, iv )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
{
    // Calculation of lower bound, depending on the model.

    // calculation of lower bound of common HMM
    tempParameters *p = (tempParameters*)gv->params;
    tempStats *s = (tempStats*)iv->stats;
    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 *avgLnPi = p->avgLnPi, **avgLnA = p->avgLnA;
    double *Nii = s->Nii, **Nij = s->Nij;
    size_t n;
    int i, j;

    double lnpPi = gsl_sf_lngamma(sumUPi);
    double lnpA = 0.0;
    double lnpPhi = 0.0;            // ln p(phi) for model-specific parameters
    double lnqPi = gsl_sf_lngamma(sumUPi + 1.0);
    double lnqA = 0.0;
    double lnqPhi = 0.0;            // ln q(phi) for model-specific parameters
    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 + lnpPhi;
    val -= lnqPi + lnqA + lnqPhi;
    val += lnpX;
    val += log(gsl_sf_fact(sNo));

    return val;
}

double varLowerBoundG__( xns, gv, ivs )
xnDataBundle *xns;
globalVars *gv;
indVarBundle *ivs;
{
    // Calculation of lower bound for global analysis, depending on the model.

    // calculation of lower bound for global analysis of common HMM
    int sNo = gv->sNo, rNo = xns->R;
    tempParameters *p = (tempParameters*)gv->params;
    double *uPiArr = p->uPiArr, sumUPi = p->sumUPi;
    double **uAMat = p->uAMat, *sumUAArr = p->sumUAArr;
    double *avgLnPi = p->avgLnPi, **avgLnA = p->avgLnA;
    tempGlobalStats *gs = (tempGlobalStats*)ivs->stats;
    double *NiiR = gs->NiiR, **NijR = gs->NijR, *z1iR = gs->z1iR, dR = (double)xns->R;
    size_t n;
    int i, j, r;

    double lnpPi = gsl_sf_lngamma(sumUPi);
    double lnpA = 0.0;
    double lnpPhi = 0.0;            // ln p(phi) for model-specific parameters
    double lnqPi = gsl_sf_lngamma(sumUPi + dR);
    double lnqA = 0.0;
    double lnqPhi = 0.0;            // ln q(phi) for model-specific parameters
    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 + lnpPhi;
    val -= lnqPi + lnqA + lnqPhi;
    val += lnpX;
    val += log(gsl_sf_fact(sNo));

    return val;
}

void reorderParameters__( xn, gv, iv )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
{
    // Reorder states in the final result, for example, in the order of intensity.

    // an example for common HMM.
    tempParameters *p = (tempParameters*)gv->params;
    tempStats *s = (tempStats*)iv->stats;
    size_t dLen = xn->N;
    int sNo = gv->sNo;
    double **gmMat = iv->gmMat, ***xiMat = iv->xiMat;
    double *avgPi = p->avgPi, *avgLnPi = p->avgLnPi, **avgA = p->avgA, **avgLnA = p->avgLnA;
    double *Ni = s->Ni;
    size_t n;
    int i, j;

    int *index = (int*)malloc( sNo * sizeof(int) );
    double *store = (double*)malloc( sNo * sizeof(double) );
    double **s2D = (double**)malloc( sNo * sizeof(double*) );
    for( i = 0 ; i < sNo ; i++ )
    {   s2D[i] = (double*)malloc( MAX(sNo,2) * sizeof(double) );   }

    // index indicates order of avgPi values (0=biggest avgMu -- sNo=smallest avgMu).
    for( i = 0 ; i < sNo ; i++ ){
        index[i] = sNo - 1;
        for( j = 0 ; j < sNo ; j++ ){
            if( j != i ){
                if( avgPi[i] < avgPi[j] ){
                    index[i]--;
                } else if( avgPi[i] == avgPi[j] ){
                    if( j > i )
                    {   index[i]--;   }
                }
            }
        }
    }

    for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = avgPi[i];   }
    for( i = 0 ; i < sNo ; i++ ){   avgPi[i] = store[i];   }

    for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = avgLnPi[i];   }
    for( i = 0 ; i < sNo ; i++ ){   avgLnPi[i] = store[i];   }

    for( j = 0 ; j < sNo ; j++ ){
        for( i = 0 ; i < sNo ; i++ ){   s2D[index[i]][index[j]] = avgA[i][j];   }
    }
    for( j = 0 ; j < sNo ; j++ ){
        for( i = 0 ; i < sNo ; i++ ){   avgA[i][j] = s2D[i][j];   }
    }

    for( j = 0 ; j < sNo ; j++ ){
        for( i = 0 ; i < sNo ; i++ ){   s2D[index[i]][index[j]] = avgLnA[i][j];   }
    }
    for( j = 0 ; j < sNo ; j++ ){
        for( i = 0 ; i < sNo ; i++ ){   avgLnA[i][j] = s2D[i][j];   }
    }

    for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = Ni[i];   }
    for( i = 0 ; i < sNo ; i++ ){   Ni[i] = store[i];   }

    for( n = 0 ; n < dLen ; n++ ){
        for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = gmMat[n][i];   }
        for( i = 0 ; i < sNo ; i++ ){   gmMat[n][i] = store[i];   }
    }

    for( n = 0 ; n < dLen ; n++ ){
        for( j = 0 ; j < sNo ; j++ ){
            for( i = 0 ; i < sNo ; i++ ){   s2D[index[i]][index[j]] = xiMat[n][i][j];   }
        }
        for( j = 0 ; j < sNo ; j++ ){
            for( i = 0 ; i < sNo ; i++ ){   xiMat[n][i][j] = s2D[i][j];   }
        }
    }

    for( i = 0 ; i < sNo ; i++ ){   free( s2D[i] );   }
    free( s2D );
    free( store );
    free( index );
}

void reorderParametersG__( xns, gv, ivs )
xnDataBundle *xns;
globalVars *gv;
indVarBundle *ivs;
{
    // Reorder states in the final result, for example, in the order of intensity.

    // an example for global analysis of common HMM.
    size_t dLen;
    int sNo = gv->sNo, rNo = xns->R;
    tempParameters *p = (tempParameters*)gv->params;
    double *avgPi = p->avgPi, *avgLnPi = p->avgLnPi, **avgA = p->avgA, **avgLnA = p->avgLnA;

    size_t n;
    int i, j, r;

    int *index = (int*)malloc( sNo * sizeof(int) );
    double *store = (double*)malloc( sNo * sizeof(double) );
    double **s2D = (double**)malloc( sNo * sizeof(double*) );
    for( i = 0 ; i < sNo ; i++ )
    {   s2D[i] = (double*)malloc( MAX(sNo,2) * sizeof(double) );   }

    // index indicates order of avgPi values (0=biggest avgMu -- sNo=smallest avgMu).
    for( i = 0 ; i < sNo ; i++ ){
        index[i] = sNo - 1;
        for( j = 0 ; j < sNo ; j++ ){
            if( j != i ){
                if( avgPi[i] < avgPi[j] ){
                    index[i]--;
                } else if( avgPi[i] == avgPi[j] ){
                    if( j > i )
                    {   index[i]--;   }
                }
            }
        }
    }

    for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = avgPi[i];   }
    for( i = 0 ; i < sNo ; i++ ){   avgPi[i] = store[i];   }

    for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = avgLnPi[i];   }
    for( i = 0 ; i < sNo ; i++ ){   avgLnPi[i] = store[i];   }

    for( j = 0 ; j < sNo ; j++ ){
        for( i = 0 ; i < sNo ; i++ ){   s2D[index[i]][index[j]] = avgA[i][j];   }
    }
    for( j = 0 ; j < sNo ; j++ ){
        for( i = 0 ; i < sNo ; i++ ){   avgA[i][j] = s2D[i][j];   }
    }

    for( j = 0 ; j < sNo ; j++ ){
        for( i = 0 ; i < sNo ; i++ ){   s2D[index[i]][index[j]] = avgLnA[i][j];   }
    }
    for( j = 0 ; j < sNo ; j++ ){
        for( i = 0 ; i < sNo ; i++ ){   avgLnA[i][j] = s2D[i][j];   }
    }

    double *NiR = ((tempGlobalStats*)ivs->stats)->NiR;
    for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = NiR[i];   }
    for( i = 0 ; i < sNo ; i++ ){   NiR[i] = store[i];   }

    for( r = 0 ; r < rNo ; r++ ){
        double **gmMat = ivs->indVars[r]->gmMat;
        dLen = xns->xn[r]->N;
        for( n = 0 ; n < dLen ; n++ ){
            for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = gmMat[n][i];   }
            for( i = 0 ; i < sNo ; i++ ){   gmMat[n][i] = store[i];   }
        }
    }

    for( r = 0 ; r < rNo ; r++ ){
        double ***xiMat = ivs->indVars[r]->xiMat;
        dLen = xns->xn[r]->N;
        for( n = 0 ; n < dLen ; n++ ){
            for( j = 0 ; j < sNo ; j++ ){
                for( i = 0 ; i < sNo ; i++ ){   s2D[index[i]][index[j]] = xiMat[n][i][j];   }
            }
            for( j = 0 ; j < sNo ; j++ ){
                for( i = 0 ; i < sNo ; i++ ){   xiMat[n][i][j] = s2D[i][j];   }
            }
        }
    }

    for( i = 0 ; i < sNo ; i++ ){   free( s2D[i] );   }
    free( s2D );
    free( store );
    free( index );
}

void output__Results( xn, gv, iv, logFP )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
FILE *logFP;
{
    // Output the final result (after reordering) in any format, e.g. stdout, stderr or files.

    // an example for common HMM parameters.
    tempParameters *p = (tempParameters*)gv->params;
    int sNo = gv->sNo;

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

    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, "pi");
        for( i = 0 ; i < sNo ; i++ )
        {   fprintf(fp, ", A%dx", i);   }
        fprintf(fp, "\n");

        for( i = 0 ; i < sNo ; i++ ){
            fprintf(fp, "%g", 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 output__ResultsG( xns, gv, ivs, logFP )
xnDataBundle *xns;
globalVars *gv;
indVarBundle *ivs;
FILE *logFP;
{
    // Output the final result (after reordering) in any format, e.g. stdout, stderr or files.

    // an example for parameters of global analysis of common HMM.
    int sNo = gv->sNo, rNo = xns->R;
    tempParameters *p = (tempParameters*)gv->params;
    int i, j, r;

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

    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, "pi");
        for( i = 0 ; i < sNo ; i++ )
        {   fprintf(fp, ", A%dx", i);   }
        fprintf(fp, "\n");

        for( i = 0 ; i < sNo ; i++ ){
            fprintf(fp, "%g", 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);
    }

    for( r = 0 ; r < rNo ; r++ ){
        xnDataSet *xn = xns->xn[r];
        indVars *iv = ivs->indVars[r];

        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);
        }
    }
}

//

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:
/*
 * Copyright (c) 2016-2021 lymastee, All rights reserved.
 * Contact: lymastee@hotmail.com
 *
 * This file is part of the gslib project.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 * SOFTWARE.
 */

#pragma once

#ifndef base_type_2aa01495_7feb_4c5e_ae38_9c775c29c9f1_h
#define base_type_2aa01495_7feb_4c5e_ae38_9c775c29c9f1_h

#include <assert.h>
#include <stdint.h>
#include <gslib/config.h>

__gslib_begin__

#ifdef _GS_X86

typedef int8_t int8;
typedef int16_t int16;
typedef int32_t int32;
typedef int64_t int64;
typedef wchar_t wchar;
typedef float real32;
typedef double real;
typedef uint8_t byte;
typedef uint16_t word;
typedef uint32_t dword;
typedef uint64_t qword;
typedef uint32_t uint;
typedef uint16_t uint16;
typedef uint32_t uint32;
typedef uint64_t uint64;
typedef void* addrptr;

#ifdef _UNICODE
typedef wchar gchar;
#else
typedef char gchar;
#endif

template<class t>
struct typeof { typedef void type; };

/* virtual class ptr convert tools */
template<class ccp, class bcp>
inline int virtual_bias()
{
    static typeof<ccp>::type c;
    static const int bias = (int)&c - (int)static_cast<bcp>(&c);
    return bias;
}

template<class ccp, class bcp>
inline ccp virtual_cast(bcp p)
{
    byte* ptr = (byte*)p;
    ptr += virtual_bias<ccp, bcp>();
    return reinterpret_cast<ccp>(ptr);
}

#endif  /* end of _GS_X86 */

template<class _ty>
inline _ty gs_min(_ty a, _ty b) { return a < b ? a : b; }
template<class _ty>
inline _ty gs_max(_ty a, _ty b) { return a > b ? a : b; }
template<class _ty>
inline _ty gs_clamp(_ty v, _ty a, _ty b)
{
    assert(a <= b);
    return gs_min(b, gs_max(a, v));
}

template<class _ty>
inline void gs_swap(_ty& a, _ty& b)
{
    auto t = a;
    a = b;
    b = t;
}

__gslib_end__

#endif

ext: h

max_stars_repo_name: lymastee/gslib

max_issues_repo_name: lymastee/gslib

max_forks_repo_name: lymastee/gslib


text:
// Copyright András Vukics 2021. Distributed under the Boost Software License, Version 1.0. (See accompanying file LICENSE.txt)
/// \briefFileDefault
#ifndef CPPQEDCORE_UTILS_ODE_GSL_H_INCLUDED
#define CPPQEDCORE_UTILS_ODE_GSL_H_INCLUDED

#include "ArrayTraits.h"
#include "ODE.h"

#include <gsl/gsl_errno.h>
#include <gsl/gsl_odeiv2.h>

#include <cstddef>

namespace cppqedutils {

namespace ode_engine {

template <typename StateType>
class ControlledErrorStepperGSL
{
public:
  using state_type=StateType;
  using deriv_type=state_type;
  using time_type=double;
  using value_type=ElementType_t<StateType>;
  using stepper_category=bno::controlled_stepper_tag;

  ControlledErrorStepperGSL(double epsAbs, double epsRel)
    : control_{gsl_odeiv2_control_standard_new(epsAbs,epsRel,1,1), [] (auto* ptr) {if (ptr) gsl_odeiv2_control_free(ptr);}} {}

private:
  using ValueVector=std::vector<double>;
  using Aux_t=std::tuple<SystemFunctional_t<double,StateType>&,Extents_t<Rank_v<StateType>>>;

public:
  /// Realizes part of the logic of [gsl_odeiv2_evolve_apply](http://www.gnu.org/software/gsl/doc/html/ode-initval.html#c.gsl_odeiv2_evolve_apply)
  /**
   * Differences:
   * - First-same-as-last steppers not supported (at the moment only RKCK is supported)
   * - There is no guard against the stepsize becoming infinitesimal
   */
  template <typename System, typename S>
  auto try_step(System sys, S&& stateInOut, double& time, double& dtTry)
  {
    auto totalExtent=Size<StateType>::_(stateInOut);

    if (!step_) { // this means that tryStep is called for the first time
      step_.reset(gsl_odeiv2_step_alloc(gsl_odeiv2_step_rkck,totalExtent), [] (auto* ptr) {if (ptr) gsl_odeiv2_step_free(ptr);});
      yerr_=std::make_shared<ValueVector>(totalExtent);
      dydt_out_=std::make_shared<ValueVector>(totalExtent);
    }
#ifndef NDEBUG
    else if (!IsStorageContiguous<StateType>::_(stateInOut)) {
      throw NonContiguousStorageException{"ControlledErrorStepperGSL::tryStep"};
    }
    else if (yerr_->size()!=Size<StateType>::_(stateInOut) || dydt_out_->size()!=Size<StateType>::_(stateInOut)) {
      throw std::runtime_error("Dimensionality mismatch in ControlledErrorStepperGSL::tryStep");
    }
#endif // NDEBUG

    SystemFunctional_t<double,StateType> sysVariable=sys; // it’s important to convert `sys` (which can be a lambda) to a variable of known type

    Aux_t aux{sysVariable,Extents<StateType>::_(stateInOut)};

    gsl_odeiv2_system dydt{ControlledErrorStepperGSL::lowLevelSystemFunction,nullptr,totalExtent,&aux};

    auto step_status = gsl_odeiv2_step_apply (step_.get(), time, dtTry, Data<StateType>::_(stateInOut), yerr_->data(), nullptr, dydt_out_->data(), &dydt);

    if (step_status == GSL_EFAULT || step_status == GSL_EBADFUNC) throw std::runtime_error("GSL bad step_status");

    if (step_status != GSL_SUCCESS) {
      dtTry *= 0.5;
      return bno::fail;
    }
    else {
      time+=dtTry;
      /* const auto hadjust_status = */
      gsl_odeiv2_control_hadjust (control_.get(), step_.get(), Data<StateType>::_(stateInOut), yerr_->data(), dydt_out_->data(), &dtTry);
      // if (hadjust_status == GSL_ODEIV_HADJ_DEC) throw std::runtime_error{"Unexpected change of dtTry in ControlledErrorStepperGSL::tryStep"};
      return bno::success;
    }
  }

private:
  static int lowLevelSystemFunction(double t, const double* y, double* dydt, void* aux)
  {
    auto [highLevelSystemFunction,arrayExtents] = *static_cast<Aux_t*>(aux);

    const StateType yInterface(Create_c<StateType>::_(y,arrayExtents));

    StateType dydtInterface(Create<StateType>::_(dydt,arrayExtents));

    highLevelSystemFunction(yInterface,dydtInterface,t);

    return GSL_SUCCESS;
  }

  std::shared_ptr<gsl_odeiv2_step> step_;

  const std::shared_ptr<gsl_odeiv2_control> control_;

  // Needs to be easily copyable, as Steppers are carried around by value
  std::shared_ptr<ValueVector> yerr_, dydt_out_;

};

template <typename StateType>
struct MakeControlledErrorStepper<ControlledErrorStepperGSL<StateType>>
{
  static auto _(double epsRel, double epsAbs)
  {
    return ControlledErrorStepperWithParameters<ControlledErrorStepperGSL<StateType>>{
      {epsRel,epsAbs},epsRel,epsAbs};
  }

  template <typename ParsBase>
  static auto _(const Pars<ParsBase>& p)
  {
    return _(p.epsRel,p.epsAbs);
  }

};

template <typename StateType>
inline const std::string StepperDescriptor<ControlledErrorStepperGSL<StateType>> = "GSL RKCK controlled stepper";

} // ode_engine

template <typename StateType>
using ODE_EngineGSL = ode_engine::Base<ode_engine::ControlledErrorStepperGSL<StateType>>;

} // cppqedutils

#endif // CPPQEDCORE_UTILS_ODE_GSL_H_INCLUDED

max_stars_repo_name: vukics/cppqed

max_issues_repo_name: vukics/cppqed

max_forks_repo_name: vukics/cppqed


text:
#include "cconfigspace_internal.h"
#include <stdlib.h>
#include <gsl/gsl_rng.h>

const ccs_datum_t ccs_none = CCS_NONE_VAL;
const ccs_datum_t ccs_inactive = CCS_INACTIVE_VAL;
const ccs_datum_t ccs_true = CCS_TRUE_VAL;
const ccs_datum_t ccs_false = CCS_FALSE_VAL;
const ccs_version_t ccs_version = { 0, 1, 0, 0 };

ccs_result_t
ccs_init() {
	gsl_rng_env_setup();
	return CCS_SUCCESS;
}

ccs_result_t
ccs_fini() {
	return CCS_SUCCESS;
}

ccs_version_t
ccs_get_version() {
	return ccs_version;
}

ccs_result_t
ccs_retain_object(ccs_object_t object) {
	_ccs_object_internal_t *obj = (_ccs_object_internal_t *)object;
        if (!obj || obj->refcount <= 0)
		return -CCS_INVALID_OBJECT;
	obj->refcount += 1;
	return CCS_SUCCESS;
}

ccs_result_t
ccs_release_object(ccs_object_t object) {
	_ccs_object_internal_t *obj = (_ccs_object_internal_t *)object;
	if (!obj || obj->refcount <= 0)
		return -CCS_INVALID_OBJECT;
	obj->refcount -= 1;
	if (obj->refcount == 0) {
		if (obj->callbacks) {
			_ccs_object_callback_t *cb = NULL;
			while ( (cb = (_ccs_object_callback_t *)
			              utarray_prev(obj->callbacks, cb)) ) {
				cb->callback(object, cb->user_data);
			}
			utarray_free(obj->callbacks);
		}
		CCS_VALIDATE(obj->ops->del(object));
		free(object);
	}
	return CCS_SUCCESS;
}

ccs_result_t
ccs_object_get_type(ccs_object_t       object,
                     ccs_object_type_t *type_ret) {
	_ccs_object_internal_t *obj = (_ccs_object_internal_t *)object;
	if (!obj)
		return -CCS_INVALID_OBJECT;
	CCS_CHECK_PTR(type_ret);
	*type_ret = obj->type;
	return CCS_SUCCESS;
}

ccs_result_t
ccs_object_get_refcount(ccs_object_t  object,
                         int32_t      *refcount_ret) {
	_ccs_object_internal_t *obj = (_ccs_object_internal_t *)object;
	if (!obj)
		return -CCS_INVALID_OBJECT;
	CCS_CHECK_PTR(refcount_ret);
	*refcount_ret = obj->refcount;
	return CCS_SUCCESS;
}

static const UT_icd _object_callback_icd = {
	sizeof(_ccs_object_callback_t),
	NULL,
	NULL,
	NULL
};

ccs_result_t
ccs_object_set_destroy_callback(ccs_object_t  object,
                                void (*callback)(
                                  ccs_object_t object,
                                  void *user_data),
                                void *user_data) {
	_ccs_object_internal_t *obj = (_ccs_object_internal_t *)object;
	if (!obj)
		return -CCS_INVALID_OBJECT;
	if (!callback)
		return -CCS_INVALID_VALUE;
	if (!obj->callbacks)
		utarray_new(obj->callbacks, &_object_callback_icd);

	_ccs_object_callback_t cb = { callback, user_data };
	utarray_push_back(obj->callbacks, &cb);
	return CCS_SUCCESS;
}

ext: c

max_stars_repo_name: deephyper/CCS

max_issues_repo_name: deephyper/CCS

max_forks_repo_name: deephyper/CCS


text:
static char help[] = "1D reaction-diffusion problem with DMDA and SNES.  Option prefix -rct_.\n\n";

#include <petsc.h>

typedef struct {
    PetscReal  rho, M, alpha, beta;
    PetscBool  noRinJ;
} AppCtx;

extern PetscReal f_source(PetscReal);
extern PetscErrorCode InitialAndExact(DMDALocalInfo*, PetscReal*, PetscReal*, AppCtx*);
extern PetscErrorCode FormFunctionLocal(DMDALocalInfo*, PetscReal*, PetscReal*, AppCtx*);
extern PetscErrorCode FormJacobianLocal(DMDALocalInfo*, PetscReal*, Mat, Mat, AppCtx*);

//STARTMAIN
int main(int argc,char **args) {
  PetscErrorCode ierr;
  DM            da;
  SNES          snes;
  AppCtx        user;
  Vec           u, uexact;
  PetscReal     errnorm, *au, *auex;
  DMDALocalInfo info;

  ierr = PetscInitialize(&argc,&args,NULL,help); if (ierr) return ierr;
  user.rho   = 10.0;
  user.M     = PetscSqr(user.rho / 12.0);
  user.alpha = user.M;
  user.beta  = 16.0 * user.M;
  user.noRinJ = PETSC_FALSE;

  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"rct_","options for reaction",""); CHKERRQ(ierr);
  ierr = PetscOptionsBool("-noRinJ","do not include R(u) term in Jacobian",
      "reaction.c",user.noRinJ,&(user.noRinJ),NULL); CHKERRQ(ierr);
  ierr = PetscOptionsEnd(); CHKERRQ(ierr);

  ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,9,1,1,NULL,&da); CHKERRQ(ierr);
  ierr = DMSetFromOptions(da); CHKERRQ(ierr);
  ierr = DMSetUp(da); CHKERRQ(ierr);
  ierr = DMSetApplicationContext(da,&user); CHKERRQ(ierr);

  ierr = DMCreateGlobalVector(da,&u); CHKERRQ(ierr);
  ierr = VecDuplicate(u,&uexact); CHKERRQ(ierr);
  ierr = DMDAVecGetArray(da,u,&au); CHKERRQ(ierr);

  ierr = DMDAGetLocalInfo(da,&info); CHKERRQ(ierr);
  ierr = DMDAVecGetArray(da,uexact,&auex); CHKERRQ(ierr);
  ierr = InitialAndExact(&info,au,auex,&user); CHKERRQ(ierr);
  ierr = DMDAVecRestoreArray(da,u,&au); CHKERRQ(ierr);
  ierr = DMDAVecRestoreArray(da,uexact,&auex); CHKERRQ(ierr);

  ierr = SNESCreate(PETSC_COMM_WORLD,&snes); CHKERRQ(ierr);
  ierr = SNESSetDM(snes,da); CHKERRQ(ierr);
  ierr = DMDASNESSetFunctionLocal(da,INSERT_VALUES,
             (DMDASNESFunction)FormFunctionLocal,&user); CHKERRQ(ierr);
  ierr = DMDASNESSetJacobianLocal(da,
             (DMDASNESJacobian)FormJacobianLocal,&user); CHKERRQ(ierr);
  ierr = SNESSetFromOptions(snes); CHKERRQ(ierr);

  ierr = SNESSolve(snes,NULL,u); CHKERRQ(ierr);

  ierr = VecAXPY(u,-1.0,uexact); CHKERRQ(ierr);    // u <- u + (-1.0) uexact
  ierr = VecNorm(u,NORM_INFINITY,&errnorm); CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD,
      "on %d point grid:  |u-u_exact|_inf = %g\n",info.mx,errnorm); CHKERRQ(ierr);

  VecDestroy(&u);  VecDestroy(&uexact);
  SNESDestroy(&snes);  DMDestroy(&da);
  return PetscFinalize();
}
//ENDMAIN

PetscReal f_source(PetscReal x) {
    return 0.0;
}

PetscErrorCode InitialAndExact(DMDALocalInfo *info, PetscReal *u0,
                               PetscReal *uex, AppCtx *user) {
    PetscInt   i;
    PetscReal  h = 1.0 / (info->mx-1), x;
    for (i=info->xs; i<info->xs+info->xm; i++) {
        x = h * i;
        u0[i]  = user->alpha * (1.0 - x) + user->beta * x;
        uex[i] = user->M * PetscPowReal(x + 1.0,4.0);
    }
    return 0;
}

//STARTFUNCTIONS
PetscErrorCode FormFunctionLocal(DMDALocalInfo *info, PetscReal *u,
                                 PetscReal *FF, AppCtx *user) {
    PetscInt   i;
    PetscReal  h = 1.0 / (info->mx-1), x, R;
    for (i=info->xs; i<info->xs+info->xm; i++) {
        if (i == 0) {
            FF[i] = u[i] - user->alpha;
        } else if (i == info->mx-1) {
            FF[i] = u[i] - user->beta;
        } else {  // interior location
            if (i == 1) {
                FF[i] = - u[i+1] + 2.0 * u[i] - user->alpha;
            } else if (i == info->mx-2) {
                FF[i] = - user->beta + 2.0 * u[i] - u[i-1];
            } else {
                FF[i] = - u[i+1] + 2.0 * u[i] - u[i-1];
            }
            R = - user->rho * PetscSqrtReal(u[i]);
            x = i * h;
            FF[i] -= h*h * (R + f_source(x));
        }
    }
    return 0;
}

PetscErrorCode FormJacobianLocal(DMDALocalInfo *info, PetscReal *u,
                                 Mat J, Mat P, AppCtx *user) {
    PetscErrorCode ierr;
    PetscInt   i, col[3];
    PetscReal  h = 1.0 / (info->mx-1), dRdu, v[3];
    for (i=info->xs; i<info->xs+info->xm; i++) {
        if ((i == 0) | (i == info->mx-1)) {
            v[0] = 1.0;
            ierr = MatSetValues(P,1,&i,1,&i,v,INSERT_VALUES); CHKERRQ(ierr);
        } else {
            col[0] = i;
            v[0] = 2.0;
            if (!user->noRinJ) {
                dRdu = - (user->rho / 2.0) / PetscSqrtReal(u[i]);
                v[0] -= h*h * dRdu;
            }
            col[1] = i-1;   v[1] = (i > 1) ? - 1.0 : 0.0;
            col[2] = i+1;   v[2] = (i < info->mx-2) ? - 1.0 : 0.0;
            ierr = MatSetValues(P,1,&i,3,col,v,INSERT_VALUES); CHKERRQ(ierr);
        }
    }
    ierr = MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
    ierr = MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
    if (J != P) {
        ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
        ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    }
    return 0;
}
//ENDFUNCTIONS


text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <string.h>

double likelihood(double *y_obs, double* y_model);
void my_model(double *y_init, double *x_obs, double a, double b, double c, double d);
int minimo (double *arr, int nsteps, int nburn);
double std_dev(double *arr, int n_steps, int n_burn);

int main(int argc, char **argv){

  int n_steps = atoi(argv[1]);
  int n_burn = atoi(argv[2]);
  int npoints = n_steps;

  double *x_obs;
  double *y_obs;
  double *a_walk;
  double *b_walk;
  double *c_walk;
  double *d_walk;
  double *l_walk;
  double *y_init;
  double *y_prime;
  FILE *in;
  double var0;
  double var1;
  double var2;
  double var3;
  double var4;

  x_obs  = malloc(npoints * sizeof(double));
  y_obs  = malloc(npoints * sizeof(double));
  a_walk = malloc(npoints * sizeof(double));
  b_walk = malloc(npoints * sizeof(double));
  c_walk = malloc(npoints * sizeof(double));
  d_walk = malloc(npoints * sizeof(double));
  l_walk = malloc(npoints * sizeof(double));
  y_init = malloc(npoints * sizeof(double));
  y_prime = malloc(npoints * sizeof(double));

  in = fopen("monthrg.dat", "r");
  if(!in){
    printf("problems opening the file %s\n", "monthrg.dat");
    exit(1);
  }

  int i;
  for(i=0;i<4141;i++){
    fscanf(in,"%lf %lf %lf %lf %lf\n", &var0, &var1, &var2, &var3, &var4);
    if (var3!=-99){
      x_obs[i] = var0 + var1/12.0;
      y_obs[i] = var3;
      // printf("%f %f\n",x_obs[i], var3);
    }

  }
  fclose(in);

  char output1[512];
  char output2[512];
  char* line = "_";
  char* mcmc_solar = "mcmc_solar";
  char* dat = ".dat";
  char name[512] = "";
  strcat(name,mcmc_solar);
  strcat(name,line);
  strcat(name,argv[1]);
  strcat(name,line);
  strcat(name,argv[2]);
  strcat(name,dat);
  //printf("%s\n", name);

  srand((unsigned)time(NULL));

  FILE *fp;
  fp = fopen(name, "ab+");

  a_walk[0]=((double)rand() / (double) (RAND_MAX/50.0));
  b_walk[0]=((double)rand() / (double) (RAND_MAX/50.0));
  c_walk[0]=((double)rand() / (double) (RAND_MAX/50.0));
  d_walk[0]=((double)rand() / (double) (RAND_MAX/50.0));

  // printf("%f %f %f %f\n",a_walk[0], b_walk[0], c_walk[0], d_walk[0]);
  my_model (y_init, x_obs, a_walk[0], b_walk[0], c_walk[0], d_walk[0]);
  l_walk[0] = likelihood(y_obs, y_init);
  // printf("%f\n",l_walk[0]);

  int j;

   const gsl_rng_type * T;
   gsl_rng * r;
   gsl_rng_env_setup();
   T = gsl_rng_default;
   r = gsl_rng_alloc(T);

   double a_prime;
   double b_prime;
   double c_prime;
   double d_prime;
   double l_prime;
   double l_init;
   double gamma;
   double beta;

   for (j = 0; j<n_steps; j++){

    a_prime =(double)gsl_ran_gaussian(r, 0.1)+a_walk[j];
    b_prime =(double)gsl_ran_gaussian(r, 0.1)+b_walk[j];
    c_prime =(double)gsl_ran_gaussian(r, 0.1)+c_walk[j];
    d_prime =(double)gsl_ran_gaussian(r, 0.1)+d_walk[j];

    // printf("%f %f %f %f\n",a_walk[j], b_walk[j], c_walk[j], d_walk[j]);

    my_model(y_init, x_obs, a_walk[j], b_walk[j], c_walk[j], d_walk[j]);
    my_model(y_prime, x_obs, a_prime, b_prime, c_prime, d_prime);

    l_prime = likelihood(y_obs, y_prime);
    l_init = likelihood(y_obs, y_init);

    gamma = -l_prime + l_init;
    // printf("%f\n",gamma);

    if(gamma>=0.0){

      a_walk[j+1] = a_prime;
      b_walk[j+1] = b_prime;
      c_walk[j+1] = c_prime;
      d_walk[j+1] = d_prime;
      l_walk[j+1] = l_prime;

    }
    else{

      beta = ((double)rand() / (double) (RAND_MAX/50.0));

      if(beta<=exp(gamma)){

      a_walk[j+1] = a_prime;
      b_walk[j+1] = b_prime;
      c_walk[j+1] = c_prime;
      d_walk[j+1] = d_prime;
      l_walk[j+1] = l_prime;

      }
      else{

      a_walk[j+1] = a_walk[j];
      b_walk[j+1] = b_walk[j];
      c_walk[j+1] = c_walk[j];
      d_walk[j+1] = d_walk[j];
      l_walk[j+1] = l_init;

      }

    }
    // printf("%f %f %f %f\n",a_walk[j], b_walk[j], c_walk[j], d_walk[j]);

  }
  double *a_burned;
  double *b_burned;
  double *c_burned;
  double *d_burned;
  double *l_burned;
  double *x_burned;

  a_burned = malloc(npoints * sizeof(double));
  b_burned = malloc(npoints * sizeof(double));
  c_burned = malloc(npoints * sizeof(double));
  d_burned = malloc(npoints * sizeof(double));
  l_burned = malloc(npoints * sizeof(double));

  int co;

  for (co = 0; co < n_steps-n_burn; co++){

    a_burned[co] = a_walk[n_burn + co];
    b_burned[co] = b_walk[n_burn + co];
    c_burned[co] = c_walk[n_burn + co];
    d_burned[co] = d_walk[n_burn + co];
    l_burned[co] = l_walk[n_burn + co];
    fprintf(fp, "%f %f %f %f %f\n",a_burned[co], b_burned[co], c_burned[co], d_burned[co], l_burned[co]);
    //printf("%f %f %f %f %f\n",a_walk[co], b_walk[co], c_walk[co], d_walk[co], l_walk[co]);

  }

  int max_likelihood_id = minimo(l_burned, n_steps, n_burn);
  // printf("%i\n", max_likelihood_id);
  double best_a = a_burned[max_likelihood_id];
  double best_b = b_burned[max_likelihood_id];
  double best_c = c_burned[max_likelihood_id];
  double best_d = d_burned[max_likelihood_id];
  // printf("%f %f %f %f\n",a_burned[max_likelihood_id], b_burned[max_likelihood_id], c_burned[max_likelihood_id], d_burned[max_likelihood_id]);

  printf("a = %f +/- %f\n",best_a, std_dev(a_burned, n_steps, n_burn));
  printf("b = %f +/- %f\n",best_b, std_dev(b_burned, n_steps, n_burn));
  printf("c = %f +/- %f\n",best_c, std_dev(c_burned, n_steps, n_burn));
  printf("d = %f +/- %f\n",best_d, std_dev(d_burned, n_steps, n_burn));

 fclose(fp);

  gsl_rng_free(r);
  return 0;
}

double likelihood(double *y_obs, double* y_model){

  double chi_squared = 0.0;
  int counter;

  for (counter = 0; counter < 4141; counter++){
    chi_squared = chi_squared + pow((y_obs[counter]-y_model[counter]),2);
  }
  return (1.0/2.0)*chi_squared;
}

void my_model(double *y_init, double *x_obs, double a, double b, double c, double d){

  int counter2;
  double pi = 3.14159265359;

  for (counter2 = 0; counter2<4141; counter2++){

    y_init[counter2] = a*cos((2*pi/d)*x_obs[counter2] + b) + c;

  }
}

int minimo (double *arr, int nsteps, int nburn){

    int c = 0;
    double min = arr[0];
    int minimo = 0;

    for (c = 1; c<nsteps-nburn; c++){

      if (arr[c] < min){

	min = arr[c];
	minimo = c;

      }

    }
    return minimo;

}

double std_dev(double *arr, int n_steps, int n_burn){

  double mean = 0.0, sum_dev = 0.0;

  int k;
  int n = n_steps-n_burn;

  for (k=0; k<n; k++){

    mean+=arr[k];

  }
  mean = mean/n;
  for (k = 0; k<n; k++){

    sum_dev+=(arr[k]-mean)*(arr[k]-mean);

  }
  return pow(sum_dev/n,0.5);

}


text:
#ifndef PETSCWRAPPER_H
#define PETSCWRAPPER_H

#include <petsc.h>
#include "../TPLs/eigen-git-mirror/Eigen/Eigen"
#include <iostream>

using namespace std;

//==============================================================================

int initPETScMat(Mat * A, int squareSize, int nonZeros);
int initPETScRectMat(Mat * A, int rows, int cols, int nonZeros);
int initPETScVec(Vec * A, int size);
int eigenVecToPETScVec(Eigen::VectorXd * x_e,Vec * x_p);
int petscVecToEigenVec(Vec * x_p,Eigen::VectorXd * x_e);

//==============================================================================

#endif

max_stars_repo_name: khurrumsaleem/QuasiMolto

max_issues_repo_name: khurrumsaleem/QuasiMolto

max_forks_repo_name: khurrumsaleem/QuasiMolto


text:
#ifndef __DREAM_H__
#define __DREAM_H__

#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <string>
#include <map>
#include <vector>
#include <algorithm>
using namespace std;

#include <gsl/gsl_statistics.h>
#include <gsl/gsl_math.h>

#include <rapidjson/document.h>
#include <rng/RngStream.h>
#include "array.h"

typedef double (*LikFun)(int chain_id, int gen, const double* state,
                         const void* pars, bool recalc);
typedef void (*ReportFun)(int gen, int chain_id, int n,
                          const double* state, double lik,
                          int burnIn, const void* pars);

typedef struct t_dream_pars {
  int vflag;                 /* vebose flag */
  int maxEvals;              /* max number of function evaluations */
  int optimAlg;              /* lock variables */
  int numChains;
  string fn;                 /* input filename */
  string out_fn;             /* output filename */
  int appendFile;            /* continue from previous state */
  int report_interval;       /* report interval for state */
  int diagnostics;           /* report diagnostics at the end of the run */
  int burnIn;                /* number of steps for which to run an adaptive proposal size */
  int recalcLik;             /* recalculate likelihood of previously evaluated states */

  // DREAM variables
  int collapseOutliers;
  int gelmanEvals;
  int loopSteps;
  int deltaMax;
  int pCR_update;
  int nCR;

  double noise;
  double bstar_zero;
  double scaleReductionCrit;
  double reenterBurnin;

  int nfree;
  size_t nvar;
  double* varLo;
  double* varHi;
  double* varInit;
  int* varLock;
  string* varName;
  char* scale;

  LikFun fun;
  void* funPars;
  ReportFun rfun;
} dream_pars;

void dream_pars_default(dream_pars* p);
void dream_pars_init_vars(dream_pars* p, size_t n);
int  dream_pars_free_vars(dream_pars* p);

size_t dream_par_by_name(const dream_pars* p, string name);

void dream_pars_read_json(dream_pars* p, rapidjson::Value& jpars);
void dream_pars_from_json(dream_pars* p, rapidjson::Value& jpars);
void dream_set_init(dream_pars* p, int n,
                    const double* init, const string* name, const int* lock,
                    const double* lo, const double* hi, const char* scale);

int dream_restore_state(const dream_pars* p, Array3D<double>& state, Array2D<double>& lik, vector<double>& pCR, int& inBurnIn);
void dream_initialize(const dream_pars* p, rng::RngStream* rng, Array2DView<double>& state, ArrayView<double>& lik);
int dream(const dream_pars* p, rng::RngStream* rng);

void check_outliers(int t, Array2D<double>& lik, vector<double>& meanlik,
                    vector<bool>& outliers);
void gen_CR(rng::RngStream* rng, const vector<double>& pCR,
            Array2D<int>& CRm, vector<unsigned>& L);
void gelman_rubin(Array3D<double>& state, vector<double>& scaleReduction,
                  const int* lockVar, int first = 0, int numIter = -1,
                  int adjustDF = 0);

#endif

max_stars_repo_name: gaberoo/cdream

max_issues_repo_name: gaberoo/cdream

max_forks_repo_name: gaberoo/cdream


text:
/* const/test.c
 *
 * Copyright (C) 2003, 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.
 */

#include <config.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_const.h>
#include <gsl/gsl_test.h>

#include <gsl/gsl_ieee_utils.h>

int
main (void)
{
  gsl_ieee_env_setup ();

  /* Basic check to make sure the header files are functioning */

  {
    double c = GSL_CONST_MKS_SPEED_OF_LIGHT;
    double eps = GSL_CONST_MKS_VACUUM_PERMITTIVITY;
    double mu = GSL_CONST_MKS_VACUUM_PERMEABILITY;

    gsl_test_rel (c, 1.0/sqrt(eps*mu), 1e-6, "speed of light (mks)");
  }

  {
    double ly = GSL_CONST_CGS_LIGHT_YEAR;
    double c = GSL_CONST_CGS_SPEED_OF_LIGHT;
    double y = 365.2425 * GSL_CONST_CGS_DAY;

    gsl_test_rel (ly, c * y, 1e-6, "light year (cgs)");
  }

  {
    double c = GSL_CONST_MKSA_SPEED_OF_LIGHT;
    double eps = GSL_CONST_MKSA_VACUUM_PERMITTIVITY;
    double mu = GSL_CONST_MKSA_VACUUM_PERMEABILITY;

    gsl_test_rel (c, 1.0/sqrt(eps*mu), 1e-6, "speed of light (mksa)");
  }

  {
    double ly = GSL_CONST_CGSM_LIGHT_YEAR;
    double c = GSL_CONST_CGSM_SPEED_OF_LIGHT;
    double y = 365.2425 * GSL_CONST_CGSM_DAY;

    gsl_test_rel (ly, c * y, 1e-6, "light year (cgsm)");
  }

  {
    double micro = GSL_CONST_NUM_MICRO;
    double mega = GSL_CONST_NUM_MEGA;
    double kilo = GSL_CONST_NUM_KILO;

    gsl_test_rel (mega/kilo, 1/(micro*kilo), 1e-10, "kilo (mega/kilo, 1/(micro*kilo))");
  }

  {
    double d = GSL_CONST_MKSA_DEBYE;
    double c = GSL_CONST_MKSA_SPEED_OF_LIGHT;
    double desu = d * c * 1000.0;

    gsl_test_rel (desu, 1e-18, 1e-10, "debye (esu)");
  }

  {
    double k = GSL_CONST_MKSA_BOLTZMANN;
    double c = GSL_CONST_MKSA_SPEED_OF_LIGHT;
    double h = GSL_CONST_MKSA_PLANCKS_CONSTANT_H;
    double s = 2 * pow(M_PI, 5.0) * pow(k, 4.0) / (15 * pow(c, 2.0) * pow(h, 3.0));
    double sigma = GSL_CONST_MKSA_STEFAN_BOLTZMANN_CONSTANT;

    gsl_test_rel(s, sigma, 1e-10, "stefan boltzmann constant");
  }

  exit (gsl_test_summary ());
}

max_stars_repo_name: utdsimmons/ohpc

max_issues_repo_name: utdsimmons/ohpc

max_forks_repo_name: utdsimmons/ohpc


text:
/**
 * APER_CHECK - create a FITS image showing the locations of
 *               apertures.
*/

#ifndef _APER_CHECK_H
#define _APER_CHECK_H
#include <gsl/gsl_matrix.h>
#include "aXe_grism.h"

typedef struct
{
  gsl_matrix *img;
}
aXe_mask;

extern aXe_mask *
aXe_mask_init (observation * ob);

extern void
add_ap_p_to_aXe_mask (ap_pixel * ap_p, aXe_mask * mask);

extern void
mark_trace_in_aXe_mask(ap_pixel * ap_p, aXe_mask *mask);
#endif

max_stars_repo_name: sosey/pyaxe

max_issues_repo_name: sosey/pyaxe

max_forks_repo_name: sosey/pyaxe


text:
/**
 * original author: Jochen K"upper
 * author: Pierre Schnizer
 * created: Jan 2002
 * modified: May 2009
 * file: pygsl/src/statistics/longmodule.c
 * $Id: charmodule.c,v 1.9 2009/05/11 08:53:03 schnizer Exp $
 *
 *
 *"
 */

#include <Python.h>
#include <pygsl/error_helpers.h>
#include <pygsl/block_helpers.h>
#include <gsl/gsl_statistics.h>

/* include real functions for different data-types */

#define STATMOD_APPEND_PY_TYPE(X) X ## Int

#if (defined PyGSL_NUMARRAY) || (defined PyGSL_NUMERIC)
#define STATMOD_APPEND_PYC_TYPE(X) X ## CHAR
#else /* PyGSL_NUMARRAY */
#define STATMOD_APPEND_PYC_TYPE(X) X ## BYTE
#endif /* PyGSL_NUMARRAY */

#define STATMOD_FUNC_EXT(X, Y) X ## _char ## Y
#define STATMOD_PY_AS_C PyInt_AsLong
#define STATMOD_C_TYPE char
#include "functions.c"

PyGSL_STATISTICS_INIT(char, "char")

/*
 * Local Variables:
 * mode: c
 * c-file-style: "python"
 * End:
 */

max_stars_repo_name: juhnowski/FishingRod

max_issues_repo_name: juhnowski/FishingRod

max_forks_repo_name: juhnowski/FishingRod


text:
#include "AFS.h"
#include "cs.h"
#include "AFS_ctmc.h"
#include "sparseExp.h"
#include "adkCSparse.h"
#include "adkGSL.h"
//#include "cuba.h"
#include <gsl/gsl_blas.h>

void fillExpectedAFS_unnorm(afsStateSpace *S, double *visitMat ,gsl_vector *rates, gsl_matrix *expAFS);

//fills a sparse transition matrix for Continuous Time model
struct cs_di_sparse *fillTransitionMatrix_ctmc(afsStateSpace *S, double *topol, int *moveType, int *nzCount,
	int *dim1, int *dim2, double theta1, double theta2, double mig1, double mig2){
	int i, j,k;
	int N = S->nstates;
	double  c0r[N],c1r[N],m0r[N],m1r[N],totR[N],rowSums[N], tmp;;

	struct cs_di_sparse *triplet;

	//now allocate new cs_sparse obj
	triplet = cs_spalloc(N, N, *nzCount + N , 1, 1); //alloc sparse mat with extra space for nonzero identity mats
	for(i=0;i<N;i++){
		c0r[i] = S->states[i]->aCounts[0] * (S->states[i]->aCounts[0]-1) / theta1 ;
		c1r[i] = S->states[i]->aCounts[1] * (S->states[i]->aCounts[1]-1) / theta2;
		m0r[i] = S->states[i]->aCounts[0] * mig1 ;
		m1r[i] = S->states[i]->aCounts[1] * mig2 ;
		totR[i] = c0r[i] + c1r[i] + m0r[i] + m1r[i];
		rowSums[i] = 0.0;
	}

	for(k=0;k<*nzCount;k++){
		i = dim1[k];
		j= dim2[k];
		switch(moveType[k]){
			case 0: //coal pop0
			tmp = c0r[i] / totR[i] * topol[k];
			cs_entry(triplet,i,j,tmp);
			rowSums[i] += tmp;
			break;
			case 1: //coal pop1
			tmp = c1r[i] / totR[i] * topol[k];
			cs_entry(triplet,i,j,tmp);
			rowSums[i] += tmp;
			break;
			case 2: //mig pop0
			tmp = m0r[i] / totR[i] * topol[k];
			cs_entry(triplet,i,j,tmp);
			rowSums[i] += tmp;
			break;
			case 3: //mig pop1
			tmp = m1r[i] / totR[i] * topol[k];
			cs_entry(triplet,i,j,tmp);
			rowSums[i] += tmp;
			break;
		}
	}
	//now add diagonal elements; topol pre-set to -1
	for(k=0;k<N;k++){
		if((rowSums[k]) != 0.0){
	 		cs_entry(triplet,k,k,0-rowSums[k]);
			*nzCount += 1;
		}
	}
//	cs_dropzeros(triplet);
//	tmat = cs_compress(triplet);
//	cs_spfree(triplet);
	return(triplet);

}

//this only fills a new set of arrays
void fillTransitionMatrixArray_ctmc(afsStateSpace *S, double *topol, int *moveType, int *nzCount,
	int *dim1, int *dim2, double theta1, double theta2, double mig1, double mig2,double *newArray,
	int *newDim1, int *newDim2){
	int i, j,k, newnz;
	int N = S->nstates;
	double  c0r[N],c1r[N],m0r[N],m1r[N],totR[N],rowSums[N], tmp;;

	for(i=0;i<N;i++){
		c0r[i] = S->states[i]->aCounts[0] * (S->states[i]->aCounts[0]-1) / theta1 ;
		c1r[i] = S->states[i]->aCounts[1] * (S->states[i]->aCounts[1]-1) / theta2;
		m0r[i] = S->states[i]->aCounts[0] * mig1 ;
		m1r[i] = S->states[i]->aCounts[1] * mig2 ;
		totR[i] = c0r[i] + c1r[i] + m0r[i] + m1r[i];
		rowSums[i] = 0.0;
	}
	newnz=0;
	for(k=0;k<*nzCount;++k){
		i = dim1[k];
		j= dim2[k];
		switch(moveType[k]){
			case 0: //coal pop0
			tmp = c0r[i] / totR[i] * topol[k];
			newArray[k]=tmp;
			newDim1[k] = i+1;
			newDim2[k] = j+1;
			rowSums[i] += tmp;
			newnz++;
			break;
			case 1: //coal pop1
			tmp = c1r[i] / totR[i] * topol[k];
			newArray[k]=tmp;
			newDim1[k] = i+1;
			newDim2[k] = j+1;
			rowSums[i] += tmp;
			newnz++;
			break;
			case 2: //mig pop0
			tmp = m0r[i] / totR[i] * topol[k];
			newArray[k]=tmp;
			newDim1[k] = i+1;
			newDim2[k] = j+1;
			rowSums[i] += tmp;
			newnz++;
			break;
			case 3: //mig pop1
			tmp = m1r[i] / totR[i] * topol[k];
			newArray[k]=tmp;
			newDim1[k] = i+1;
			newDim2[k] = j+1;
			rowSums[i] += tmp;
			newnz++;
			break;
		}
	}

	//now add diagonal elements; topol pre-set to -1
	for(k=0;k<N;k++){
		if((rowSums[k]) != 0.0){
			newArray[newnz]=0-rowSums[k];
			newDim1[newnz] = k+1;
			newDim2[newnz] = k+1;
			newnz++;
		}
	}

	//sanity check
//	for(k=0;k<newnz;k++){
//		printf("row: %d col: %d val: %g\n",newDim1[k],newDim2[k],newArray[k]);
//	}
	*nzCount = newnz;
}

//Functions for jointly computing CTMC transitions and embedded chain

//this function returns pointer to embedded Discrete MC trans Mat in sparse form, while filling rates vector and filling
// arrays for CTMC representation
struct cs_di_sparse * fillTransitionMatrixArray_embed_ctmc(afsStateSpace *S, double *topol, int *moveType, int *nzCount,
	int *dim1, int *dim2, double theta1, double theta2, double mig1, double mig2,double *newArray,
	int *newDim1, int *newDim2, gsl_vector *rates){
	int i, j,k, newnz, abcount=0;
	int N = S->nstates,abFlag;
	double  c0r[N],c1r[N],m0r[N],m1r[N],totR[N],rowSums[N], tmp;;
	struct cs_di_sparse *triplet, *tmat;

	gsl_vector_set_zero(rates);
	//allocate new cs_sparse obj
	triplet = cs_spalloc(N, N, *nzCount + N , 1, 1); //alloc sparse mat with extra space for nonzero identity mats

	for(i=0;i<N;i++){
		c0r[i] = S->states[i]->aCounts[0] * (S->states[i]->aCounts[0]-1)  / theta1 ;
	//	printf("co check:%f\n",c0r[i]);
		if(theta2 == 0)
			c1r[i] = 0.0;
		else
			c1r[i] = S->states[i]->aCounts[1] * (S->states[i]->aCounts[1]-1)  / theta2;
		m0r[i] = S->states[i]->aCounts[0] * mig1 ;
		m1r[i] = S->states[i]->aCounts[1] * mig2 ;
		totR[i] = c0r[i] + c1r[i] + m0r[i] + m1r[i];

		rowSums[i] = 0.0;
		if(S->states[i]->nalleles > 1  && c0r[i] >= 0 && c1r[i] >= 0 && totR[i] != 0.0)
			gsl_vector_set(rates,i,1.0/totR[i]);
		else
			gsl_vector_set(rates,i,0);
	}

	newnz=0;
	for(k=0;k<*nzCount;k++){
		i = dim1[k];
		j= dim2[k];
//		printf("i: %d, j:%d, movetype: %d\n",i,j,moveType[k]);
		if(S->states[i]->nalleles ==1 || S->states[j]->nalleles ==1 )abFlag=1;
	//	if(S->states[i]->nalleles ==1 )abFlag=1;
		else abFlag = 0;
		switch(moveType[k]){
			case 0: //coal pop0
			if(totR[i]==0)tmp=0;
			else tmp = c0r[i] / totR[i] * topol[k];
			newArray[k]=tmp* totR[i];
			newDim1[k] = i;
			newDim2[k] = j;
			rowSums[i] += tmp* totR[i];
			if(abFlag==0)cs_entry(triplet,i,j,tmp);
			break;
			case 1: //coal pop1
			if(totR[i]==0)tmp=0;
			else tmp = c1r[i] / totR[i] * topol[k];
			newArray[k]=tmp* totR[i];
			newDim1[k] = i;
			newDim2[k] = j;
			rowSums[i] += tmp* totR[i];
			if(abFlag==0)cs_entry(triplet,i,j,tmp);
			break;
			case 2: //mig pop0
			if(totR[i]==0)tmp=0;
			else tmp = m0r[i] / totR[i] * topol[k];
			newArray[k]=tmp* totR[i];
			newDim1[k] = i;
			newDim2[k] = j;
			rowSums[i] += tmp* totR[i];
			if(abFlag==0)cs_entry(triplet,i,j,tmp);
			break;
			case 3: //mig pop1
			if(totR[i]==0)tmp=0;
			else tmp = m1r[i] / totR[i] * topol[k];
			newArray[k]=tmp * totR[i];
			newDim1[k] = i;
			newDim2[k] = j;
			rowSums[i] += tmp* totR[i];
			if(abFlag==0)cs_entry(triplet,i,j,tmp);
			break;
			case -1:
			abcount+=1;
			break;

		}
	}
	newnz=*nzCount ;
	//now add diagonal elements; topol pre-set to -1
	for(k=0;k<N;k++){
		if((rowSums[k]) != 0.0){
			newArray[newnz]=0.0-rowSums[k];
			newDim1[newnz] = k;
			newDim2[newnz] = k;
			newnz++;
		}
	}

	tmat = cs_compress(triplet);
	cs_spfree(triplet);
	*nzCount = newnz;
	return(tmat);
}

//fills up an AFS matrix using expected times and a precalculated visitMat
void fillExpectedAFS_ctmc(afsStateSpace *S, double *visitMat, gsl_matrix *expAFS){
	int  j, k, l, m1, m2, count;

	double tmpRes;

	m1=expAFS->size1;
	m2=expAFS->size2;

	gsl_matrix_set_zero(expAFS);

	for(j=0;j<S->nstates;j++){
		if(S->states[j]->nalleles != 1){
			for(k=0;k<m1;k++){
			for(l=0;l<m2;l++){
				count = gsl_matrix_int_get(S->states[j]->popMats[0],k,l) +  gsl_matrix_int_get(S->states[j]->popMats[1],k,l);
				tmpRes = gsl_matrix_get(expAFS,k,l) + (count * visitMat[j]);
				gsl_matrix_set(expAFS,k,l, tmpRes);
				}
			}
		}
	}
//	tmpRes = matrixSumDouble(expAFS);
//	gsl_matrix_scale(expAFS, 1.0 / tmpRes);
}

//fills up an AFS matrix using expected times and a precalculated visitMat
void fillExpectedAFS_unnorm(afsStateSpace *S, double *visitMat ,gsl_vector *rates, gsl_matrix *expAFS){
	int  j, k, l, m1, m2, count;

	double tmpRes;

	m1=expAFS->size1;
	m2=expAFS->size2;

	gsl_matrix_set_zero(expAFS);
	for(j=0;j<S->nstates;j++){
		for(k=0;k<m1;k++){
			for(l=0;l<m2;l++){
				count = gsl_matrix_int_get(S->states[j]->popMats[0],k,l) +  gsl_matrix_int_get(S->states[j]->popMats[1],k,l);
				tmpRes = gsl_matrix_get(expAFS,k,l) + (count * gsl_vector_get(rates,j) * visitMat[j]);
//				printf("tmpRes[%d][%d]: %f rates[j]: %f visits[j]: %f count: %d\n",tmpRes,k,l,gsl_vector_get(rates,j),visitMat[j],count);
				gsl_matrix_set(expAFS,k,l, tmpRes);
			}
		}
	}
}

//fills up an AFS matrix using expected times and a precalculated visitMat
void fillExpectedAFS_unnorm_mat(afsStateSpace *S, double **visitMat , double *initialStateProbs, gsl_vector *rates, gsl_matrix *expAFS){
	int  j, k, l,i, m1, m2, count;

	double tmpRes,sum;

	m1=expAFS->size1;
	m2=expAFS->size2;

	gsl_matrix_set_zero(expAFS);
	for(j=0;j<S->nstates;j++){
		for(k=0;k<m1;k++){
			for(l=0;l<m2;l++){
				count = gsl_matrix_int_get(S->states[j]->popMats[0],k,l) +  gsl_matrix_int_get(S->states[j]->popMats[1],k,l);
				sum=0.0;
				for(i=0;i<S->nstates;i++){
					sum += visitMat[i][j] * initialStateProbs[i];
				}
				tmpRes	= gsl_matrix_get(expAFS,k,l) + (count * gsl_vector_get(rates,j) * sum);
			//	printf("tmpRes[%d][%d]: %f rates[j]: %f visits[j]: %f count: %d\n",k,l,tmpRes,gsl_vector_get(rates,j),sum,count);
				gsl_matrix_set(expAFS,k,l, tmpRes);
			}
		}
	}
}

/// I believe the below is now completely unused
//fillLogAFS -- uses cxsparse library for inversion -- will replace above
// void fillLogAFS(void * p){
// 	struct im_lik_params * params = (struct im_lik_params *) p;
// 	struct expoMatObj anExpoMatObj;
// 	int i,j, N = params->stateSpace->nstates;
//  	cs *spMat, *mt, *ident, *eye, *tmpMat ,*tmpMat2;
// 	double timeV, sum, thetaA, theta2, mig1, mig2;
// 	int Na;
// 	gsl_vector *tmpStates;
//
//
// 	//initialize some vectors
// 	gsl_vector_set_zero(params->rates);
// 	tmpStates = gsl_vector_alloc(N);
//
// 	//For straight MLE the paramVector takes the form [N2,NA,m1,m2,t]
// 	theta2 = gsl_vector_get(params->paramVector,0);
// 	thetaA = gsl_vector_get(params->paramVector,1);
// 	mig1 = gsl_vector_get(params->paramVector,2);
// 	mig2 = gsl_vector_get(params->paramVector,3);
// 	timeV = gsl_vector_get(params->paramVector,4);
// //	printf("params-> %f %f %f %f %f\n",theta2,thetaA,mig1,mig2,timeV);
// 	//fill transMat
// 	tmpMat= fillTransitionMatrixArray_embed_ctmc(params->stateSpace, params->topA, params->moveA,
// 	 	&(params->nnz),params->dim1, params->dim2,1.0, theta2, mig1, mig2,
// 	 	params->new, params->newdim1, params->newdim2,params->rates);
//
// 	//using CSparse
// 	// S = (I-P)^-1
//
// 	//add negative ident
// 	ident = cs_spalloc(N,N,N,1,1);
// 	for(i=0;i<N;i++) cs_entry(ident,i,i,1);
// 	eye = cs_compress(ident);
// 	spMat = cs_add(eye,tmpMat,1.0,-1.0);
// 	mt = cs_transpose(spMat,1);
//
// 	////Compute Entire Inverse Mat
// 	for(j=0;j<N;j++){
// 		//create unit array for solve
// 		for(i=0; i<N; i++)params->b[i] = 0.0;
// 		params->b[j]=1.0;
// 		cs_lusol(0, mt, params->b, 1e-12);
// 		for(i=0; i<N; i++){
// 			params->invMat[j][i]=params->b[i];
// 			gsl_matrix_set(params->invMatGSL,j,i,params->b[i]);
// 		}
// 		printf("inv mat row %d complete\n",j);
// 	}
//
// 	//Get Island Time 0-INF Unnormal
// 	gsl_matrix_set_zero(params->expAFS);
// 	for(i=0; i<N; i++)params->b[i] = params->invMat[0][i];
// 	fillExpectedAFS_unnorm(params->stateSpace, params->b,params->rates,params->expAFS);
// 	// printf("////////////////////Island Time 0 - INF unnormalized\n");
// 	// for(i=0;i< (params->expAFS->size1) ;i++){
// 	// 	for(j=0;j< (params->expAFS->size2);j++){
// 	// 		printf("%.5f ",gsl_matrix_get(params->expAFS,i,j));
// 	// 	}
// 	// 	printf("\n");
// 	// }
//
// 	//Matrix Exponentiation to get state vector at time t
// 	//
//
// 	//push indices
// 	for(i=0;i<params->nnz;i++){
// 		params->expodim1[i]=params->newdim1[i]+1;
// 		params->expodim2[i]=params->newdim2[i]+1;
//
// 	}
// 	//reset expoInit
// 	for(i = 0;i<N;i++){
// 		params->expoInit[i]=0.0;
// 		params->expoResult[i]=0.0;
// 	}
// 	params->expoInit[0] = 1.0;
//
// 	//sanity check ///////
// //	for(i=0;i<params->nnz;i++){
// //		printf("a[%d]=%f expodim1=%d expodim2=%d \n", i,params->anExpoMatObj->a[i],params->anExpoMatObj->ai[i],params->anExpoMatObj->aj[i]);
// //	}
// //	printf("rank=%d nnz=%d\n", params->anExpoMatObj->rank, params->anExpoMatObj->nnz);
//
// //	for(i=0;i<N;i++)printf("params->anExpoMatObj->resultSpace[%d]=%f\n",i,params->anExpoMatObj->resultSpace[i]);
// //	for(i=0;i<N;i++)printf("params->anExpoMatObj->initVec[%d]=%f\n",i,params->anExpoMatObj->initVec[i]);
//
// 	anExpoMatObj.ai = params->expodim1;
// 	anExpoMatObj.aj = params->expodim2;
// 	anExpoMatObj.a = params->new;
// 	anExpoMatObj.initVec = params->expoInit;
// 	anExpoMatObj.resultSpace = params->expoResult;
// //	sparse_exponential_single_row(N, params->nnz, anExpoMatObj.ai, anExpoMatObj.aj,
// //		anExpoMatObj.a, timeV,anExpoMatObj.initVec, params->expoResult, 0);
// 	//state vector at time t stored in b[i]
// 	for(i=0; i<N; i++){
// 		params->b[i] = params->expoResult[i];
// //		printf("params->expoResult[%d]=%f\n",i,params->expoResult[i] );
// 		gsl_vector_set(tmpStates,i,params->expoResult[i]);
// 	}
//
// 	gsl_blas_dgemv(CblasTrans, 1.0, params->invMatGSL,tmpStates, 0.0, params->resVec);
//
//
// 	for(i=0; i<N; i++){
// 		params->b[i] = gsl_vector_get(params->resVec,i);
// 	}
// 	gsl_matrix_set_zero(params->expAFS2);
// 	fillExpectedAFS_unnorm(params->stateSpace, params->b,params->rates,params->expAFS2);
//
// 	//now subtract older AFS (t_t - t_INF) from total AFS (t_0 - t_INF)
// 	/////////////////
// 	///////
// //	printf("////////////////////Island time0-t=%f unnormalized\n",timeV);
//
// 	gsl_matrix_sub(params->expAFS,params->expAFS2);
// 	// for(i=0;i< (params->expAFS->size1) ;i++){
// 	// 	for(j=0;j< (params->expAFS->size2);j++){
// 	// 		printf("%.5f ",gsl_matrix_get(params->expAFS,i,j));
// 	// 	}
// 	// 	printf("\n");
// 	// }
// 	// printf("////////////////////\n");
//
//
// 	//mapping of states already complete when params initialized
// 	//state vector at time t stored in st[] for largerStateSpace; use reverseMap
// 	for(i=0;i<N;i++)params->st[i]=0.0;
// 	for(i=0;i<N;i++)params->st[params->map[i]]+=params->expoResult[i];
// 	for(i=0; i<params->reducedStateSpace->nstates; i++){
// 		gsl_vector_set(params->ancStateVec,i,params->st[params->reverseMap[i]]);
// 	}
// 	Na = params->reducedStateSpace->nstates;
// 	//fill up ancestral transition matrix
// 	gsl_vector_set_zero(params->rates);
// 	tmpMat2= fillTransitionMatrixArray_embed_ctmc(params->reducedStateSpace, params->topA2, params->moveA2, &(params->nnzA),
// 		params->dim1A, params->dim2A,thetaA, 0, 0, 0,params->new, params->newdim1, params->newdim2,params->rates);
//
// 	ident = cs_spalloc(Na,Na,Na,1,1);
// 	for(i=0;i<Na;i++) cs_entry(ident,i,i,1);
// 	eye = cs_compress(ident);
// 	spMat = cs_add(eye,tmpMat2,1.0,-1.0);
// 	mt = cs_transpose(spMat,1);
//
// //	gsl_matrix_free(invMatGSL);
// //	invMatGSLA = gsl_matrix_calloc(Na,Na);
//
// 	for(j=0;j<Na;j++){
// 		//create unit array for solve
// 		for(i=0; i<Na; i++)params->b[i] = 0.0;
// 		params->b[j]=1.0;
// 		cs_lusol(0, mt, params->b, 1e-12);
// 		for(i=0; i<Na; i++){
// 			gsl_matrix_set(params->invMatGSLA,j,i,params->b[i]);
// 		}
// 	}
//
// 	//get contribution starting at st
// 	gsl_blas_dgemv(CblasTrans, 1.0, params->invMatGSLA, params->ancStateVec, 0.0, params->ancResVec);
// 	for(i=0; i<Na; i++)params->b[i] = gsl_vector_get(params->ancResVec,i);
// 	gsl_matrix_set_zero(params->expAFS2);
// 	fillExpectedAFS_unnorm(params->reducedStateSpace, params->b,params->rates,params->expAFS2);
//
// 	////////////////////normalized IM AFS
// 	gsl_matrix_add(params->expAFS,params->expAFS2);
// 	sum = matrixSumDouble(params->expAFS);
// 	gsl_matrix_scale(params->expAFS, 1.0 / sum);
//
//
// 	///////////// Replace AFS with log(AFS)
// 	// for(i=0;i<params->expAFS->size1;i++){
// 	// 	for(j=0;j<params->expAFS->size2;j++){
// 	// 		gsl_matrix_set(params->expAFS,i,j,log(gsl_matrix_get(params->expAFS,i,j)));
// 	// 	}
// 	// }
//
// 	//clean up
// 	cs_spfree(spMat);
// 	cs_spfree(mt);
// 	cs_spfree(ident);
// 	cs_spfree(eye);
// 	cs_spfree(tmpMat);
// 	cs_spfree(tmpMat2);
// 	gsl_vector_free(tmpStates);
// }

//////////////////// Uniformization Algorithm for Exponentiation
/// P. 413 from Stewart Book
void uniformizExpSparseMat(int N, int nnz, double *vals, int *dim1, int *dim2, double tSoln, double *result){
	int i, k, ii;
	double cosi, sigma, nu, gamma, eps,y[N],tmpY[N],scale;
	struct cs_di_sparse *tmpMat, *ctmc, *ident, *eye, *new;

	//stuff values into more convenient matrix
	//uniformization of matrix
	tmpMat = cs_spalloc(N,N,N,1,1);
	ident = cs_spalloc(N,N,N,1,1);
	for(i=0;i<nnz;i++)cs_entry(tmpMat,dim1[i],dim2[i],vals[i]);
	ctmc=cs_compress(tmpMat);
	gamma = cs_max(ctmc);
	printf("gamma=%f\n------------------\n",gamma);
	cs_print_adk(ctmc);
	cs_scale(ctmc,1.0/gamma);
	cs_print_adk(ctmc);

	for(i=0;i<N;i++) cs_entry(ident,i,i,1);
	eye = cs_compress(ident);
	new = cs_add(eye,ctmc,1.0,1.0);
	cs_print_adk(new);

	//tmpMat now holds the subjugated dtmc

	//step one of uniformization is to compute K, number of terms in the summation
	eps=0.001;
	k=0; cosi=1.0; sigma=1.0; nu = (1-eps)/exp(-1.0*gamma*tSoln);
	while(sigma < nu){
		k+=1;
		cosi *= (gamma * tSoln)/k;
		sigma += cosi;
	}
	printf("here is k: %d\n",k);

	//initialize result;
	for(i=1;i<N;i++){
		result[i]=0.0;
		y[i]=0.0;
		tmpY[i]=0.0;
	}
	result[0]=1.0;y[0]=1.0;tmpY[0]=0.0;
	printf("initialization done\n");
	for(ii=0;ii<N;ii++) printf("tmpY[%d]=%f\n",ii,tmpY[ii]);
	printf("---------------------\n");
	for(ii=0;ii<N;ii++)printf("y[%d]=%f\n",ii,y[ii]);
	printf("---------------------\n");

	//now approximate pi(t)
	int testVal;
	for(i=1;i<=k;i++){
		printf("k loop k=%d\n",i);
		testVal=cs_gaxpy(new, y, tmpY); //tmpY now contains the answer
		printf("test of gaxpy: %d\n",testVal);
		for(ii=0;ii<N;ii++) printf("tmpY[%d]=%f\n",ii,tmpY[ii]);
		printf("---------------------\n");
		for(ii=0;ii<N;ii++)printf("y[%d]=%f\n",ii,y[ii]);
		printf("---------------------\n");
		cs_print_adk(tmpMat);
		printf("---------------------\n");

		for(ii=0;ii<N;ii++){
	//		printf("tmpY[%d]=%f\n",ii,tmpY[ii]);
	//		printf("y[%d]=%f\n",ii,y[ii]);
			y[ii]=tmpY[ii] * ((gamma*tSoln)/i);
			result[ii]+=y[ii];
			tmpY[ii]=0;
		}

	}
	//final scaling
	scale = exp(-1.0*gamma*tSoln);
	for(i=0;i<N;i++){
		result[i] *= scale;
		printf("result[%d]=%f\n",i,result[i]);
	}
}

max_stars_repo_name: dortegadelv/IMaDNA

max_issues_repo_name: dortegadelv/IMaDNA

max_forks_repo_name: dortegadelv/IMaDNA


text:
/*****************************************************************************
*
* Rokko: Integrated Interface for libraries of eigenvalue decomposition
*
* Copyright (C) 2012-2017 by Synge Todo <wistaria@phys.s.u-tokyo.ac.jp>
*
* Distributed under the Boost Software License, Version 1.0. (See accompanying
* file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
*
*****************************************************************************/

#include <math.h>
#include <cblas.h>
#include <lapacke.h>
#include <rokko/cmatrix.h>

int imax(int x, int y) { return (x > y) ? x : y; }

int main(int argc, char** argv) {
  int n = 5;
  int i, j;
  double norm, check;
  double **a;

  if (argc > 1) n = atoi(argv[1]);

  // generate matrix and rhs vector
  a = alloc_dmatrix(n, n);
  for (j = 0; j < n; ++j) {
    for (i = 0; i < n; ++i) {
      mat_elem(a, i, j) = n - 0.253 * imax(i, j);
    }
  }
  printf("Matrix A: ");
  fprint_dmatrix(stdout, n, n, a);

  /* calculate various norms */
  norm = LAPACKE_dlange(LAPACK_COL_MAJOR, '1', n, n, mat_ptr(a), n);
  printf("norm1(A) = %e\n", norm);

  norm = LAPACKE_dlange(LAPACK_COL_MAJOR, 'I', n, n, mat_ptr(a), n);
  printf("normI(A) = %e\n", norm);

  norm = LAPACKE_dlange(LAPACK_COL_MAJOR, 'F', n, n, mat_ptr(a), n);
  printf("normF(A) = %e\n", norm);

  /* check of result */
  check = 0;
  for (j = 0; j < n; ++j) {
    for (i = 0; i < n; ++i) {
      check += mat_elem(a, i, j) * mat_elem(a, i, j);
    }
  }
  if (fabs(check - norm * norm) > 1e-12 * check) {
    fprintf(stderr, "Error: check error %e != %e\n", norm * norm, check);
    exit(255);
  }

  free_dmatrix(a);
  return 0;
}

max_stars_repo_name: t-sakashita/rokko

max_issues_repo_name: t-sakashita/rokko

max_forks_repo_name: t-sakashita/rokko


text:
/*******************************************************************************  *
*
*       This file is part of the General Hidden Markov Model Library,
*       GHMM version __VERSION__, see http://ghmm.org
*
*       Filename: ghmm/ghmm/root_finder.c
*       Authors:  Achim Gaedke
*
*       Copyright (C) 1998-2004 Alexander Schliep
*       Copyright (C) 1998-2001 ZAIK/ZPR, Universitaet zu Koeln
*       Copyright (C) 2002-2004 Max-Planck-Institut fuer Molekulare Genetik,
*                               Berlin
*
*       Contact: schliep@ghmm.org
*
*       This library is free software; you can redistribute it and/or
*       modify it under the terms of the GNU Library General Public
*       License as published by the Free Software Foundation; either
*       version 2 of the License, or (at your option) any later version.
*
*       This library 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
*       Library General Public License for more details.
*
*       You should have received a copy of the GNU Library General Public
*       License along with this library; if not, write to the Free
*       Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*
*       This file is version $Revision: 2267 $
*                       from $Date: 2009-04-24 11:01:58 -0400 (Fri, 24 Apr 2009) $
*             last change by $Author: grunau $.
*
*******************************************************************************/

#ifdef HAVE_CONFIG_H
#  include "../config.h"
#endif  /*  */

#ifndef DO_WITH_GSL

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

#include "ghmm_internals.h"

double ghmm_zbrent_AB (double (*func) (double, double, double, double),
                   double x1, double x2, double tol, double A, double B,
                   double eps)
{
  fprintf (stderr, "Function ghmm_zbrent_AB() not implemented!\n");
  exit (1);

    /*
       double a, b, c;
       double fa, fb, fc;

       a = min(x1, x2);
       fa = (*func)(a);
       c = max(x1, x2);
       fc = (*func)(c);
       b = (c - a)/2.0;
       fb = (*func)(b);

       while (fabs(c - a) > tol + (tol * min(fabs(a), fabs(c))))
       {
       r = fb/fc;
       s = fb/fa;
       t = fa/fc;
       p = s * (t * (r - t) * (c - b) - (1.0 - r) * (b - a));
       q = (t - 1.0) * (r - 1.0) * (s - 1.0);

       x = b + p/q;

       if (x > a && x < c)
       {
       / * Accept interpolating point * /
       if (x < b)
       {
       c = b;
       fc = fb;
       b = x;
       fb = (*func)(b);
       }
       else if (x > b)
       {
       a = b;
       fa = fb;
       b = x;
       fb = (*func)(b);
       }
       }
       else
       {
       / * Use bisection * /
       }
       }
     */
}
#else   /*  */

#include <gsl/gsl_errno.h>
#include <gsl/gsl_roots.h>

/* struct for function pointer and parameters except 1st one */
  struct parameter_wrapper {
  double (*func) (double, double, double, double);
   double x2;
    double x3;
    double x4;
 } parameter;

/* calls the given function during gsl root solving iteration
   the first parameter variates, all other are kept constant
*/
double function_wrapper (double x, void *p)
{
  struct parameter_wrapper *param = (struct parameter_wrapper *) p;
  return param->func (x, param->x2, param->x3, param->x4);
}

/*
  this interface is used in sreestimate.c
 */
double ghmm_zbrent_AB (double (*func) (double, double, double, double), double x1,
                  double x2, double tol, double A, double B, double eps)
{

    /* initialisation of wrapper structure */
  struct parameter_wrapper param;
  gsl_function f;

#ifdef HAVE_GSL_INTERVAL /* gsl_interval vanished with version 0.9 */
    gsl_interval x;

#endif  /*  */
    gsl_root_fsolver * s;
  double tolerance;
  int success = 0;
  double result = 0;
  param.func = func;
  param.x2 = A;
  param.x3 = B;
  param.x4 = eps;
  f.function = &function_wrapper;
  f.params = (void *) &param;
  tolerance = tol;

    /* initialisation */
#ifdef HAVE_GSL_INTERVAL
# ifdef GSL_ROOT_FSLOVER_ALLOC_WITH_ONE_ARG
    x.lower = x1;
  x.upper = x2;
  s = gsl_root_fsolver_alloc (gsl_root_fsolver_brent);
  gsl_root_fsolver_set (s, &f, x);

# else
    s = gsl_root_fsolver_alloc (gsl_root_fsolver_brent, &f, x);

# endif
#else   /* gsl_interval vanished with version 0.9 */
    s = gsl_root_fsolver_alloc (gsl_root_fsolver_brent);
  gsl_root_fsolver_set (s, &f, x1, x2);

#endif  /*  */

    /* iteration */
    do {
    success = gsl_root_fsolver_iterate (s);
    if (success == GSL_SUCCESS)
       {

#ifdef HAVE_GSL_INTERVAL
        gsl_interval new_x;
      new_x = gsl_root_fsolver_interval (s);
      success = gsl_root_test_interval (new_x, tolerance, tolerance);

#else   /* gsl_interval vanished with version 0.9 */
      double x_up;
      double x_low;
      (void) gsl_root_fsolver_iterate (s);
      x_up = gsl_root_fsolver_x_upper (s);
      x_low = gsl_root_fsolver_x_lower (s);
      success = gsl_root_test_interval (x_low, x_up, tolerance, tolerance);

#endif  /*  */
      }
  } while (success == GSL_CONTINUE);

    /* result */
    if (success != GSL_SUCCESS)
     {
    gsl_error ("solver failed", __FILE__, __LINE__, success);
    }

  else
     {
    result = gsl_root_fsolver_root (s);
    }

    /* destruction */
    gsl_root_fsolver_free (s);
  return result;
}

#endif  /*  */

max_stars_repo_name: ruslankuzmin/julia

max_issues_repo_name: ruslankuzmin/julia

max_forks_repo_name: ruslankuzmin/julia


text:
//   Copyright (c) 2018 PaddlePaddle Authors. All Rights Reserved.
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
//     http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.

#pragma once

#include "paddle/fluid/framework/operator.h"
#include "paddle/fluid/framework/tensor.h"

#ifdef PADDLE_WITH_MKLML
#include <mkl_cblas.h>
#include <mkl_lapacke.h>
#include <mkl_vml_functions.h>
#endif

#ifdef PADDLE_USE_OPENBLAS
#include <cblas.h>
#include <lapacke.h>
#endif

#ifndef LAPACK_FOUND
extern "C" {
#include <cblas.h>  // NOLINT
int LAPACKE_sgetrf(int matrix_layout, int m, int n, float* a, int lda,
                   int* ipiv);
int LAPACKE_dgetrf(int matrix_layout, int m, int n, double* a, int lda,
                   int* ipiv);
int LAPACKE_sgetri(int matrix_layout, int n, float* a, int lda,
                   const int* ipiv);
int LAPACKE_dgetri(int matrix_layout, int n, double* a, int lda,
                   const int* ipiv);
}
#endif

namespace paddle {
namespace operators {
namespace math {

template <typename DeviceContext>
class Blas {
 public:
  explicit Blas(const DeviceContext& context) : context_(context) {}

  template <typename T>
  void GEMM(CBLAS_TRANSPOSE transA, CBLAS_TRANSPOSE transB, int M, int N, int K,
            T alpha, const T* A, const T* B, T beta, T* C) const;

  template <typename T>
  void GEMM(bool transA, bool transB, int M, int N, int K, T alpha, const T* A,
            int lda, const T* B, int ldb, T beta, T* C, int ldc) const;

  template <typename T>
  void MatMul(const framework::Tensor& mat_a, bool trans_a,
              const framework::Tensor& mat_b, bool trans_b, T alpha,
              framework::Tensor* mat_out, T beta) const;

  template <typename T>
  void MatMul(const framework::Tensor& mat_a, bool trans_a,
              const framework::Tensor& mat_b, bool trans_b,
              framework::Tensor* mat_out) const {
    MatMul(mat_a, trans_a, mat_b, trans_b, static_cast<T>(1.0), mat_out,
           static_cast<T>(0.0));
  }

  template <typename T>
  void MatMul(const framework::Tensor& mat_a, const framework::Tensor& mat_b,
              framework::Tensor* mat_out) const {
    this->template MatMul<T>(mat_a, false, mat_b, false, mat_out);
  }

  template <typename T>
  void AXPY(int n, T alpha, const T* x, T* y) const;

  template <typename T>
  void GEMV(bool trans_a, int M, int N, T alpha, const T* A, const T* B, T beta,
            T* C) const;

  template <typename T>
  void BatchedGEMM(CBLAS_TRANSPOSE transA, CBLAS_TRANSPOSE transB, int M, int N,
                   int K, T alpha, const T* A, const T* B, T beta, T* C,
                   int batchCount, int64_t strideA, int64_t strideB) const;

 private:
  const DeviceContext& context_;
};

template <typename DeviceContext, typename T>
class BlasT : private Blas<DeviceContext> {
 public:
  using Blas<DeviceContext>::Blas;

  template <typename... ARGS>
  void GEMM(ARGS... args) const {
    Base()->template GEMM<T>(args...);
  }

  template <typename... ARGS>
  void MatMul(ARGS... args) const {
    Base()->template MatMul<T>(args...);
  }

  template <typename... ARGS>
  void AXPY(ARGS... args) const {
    Base()->template AXPY<T>(args...);
  }

  template <typename... ARGS>
  void GEMV(ARGS... args) const {
    Base()->template GEMV<T>(args...);
  }

  template <typename... ARGS>
  void BatchedGEMM(ARGS... args) const {
    Base()->template BatchedGEMM<T>(args...);
  }

 private:
  const Blas<DeviceContext>* Base() const {
    return static_cast<const Blas<DeviceContext>*>(this);
  }
};

template <typename DeviceContext, typename T>
inline BlasT<DeviceContext, T> GetBlas(
    const framework::ExecutionContext& exe_ctx) {
  return BlasT<DeviceContext, T>(
      exe_ctx.template device_context<DeviceContext>());
}

template <typename DeviceContext, typename T>
inline BlasT<DeviceContext, T> GetBlas(const DeviceContext& dev_ctx) {
  return BlasT<DeviceContext, T>(dev_ctx);
}

}  // namespace math
}  // namespace operators
}  // namespace paddle

#include "paddle/fluid/operators/math/blas_impl.h"
#ifdef PADDLE_WITH_CUDA
#include "paddle/fluid/operators/math/blas_impl.cu.h"
#endif

max_stars_repo_name: abhinavarora/Paddle

max_issues_repo_name: abhinavarora/Paddle

max_forks_repo_name: abhinavarora/Paddle


text:
#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_eigen.h>

int
main (void)
{
  double data[] = { -1.0, 1.0, -1.0, 1.0,
                    -8.0, 4.0, -2.0, 1.0,
                    27.0, 9.0, 3.0, 1.0,
                    64.0, 16.0, 4.0, 1.0 };

  gsl_matrix_view m
    = gsl_matrix_view_array (data, 4, 4);

  gsl_vector_complex *eval = gsl_vector_complex_alloc (4);
  gsl_matrix_complex *evec = gsl_matrix_complex_alloc (4, 4);

  gsl_eigen_nonsymmv_workspace * w =
    gsl_eigen_nonsymmv_alloc (4);

  gsl_eigen_nonsymmv (&m.matrix, eval, evec, w);

  gsl_eigen_nonsymmv_free (w);

  gsl_eigen_nonsymmv_sort (eval, evec,
                           GSL_EIGEN_SORT_ABS_DESC);

  {
    int i, j;

    for (i = 0; i < 4; i++)
      {
        gsl_complex eval_i
           = gsl_vector_complex_get (eval, i);
        gsl_vector_complex_view evec_i
           = gsl_matrix_complex_column (evec, i);

        printf ("eigenvalue = %g + %gi\n",
                GSL_REAL(eval_i), GSL_IMAG(eval_i));
        printf ("eigenvector = \n");
        for (j = 0; j < 4; ++j)
          {
            gsl_complex z =
              gsl_vector_complex_get(&evec_i.vector, j);
            printf("%g + %gi\n", GSL_REAL(z), GSL_IMAG(z));
          }
      }
  }

  gsl_vector_complex_free(eval);
  gsl_matrix_complex_free(evec);

  return 0;
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#if !defined(PETIGABL_H)
#define PETIGABL_H

#include <petsc.h>
#include <petscblaslapack.h>
#if defined(PETSC_BLASLAPACK_UNDERSCORE)
   #define sgetri_ sgetri_
   #define dgetri_ dgetri_
   #define qgetri_ qgetri_
   #define cgetri_ cgetri_
   #define zgetri_ zgetri_
#elif defined(PETSC_BLASLAPACK_CAPS)
   #define sgetri_ SGETRI
   #define dgetri_ DGETRI
   #define qgetri_ QGETRI
   #define cgetri_ CGETRI
   #define zgetri_ ZGETRI
#else /* (PETSC_BLASLAPACK_C) */
   #define sgetri_ sgetri
   #define dgetri_ dgetri
   #define qgetri_ qgetri
   #define cgetri_ cgetri
   #define zgetri_ zgetri
#endif
#if !defined(PETSC_USE_COMPLEX)
  #if defined(PETSC_USE_REAL_SINGLE)
    #define LAPACKgetri_ sgetri_
  #elif defined(PETSC_USE_REAL_DOUBLE)
    #define LAPACKgetri_ dgetri_
  #else /* (PETSC_USE_REAL_QUAD) */
    #define LAPACKgetri_ qgetri_
  #endif
#else
  #if defined(PETSC_USE_REAL_SINGLE)
    #define LAPACKgetri_ cgetri_
  #elif defined(PETSC_USE_REAL_DOUBLE)
    #define LAPACKgetri_ zgetri_
  #else /* (PETSC_USE_REAL_QUAD) */
    #define LAPACKgetri_ wgetri_
  #endif
#endif
EXTERN_C_BEGIN
extern void LAPACKgetri_(PetscBLASInt*,PetscScalar*,PetscBLASInt*,
                         PetscBLASInt*,PetscScalar*,PetscBLASInt*,
                         PetscBLASInt*);
EXTERN_C_END

#if PETSC_VERSION_LE(3,3,0)
#undef PetscBLASIntCast
#undef  __FUNCT__
#define __FUNCT__ "PetscBLASIntCast"
PETSC_STATIC_INLINE PetscErrorCode PetscBLASIntCast(PetscInt a,PetscBLASInt *b)
{
  PetscFunctionBegin;
#if defined(PETSC_USE_64BIT_INDICES) && !defined(PETSC_HAVE_64BIT_BLAS_INDICES)
  if ((a) > PETSC_BLAS_INT_MAX) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Array too long for BLAS/LAPACK");
#endif
  *b =  (PetscBLASInt)(a);
  PetscFunctionReturn(0);
}
#endif

#endif/*PETIGABL_H*/

max_stars_repo_name: otherlab/petiga

max_issues_repo_name: otherlab/petiga

max_forks_repo_name: otherlab/petiga