text:
/*
 *  vbHmmTsFret.c
 *  Model-specific core functions for VB-HMM-TS.
 *
 *  Created by OKAMOTO Kenji, SAKO Yasushi and RIKEN
 *  Copyright 2011-2015
 *  Cellular Informatics Laboratory, Advance Science Institute, RIKEN, Japan.
 *  All rights reserved.
 *
 *  Ver. 1.1.0
 *  Last modified on 2015.09.17
 */

#include "vbHmmTsFret.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

//// Uncomment one/both of the following defenitions to activate constraint on I or/and K.
//#define  INTENSITY_CAP
//#define  TRANSITION_RATE_CAP
#ifdef INTENSITY_CAP
#define  maxIntensityRatio  10.0
#endif
#ifdef TRANSITION_RATE_CAP
#define  minPhotonNumPerState  10.0
#endif

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

//static int isGlobalAnalysis = 0;

void setFunctions_tsFret(){
    commonFunctions funcs;
    funcs.newModelParameters    = newModelParameters_tsFret;
    funcs.freeModelParameters   = freeModelParameters_tsFret;
    funcs.newModelStats         = newModelStats_tsFret;
    funcs.freeModelStats        = freeModelStats_tsFret;
    funcs.initializeVbHmm       = initializeVbHmm_tsFret;
    funcs.pTilde_z1             = pTilde_z1_tsFret;
    funcs.pTilde_zn_zn1         = pTilde_zn_zn1_tsFret;
    funcs.pTilde_xn_zn          = pTilde_xn_zn_tsFret;
    funcs.calcStatsVars         = calcStatsVars_tsFret;
    funcs.maximization          = maximization_tsFret;
    funcs.varLowerBound         = varLowerBound_tsFret;
    funcs.reorderParameters     = reorderParameters_tsFret;
    funcs.outputResults         = outputResults_tsFret;
    setFunctions( funcs );
}

//void setGFunctions_tsFret(){
//    gCommonFunctions funcs;
//    funcs.newModelParameters    = newModelParameters_tsFret;
//    funcs.freeModelParameters   = freeModelParameters_tsFret;
//    funcs.newModelStats         = newModelStats_tsFret;
//    funcs.freeModelStats        = freeModelStats_tsFret;
//    funcs.newModelStatsG        = newModelStatsG_tsFret;
//    funcs.freeModelStatsG       = freeModelStatsG_tsFret;
//    funcs.initializeVbHmmG      = initializeVbHmmG_tsFret;
//    funcs.pTilde_z1             = pTilde_z1_tsFret;
//    funcs.pTilde_zn_zn1         = pTilde_zn_zn1_tsFret;
//    funcs.pTilde_xn_zn          = pTilde_xn_zn_tsFret;
//    funcs.calcStatsVarsG        = calcStatsVarsG_tsFret;
//    funcs.maximizationG         = maximizationG_tsFret;
//    funcs.varLowerBoundG        = varLowerBoundG_tsFret;
//    funcs.reorderParametersG    = reorderParametersG_tsFret;
//    funcs.outputResultsG        = outputResultsG_tsFret;
//    setGFunctions( funcs );
//    isGlobalAnalysis = 1;
//}

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

//void outputResultsG_tsFret( xns, gv, ivs, logFP )
//xnDataBundle *xns;
//globalVars *gv;
//indVarBundle *ivs;
//FILE *logFP;
//{
//    outputTsFretResultsG( xns, gv, ivs, logFP );
//}

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

    p->uPiArr = (double*)malloc( sNo * sizeof(double) );
    p->sumUPi = 0.0;
    // hyperparameter for p( k(i,j) ) (i != j)
    p->uKMat = (double**)malloc( sNo * sizeof(double*) );
    for( i = 0 ; i < sNo ; i++ ){
        p->uKMat[i] = (double*)malloc( sNo * sizeof(double) );
    }
    p->sumUKArr = (double*)malloc( sNo * sizeof(double) );
    // hyperparameter for p( I(k) )
    p->aIArr = (double*)malloc( sNo * sizeof(double) );
    p->bIArr = (double*)malloc( sNo * sizeof(double) );
    // hyperparameter for p( E(i) )
    p->uEArr = (double*)malloc( sNo * sizeof(double) );
    p->vEArr = (double*)malloc( sNo * sizeof(double) );

    // parameters
    p->avgPi = (double *)malloc( sNo * sizeof(double) );
    p->avgLnPi = (double *)malloc( sNo * sizeof(double) );
    p->avgI = (double *)malloc( sNo * sizeof(double) );
    p->avgLnI = (double *)malloc( sNo * sizeof(double) );
    p->avgK = (double **)malloc( sNo * sizeof(double*) );
    p->avgLnK = (double **)malloc( sNo * sizeof(double*) );
    for( i = 0 ; i < sNo ; i++ ){
        p->avgK[i] = (double *)malloc( sNo * sizeof(double) );
        p->avgLnK[i] = (double *)malloc( sNo * sizeof(double) );
    }
    p->avgLnKI = (double *)malloc( sNo * sizeof(double) );
    p->avgE = (double *)malloc( sNo * sizeof(double) );
    p->avgLnE = (double **)malloc( sNo * sizeof(double*) );
    for( i = 0 ; i < sNo ; i++ )
    {   p->avgLnE[i] = (double *)malloc( 2 * sizeof(double) );   }

    return p;
}

void freeModelParameters_tsFret( p, xn, sNo )
void **p;
xnDataSet *xn;
int sNo;
{
    tsFretParameters *gp = *p;
    int i;

    free( gp->uPiArr );
    for( i = 0 ; i < sNo ; i++ ){
        free( gp->uKMat[i] );
    }
    free( gp->uKMat );
    free( gp->sumUKArr );
    free( gp->aIArr );
    free( gp->bIArr );
    free( gp->uEArr );
    free( gp->vEArr );

    free( gp->avgPi );
    free( gp->avgLnPi );
    free( gp->avgI );
    free( gp->avgLnI );
    for( i = 0 ; i < sNo ; i++ ){
        free( gp->avgK[i] );
        free( gp->avgLnK[i] );
    }
    free( gp->avgK );
    free( gp->avgLnK );
    free( gp->avgLnKI );
    free( gp->avgE );
    for( i = 0 ; i < sNo ; i++ ){
        free( gp->avgLnE[i] );
    }
    free( gp->avgLnE );

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

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

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

    return s;

//    } else {
//
//        return NULL;
//
//    }
}

void freeModelStats_tsFret( s, xn, gv, iv )
void **s;
xnDataSet *xn;
globalVars *gv;
indVars *iv;
{
//    if( isGlobalAnalysis == 0 ){
        int sNo = gv->sNo;
        tsFretStats *gs = *s;
        int i;
        free( gs->Ni );
        free( gs->Ti );
        free( gs->Nii );
        free( gs->Nij );
        for( i = 0 ; i < sNo ; i++ )
        {   free( gs->Mij[i] );   }
        free( gs->Mij );
        free( gs->eps );

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

//void *newModelStatsG_tsFret( xns, gv, ivs)
//xnDataBundle *xns;
//globalVars *gv;
//indVarBundle *ivs;
//{
//    int sNo = gv->sNo;
//    tsGlobalStats *gs = (tsGlobalStats*)malloc( sizeof(tsGlobalStats) );
//
//    return gs;
//}

//void freeModelStatsG_tsFret( gs, xns, gv, ivs )
//void **gs;
//xnDataBundle *xns;
//globalVars *gv;
//indVarBundle *ivs;
//{
//
//    free( *gs );
//    *gs = NULL;
//}

void initializeVbHmm_tsFret( xn, gv, iv )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
{
    tsFretData *d = xn->data;
    int sNo = gv->sNo;
    tsFretParameters *p = gv->params;
    int i, j;

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

    // hyperparameter for p( k(i,j) ) (i != j)
    for( i = 0 ; i < sNo ; i++ ){
        p->uKMat[i] = (double*)malloc( sNo * sizeof(double) );
        p->sumUKArr[i] = 0.0;
        for( j = 0 ; j < sNo ; j++ ){
            p->uKMat[i][j] = 1.0;
            if( j != i ){
                p->sumUKArr[i] += p->uKMat[i][j];
            }
        }
    }

    double meanI = (double)xn->N / d->T;

    // hyperparameter for p( I(k) )
    for( i = 0 ; i < sNo ; i++ ){
        p->aIArr[i] = 1.0;
        p->bIArr[i] = 1.0 / meanI;
    }

    // hyperparameter for p( E(i) )
    for( i = 0 ; i < sNo ; i++ ){
        p->uEArr[i] = 1.0;
        p->vEArr[i] = 1.0;
    }

    initialize_indVars_tsFret( xn, gv, iv );

    calcStatsVars_tsFret( xn, gv, iv );
    maximization_tsFret( xn, gv, iv );
}

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

void initialize_indVars_tsFret( 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_tsFret( 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 = (tsFretData*)malloc( sizeof(tsFretData) );
    tsFretData *d = (tsFretData*)xn->data;
    d->T = 0.0;
    d->dt = NULL;
    d->time = NULL;
    d->ch = NULL;
    return xn;
}

void freeXnDataSet_tsFret( xn )
xnDataSet **xn;
{
    tsFretData *d = (tsFretData*)(*xn)->data;
    free( d->dt );
    free( d->time );
    free( d->ch );
    free( (*xn)->data );
    free( (*xn)->name );
    free( *xn );
    *xn = NULL;
}

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

double pTilde_zn_zn1_tsFret( i, j, params )
int i, j;
void *params;
{
    tsFretParameters *p = (tsFretParameters*)params;
    if( i == j ){
        return exp( p->avgLnI[i] - p->avgLnKI[i] );
    } else {
        return exp( p->avgLnK[i][i] + p->avgLnK[i][j] - p->avgLnKI[i] );
    }
}

double pTilde_xn_zn_tsFret( xn, n, i, params )
xnDataSet *xn;
size_t n;
int i;
void *params;
{
    tsFretParameters *p = (tsFretParameters*)params;
    tsFretData *d = (tsFretData*)xn->data;
    return exp( p->avgLnI[i] - (p->avgI[i] * d->dt[n]) + p->avgLnE[i][d->ch[n]] );
}

void calcStatsVars_tsFret( xn, gv, iv )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
{
    tsFretData *d = (tsFretData*)xn->data;
    tsFretStats *s = (tsFretStats*)iv->stats;
    size_t dLen = xn->N;
    int sNo = gv->sNo;
    double **gmMat = iv->gmMat, ***xiMat = iv->xiMat;
    double *Ni = s->Ni, *Ti = s->Ti, *eps = s->eps;
    double *Nii = s->Nii, *Nij = s->Nij, **Mij = s->Mij;
    size_t n;
    int i, j;

    for( i = 0 ; i < sNo ; i++ ){
        Ni[i] = 1e-10;
        Ti[i] = 1e-10;
        eps[i] = 1e-10;
        Nii[i] = 0.0;
        Nij[i] = 0.0;
        for( j = 0 ; j < sNo ; j++ ){
            Mij[i][j] = 1e-10;
        }
    }
    for( n = 0 ; n < dLen ; n++ ){
        for( i = 0 ; i < sNo ; i++ ){
            Ni[i]  += gmMat[n][i];
            Ti[i]  += gmMat[n][i] * d->dt[n];
            eps[i] += gmMat[n][i] * d->ch[n];
            Nii[i] += xiMat[n][i][i];
            for( j = 0 ; j < sNo ; j++ ){
                if( j != i ){
                    Mij[i][j] += xiMat[n][i][j];
                    Nij[i]    += xiMat[n][i][j];
                }
            }
        }
    }
    for( i = 0 ; i < sNo ; i++ ){
        Nii[i] = MAX( Nii[i], 1.0 );
        Nij[i] = MAX( Nij[i], 1.0 );
    }
}

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

void maximization_tsFret( xn, gv, iv )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
{
    tsFretParameters *p = (tsFretParameters*)gv->params;
    tsFretStats *s = (tsFretStats*)iv->stats;
    int sNo = gv->sNo;
    double **gmMat = iv->gmMat;
    double *avgPi = p->avgPi, *avgLnPi = p->avgLnPi, **avgK = p->avgK, **avgLnK = p->avgLnK;
    double *avgLnKI = p->avgLnKI, *avgI = p->avgI, *avgLnI = p->avgLnI;
    double *avgE = p->avgE, **avgLnE = p->avgLnE;
    double *uPiArr = p->uPiArr, sumUPi = p->sumUPi, *aIArr = p->aIArr, *bIArr = p->bIArr;
    double **uKMat = p->uKMat, *sumUKArr = p->sumUKArr;
    double *uEArr = p->uEArr, *vEArr = p->vEArr;
    double *Ni = s->Ni, *Ti = s->Ti, *eps = s->eps;
    double *Nii = s->Nii, *Nij = s->Nij, **Mij = s->Mij;
    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 );

        avgK[i][i] = (Ni[i] + aIArr[i]) * Nij[i] / (Nii[i] - 1.0) / (Ti[i] + bIArr[i]);
        avgLnK[i][i]  = gsl_sf_psi(Ni[i] + aIArr[i]) + gsl_sf_psi(Nij[i]);
        avgLnK[i][i] += -gsl_sf_psi(Nii[i]) - log(Ti[i] + bIArr[i]);

        avgI[i] = (Ni[i] + aIArr[i]) / (Ti[i] + bIArr[i]);
        avgLnI[i] = gsl_sf_psi(Ni[i] + aIArr[i]) - log(Ti[i] + bIArr[i]);

        avgLnKI[i]  = gsl_sf_psi(Ni[i] + aIArr[i]) + gsl_sf_psi(Nii[i] + Nij[i]);
        avgLnKI[i] += -gsl_sf_psi(Nii[i]) - log(Ti[i] + bIArr[i]);

#ifdef INTENSITY_CAP
        double meanI = (double)xnWv->N / xnWv->T;
        avgI[i]  = MIN( avgI[i], maxIntensityRatio * meanI );
        avgLnI[i]  = MIN( avgLnI[i], log(maxIntensityRatio * meanI) );
#endif
#ifdef TRANSITION_RATE_CAP
        avgK[i][i] = MIN( avgK[i][i], avgI[i]/minPhotonNumPerState );
        avgLnK[i][i]  = MIN( avgLnK[i][i], avgLnI[i] - log(minPhotonNumPerState) );
        avgLnKI[i]  = MIN( avgLnKI[i], avgLnI[i] + log(1.0 + 1.0/minPhotonNumPerState) );
#endif

        for( j = 0 ; j < sNo ; j++ ){
            if( i != j ){
                avgK[i][j] = ( uKMat[i][j] + Mij[i][j] ) / ( sumUKArr[i] + Nij[i] );
                avgLnK[i][j] = gsl_sf_psi( uKMat[i][j] + Mij[i][j] ) - gsl_sf_psi( sumUKArr[i] + Nij[i] );
            }
        }

        avgE[i] = ( uEArr[i] + eps[i] ) / ( uEArr[i] + vEArr[i] + Ni[i] );
        avgLnE[i][0]  = gsl_sf_psi( vEArr[i] + Ni[i] - eps[i] );    // ln(1-E) for donor
        avgLnE[i][0] -= gsl_sf_psi( uEArr[i] + vEArr[i] + Ni[i] );
        avgLnE[i][1]  = gsl_sf_psi( uEArr[i] + eps[i] );            // ln(E) for acceptor
        avgLnE[i][1] -= gsl_sf_psi( uEArr[i] + vEArr[i] + Ni[i] );
    }
}

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

double varLowerBound_tsFret( xn, gv, iv )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
{
    tsFretParameters *p = (tsFretParameters*)gv->params;
    tsFretStats *s = (tsFretStats*)iv->stats;
    size_t dLen = xn->N;
    int sNo = gv->sNo;
    double **gmMat = iv->gmMat, *cn = iv->cn;
    double *avgLnPi = p->avgLnPi, **avgLnK = p->avgLnK, *avgI = p->avgI, *avgLnI = p->avgLnI;
    double *avgLnKI = p->avgLnKI, **avgLnE = p->avgLnE;
    double *uPiArr = p->uPiArr, sumUPi = p->sumUPi, *aIArr = p->aIArr, *bIArr = p->bIArr;
    double **uKMat = p->uKMat, *sumUKArr = p->sumUKArr;
    double *uEArr = p->uEArr, *vEArr = p->vEArr;
    double *Ni = s->Ni, *Ti = s->Ti, *eps = s->eps;
    double *Nii = s->Nii, *Nij = s->Nij, **Mij = s->Mij;
    size_t n;
    int i, j;

    double lnpPi = gsl_sf_lngamma(sumUPi);
    double Ck = 1.0, lnpKii = sNo * log(Ck);
    double lnpKij = 0.0;
    double lnpI = 0.0;
    double lnpE = 0.0;
    double lnqPi = gsl_sf_lngamma(sumUPi + 1);
    double lnqKiiI = 0.0;
    double lnqKij = 0.0;
    double lnqE = 0.0;
    for( i = 0 ; i < sNo ; i++ ){
        lnpPi += (uPiArr[i]-1.0) * avgLnPi[i] - gsl_sf_lngamma(uPiArr[i]);

        lnpKii -= avgLnK[i][i];

        lnpI += aIArr[i] * log(bIArr[i]) - gsl_sf_lngamma(aIArr[i]);
        lnpI += (aIArr[i]-1.0)*avgLnI[i] - bIArr[i]*avgI[i];

        lnpE += gsl_sf_lngamma(uEArr[i]+vEArr[i]) - gsl_sf_lngamma(uEArr[i]);
        lnpE += -gsl_sf_lngamma(vEArr[i]);
        lnpE += (uEArr[i]-1.0)*avgLnE[i][1] + (vEArr[i]-1.0)*avgLnE[i][0];

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

        lnqKiiI += gsl_sf_lngamma(Nii[i] + Nij[i]) + (Ni[i] + aIArr[i]) * log(Ti[i] + bIArr[i]);
        lnqKiiI += -gsl_sf_lngamma(Ni[i] + aIArr[i]) -gsl_sf_lngamma(Nii[i]);
        lnqKiiI += -gsl_sf_lngamma(Nij[i]) + (Nij[i]-1.0)*avgLnK[i][i];
        lnqKiiI += (Ni[i] + Nii[i] + aIArr[i] - 1.0)*avgLnI[i];
        lnqKiiI += -(Nii[i] + Nij[i])*avgLnKI[i] - (Ti[i] + bIArr[i])*avgI[i];

        lnqE += gsl_sf_lngamma(uEArr[i] + vEArr[i] + Ni[i]);
        lnqE -= gsl_sf_lngamma(uEArr[i] + eps[i]);
        lnqE -= gsl_sf_lngamma(vEArr[i] + Ni[i] - eps[i]);
        lnqE += (uEArr[i] + eps[i] - 1.0)*avgLnE[i][1];
        lnqE += (vEArr[i] + Ni[i] - eps[i] - 1.0)*avgLnE[i][0];

        lnpKij += gsl_sf_lngamma(sumUKArr[i]);
        lnqKij += gsl_sf_lngamma(sumUKArr[i] + Nij[i]);
        for( j = 0 ; j < sNo ; j++ ){
            if( j != i ){
                lnpKij += (uKMat[i][j]-1.0)*avgLnK[i][j] - gsl_sf_lngamma(uKMat[i][j]);
                lnqKij += (uKMat[i][j]+Mij[i][j]-1)*(gsl_sf_psi(uKMat[i][j]+Mij[i][j])-gsl_sf_psi(sumUKArr[i]+Nij[i]));
                lnqKij -= gsl_sf_lngamma( uKMat[i][j] + Mij[i][j] );
            }
        }
    }

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

    double val;
    val  = lnpPi + lnpKii + lnpKij + lnpI + lnpE;
    val -= lnqPi + lnqKiiI + lnqKij + lnqE;
    val += lnpX;

    return val;
}

//double varLowerBoundG_tsFret( xns, gv, ivs )
//xnDataBundle *xns;
//globalVars *gv;
//indVarBundle *ivs;
//{
//}

void reorderParameters_tsFret( xn, gv, iv )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
{
    tsFretParameters *p = (tsFretParameters*)gv->params;
    tsFretStats *s = (tsFretStats*)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, **avgK = p->avgK;
    double **avgLnK = p->avgLnK, *avgLnKI = p->avgLnKI;
    double *avgI = p->avgI, *avgLnI = p->avgLnI, *avgE = p->avgE, **avgLnE = p->avgLnE;
    double *Ni = s->Ni, *Ti = s->Ti, *eps = s->eps;
    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 avgE values (0=biggest avgE -- sNo=smallest avgE).
    for( i = 0 ; i < sNo ; i++ ){
        index[i] = sNo - 1;
        for( j = 0 ; j < sNo ; j++ ){
            if( j != i ){
                if( avgE[i] < avgE[j] ){
                    index[i]--;
                } else if( avgE[i] == avgE[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( i = 0 ; i < sNo ; i++ ){   store[index[i]] = avgI[i];   }
    for( i = 0 ; i < sNo ; i++ ){   avgI[i] = store[i];   }

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

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

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

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

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

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

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

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

    for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = eps[i];   }
    for( i = 0 ; i < sNo ; i++ ){   eps[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_tsFret( xns, gv, ivs )
//xnDataBundle *xns;
//globalVars *gv;
//indVarBundle *ivs;
//{
//}

void outputTsFretResults( xn, gv, iv, logFP )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
FILE *logFP;
{
    tsFretParameters *p = (tsFretParameters*)gv->params;
    int sNo = gv->sNo;
    int i, j;

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

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

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

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

    fprintf(logFP, "   k_matrix: [");
    for( i = 0 ; i < sNo ; i++ ){
            fprintf(logFP, " ( %g", p->avgK[i][0]);
            for( j = 1 ; j < sNo ; j++ )
            {   fprintf(logFP, ", %g", p->avgK[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, "I, E, pi");
        for( i = 0 ; i < sNo ; i++ )
        {   fprintf(fp, ", K%dx", i);   }
        fprintf(fp, "\n");

        for( i = 0 ; i < sNo ; i++ ){
            fprintf(fp, "%g, %g, %g", p->avgI[i], p->avgE[i], p->avgPi[i]);
            for( j = 0 ; j < sNo ; j++ )
            {   fprintf(fp, ", %g", p->avgK[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 outputTsFretResultsG( xns, gv, ivs, logFP )
//xnDataBundle *xns;
//globalVars *gv;
//indVarBundle *ivs;
//FILE *logFP;
//{
//}

//

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) 1997-1999 Massachusetts Institute of Technology
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 */

/* This file was automatically generated --- DO NOT EDIT */
/* Generated on Sun Nov  7 20:45:16 EST 1999 */

#include <fftw-int.h>
#include <fftw.h>

/* Generated by: ./genfft -magic-alignment-check -magic-twiddle-load-all -magic-variables 4 -magic-loopi -hc2hc-backward 16 */

/*
 * This function contains 298 FP additions, 134 FP multiplications,
 * (or, 244 additions, 80 multiplications, 54 fused multiply/add),
 * 49 stack variables, and 128 memory accesses
 */
static const fftw_real K1_961570560 = FFTW_KONST(+1.961570560806460898252364472268478073947867462);
static const fftw_real K390180644 = FFTW_KONST(+0.390180644032256535696569736954044481855383236);
static const fftw_real K1_111140466 = FFTW_KONST(+1.111140466039204449485661627897065748749874382);
static const fftw_real K1_662939224 = FFTW_KONST(+1.662939224605090474157576755235811513477121624);
static const fftw_real K707106781 = FFTW_KONST(+0.707106781186547524400844362104849039284835938);
static const fftw_real K382683432 = FFTW_KONST(+0.382683432365089771728459984030398866761344562);
static const fftw_real K923879532 = FFTW_KONST(+0.923879532511286756128183189396788286822416626);
static const fftw_real K765366864 = FFTW_KONST(+0.765366864730179543456919968060797733522689125);
static const fftw_real K1_847759065 = FFTW_KONST(+1.847759065022573512256366378793576573644833252);
static const fftw_real K1_414213562 = FFTW_KONST(+1.414213562373095048801688724209698078569671875);
static const fftw_real K2_000000000 = FFTW_KONST(+2.000000000000000000000000000000000000000000000);

/*
 * Generator Id's :
 * $Id: exprdag.ml,v 1.41 1999/05/26 15:44:14 fftw Exp $
 * $Id: fft.ml,v 1.43 1999/05/17 19:44:18 fftw Exp $
 * $Id: to_c.ml,v 1.25 1999/10/26 21:41:32 stevenj Exp $
 */

void fftw_hc2hc_backward_16(fftw_real *A, const fftw_complex *W, int iostride, int m, int dist)
{
     int i;
     fftw_real *X;
     fftw_real *Y;
     X = A;
     Y = A + (16 * iostride);
     {
	  fftw_real tmp279;
	  fftw_real tmp324;
	  fftw_real tmp312;
	  fftw_real tmp299;
	  fftw_real tmp276;
	  fftw_real tmp296;
	  fftw_real tmp309;
	  fftw_real tmp323;
	  fftw_real tmp283;
	  fftw_real tmp291;
	  fftw_real tmp286;
	  fftw_real tmp294;
	  fftw_real tmp301;
	  fftw_real tmp319;
	  fftw_real tmp327;
	  fftw_real tmp326;
	  fftw_real tmp316;
	  fftw_real tmp302;
	  ASSERT_ALIGNED_DOUBLE;
	  {
	       fftw_real tmp277;
	       fftw_real tmp278;
	       fftw_real tmp310;
	       fftw_real tmp297;
	       fftw_real tmp298;
	       fftw_real tmp311;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp277 = X[2 * iostride];
	       tmp278 = X[6 * iostride];
	       tmp310 = tmp277 - tmp278;
	       tmp297 = Y[-2 * iostride];
	       tmp298 = Y[-6 * iostride];
	       tmp311 = tmp298 + tmp297;
	       tmp279 = K2_000000000 * (tmp277 + tmp278);
	       tmp324 = K1_414213562 * (tmp310 + tmp311);
	       tmp312 = K1_414213562 * (tmp310 - tmp311);
	       tmp299 = K2_000000000 * (tmp297 - tmp298);
	  }
	  {
	       fftw_real tmp275;
	       fftw_real tmp308;
	       fftw_real tmp273;
	       fftw_real tmp306;
	       ASSERT_ALIGNED_DOUBLE;
	       {
		    fftw_real tmp274;
		    fftw_real tmp307;
		    fftw_real tmp271;
		    fftw_real tmp272;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp274 = X[4 * iostride];
		    tmp275 = K2_000000000 * tmp274;
		    tmp307 = Y[-4 * iostride];
		    tmp308 = K2_000000000 * tmp307;
		    tmp271 = X[0];
		    tmp272 = X[8 * iostride];
		    tmp273 = tmp271 + tmp272;
		    tmp306 = tmp271 - tmp272;
	       }
	       tmp276 = tmp273 + tmp275;
	       tmp296 = tmp273 - tmp275;
	       tmp309 = tmp306 - tmp308;
	       tmp323 = tmp306 + tmp308;
	  }
	  {
	       fftw_real tmp314;
	       fftw_real tmp318;
	       fftw_real tmp317;
	       fftw_real tmp315;
	       ASSERT_ALIGNED_DOUBLE;
	       {
		    fftw_real tmp281;
		    fftw_real tmp282;
		    fftw_real tmp289;
		    fftw_real tmp290;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp281 = X[iostride];
		    tmp282 = X[7 * iostride];
		    tmp283 = tmp281 + tmp282;
		    tmp314 = tmp281 - tmp282;
		    tmp289 = Y[-iostride];
		    tmp290 = Y[-7 * iostride];
		    tmp291 = tmp289 - tmp290;
		    tmp318 = tmp289 + tmp290;
	       }
	       {
		    fftw_real tmp284;
		    fftw_real tmp285;
		    fftw_real tmp292;
		    fftw_real tmp293;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp284 = X[3 * iostride];
		    tmp285 = X[5 * iostride];
		    tmp286 = tmp284 + tmp285;
		    tmp317 = tmp285 - tmp284;
		    tmp292 = Y[-3 * iostride];
		    tmp293 = Y[-5 * iostride];
		    tmp294 = tmp292 - tmp293;
		    tmp315 = tmp293 + tmp292;
	       }
	       tmp301 = tmp283 - tmp286;
	       tmp319 = tmp317 + tmp318;
	       tmp327 = tmp318 - tmp317;
	       tmp326 = tmp314 + tmp315;
	       tmp316 = tmp314 - tmp315;
	       tmp302 = tmp294 + tmp291;
	  }
	  {
	       fftw_real tmp280;
	       fftw_real tmp287;
	       fftw_real tmp288;
	       fftw_real tmp295;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp280 = tmp276 + tmp279;
	       tmp287 = K2_000000000 * (tmp283 + tmp286);
	       X[8 * iostride] = tmp280 - tmp287;
	       X[0] = tmp280 + tmp287;
	       tmp288 = tmp276 - tmp279;
	       tmp295 = K2_000000000 * (tmp291 - tmp294);
	       X[12 * iostride] = tmp288 + tmp295;
	       X[4 * iostride] = tmp288 - tmp295;
	  }
	  {
	       fftw_real tmp300;
	       fftw_real tmp303;
	       fftw_real tmp304;
	       fftw_real tmp305;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp300 = tmp296 - tmp299;
	       tmp303 = K1_414213562 * (tmp301 - tmp302);
	       X[10 * iostride] = tmp300 - tmp303;
	       X[2 * iostride] = tmp300 + tmp303;
	       tmp304 = tmp296 + tmp299;
	       tmp305 = K1_414213562 * (tmp301 + tmp302);
	       X[6 * iostride] = tmp304 - tmp305;
	       X[14 * iostride] = tmp304 + tmp305;
	  }
	  {
	       fftw_real tmp313;
	       fftw_real tmp320;
	       fftw_real tmp321;
	       fftw_real tmp322;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp313 = tmp309 + tmp312;
	       tmp320 = (K1_847759065 * tmp316) - (K765366864 * tmp319);
	       X[9 * iostride] = tmp313 - tmp320;
	       X[iostride] = tmp313 + tmp320;
	       tmp321 = tmp309 - tmp312;
	       tmp322 = (K765366864 * tmp316) + (K1_847759065 * tmp319);
	       X[5 * iostride] = tmp321 - tmp322;
	       X[13 * iostride] = tmp321 + tmp322;
	  }
	  {
	       fftw_real tmp325;
	       fftw_real tmp328;
	       fftw_real tmp329;
	       fftw_real tmp330;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp325 = tmp323 - tmp324;
	       tmp328 = (K765366864 * tmp326) - (K1_847759065 * tmp327);
	       X[11 * iostride] = tmp325 - tmp328;
	       X[3 * iostride] = tmp325 + tmp328;
	       tmp329 = tmp323 + tmp324;
	       tmp330 = (K1_847759065 * tmp326) + (K765366864 * tmp327);
	       X[7 * iostride] = tmp329 - tmp330;
	       X[15 * iostride] = tmp329 + tmp330;
	  }
     }
     X = X + dist;
     Y = Y - dist;
     for (i = 2; i < m; i = i + 2, X = X + dist, Y = Y - dist, W = W + 15) {
	  fftw_real tmp73;
	  fftw_real tmp98;
	  fftw_real tmp135;
	  fftw_real tmp160;
	  fftw_real tmp182;
	  fftw_real tmp236;
	  fftw_real tmp210;
	  fftw_real tmp248;
	  fftw_real tmp95;
	  fftw_real tmp124;
	  fftw_real tmp138;
	  fftw_real tmp164;
	  fftw_real tmp197;
	  fftw_real tmp216;
	  fftw_real tmp244;
	  fftw_real tmp252;
	  fftw_real tmp80;
	  fftw_real tmp128;
	  fftw_real tmp105;
	  fftw_real tmp161;
	  fftw_real tmp213;
	  fftw_real tmp237;
	  fftw_real tmp189;
	  fftw_real tmp249;
	  fftw_real tmp88;
	  fftw_real tmp115;
	  fftw_real tmp137;
	  fftw_real tmp163;
	  fftw_real tmp204;
	  fftw_real tmp215;
	  fftw_real tmp241;
	  fftw_real tmp251;
	  ASSERT_ALIGNED_DOUBLE;
	  {
	       fftw_real tmp69;
	       fftw_real tmp180;
	       fftw_real tmp131;
	       fftw_real tmp209;
	       fftw_real tmp72;
	       fftw_real tmp208;
	       fftw_real tmp134;
	       fftw_real tmp181;
	       ASSERT_ALIGNED_DOUBLE;
	       {
		    fftw_real tmp67;
		    fftw_real tmp68;
		    fftw_real tmp129;
		    fftw_real tmp130;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp67 = X[0];
		    tmp68 = Y[-8 * iostride];
		    tmp69 = tmp67 + tmp68;
		    tmp180 = tmp67 - tmp68;
		    tmp129 = Y[0];
		    tmp130 = X[8 * iostride];
		    tmp131 = tmp129 - tmp130;
		    tmp209 = tmp129 + tmp130;
	       }
	       {
		    fftw_real tmp70;
		    fftw_real tmp71;
		    fftw_real tmp132;
		    fftw_real tmp133;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp70 = X[4 * iostride];
		    tmp71 = Y[-12 * iostride];
		    tmp72 = tmp70 + tmp71;
		    tmp208 = tmp70 - tmp71;
		    tmp132 = Y[-4 * iostride];
		    tmp133 = X[12 * iostride];
		    tmp134 = tmp132 - tmp133;
		    tmp181 = tmp132 + tmp133;
	       }
	       tmp73 = tmp69 + tmp72;
	       tmp98 = tmp69 - tmp72;
	       tmp135 = tmp131 - tmp134;
	       tmp160 = tmp131 + tmp134;
	       tmp182 = tmp180 - tmp181;
	       tmp236 = tmp180 + tmp181;
	       tmp210 = tmp208 + tmp209;
	       tmp248 = tmp209 - tmp208;
	  }
	  {
	       fftw_real tmp91;
	       fftw_real tmp194;
	       fftw_real tmp119;
	       fftw_real tmp192;
	       fftw_real tmp94;
	       fftw_real tmp191;
	       fftw_real tmp122;
	       fftw_real tmp195;
	       fftw_real tmp116;
	       fftw_real tmp123;
	       ASSERT_ALIGNED_DOUBLE;
	       {
		    fftw_real tmp89;
		    fftw_real tmp90;
		    fftw_real tmp117;
		    fftw_real tmp118;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp89 = Y[-15 * iostride];
		    tmp90 = X[7 * iostride];
		    tmp91 = tmp89 + tmp90;
		    tmp194 = tmp89 - tmp90;
		    tmp117 = Y[-7 * iostride];
		    tmp118 = X[15 * iostride];
		    tmp119 = tmp117 - tmp118;
		    tmp192 = tmp117 + tmp118;
	       }
	       {
		    fftw_real tmp92;
		    fftw_real tmp93;
		    fftw_real tmp120;
		    fftw_real tmp121;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp92 = X[3 * iostride];
		    tmp93 = Y[-11 * iostride];
		    tmp94 = tmp92 + tmp93;
		    tmp191 = tmp92 - tmp93;
		    tmp120 = Y[-3 * iostride];
		    tmp121 = X[11 * iostride];
		    tmp122 = tmp120 - tmp121;
		    tmp195 = tmp120 + tmp121;
	       }
	       tmp95 = tmp91 + tmp94;
	       tmp116 = tmp91 - tmp94;
	       tmp123 = tmp119 - tmp122;
	       tmp124 = tmp116 + tmp123;
	       tmp138 = tmp123 - tmp116;
	       tmp164 = tmp119 + tmp122;
	       {
		    fftw_real tmp193;
		    fftw_real tmp196;
		    fftw_real tmp242;
		    fftw_real tmp243;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp193 = tmp191 - tmp192;
		    tmp196 = tmp194 - tmp195;
		    tmp197 = (K923879532 * tmp193) - (K382683432 * tmp196);
		    tmp216 = (K382683432 * tmp193) + (K923879532 * tmp196);
		    tmp242 = tmp194 + tmp195;
		    tmp243 = tmp191 + tmp192;
		    tmp244 = (K382683432 * tmp242) - (K923879532 * tmp243);
		    tmp252 = (K382683432 * tmp243) + (K923879532 * tmp242);
	       }
	  }
	  {
	       fftw_real tmp76;
	       fftw_real tmp183;
	       fftw_real tmp104;
	       fftw_real tmp184;
	       fftw_real tmp79;
	       fftw_real tmp186;
	       fftw_real tmp101;
	       fftw_real tmp187;
	       ASSERT_ALIGNED_DOUBLE;
	       {
		    fftw_real tmp74;
		    fftw_real tmp75;
		    fftw_real tmp102;
		    fftw_real tmp103;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp74 = X[2 * iostride];
		    tmp75 = Y[-10 * iostride];
		    tmp76 = tmp74 + tmp75;
		    tmp183 = tmp74 - tmp75;
		    tmp102 = Y[-2 * iostride];
		    tmp103 = X[10 * iostride];
		    tmp104 = tmp102 - tmp103;
		    tmp184 = tmp102 + tmp103;
	       }
	       {
		    fftw_real tmp77;
		    fftw_real tmp78;
		    fftw_real tmp99;
		    fftw_real tmp100;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp77 = Y[-14 * iostride];
		    tmp78 = X[6 * iostride];
		    tmp79 = tmp77 + tmp78;
		    tmp186 = tmp77 - tmp78;
		    tmp99 = Y[-6 * iostride];
		    tmp100 = X[14 * iostride];
		    tmp101 = tmp99 - tmp100;
		    tmp187 = tmp99 + tmp100;
	       }
	       tmp80 = tmp76 + tmp79;
	       tmp128 = tmp76 - tmp79;
	       tmp105 = tmp101 - tmp104;
	       tmp161 = tmp104 + tmp101;
	       {
		    fftw_real tmp211;
		    fftw_real tmp212;
		    fftw_real tmp185;
		    fftw_real tmp188;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp211 = tmp183 + tmp184;
		    tmp212 = tmp186 + tmp187;
		    tmp213 = K707106781 * (tmp211 - tmp212);
		    tmp237 = K707106781 * (tmp211 + tmp212);
		    tmp185 = tmp183 - tmp184;
		    tmp188 = tmp186 - tmp187;
		    tmp189 = K707106781 * (tmp185 + tmp188);
		    tmp249 = K707106781 * (tmp185 - tmp188);
	       }
	  }
	  {
	       fftw_real tmp84;
	       fftw_real tmp201;
	       fftw_real tmp110;
	       fftw_real tmp199;
	       fftw_real tmp87;
	       fftw_real tmp198;
	       fftw_real tmp113;
	       fftw_real tmp202;
	       fftw_real tmp107;
	       fftw_real tmp114;
	       ASSERT_ALIGNED_DOUBLE;
	       {
		    fftw_real tmp82;
		    fftw_real tmp83;
		    fftw_real tmp108;
		    fftw_real tmp109;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp82 = X[iostride];
		    tmp83 = Y[-9 * iostride];
		    tmp84 = tmp82 + tmp83;
		    tmp201 = tmp82 - tmp83;
		    tmp108 = Y[-iostride];
		    tmp109 = X[9 * iostride];
		    tmp110 = tmp108 - tmp109;
		    tmp199 = tmp108 + tmp109;
	       }
	       {
		    fftw_real tmp85;
		    fftw_real tmp86;
		    fftw_real tmp111;
		    fftw_real tmp112;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp85 = X[5 * iostride];
		    tmp86 = Y[-13 * iostride];
		    tmp87 = tmp85 + tmp86;
		    tmp198 = tmp85 - tmp86;
		    tmp111 = Y[-5 * iostride];
		    tmp112 = X[13 * iostride];
		    tmp113 = tmp111 - tmp112;
		    tmp202 = tmp111 + tmp112;
	       }
	       tmp88 = tmp84 + tmp87;
	       tmp107 = tmp84 - tmp87;
	       tmp114 = tmp110 - tmp113;
	       tmp115 = tmp107 - tmp114;
	       tmp137 = tmp107 + tmp114;
	       tmp163 = tmp110 + tmp113;
	       {
		    fftw_real tmp200;
		    fftw_real tmp203;
		    fftw_real tmp239;
		    fftw_real tmp240;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp200 = tmp198 + tmp199;
		    tmp203 = tmp201 - tmp202;
		    tmp204 = (K923879532 * tmp200) + (K382683432 * tmp203);
		    tmp215 = (K923879532 * tmp203) - (K382683432 * tmp200);
		    tmp239 = tmp201 + tmp202;
		    tmp240 = tmp199 - tmp198;
		    tmp241 = (K382683432 * tmp239) - (K923879532 * tmp240);
		    tmp251 = (K382683432 * tmp240) + (K923879532 * tmp239);
	       }
	  }
	  {
	       fftw_real tmp81;
	       fftw_real tmp96;
	       fftw_real tmp158;
	       fftw_real tmp162;
	       fftw_real tmp165;
	       fftw_real tmp166;
	       fftw_real tmp157;
	       fftw_real tmp159;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp81 = tmp73 + tmp80;
	       tmp96 = tmp88 + tmp95;
	       tmp158 = tmp81 - tmp96;
	       tmp162 = tmp160 + tmp161;
	       tmp165 = tmp163 + tmp164;
	       tmp166 = tmp162 - tmp165;
	       X[0] = tmp81 + tmp96;
	       Y[-15 * iostride] = tmp162 + tmp165;
	       tmp157 = c_re(W[7]);
	       tmp159 = c_im(W[7]);
	       X[8 * iostride] = (tmp157 * tmp158) + (tmp159 * tmp166);
	       Y[-7 * iostride] = (tmp157 * tmp166) - (tmp159 * tmp158);
	  }
	  {
	       fftw_real tmp170;
	       fftw_real tmp176;
	       fftw_real tmp174;
	       fftw_real tmp178;
	       ASSERT_ALIGNED_DOUBLE;
	       {
		    fftw_real tmp168;
		    fftw_real tmp169;
		    fftw_real tmp172;
		    fftw_real tmp173;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp168 = tmp160 - tmp161;
		    tmp169 = tmp88 - tmp95;
		    tmp170 = tmp168 - tmp169;
		    tmp176 = tmp169 + tmp168;
		    tmp172 = tmp73 - tmp80;
		    tmp173 = tmp164 - tmp163;
		    tmp174 = tmp172 - tmp173;
		    tmp178 = tmp172 + tmp173;
	       }
	       {
		    fftw_real tmp167;
		    fftw_real tmp171;
		    fftw_real tmp175;
		    fftw_real tmp177;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp167 = c_re(W[11]);
		    tmp171 = c_im(W[11]);
		    Y[-3 * iostride] = (tmp167 * tmp170) - (tmp171 * tmp174);
		    X[12 * iostride] = (tmp171 * tmp170) + (tmp167 * tmp174);
		    tmp175 = c_re(W[3]);
		    tmp177 = c_im(W[3]);
		    Y[-11 * iostride] = (tmp175 * tmp176) - (tmp177 * tmp178);
		    X[4 * iostride] = (tmp177 * tmp176) + (tmp175 * tmp178);
	       }
	  }
	  {
	       fftw_real tmp126;
	       fftw_real tmp142;
	       fftw_real tmp140;
	       fftw_real tmp144;
	       ASSERT_ALIGNED_DOUBLE;
	       {
		    fftw_real tmp106;
		    fftw_real tmp125;
		    fftw_real tmp136;
		    fftw_real tmp139;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp106 = tmp98 + tmp105;
		    tmp125 = K707106781 * (tmp115 + tmp124);
		    tmp126 = tmp106 - tmp125;
		    tmp142 = tmp106 + tmp125;
		    tmp136 = tmp128 + tmp135;
		    tmp139 = K707106781 * (tmp137 + tmp138);
		    tmp140 = tmp136 - tmp139;
		    tmp144 = tmp136 + tmp139;
	       }
	       {
		    fftw_real tmp97;
		    fftw_real tmp127;
		    fftw_real tmp141;
		    fftw_real tmp143;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp97 = c_re(W[9]);
		    tmp127 = c_im(W[9]);
		    X[10 * iostride] = (tmp97 * tmp126) + (tmp127 * tmp140);
		    Y[-5 * iostride] = (tmp97 * tmp140) - (tmp127 * tmp126);
		    tmp141 = c_re(W[1]);
		    tmp143 = c_im(W[1]);
		    X[2 * iostride] = (tmp141 * tmp142) + (tmp143 * tmp144);
		    Y[-13 * iostride] = (tmp141 * tmp144) - (tmp143 * tmp142);
	       }
	  }
	  {
	       fftw_real tmp148;
	       fftw_real tmp154;
	       fftw_real tmp152;
	       fftw_real tmp156;
	       ASSERT_ALIGNED_DOUBLE;
	       {
		    fftw_real tmp146;
		    fftw_real tmp147;
		    fftw_real tmp150;
		    fftw_real tmp151;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp146 = tmp135 - tmp128;
		    tmp147 = K707106781 * (tmp115 - tmp124);
		    tmp148 = tmp146 - tmp147;
		    tmp154 = tmp146 + tmp147;
		    tmp150 = tmp98 - tmp105;
		    tmp151 = K707106781 * (tmp138 - tmp137);
		    tmp152 = tmp150 - tmp151;
		    tmp156 = tmp150 + tmp151;
	       }
	       {
		    fftw_real tmp145;
		    fftw_real tmp149;
		    fftw_real tmp153;
		    fftw_real tmp155;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp145 = c_re(W[13]);
		    tmp149 = c_im(W[13]);
		    Y[-iostride] = (tmp145 * tmp148) - (tmp149 * tmp152);
		    X[14 * iostride] = (tmp149 * tmp148) + (tmp145 * tmp152);
		    tmp153 = c_re(W[5]);
		    tmp155 = c_im(W[5]);
		    Y[-9 * iostride] = (tmp153 * tmp154) - (tmp155 * tmp156);
		    X[6 * iostride] = (tmp155 * tmp154) + (tmp153 * tmp156);
	       }
	  }
	  {
	       fftw_real tmp206;
	       fftw_real tmp220;
	       fftw_real tmp218;
	       fftw_real tmp222;
	       ASSERT_ALIGNED_DOUBLE;
	       {
		    fftw_real tmp190;
		    fftw_real tmp205;
		    fftw_real tmp214;
		    fftw_real tmp217;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp190 = tmp182 - tmp189;
		    tmp205 = tmp197 - tmp204;
		    tmp206 = tmp190 - tmp205;
		    tmp220 = tmp190 + tmp205;
		    tmp214 = tmp210 - tmp213;
		    tmp217 = tmp215 - tmp216;
		    tmp218 = tmp214 - tmp217;
		    tmp222 = tmp214 + tmp217;
	       }
	       {
		    fftw_real tmp179;
		    fftw_real tmp207;
		    fftw_real tmp219;
		    fftw_real tmp221;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp179 = c_re(W[12]);
		    tmp207 = c_im(W[12]);
		    X[13 * iostride] = (tmp179 * tmp206) + (tmp207 * tmp218);
		    Y[-2 * iostride] = (tmp179 * tmp218) - (tmp207 * tmp206);
		    tmp219 = c_re(W[4]);
		    tmp221 = c_im(W[4]);
		    X[5 * iostride] = (tmp219 * tmp220) + (tmp221 * tmp222);
		    Y[-10 * iostride] = (tmp219 * tmp222) - (tmp221 * tmp220);
	       }
	  }
	  {
	       fftw_real tmp226;
	       fftw_real tmp232;
	       fftw_real tmp230;
	       fftw_real tmp234;
	       ASSERT_ALIGNED_DOUBLE;
	       {
		    fftw_real tmp224;
		    fftw_real tmp225;
		    fftw_real tmp228;
		    fftw_real tmp229;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp224 = tmp210 + tmp213;
		    tmp225 = tmp204 + tmp197;
		    tmp226 = tmp224 - tmp225;
		    tmp232 = tmp224 + tmp225;
		    tmp228 = tmp182 + tmp189;
		    tmp229 = tmp215 + tmp216;
		    tmp230 = tmp228 - tmp229;
		    tmp234 = tmp228 + tmp229;
	       }
	       {
		    fftw_real tmp223;
		    fftw_real tmp227;
		    fftw_real tmp231;
		    fftw_real tmp233;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp223 = c_re(W[8]);
		    tmp227 = c_im(W[8]);
		    Y[-6 * iostride] = (tmp223 * tmp226) - (tmp227 * tmp230);
		    X[9 * iostride] = (tmp227 * tmp226) + (tmp223 * tmp230);
		    tmp231 = c_re(W[0]);
		    tmp233 = c_im(W[0]);
		    Y[-14 * iostride] = (tmp231 * tmp232) - (tmp233 * tmp234);
		    X[iostride] = (tmp233 * tmp232) + (tmp231 * tmp234);
	       }
	  }
	  {
	       fftw_real tmp246;
	       fftw_real tmp256;
	       fftw_real tmp254;
	       fftw_real tmp258;
	       ASSERT_ALIGNED_DOUBLE;
	       {
		    fftw_real tmp238;
		    fftw_real tmp245;
		    fftw_real tmp250;
		    fftw_real tmp253;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp238 = tmp236 - tmp237;
		    tmp245 = tmp241 + tmp244;
		    tmp246 = tmp238 - tmp245;
		    tmp256 = tmp238 + tmp245;
		    tmp250 = tmp248 + tmp249;
		    tmp253 = tmp251 - tmp252;
		    tmp254 = tmp250 - tmp253;
		    tmp258 = tmp250 + tmp253;
	       }
	       {
		    fftw_real tmp235;
		    fftw_real tmp247;
		    fftw_real tmp255;
		    fftw_real tmp257;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp235 = c_re(W[10]);
		    tmp247 = c_im(W[10]);
		    X[11 * iostride] = (tmp235 * tmp246) + (tmp247 * tmp254);
		    Y[-4 * iostride] = (tmp235 * tmp254) - (tmp247 * tmp246);
		    tmp255 = c_re(W[2]);
		    tmp257 = c_im(W[2]);
		    X[3 * iostride] = (tmp255 * tmp256) + (tmp257 * tmp258);
		    Y[-12 * iostride] = (tmp255 * tmp258) - (tmp257 * tmp256);
	       }
	  }
	  {
	       fftw_real tmp262;
	       fftw_real tmp268;
	       fftw_real tmp266;
	       fftw_real tmp270;
	       ASSERT_ALIGNED_DOUBLE;
	       {
		    fftw_real tmp260;
		    fftw_real tmp261;
		    fftw_real tmp264;
		    fftw_real tmp265;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp260 = tmp248 - tmp249;
		    tmp261 = tmp241 - tmp244;
		    tmp262 = tmp260 + tmp261;
		    tmp268 = tmp260 - tmp261;
		    tmp264 = tmp236 + tmp237;
		    tmp265 = tmp251 + tmp252;
		    tmp266 = tmp264 - tmp265;
		    tmp270 = tmp264 + tmp265;
	       }
	       {
		    fftw_real tmp259;
		    fftw_real tmp263;
		    fftw_real tmp267;
		    fftw_real tmp269;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp259 = c_re(W[6]);
		    tmp263 = c_im(W[6]);
		    Y[-8 * iostride] = (tmp259 * tmp262) - (tmp263 * tmp266);
		    X[7 * iostride] = (tmp263 * tmp262) + (tmp259 * tmp266);
		    tmp267 = c_re(W[14]);
		    tmp269 = c_im(W[14]);
		    Y[0] = (tmp267 * tmp268) - (tmp269 * tmp270);
		    X[15 * iostride] = (tmp269 * tmp268) + (tmp267 * tmp270);
	       }
	  }
     }
     if (i == m) {
	  fftw_real tmp7;
	  fftw_real tmp51;
	  fftw_real tmp19;
	  fftw_real tmp43;
	  fftw_real tmp39;
	  fftw_real tmp47;
	  fftw_real tmp59;
	  fftw_real tmp64;
	  fftw_real tmp14;
	  fftw_real tmp56;
	  fftw_real tmp24;
	  fftw_real tmp32;
	  fftw_real tmp29;
	  fftw_real tmp33;
	  fftw_real tmp54;
	  fftw_real tmp65;
	  fftw_real tmp63;
	  fftw_real tmp66;
	  ASSERT_ALIGNED_DOUBLE;
	  {
	       fftw_real tmp3;
	       fftw_real tmp15;
	       fftw_real tmp38;
	       fftw_real tmp57;
	       fftw_real tmp6;
	       fftw_real tmp35;
	       fftw_real tmp18;
	       fftw_real tmp58;
	       ASSERT_ALIGNED_DOUBLE;
	       {
		    fftw_real tmp1;
		    fftw_real tmp2;
		    fftw_real tmp36;
		    fftw_real tmp37;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp1 = X[0];
		    tmp2 = X[7 * iostride];
		    tmp3 = tmp1 + tmp2;
		    tmp15 = tmp1 - tmp2;
		    tmp36 = Y[0];
		    tmp37 = Y[-7 * iostride];
		    tmp38 = tmp36 + tmp37;
		    tmp57 = tmp36 - tmp37;
	       }
	       {
		    fftw_real tmp4;
		    fftw_real tmp5;
		    fftw_real tmp16;
		    fftw_real tmp17;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp4 = X[4 * iostride];
		    tmp5 = X[3 * iostride];
		    tmp6 = tmp4 + tmp5;
		    tmp35 = tmp4 - tmp5;
		    tmp16 = Y[-4 * iostride];
		    tmp17 = Y[-3 * iostride];
		    tmp18 = tmp16 + tmp17;
		    tmp58 = tmp16 - tmp17;
	       }
	       tmp7 = tmp3 + tmp6;
	       tmp51 = tmp3 - tmp6;
	       tmp19 = tmp15 - tmp18;
	       tmp43 = tmp15 + tmp18;
	       tmp39 = tmp35 + tmp38;
	       tmp47 = tmp38 - tmp35;
	       tmp59 = tmp57 - tmp58;
	       tmp64 = tmp58 + tmp57;
	  }
	  {
	       fftw_real tmp10;
	       fftw_real tmp20;
	       fftw_real tmp23;
	       fftw_real tmp53;
	       fftw_real tmp13;
	       fftw_real tmp25;
	       fftw_real tmp28;
	       fftw_real tmp52;
	       ASSERT_ALIGNED_DOUBLE;
	       {
		    fftw_real tmp8;
		    fftw_real tmp9;
		    fftw_real tmp21;
		    fftw_real tmp22;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp8 = X[2 * iostride];
		    tmp9 = X[5 * iostride];
		    tmp10 = tmp8 + tmp9;
		    tmp20 = tmp8 - tmp9;
		    tmp21 = Y[-2 * iostride];
		    tmp22 = Y[-5 * iostride];
		    tmp23 = tmp21 + tmp22;
		    tmp53 = tmp21 - tmp22;
	       }
	       {
		    fftw_real tmp11;
		    fftw_real tmp12;
		    fftw_real tmp26;
		    fftw_real tmp27;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp11 = X[iostride];
		    tmp12 = X[6 * iostride];
		    tmp13 = tmp11 + tmp12;
		    tmp25 = tmp11 - tmp12;
		    tmp26 = Y[-iostride];
		    tmp27 = Y[-6 * iostride];
		    tmp28 = tmp26 + tmp27;
		    tmp52 = tmp27 - tmp26;
	       }
	       tmp14 = tmp10 + tmp13;
	       tmp56 = tmp10 - tmp13;
	       tmp24 = tmp20 - tmp23;
	       tmp32 = tmp20 + tmp23;
	       tmp29 = tmp25 - tmp28;
	       tmp33 = tmp25 + tmp28;
	       tmp54 = tmp52 - tmp53;
	       tmp65 = tmp53 + tmp52;
	  }
	  X[0] = K2_000000000 * (tmp7 + tmp14);
	  X[8 * iostride] = -(K2_000000000 * (tmp65 + tmp64));
	  tmp63 = tmp7 - tmp14;
	  tmp66 = tmp64 - tmp65;
	  X[4 * iostride] = K1_414213562 * (tmp63 - tmp66);
	  X[12 * iostride] = -(K1_414213562 * (tmp63 + tmp66));
	  {
	       fftw_real tmp61;
	       fftw_real tmp62;
	       fftw_real tmp55;
	       fftw_real tmp60;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp61 = tmp51 - tmp54;
	       tmp62 = tmp59 - tmp56;
	       X[6 * iostride] = (K765366864 * tmp61) - (K1_847759065 * tmp62);
	       X[14 * iostride] = -((K1_847759065 * tmp61) + (K765366864 * tmp62));
	       tmp55 = tmp51 + tmp54;
	       tmp60 = tmp56 + tmp59;
	       X[2 * iostride] = (K1_847759065 * tmp55) - (K765366864 * tmp60);
	       X[10 * iostride] = -((K765366864 * tmp55) + (K1_847759065 * tmp60));
	  }
	  {
	       fftw_real tmp45;
	       fftw_real tmp49;
	       fftw_real tmp48;
	       fftw_real tmp50;
	       fftw_real tmp44;
	       fftw_real tmp46;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp44 = K707106781 * (tmp32 + tmp33);
	       tmp45 = tmp43 - tmp44;
	       tmp49 = tmp43 + tmp44;
	       tmp46 = K707106781 * (tmp24 - tmp29);
	       tmp48 = tmp46 + tmp47;
	       tmp50 = tmp47 - tmp46;
	       X[3 * iostride] = (K1_662939224 * tmp45) - (K1_111140466 * tmp48);
	       X[11 * iostride] = -((K1_111140466 * tmp45) + (K1_662939224 * tmp48));
	       X[7 * iostride] = (K390180644 * tmp49) - (K1_961570560 * tmp50);
	       X[15 * iostride] = -((K1_961570560 * tmp49) + (K390180644 * tmp50));
	  }
	  {
	       fftw_real tmp31;
	       fftw_real tmp41;
	       fftw_real tmp40;
	       fftw_real tmp42;
	       fftw_real tmp30;
	       fftw_real tmp34;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp30 = K707106781 * (tmp24 + tmp29);
	       tmp31 = tmp19 + tmp30;
	       tmp41 = tmp19 - tmp30;
	       tmp34 = K707106781 * (tmp32 - tmp33);
	       tmp40 = tmp34 + tmp39;
	       tmp42 = tmp39 - tmp34;
	       X[iostride] = (K1_961570560 * tmp31) - (K390180644 * tmp40);
	       X[9 * iostride] = -((K390180644 * tmp31) + (K1_961570560 * tmp40));
	       X[5 * iostride] = (K1_111140466 * tmp41) - (K1_662939224 * tmp42);
	       X[13 * iostride] = -((K1_662939224 * tmp41) + (K1_111140466 * tmp42));
	  }
     }
}

static const int twiddle_order[] =
{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15};
fftw_codelet_desc fftw_hc2hc_backward_16_desc =
{
     "fftw_hc2hc_backward_16",
     (void (*)()) fftw_hc2hc_backward_16,
     16,
     FFTW_BACKWARD,
     FFTW_HC2HC,
     366,
     15,
     twiddle_order,
};

ext: c

max_stars_repo_name: albertsgrc/ftdock-opt

max_issues_repo_name: albertsgrc/ftdock-opt

max_forks_repo_name: albertsgrc/ftdock-opt


text:
#pragma once

#include <initializer_list>
#include <list>
#include <vector>

#include <gsl/gsl_complex.h>

#include "Expression.h"

class TupleExpression: public Expression {

    std::vector<expression> data;

    TupleExpression();
    TupleExpression(std::vector<expression>&& tuple);
    TupleExpression(std::list<expression>& tuple);
    TupleExpression(std::initializer_list<expression> tuple);
    TupleExpression(std::initializer_list<double> tuple);
    TupleExpression(std::initializer_list<gsl_complex> tuple);

    public:
        static expression construct();
        static expression construct(std::vector<expression>&& tuple);
        static expression construct(std::list<expression>& tuple);
        static expression construct(std::initializer_list<expression> tuple);
        static expression construct(std::initializer_list<double> tuple);
        static expression construct(std::initializer_list<gsl_complex> tuple);

        expression simplify() override;
        expression derivative(const std::string& var) override;
        expression integrate(const std::string& var) override;

        expression at(const int index) override;
        size_t size() const override;

        expression apply(TransformerFunction f) override;

        EXPRESSION_OVERRIDES

};

ext: h

max_stars_repo_name: antoniojkim/CalcPlusPlus

max_issues_repo_name: antoniojkim/CalcPlusPlus

max_forks_repo_name: antoniojkim/CalcPlusPlus


text:
// @Copyright 2007 Kristjan Haule
//
#include <cstdio>
#include <cstdlib>
#include <gsl/gsl_rng.h>
#include <map>

class RanGSL{
  const gsl_rng_type *T;
  gsl_rng *r;
public:
  RanGSL(int idum)
  {
    gsl_rng_env_setup();
    //    T = gsl_rng_default;
    //    T = gsl_rng_taus;
    T = gsl_rng_ranlux389;
    r = gsl_rng_alloc (T);
    gsl_rng_set(r, idum);
  }
  double operator()()
  { return gsl_rng_uniform (r);}
  ~RanGSL()
  {gsl_rng_free (r);}
};

class dRand48{
public:
  dRand48(int idum){srand48(idum);}
  double operator()(){
    return drand48();
  }
};

class Ran0{
  long idum;
public:
  Ran0(long idum_): idum(idum_){};
  double operator()(){
    // Minimal random number generator of Park and Miller. Returns a uniform random deviate
    // between 0.0 and 1.0. Set or reset idum to any integer value (except the unlikely value MASK)
    // to initialize the sequence; idum must not be altered between calls for successive deviates in
    // a sequence.
    static const int IA=16807;
    static const int IM=2147483647;
    static const double AM=(1.0/IM);
    static const int IQ=127773;
    static const int IR=2836;
    static const int MASK=123459876;
    static long idum;
    long k;
    float ans;
    idum ^= MASK; //XORing with MASK allows use of zero and other simple bit patterns for idum.
    k=idum/IQ;
    idum = IA*(idum-k*IQ)-IR*k; // Compute idum=(IA*idum) % IM without over
    if (idum < 0) idum += IM; // flows by Schrages method.
    ans=AM*idum; // Convert idum to a floating result.
    idum ^= MASK; // Unmask before return.
    return ans;
  }
};


text:
#ifndef libceed_solids_examples_setup_dm_h
#define libceed_solids_examples_setup_dm_h

#include <ceed.h>
#include <petsc.h>
#include <petscdmplex.h>
#include <petscfe.h>
#include "../include/structs.h"

// -----------------------------------------------------------------------------
// Setup DM
// -----------------------------------------------------------------------------
PetscErrorCode CreateBCLabel(DM dm, const char name[]);

// Read mesh and distribute DM in parallel
PetscErrorCode CreateDistributedDM(MPI_Comm comm, AppCtx app_ctx, DM *dm);

// Setup DM with FE space of appropriate degree
PetscErrorCode SetupDMByDegree(DM dm, AppCtx app_ctx, PetscInt order,
                               PetscBool boundary, PetscInt num_comp_u);

#endif // libceed_solids_examples_setup_dm_h

max_stars_repo_name: wence-/libCEED

max_issues_repo_name: wence-/libCEED

max_forks_repo_name: wence-/libCEED


text:
/**
 *
 * @file dtile.c
 *
 *  PLASMA auxiliary routines
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Jakub Kurzak
 * @author Mathieu Faverge
 * @date 2010-11-15
 * @generated d Tue Jan  7 11:45:15 2014
 *
 **/
#include "common.h"
#ifdef USE_MKL
#include <mkl_lapacke.h>
#else
#include <lapacke.h>
#endif
#include "auxiliary.h"
#include "tile.h"
/*#include "quark.h"*/

#define AF77(m, n) &(Af77[ ((int64_t)A.nb*(int64_t)lda*(int64_t)(n)) + (int64_t)(A.mb*(m)) ])
#define ABDL(m, n) BLKADDR(A, double, m, n)

/***************************************************************************//**
 *  Conversion from LAPACK F77 matrix layout to tile layout - dynamic scheduling
 **/
void plasma_pdtile_to_lapack_quark(PLASMA_desc A, double *Af77, int lda)
{
    double *f77;
    double *bdl;
    int X1, Y1;
    int X2, Y2;
    int n, m, ldt;
    for (m = 0; m < A.mt; m++)
    {
        ldt = BLKLDD(A, m);
        for (n = 0; n < A.nt; n++)
        {
            X1 = n == 0 ? A.j%A.nb : 0;
            Y1 = m == 0 ? A.i%A.mb : 0;
            X2 = n == A.nt-1 ? (A.j+A.n-1)%A.nb+1 : A.nb;
            Y2 = m == A.mt-1 ? (A.i+A.m-1)%A.mb+1 : A.mb;

            f77 = AF77(m, n);
            bdl = ABDL(m, n);
            LAPACKE_dlacpy_work(
                LAPACK_COL_MAJOR,
                lapack_const(PlasmaUpperLower),
                (Y2-Y1), (X2-X1),
                &(bdl[X1*lda+Y1]), ldt,
                &(f77[X1*lda+Y1]), lda);
        }
    }
}


text:
#pragma once

#include "Cesium3DTilesReader/Library.h"

#include <Cesium3DTiles/Subtree.h>
#include <CesiumJsonReader/ExtensionReaderContext.h>

#include <gsl/span>

#include <functional>
#include <memory>
#include <optional>
#include <string>
#include <unordered_map>
#include <vector>

namespace Cesium3DTilesReader {

/**
 * @brief The result of reading a subtree with
 * {@link SubtreeReader::readSubtree}.
 */
struct CESIUM3DTILESREADER_API SubtreeReaderResult {
  /**
   * @brief The read subtree, or std::nullopt if the subtree could not be read.
   */
  std::optional<Cesium3DTiles::Subtree> subtree;

  /**
   * @brief Errors, if any, that occurred during the load process.
   */
  std::vector<std::string> errors;

  /**
   * @brief Warnings, if any, that occurred during the load process.
   */
  std::vector<std::string> warnings;
};

/**
 * @brief Reads subtrees.
 */
class CESIUM3DTILESREADER_API SubtreeReader {
public:
  /**
   * @brief Constructs a new instance.
   */
  SubtreeReader();

  /**
   * @brief Gets the context used to control how extensions are loaded from a
   * subtree.
   */
  CesiumJsonReader::ExtensionReaderContext& getExtensions();

  /**
   * @brief Gets the context used to control how extensions are loaded from a
   * subtree.
   */
  const CesiumJsonReader::ExtensionReaderContext& getExtensions() const;

  /**
   * @brief Reads a subtree.
   *
   * @param data The buffer from which to read the subtree.
   * @param options Options for how to read the subtree.
   * @return The result of reading the subtree.
   */
  SubtreeReaderResult readSubtree(const gsl::span<const std::byte>& data) const;

private:
  CesiumJsonReader::ExtensionReaderContext _context;
};

} // namespace Cesium3DTilesReader

max_stars_repo_name: yieryi/cesium-native

max_issues_repo_name: yieryi/cesium-native

max_forks_repo_name: yieryi/cesium-native


text:
// Copyright (c) Microsoft Corporation.
// Licensed under the MIT License.

//------------------------------------------------------------------------------
// KeypointSpatialIndex.h
//
// Allows to construction, query, and removal of Keypoint entries based on pixel position.
//------------------------------------------------------------------------------

#pragma once

#include "Utils\thread_memory.h"

#include <memory>
#include <vector>
#include <functional>
#include <opencv2\features2d\features2d.hpp>

#include <gsl/span>

namespace mage
{
    class KeypointSpatialIndex
    {
        struct Impl;

    public:
        KeypointSpatialIndex(gsl::span<const cv::KeyPoint> keypoints);

        void Query(const cv::Point2f& center, int octave, float radius, temp::vector<size_t>& results) const;
        void Remove(const cv::KeyPoint& keypoint, size_t value);

        ~KeypointSpatialIndex();

    private:
        static constexpr float octaveSpacing = 100;
        static constexpr float octaveQueryRange = 1;

        std::unique_ptr<Impl> m_impl;
    };
}

max_stars_repo_name: syntheticmagus/mageslam

max_issues_repo_name: syntheticmagus/mageslam

max_forks_repo_name: syntheticmagus/mageslam


text:
// Alessandro Casalino
//
// Compile with (Mac with default homebrew gsl 2.6) gcc-9 -O2 quintessence_evolve.c -o quintessence_evolve.exe -L/usr/local/Cellar/gsl/2.6/lib -I/usr/local/Cellar/gsl/2.6/include -lgsl
// Run with ./quintessence_evolve.exe

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>

// PHYSICAL PARAMETERS VALUES

// INITIAL CONDITIONS
// Initial value of the scale factor
double a_init = 1e-15;
// Final value of the scale factor
double a_end = 1.1;

// Values of fractional density for the cosmological matter today
double Omega_rad_0 = 8e-5;
double Omega_b_0 = 0.0486;
double Omega_Lambda_0 = 0.;//0.6911;
double Omega_cdm_0 = 0.2589;

// Physical constants for output conversions
double _G_ = 6.67428e-11;                 /**< Newton constant in m^3/Kg/s^2 */
double _MPc_over_m_ = 3.085677581282e22;  // Conversion factor from Mpc to m
double _Gyr_over_Mpc_ = 3.06601394e2;     // Conversion factor from Gyr to Mpc
double _c_ = 2.99792458e8;                // Speed of light in m/s
double _H0_ = 67.74;                      // H0 of LCDM in Km/s/Mpc

double fourpiG = 1.;

double T_CONF = 1.;

// COMPUTATION PARAMETERS VALUES

// Number of points used in the computation
int points = (int) 1e6;

// Raise this value to make the csv file smaller, but decreasing resolution
// The value inserted is the ratio between the number of values written in a full resolution file / decreased resolution file
int csv_resolution = 10;

double delta_bisection = 1e-2;
double mg_field_init_min = 1e-20;
double mg_field_init_max = 2e1; //1e16 max for LCDM (model 1)

// QUINTESSENCE PARAMETERS
// Initial value of the quintessence velocity (phi derivated with respect to tau)
double mg_field_p_bk_0 = 0.;

// TEST MODE
// Provides some informations on the terminal and the .csv file during the computation (1: on, others: off)
int TEST_MODE = 1;

// Definition of the POTENTIAL
// For a list of the models see above
double mg_pot(const double mg_field_bk) {

  return mg_field_bk * mg_field_bk / 2.;

}

double mg_pot_p(const double mg_field_bk) {

  return mg_field_bk;

}

// Definitions of DIFFERENTIAL EQUATION system
// Need to have first order equations to use Runge Kutta
//
double a_p_rk4(const double a, const double a_p, const double mg_field_bk, const double mg_field_p_bk)
{
	return a_p;
}
double a_pp_rk4(const double a, const double a_p, const double mg_field_bk, const double mg_field_p_bk)
{
  double a3 = a * a * a;
  double rhoa3 = (Omega_cdm_0 + Omega_b_0) + (Omega_Lambda_0 * a3) + (Omega_rad_0 / a);
  double Pa3 = - (Omega_Lambda_0 * a3) + 1./3. * (Omega_rad_0 / a );

  return fourpiG / 3. * ( (rhoa3 - 3. * Pa3) - a * mg_field_p_bk * mg_field_p_bk + 4. * a3 * mg_pot(mg_field_bk) );
}
double mg_field_bk_p_rk4(const double a, const double a_p, const double mg_field_bk, const double mg_field_p_bk)
{
	return mg_field_p_bk;
}
double mg_field_bk_pp_rk4(const double a, const double a_p, const double mg_field_bk, const double mg_field_p_bk)
{
	return - 2. * a_p/a * mg_field_p_bk - mg_pot_p(mg_field_bk) * a * a;
}
double Hconf(const double a, const double mg_field_bk, const double mg_field_p_bk)
{
	return sqrt((2. * fourpiG / 3.) * ( ((Omega_cdm_0 + Omega_b_0) / a) + (Omega_Lambda_0 * a * a) + (Omega_rad_0 / a / a) + (mg_field_p_bk * mg_field_p_bk / 2.) + (a * a * mg_pot(mg_field_bk)) ));
}

// Integrand for the particle horizon integral
double particleHorizonIntegrand(double a, double mg_field_bk, double mg_field_p_bk)
{
  //return 2. / (sqrt(a) * Hconf(a, mg_field_bk, mg_field_p_bk));
	return 1. / ( a * Hconf(a, mg_field_bk, mg_field_p_bk) );
}
// Particle horizon integral step
double particleHorizon(const int i, double * a, double * mg_field_bk, double * mg_field_p_bk) {
  double h = a[i]-a[i-1];
  double fa = particleHorizonIntegrand(a[i-1],mg_field_bk[i-1],mg_field_p_bk[i-1]);
  double fb = particleHorizonIntegrand(a[i],mg_field_bk[i],mg_field_p_bk[i]);

  return h*(fa+fb)/2.;
}

// Function used to print results stored in vectors as a csv file
void csv(double * t, double * a, double * a_p, double * mg_field_bk, double * mg_field_p_bk, double * particleHorizonVec, char * filename) {

    FILE *fp;
    fp = fopen (filename, "w+");
    fprintf(fp, "%s, %s, %s, %s, %s, %s, %s, %s", "t", "a(t)", "H(t)/H0", "H_prime(t)/H0^2", "Omega_df", "Omega_r", "Omega_b", "Omega_cdm");
    fprintf(fp, ", %s", "omega_df");
    fprintf(fp, ", %s", "PH"); // Particle horizon
    fprintf(fp, ", %s", "mg_field");
    fprintf(fp, ", %s", "mg_field_p");
    if(TEST_MODE==1) fprintf(fp, ", %s", "H_check(t)/H0");
    fprintf(fp, "\n");

    // Time conversion factor
    // double tcf = 1./_H0_/(60.*60.*24.*365.*1e9)*_MPc_over_m_/1000.;

    int i = 0;
    double particleHorizonRes = 0.;

    for(i=1;i<points;i++){
      particleHorizonRes += particleHorizon(i,a,mg_field_bk,mg_field_p_bk);
      particleHorizonVec[i] = particleHorizonRes;
    }

    i = 0;

    double H_test, H_cosmo = 0.;

    while( a[i] <= a_end){

      double H = Hconf(a[i],mg_field_bk[i],mg_field_p_bk[i]); // this is the Hubble constant with conformal time
      double H_cosmo = H / a[i];
      double H_prime = 1. / a[i] * a_pp_rk4(a[i], a_p[i], mg_field_bk[i], mg_field_p_bk[i]) - H * H; // this is Hubble prime with conformal time
      if(T_CONF==0) H_prime = ( H_prime - H * H )/a[i]/a[i]; // this is Hubble prime with cosmological time
      if(T_CONF==0) H = H_cosmo; // this it the Hubble constant with cosmological time

      fprintf(fp, "%e, %e, %e, %e, %e, %e, %e, %e", t[i], a[i], H / sqrt(2. * fourpiG / 3.), H_prime / sqrt(2. * fourpiG / 3.) / sqrt(2. * fourpiG / 3.), 2. * fourpiG / 3. * ( (mg_field_p_bk[i] * mg_field_p_bk[i] / 2.) + (a[i] * a[i] * mg_pot(mg_field_bk[i])) ) /H_cosmo /H_cosmo , 2. * fourpiG / 3. * Omega_rad_0 /a[i] /a[i] /a[i] /a[i] /H_cosmo  / H_cosmo , 2. * fourpiG / 3. * Omega_b_0 /a[i] /a[i] /a[i] /H_cosmo  / H_cosmo , 2. * fourpiG / 3. * Omega_cdm_0 /a[i] /a[i] /a[i] /H_cosmo  / H_cosmo );
      fprintf(fp, ", %e", (-mg_pot(mg_field_bk[i])+mg_field_p_bk[i]*mg_field_p_bk[i]/2./a[i]/a[i])/(mg_pot(mg_field_bk[i])+mg_field_p_bk[i]*mg_field_p_bk[i]/2./a[i]/a[i]));

      fprintf(fp, ", %e", particleHorizonVec[i]);

      fprintf(fp, ", %e", mg_field_bk[i]);
      fprintf(fp, ", %e", mg_field_p_bk[i]);

      if(TEST_MODE==1){

        H_test = a_p[i]/a[i];
        if(T_CONF==0) H_test = H_test / a[i];

        fprintf(fp, ", %e", H_test / sqrt(2. * fourpiG / 3.));

      }

      fprintf(fp, "\n");

      if(a[i+csv_resolution]<= a_end){
        i=i+csv_resolution;
      }
      else{
        i++;
      }

    }

    fclose(fp);

}

int scan_for_a0 (double * a) {

  int i=0;

  for(i=0; a[i] <= a_end; i++);

  return i;

}

// This divides an interval in a logarithmic scale of basis 10, and store results in input vector v
void logscale10 (double * v, double A, double B, int points){

    int i;

    double a = log10(A);
    double b = log10(B);

    double h = (b - a) / (points-1.0);

    for(i=0;i<points;i++){
        v[i] = pow(10.0, a + i * h);
    }

}

// This function is a temporal step of the Runge-Kutta 4 method:
// evolves the system for a time equal to dtau (h in the program)
void mg_rungekutta4bg(double * f, const double dtau)
{

  double k_F[4], k_A[4], k_rhof[4], k_rhor[4], k_rhob[4];
  double a = f[1];
  double a_p = f[2];
  double mg_field_bk = f[3];
  double mg_field_p_bk = f[4];

	double k1a, k2a, k3a, k4a;
	double k1ap, k2ap, k3ap, k4ap;
	double k1f, k2f, k3f, k4f;
	double k1fp, k2fp, k3fp, k4fp;

	k1a 	= a_p_rk4(a, a_p, mg_field_bk, mg_field_p_bk);
	k1ap 	= a_pp_rk4(a, a_p, mg_field_bk, mg_field_p_bk);
	k1f 	= mg_field_bk_p_rk4(a, a_p, mg_field_bk, mg_field_p_bk);
	k1fp 	= mg_field_bk_pp_rk4(a, a_p, mg_field_bk, mg_field_p_bk);

	k2a 	= a_p_rk4(a + k1a * dtau / 2., a_p + k1ap * dtau / 2., mg_field_bk + k1f * dtau / 2., mg_field_p_bk + k1fp * dtau / 2.);
	k2ap 	= a_pp_rk4(a + k1a * dtau / 2., a_p + k1ap * dtau / 2., mg_field_bk + k1f * dtau / 2., mg_field_p_bk + k1fp * dtau / 2.);
	k2f 	= mg_field_bk_p_rk4(a + k1a * dtau / 2., a_p + k1ap * dtau / 2., mg_field_bk + k1f * dtau / 2., mg_field_p_bk + k1fp * dtau / 2.);
	k2fp 	= mg_field_bk_pp_rk4(a + k1a * dtau / 2., a_p + k1ap * dtau / 2., mg_field_bk + k1f * dtau / 2., mg_field_p_bk + k1fp * dtau / 2.);

	k3a 	= a_p_rk4(a + k2a * dtau / 2., a_p + k2ap * dtau / 2., mg_field_bk + k2f * dtau / 2., mg_field_p_bk + k2fp * dtau / 2.);
	k3ap 	= a_pp_rk4(a + k2a * dtau / 2., a_p + k2ap * dtau / 2., mg_field_bk + k2f * dtau / 2., mg_field_p_bk + k2fp * dtau / 2.);
	k3f 	= mg_field_bk_p_rk4(a + k2a * dtau / 2., a_p + k2ap * dtau / 2., mg_field_bk + k2f * dtau / 2., mg_field_p_bk + k2fp * dtau / 2.);
	k3fp 	= mg_field_bk_pp_rk4(a + k2a * dtau / 2., a_p + k2ap * dtau / 2., mg_field_bk + k2f * dtau / 2., mg_field_p_bk + k2fp * dtau / 2.);

	k4a 	= a_p_rk4(a + k3a * dtau, a_p + k3ap * dtau, mg_field_bk + k3f * dtau, mg_field_p_bk + k3fp * dtau);
	k4ap 	= a_pp_rk4(a + k3a * dtau, a_p + k3ap * dtau, mg_field_bk + k3f * dtau, mg_field_p_bk + k3fp * dtau);
	k4f 	= mg_field_bk_p_rk4(a + k3a * dtau, a_p + k3ap * dtau, mg_field_bk + k3f * dtau, mg_field_p_bk + k3fp * dtau);
	k4fp 	= mg_field_bk_pp_rk4(a + k3a * dtau, a_p + k3ap * dtau, mg_field_bk + k3f * dtau, mg_field_p_bk + k3fp * dtau);

	f[1] += dtau * (k1a + 2. * k2a + 2. * k3a + k4a) / 6.;
	f[2] += dtau * (k1ap + 2. * k2ap + 2. * k3ap + k4ap) / 6.;
	f[3] += dtau * (k1f + 2. * k2f + 2. * k3f + k4f) / 6.;
	f[4] += dtau * (k1fp + 2. * k2fp + 2. * k3fp + k4fp) / 6.;

}

// This function evolves the system with the Runge-Kutta 4 method until t_stop
double rk4(double * t, double * a, double * a_p, double * mg_field_bk, double * mg_field_p_bk, double Omega_f_0) {

    double f[5];
    int i = 0;

    // Call the function for a logarithmic scale of time
    // We don't start from t = 0 to avoid problems with quintessence potentials
    logscale10(t,1e-9,20.,points);

    while( i < points - 1 ){

      f[0] = t[i];
      f[1] = a[i];
      f[2] = a_p[i];
      f[3] = mg_field_bk[i];
      f[4] = mg_field_p_bk[i];

      mg_rungekutta4bg(f, t[i+1]-t[i]);

      a[i+1] = f[1];
      a_p[i+1] = f[2];
      mg_field_bk[i+1] = f[3];
      mg_field_p_bk[i+1] = f[4];
      i++;

    }

    int j = scan_for_a0(a);

    double H        = Hconf(a[j], mg_field_bk[j], mg_field_p_bk[j]);
    double Omega_f  = (2. * fourpiG / 3.) * ( (mg_field_p_bk[j] * mg_field_p_bk[j] / 2.) + (a[j] * a[j] * mg_pot(mg_field_bk[j])) ) /H /H;
    return Omega_f  - Omega_f_0;

}

double bisection (double min, double max, double * t, double * a, double * a_p, double * mg_field_bk, double * mg_field_p_bk, const double Omega_f_0) {

  double C = (min+max) / 2.;

  while(fabs((max-min)/min)>delta_bisection){

    a_p[0] = a_init * Hconf(a[0], mg_field_init_min, mg_field_p_bk[0]);
    mg_field_bk[0] = min;
    double rk4_min = rk4(t, a, a_p, mg_field_bk, mg_field_p_bk, Omega_f_0);
    a_p[0] = a_init * Hconf(a[0], C, mg_field_p_bk[0]);
    mg_field_bk[0] = C;
    double rk4_C = rk4(t, a, a_p, mg_field_bk, mg_field_p_bk, Omega_f_0);

    if(rk4_min*rk4_C>=0) {
      min=C;
    }
    else {
      max=C;
    }

    C=(max+min)/2.;

    if (TEST_MODE == 1) printf("TEST_MODE ON - min: %e , max: %e, C: %e, rk4_min: %e , rk4_C: %e \n", min, max, C, rk4_min, rk4_C);

  }

  double result = (max+min)/2.;

  if(TEST_MODE==1) printf("\n");
  printf("\t-> Result of bisection method is mg_field_bk: %e (internal units).\n", result);
  //if(TEST_MODE==1) printf("\t--> Confront with LCDM value: %e (internal units).\n", (2. * fourpiG / 3.) * ( Omega_Lambda_0 + 6. * OmegaCDM_0 )/ pow(a_init,3.));

  return result;

}

int main() {

    double Omega_f_0 = 1. - Omega_Lambda_0 - Omega_rad_0 - Omega_b_0 - Omega_cdm_0;

    int i = 0;

    printf("\n\t\t----------------------------------------\n\n");

    // Definition of the vector needed for the evolution functions
    double * t; double * a; double * a_p; double * mg_field_bk; double * mg_field_p_bk; double * particleHorizonVec;
    t                   = (double *) malloc(sizeof(double) * points);
    a                   = (double *) malloc(sizeof(double) * points);
    a_p                 = (double *) malloc(sizeof(double) * points);
    mg_field_bk         = (double *) malloc(sizeof(double) * points);
    mg_field_p_bk       = (double *) malloc(sizeof(double) * points);
    particleHorizonVec  = (double *) malloc(sizeof(double) * points);

    if(!t||!a||!a_p||!mg_field_bk||!mg_field_p_bk||!particleHorizonVec){
      printf("Error! The memory cannot be allocated. The program will be terminated.\n");
      exit(1);
    }

    // Initial conditions (tau=0)
    a[0] = a_init;
    mg_field_p_bk[0] = mg_field_p_bk_0;

    printf(" Searching for best initial value for the dark fluid ... \n \n");
    mg_field_bk[0] = bisection(mg_field_init_min, mg_field_init_max, t, a, a_p, mg_field_bk, mg_field_p_bk, Omega_f_0);
    a_p[0] = a_init * Hconf(a_init, mg_field_bk[0], mg_field_p_bk[0]);

    printf("\n\n Evolving the system ...\n");

    rk4(t, a, a_p, mg_field_bk, mg_field_p_bk, Omega_f_0);

    char filename[50];
    sprintf (filename, "mg_bk.csv");
    int last_int = scan_for_a0(a);

    printf("\n RESULTS:\n");
    printf("\t-> H0: %f \n", a_p[last_int]/a[last_int]/sqrt(2. * fourpiG / 3.));
    printf("\t-> number of points %d \n", last_int);
    //printf("\t-> Age of the Universe: %f Gyr\n", t[last_int] /_H0_/(60.*60.*24.*365.*1e9)*_MPc_over_m_/1000.);

    csv(t, a, a_p, mg_field_bk, mg_field_p_bk, particleHorizonVec, filename);

    printf("\n The results are saved in '.csv' files. The name is labelled with the value of c, and the model (m).\n");

    if(TEST_MODE==1) printf("\n TEST_MODE ON: check the values of H in .csv file. They must be equal!");

    printf("\n\t\t----------------------------------------\n");

    double * a_int; double * a_p_int; double * mg_field_bk_int; double * mg_field_p_bk_int; double * particleHorizonVec_int;
    a_int                   = (double *) malloc(sizeof(double) * last_int);
    a_p_int                 = (double *) malloc(sizeof(double) * last_int);
    mg_field_bk_int         = (double *) malloc(sizeof(double) * last_int);
    mg_field_p_bk_int       = (double *) malloc(sizeof(double) * last_int);
    particleHorizonVec_int  = (double *) malloc(sizeof(double) * last_int);

    if(!a_int||!a_p_int||!mg_field_bk_int||!mg_field_p_bk_int||!particleHorizonVec_int){
      printf("Error! The memory cannot be allocated. The program will be terminated.\n");
      exit(1);
    }

    memcpy(a_int, a, last_int * sizeof(double));
    memcpy(a_p_int, a_p, last_int * sizeof(double));
    memcpy(mg_field_bk_int, mg_field_bk, last_int * sizeof(double));
    memcpy(mg_field_p_bk_int, mg_field_p_bk, last_int * sizeof(double));
    memcpy(particleHorizonVec_int, particleHorizonVec, last_int * sizeof(double));

    free(a);free(a_p);free(mg_field_bk);free(mg_field_p_bk);free(particleHorizonVec);

    // Spline interpolation with gsl
    gsl_interp_accel *acc_mg_field = gsl_interp_accel_alloc();
    gsl_spline * spline_mg_field = gsl_spline_alloc(gsl_interp_cspline,last_int);
    gsl_interp_accel *acc_mg_field_p = gsl_interp_accel_alloc();
    gsl_spline * spline_mg_field_p = gsl_spline_alloc(gsl_interp_cspline,last_int);
    gsl_interp_accel *acc_a_p = gsl_interp_accel_alloc();
    gsl_spline * spline_a_p = gsl_spline_alloc(gsl_interp_cspline,last_int);
    gsl_interp_accel *acc_particleHorizon = gsl_interp_accel_alloc();
    gsl_spline * spline_particleHorizon = gsl_spline_alloc(gsl_interp_cspline,last_int);

    gsl_spline_init(spline_mg_field,a_int,mg_field_bk_int,last_int);
    gsl_spline_init(spline_mg_field_p,a_int,mg_field_p_bk_int,last_int);
    gsl_spline_init(spline_a_p,a_int,a_p_int,last_int);
    gsl_spline_init(spline_particleHorizon,a_int,particleHorizonVec_int,last_int);

    double a_eval = 1e-2;
    printf("spline eval: %e %e \n",a_eval,gsl_spline_eval(spline_mg_field,a_eval,acc_mg_field));
    printf("spline eval: %e %e \n",a_eval,gsl_spline_eval(spline_mg_field_p,a_eval,acc_mg_field_p));
    printf("spline eval: %e %e \n",a_eval,gsl_spline_eval(spline_a_p,a_eval,acc_a_p));
    printf("spline eval: %e %e \n",a_eval,gsl_spline_eval(spline_particleHorizon,a_eval,acc_particleHorizon));
    printf("%d %e \n",last_int, particleHorizonVec_int[gsl_interp_bsearch(a_int,a_eval,0,last_int-1)]);

    gsl_spline_free(spline_mg_field);gsl_interp_accel_free(acc_mg_field);
    gsl_spline_free(spline_mg_field_p);gsl_interp_accel_free(acc_mg_field_p);
    gsl_spline_free(spline_a_p);gsl_interp_accel_free(acc_a_p);
    gsl_spline_free(spline_particleHorizon);gsl_interp_accel_free(acc_particleHorizon);
    free(t);free(a_int);free(a_p_int);free(mg_field_bk_int);free(mg_field_p_bk_int);free(particleHorizonVec_int);

    printf("\n");

    exit(0);

}

max_stars_repo_name: alessandrocasalino/mgevolution

max_issues_repo_name: alessandrocasalino/mgevolution

max_forks_repo_name: alessandrocasalino/mgevolution


text:
// Copyright (c) Microsoft Corporation.
// Licensed under the MIT License.

#pragma once

#include <memory>
#include <vector>
#include <gsl/span>
#include <map>

#include <Eigen/Geometry>

namespace mage
{
    struct BundlerParameters
    {
        bool  ArePointsFixed{ false };        // True if map points should not be optimized
    };

    class BundlerLib
    {
    public:
        BundlerLib(const BundlerParameters& bundlerParameters);
        ~BundlerLib();

        void AllocateCameras(size_t count);

        void SetCameraPose(size_t idx,
            Eigen::Map<const Eigen::Vector3f> position,
            Eigen::Map<const Eigen::Matrix3f> orientation,
            Eigen::Map<const Eigen::Vector4f> intrinsics, bool isFixed);

        void FixCameraPose(size_t idx, bool value);

        void AllocateMapPoints(size_t count);
        void SetMapPoint(size_t idx, Eigen::Map<const Eigen::Vector3f> point);

        void AllocateObservations(size_t count);
        void SetObservation(size_t idx, Eigen::Map<const Eigen::Vector2f> position, size_t cameraIndex, size_t mapPointIndex, float informationMatrixScalar);

        void AllocateFixedDistanceConstraints(size_t count);
        void SetFixedDistanceConstraint(size_t idx, size_t cameraIndex1, size_t cameraIndex2, float distance = 1.0f, float weight = 1.0f);

        void AllocateRelativeRotationConstraints(size_t count);
        void SetRelativeRotationConstraint(size_t idx, size_t cameraIndex1, size_t cameraIndex2, const Eigen::Quaternionf& deltaRotation, float weight = 1.0f);

        void AllocateRelativeTransformConstraints(size_t count);
        void SetRelativeTransformConstraint(size_t idx, size_t cameraIndex1, size_t cameraIndex2, Eigen::Map<const Eigen::Vector3f> deltaPosition, const Eigen::Quaternionf& deltaRotation, float weight);

        void SetCurrentLambda(float userLambda);
        float GetCurrentLambda() const;

        // Runs an iteration of the solver for each provided Huber width.
        // Return the average square error.
        float StepBundleAdjustment(gsl::span<const float> huberWidthPerIteration, float maxErrorSquare, std::vector<unsigned int>& outliers);

        void GetPose(size_t idx, Eigen::Map<Eigen::Vector3f> position, Eigen::Map<Eigen::Matrix3f> orientation) const;
        void GetPoint(size_t idx, Eigen::Map<Eigen::Vector3f> position) const;

    private:
        struct Impl;
        const std::unique_ptr<Impl> m_impl;

        // Description of the problem to optimize from the calling code.
        BundlerParameters m_bundlerParameters;
    };
}