text:
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_sf_bessel.h>
#include "utility.h"
#include "constants.h"
#include <gsl/gsl_cblas.h>
#include <time.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <unistd.h>
#include <omp.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_sf_pow_int.h>

#ifdef EXTERN_VARIABLES
#define EXTERN
#else
#define EXTERN extern
#endif

#define OSSZI 14000
#define	SIMULATIONEN 10
#define KUPF_BINARY	1				// set to 1 for Kupferman/ fully mobile bath
#define VIRTUAL_FLAG	0			// set to 1 to transfer force one way onto virtual particle, no bath recouoling
#define SIMPSON_BINARY 1
#define REFLECTALL_FLAG 0
#define OMMEGA_READ_BINARY 0		// if true, don't regenerate ommega. use ommega.dat
#define DIM 1
#define ORDER ((2 * OSSZI + 2) * DIM)
#define LABELPREFIX "GIT_HUB"	//a string to add to your result folder
#define T_POWERS 0 							// no function, just compiler needs to know
#define T_MULT 1.0              			// same
#define LOGSCALE_BINARY 0       			// leave at 0, toggles time setup. WIP
#define TIME_STEPS  0.05					// snapshot intervall for saved values
#define LENGTH_T 	1000  +  9 				// Number N of saved steps (format must be N +9)
#define TIME_SHIFT 500
#define TIME_END (double)((double) (LENGTH_T-9) * TIME_STEPS) //leave macro as is
//#endif
#define THREADSNR 1 						// most efficient at dimension
#define	KBOLTZ  1.380641e-8
#define	TEMP  	100.0
#define	GAMMA  1.0 							// drag
#define LOWER_CUTOFF_FLAG 0 				// 1 for lower cutoff frequency when settinh up the ommega
#define	ALPHA  0.8							// exponent free random walk msqdplcmnt t^alpha
#define	mass  1.0 							// mass tracer particle
#define	TIME_ARRIVAL 3000.0 				// Diffusion time over lattice cell
#define MAXSTEPSIZE_BINARY	0 				// set to 1 to decouple from smallest period to fixed N * DT = MAXSTEPSIEMULT
#define MAXSTEPSIZEMULT    5.0e-4			// stepsize in terms of shortest period scale
#define	NUMBER_TO_INF  (500)				// muss gerade sein, hindernisse die Symetrisch um null verteilt sind in 1-D
#define	VOL_FRAC  0.0
#define LABEL "ZWANZIG_THERMALIZED"			// "ZWANZIG" or "ZWANZIG_THERMALISED"
#define SUMMATION_METHOD Neumaier_Sum  		// Sign_, Kahan oder Neumaier_ gefolgt von sum
#define STEPPER_METHOD VVerlet_Step_deriv_Box	//bis jetzt VVerlet_Step_Target_Square oder VVerlet_Step_deriv oder VVerlet_Step_Yukawa Verlet_Step_Target_Pore_Rectangle
#define POTI_HANDLE deriv_parallel
#define NORM_INT 0  						// für Norm_Diff die art der p-norm, 0 für infty
#define MAX_NR_PLOTS 50
#define ESCAPE_CELLS  0                		// anzahl der gitterzellen ohne hinderniss -1
#define REFLEC_BINARY 1 					// toggle bewteewn reflecting new(1) or old(0) velocity, or no velcotiy adjust (2)
#define FLAG_RAND_BATH 0 					// set 1 for unifrom random bath, 0 for uniform deteriministic

#define DENSITY_SAMPLING_LENGTH 10
#define DENSITY_SAMPLING_POINTS 10
#define DIFF_COEFF_MODIFY 0.01 				// modify for VVerlet_Step_drag_sphere
#define DENSITY_SNAPSHOTS 5
#define TARGET_METHOD	Get_Lattice_Targets_Corner_Square // method to determine free cells
#define PARALLEL  //nicht definieren falls ohne -fopenmp compiliert!
#ifdef PARALLEL
	#define THREADS THREADSNR
	#define SETNUMTHREADS omp_set_num_threads(THREADSNR);
#else
	#define SETNUMTHREADS
	#define THREADS 1  //wichtig zur zeitberechnung dass dann THREADS auf 1 gesetzt wird, wenn nicht parallel operiert
#endif

//#define OCL
#ifdef OCL
	#define VVERLET Opencl_VVerlet
#else
	#define VVERLET VVerlet_parallel
#endif
EXTERN int FLAG_DUMP_VALUE;			// if error occurs set 1, dump calculation, start fresh
EXTERN int NUMBER_PORE_OBSTACLES;  	//number of spheres in Yukawa pore wall per side
EXTERN double LATTICE_SPACING;
EXTERN double DIFF_COEFF;
EXTERN double TARGET_LENGTH;
EXTERN double VIRTUAL_X[DIM];
EXTERN double VIRTUAL_V[DIM];
EXTERN int OUTPUT_FLAG;				//set 1 before first output, after that 0 to hide messages
									// from Setuo functions for ommega and gamma
EXTERN int TARGET_CONTACT;			// counts number of hits on lattice
EXTERN int TARGET_FLAG;				// Toggles after first hit per run
EXTERN int lattice_position[DIM];	// aktuelle Gitterzelle
EXTERN double massq[OSSZI];
EXTERN double time_first_contact;
EXTERN double coupling[OSSZI];  //Koppl. Konst. Gamma im Zwanzig
EXTERN double ommega[OSSZI]; 	// Osszi Kreisfreq
EXTERN double y[ORDER];			// Startvec
EXTERN double MAXSTEPSIZE;		//Schrittweite Integrator
EXTERN double LASTSTEPSIZE;
EXTERN int **POSITIONS; 	//enthält Koordinated der Gitterpunkte
EXTERN double L_BALL;		// ballistische Länge der Diffusion
EXTERN int (*Sum_Method)(int N, double input[], double *ans);
EXTERN void (*Poti_Handle)(const double *y, double *ans, double t,int N); 			// externe oder Hindernisskraft
EXTERN void (*Stepper_Method)(int N,  double  x[N/2] , double v[N/2], double a[N/2], double *t,
					void (*deriv) (double *y, double *ans, double t,int N));
					// Methode mit oder ohne Hindernisse für Verlet_schritt
EXTERN gsl_rng * RAND_GLOBAL;	// globally accessible random variable

ext: h


text:
#ifndef UTILITIES_H
#define UTILITIES_H

#include <gsl/gsl_rng.h>
#include <gsl/gsl_histogram.h>
#include <fftw3.h>

gsl_rng * init_rng(void);

void allocate_arrays(double ** signal, double ** spectrum,
    const size_t length);

gsl_histogram * init_histogram(const double r_min, const double r_max,
    const size_t bin_n);

gsl_histogram * init_histogram_uniform(const double r_min, const double r_max,
    const size_t bin_n);

void load_wisdom(const char * filename);

fftw_plan make_plan(double * signal, const size_t length);

void make_fft(fftw_plan plan);

void remove_mean(double * signal, const size_t length);

void compute_spectrum(const double * signal, double * spectrum,
    const size_t length);

void smooth(const unsigned neighbors, double *spectrum, const size_t length);

void final_spectrum(double *spectrum, const size_t length,
    const unsigned count, const double delta, const unsigned neighbors);

double frequency(const size_t i, const size_t length, const double delta);

void output_spectrum(const char * filename, const double * spectrum,
    const size_t length, const double delta, const size_t n_points);

void output_signal(const char * filename, const double * signal,
    const size_t length, const double delta, const size_t n_points);

void output_signal_full(const char * filename, const double * signal,
    const size_t length, const double delta);

void output_distribution(const char * filename, const gsl_histogram * hist);

void cleanup(double * signal, double * spectrum, fftw_plan plan, gsl_rng * r,
    gsl_histogram * hist);

#endif


text:
#ifndef _CHAIN_H_
#define _CHAIN_H_

/* ---------------------------------------------------------
   headers
   --------------------------------------------------------- */
#include <stdio.h>
#include <string>
#include <stdlib.h>
#include <iostream>
#include <vector>

// my headers
#include <md.h>
#include <debug.h>

#include "dcdio.h"

// headers 3rd
/* #include <armadillo> */
/* #include <gsl/gsl_math.h> */
/* #include <gsl/gsl_eigen.h> */

/* #include "boost/multi_array.hpp" */
#include <cassert>

/* SYMLINK: Eigen /usr/local/include/. */
/* [dale~/sop_dev/contacts/segment]$lt /usr/local/include/eigen                        */
/* lrwxrwxrwx 1 root root 25 09.25.2015 21:04 /usr/local/include/eigen -> /usr/include/eigen3/Eigen/ */
/* #include<Eigen/Eigenvalues> */
/* #include <Eigen/Dense> */

#ifdef INERTIA
#include <Eigen/Core>
#include <Eigen/Eigenvalues>
using namespace Eigen;
#endif // INERTIA

/* using namespace Eigen; */

 // namespaces
/* using namespace arma; */

/* ---------------------------------------------------------
   definitions
   --------------------------------------------------------- */
/* 2000, 100;  700,20; 500,50 */
#if defined ONEMOL || defined ALLATOM
#define MOLECULE_INITIAL_SIZE 6900
#define MOLECULE_INITIAL_SIZE_MT 6900
#else
#define MOLECULE_INITIAL_SIZE 700
#define MOLECULE_INITIAL_SIZE_MT 700
#endif // ONEMOL
#define MAX_CONTACTS 24 /* per atom */

class Contact {
public:
    /* int resid; */
    int cresid;

    int total;

    int index; // index
    int cindex; // contacted index

    double eh; // energy scaling (well depth)
    double distance;

    // Constructor
    Contact (); // Constructor declared.

    // Destructors:
    ~Contact() {
        /* delete [] pos; */
    };

};
inline Contact::Contact() {
    cresid = -1;
    index = -1;
    cindex = -1;
    total = -1;
    eh = 0.0; // keep at 0.0, used in assign_contact_eh as default
    distance = 0.0;
}

typedef struct Chain Molecule;

class Chain {
public:
    char filename[46];

    int chainid;
    int num_atoms;
    int num_atoms_ca;

    int index;
    int findex;
    int pf_pos_low; // for selecting index near midpoint of centroids, curvature of PF
    int pf_pos_high;
    int file_line_begin;
    int file_line_end;

    char atomtype[MOLECULE_INITIAL_SIZE][4];
    char resname[MOLECULE_INITIAL_SIZE][4];
    int resid[MOLECULE_INITIAL_SIZE];
    int indices[MOLECULE_INITIAL_SIZE]; // not used yet.

    Vector pos[MOLECULE_INITIAL_SIZE];
    Contact contacts[MOLECULE_INITIAL_SIZE][MAX_CONTACTS];

    int total_num_contacts_2x;
    int tally_num_contacts_2x[MOLECULE_INITIAL_SIZE];
    int num_contacts_2x_persist_ref[MOLECULE_INITIAL_SIZE];

    Vector centroid;
    int nearest_targets[32];
    double nearest_targets_dist[32];
    int nearestneighbors_chainid[8];
    double nearestneighbors_dist[8];

    int LonN; // contacted monomer
    int LatE;
    int LonS;
    int LatW;
    int LonN_total2x; // total
    int LatE_total2x;
    int LonS_total2x;
    int LatW_total2x;

    // general inter
    Contact intercontacts[MOLECULE_INITIAL_SIZE][MAX_CONTACTS];

    // latlon contacts
    Contact contactsLonN[MOLECULE_INITIAL_SIZE][MAX_CONTACTS];
    Contact contactsLatE[MOLECULE_INITIAL_SIZE][MAX_CONTACTS];
    Contact contactsLonS[MOLECULE_INITIAL_SIZE][MAX_CONTACTS];
    Contact contactsLatW[MOLECULE_INITIAL_SIZE][MAX_CONTACTS];

    // Constructors:
    Chain(int a=-1, int b = 0, double p = 0.0, double q = -1.0) {
        chainid = a;
        num_atoms = a;
        num_atoms_ca = a;
        index = a;
        findex = a;
        pf_pos_low = a;
        pf_pos_high = a;
        file_line_begin = a;
        file_line_end = a;
        total_num_contacts_2x = a;

        centroid.x = p;
        centroid.y = p;
        centroid.z = p;

        LonN = a;
        LonS = a;
        LatW = a;
        LatE = a;

        // initialize to 0.0 p
        for (int i=0; i < MOLECULE_INITIAL_SIZE; i++ ) {
            pos[i].x = p;
            pos[i].y = p;
            pos[i].z = p;

            // initialize to 0 b
            tally_num_contacts_2x[i] = b;
            num_contacts_2x_persist_ref[i] = b;

            for ( int j=0; j < MAX_CONTACTS; j++ ) {
                contacts[i][j].cresid = a;
                contacts[i][j].distance = p;
            }
            for ( int j=0; j < 8; j++ ) {
                nearestneighbors_chainid[j] = -1;
                nearestneighbors_dist[j] = p;
            }
            for ( int j=0; j < 32; j++ ) {
                nearest_targets[j] = a;
            }
        }
    }

    // Destructors:
    ~Chain() {
        /* delete [] pos; */
        /* std::cout << "Object is being deleted" << std::endl; */
    };

    // Declarations:
    void print_prop();
    void assign_indices();
    void fprintf_indices();
    void GetMomentofInertiaTensor();
    void ComputeCentroid();
    void print_centroid();

    int countcontacts(Contact (*map)[MAX_CONTACTS]);

};
inline void Chain::ComputeCentroid() {
    int i;
    double xtot = 0.0, ytot = 0.0, ztot = 0.0;
    /* printf("computing centroid ca: %d\n",num_atoms_ca); */
    for ( int i=0; i < num_atoms_ca; i++ ) {
        /* printf("%d\n",i); */
        xtot += pos[i].x;
        ytot += pos[i].y;
        ztot += pos[i].z;
    }
    centroid.x = xtot / num_atoms_ca;
    centroid.y = ytot / num_atoms_ca;
    centroid.z = ztot / num_atoms_ca;
    // PRINT HERE
    /* debug("c(%d): %f %f %f\n",chainid,centroid.x,centroid.y,centroid.z); */
}

// amazing!!
/* By declared that method with inline keyword, compiler will either inline */
/*     the whole method or, if it decides not to, it will generate anonymous */
/*     method (same method but with some unique name for a given object file), */
/*     so there will be no conflicts in object files. For example: */
/* class foo { */
/* public: */
/*   void bar (); */
/* }; */
/* inline void foo::bar () */
/* { */
/* } */
inline void Chain::print_prop() {
    /* printf("id: %d ",chainid); */
    /* printf(" filename: %s\n",filename); */
    /* printf("CA: %d ",num_atoms_ca); */
    /* printf(" index(%d) ",index); */
    /* printf(" final_index(%d)\n",findex); */
    /* printf("begin_end: [ %d %d ]\n\n",file_line_begin,file_line_end); */

    /* debug("id: %d   CA: %d\n",chainid,num_atoms_ca); */
    /* debug(" filename: %s\n",filename); */
    /* debug("CA: %d \n",num_atoms_ca); */
    /* debug(" index(%d) \n",index); */
    /* debug(" final_index(%d)\n",findex); */
    /* debug("begin_end: [ %d %d ]\n",file_line_begin,file_line_end); */

    printf("centroid: %f %f %f\n",centroid.x,centroid.y,centroid.z);
}
inline void Chain::print_centroid() {
    printf("centroid: %f %f %f\n",centroid.x,centroid.y,centroid.z);
}
inline void Chain::assign_indices() {
    /* printf("id: %d ",chainid); */
    /* printf(" filename: %s\n",filename); */
    /* printf("CA: %d ",num_atoms_ca); */
    /* printf(" index(%d) ",index); */
    /* printf(" final_index(%d)\n",findex); */
    /* printf("begin_end: [ %d %d ]\n\n",file_line_begin,file_line_end); */

    debug("id: %d   atoms:%d   CA:%d\n",chainid,num_atoms,num_atoms_ca);
    /* debug(" filename: %s\n",filename); */
    /* debug("CA: %d \n",num_atoms_ca); */
    /* debug(" index(%d) \n",index); */
    /* debug(" final_index(%d)\n",findex); */
    debug("begin_end: [ %d %d ]\n",file_line_begin,file_line_end);

    for( int i=0; i < num_atoms_ca; i++ ) {
        indices[i] = index + i;
    }
}
inline void Chain::fprintf_indices() {
    /* printf(" index(%d) ",index); */
    /* printf(" final_index(%d)\n",findex); */
    printf(" %d__[%d|%d]\n",chainid,index,findex);

    FILE * fp_indices;
    fp_indices = fopen("indices_foreach_segments.dat","a+");
    for ( int i=0; i < num_atoms_ca; i++ ) {
        fprintf(fp_indices,"%6d  %6d  %6d\n",chainid,index,findex);
    }
    fclose(fp_indices);

    /* for( int i=0; i < num_atoms_ca; i++ ) { */
    /*     /\* indices[i] = index + i; *\/ */
    /*     printf(" %d",indices[i]); */
    /*     printf(" ") */
    /* } */
    /* printf("\n"); */
}

inline int Chain::countcontacts(Contact (*map)[MAX_CONTACTS]) {

    /* printf("declare countcontacts!\n"); */
    int count = 0;

    // using 0,3 indices. 1,4, 2,5, ... 430,433
    for(int i=0; i<MOLECULE_INITIAL_SIZE; i++){
        // for(int j=0; j<MAX_CONTACTS; j++){
        for(int k=0; k<MAX_CONTACTS; k++){

            if(map[i][k].cresid == -1){
                break;
            } else {
                count ++;
                // printf("%d-%d (%d)  %d/50\n",i,map[i][k].cresid,count,k);
            }
            // printf(" %d",map[i][k].index);
            // map[num1][k].cindex = chain[num1].indices[j];
            // double eh; // energy scaling (well depth)
            // map[num1][k].distance = dist;

            // if(map[i][k].cresid != -1){
            //     count ++;
            // } // if
            // } // for 3
            // } // if
        } // for 2
    } // for 1

    if(count != 0){
        debug("declare total contacts in this map are: %d\n",count);
    } else {
        debug("empty map.\n");
    }
    // // map[chain_num].total = count;

    return count;
}

#ifdef INERTIA
inline void Chain::GetMomentofInertiaTensor() {

    /* printf("chain: %d;  atoms: %d;\n",chainid,num_atoms_ca); */
    double x2=0.0,y2=0.0,z2=0.0;
    double xy=0.0,xz=0.0,yz=0.0;

    double Ixx=0.0,Iyy=0.0,Izz=0.0;
    double Ixy=0.0,Ixz=0.0,Iyz=0.0;

    /* http://farside.ph.utexas.edu/teaching/336k/Newtonhtml/node64.html moment of Inertia Tensor */
    /* http://farside.ph.utexas.edu/teaching/336k/Newtonhtml/node67.html get eigenvectors */
    for(int a=0; a<num_atoms_ca; a++) {
    /* for(int a=0; a<5; a++) { */
        /* printf("%d",a); */
        /* printf("%f %f %f\n",pos[a].x,pos[a].y,pos[a].z); */
        x2 += pos[a].x * pos[a].x;
        y2 += pos[a].y * pos[a].y;
        z2 += pos[a].z * pos[a].z;

        xy += pos[a].x * pos[a].y;
        xz += pos[a].x * pos[a].z;
        yz += pos[a].y * pos[a].z;

        Ixx += y2 + z2;
        Iyy += x2 + z2;
        Izz += x2 + y2;

        Ixy += pos[a].x * pos[a].y;
        Ixz += pos[a].x * pos[a].z;
        Iyz += pos[a].y * pos[a].z;
    }

    Ixy = Ixy * -1.0;
    Ixz = Ixz * -1.0;
    Iyz = Iyz * -1.0;

    printf("xx_yy_zz: %f %f %f\n",Ixx,Iyy,Izz);
    printf("xy_xz_yz: %f %f %f\n",Ixy,Ixz,Iyz);

    // begin eigen library

    Matrix3d A(3,3);
    A(0,0) = Ixx;
    A(0,1) = Ixy;
    A(0,2) = Ixz;

    A(1,0) = Ixy;
    A(1,1) = Iyy;
    A(1,2) = Iyz;

    A(2,0) = Ixz;
    A(2,1) = Iyz;
    A(2,2) = Izz;
    /* std::cout << "Here is the matrix A:\n" << A << std::endl; */

    // build matrix
    double InertiaTensor[3][3] = {
        {Ixx, Ixy, Ixz},
        {Ixy, Iyy, Iyz},
        {Ixz, Iyz, Izz},
    };

    // Ixy xz yz are altered beyond this point.
    double lowest_value_diag = Ixx;
    if ( Iyy < Ixx ) {
        if ( Izz < Iyy ) {
            lowest_value_diag = Izz;
        } else {
            lowest_value_diag = Iyy;
        }
    } else if ( Izz < Ixx ) {
        lowest_value_diag = Izz;
    }

    double lowest_value_offdiag;
    if ( Ixy < 0.0 ) {
        Ixy = -1.0 * Ixy;
    }
    if ( Ixz < 0.0 ) {
        Ixz = -1.0 * Ixz;
    }
    if ( Iyz < 0.0 ) {
        Iyz = -1.0 * Iyz;
    }

    lowest_value_offdiag = Ixy;
    if ( Ixz < Ixy ) {
        if ( Iyz < Ixz ) {
            lowest_value_offdiag = Iyz;
        } else {
            lowest_value_offdiag = Ixz;
        }
    } else if ( Iyz < Ixy ) {
        lowest_value_offdiag = Iyz;
    }

    double lowest_value;
    if ( lowest_value_offdiag < lowest_value_diag ) {
        lowest_value = lowest_value_offdiag;
    } else {
        lowest_value = lowest_value_offdiag;
    }

    printf("lv: %f\n",lowest_value);
    lowest_value = lowest_value * 0.95;
    /* printf("lv: %f\n",lowest_value); */
    // Matrix entries * lowest_value;

    /* // ORIGINAL EXAMPLE */
    /* // Create a 3D array that is 3 x 4 x 2 */
    /* typedef boost::multi_array<double, 3> array_type; */
    /* typedef array_type::index index; */
    /* array_type A(boost::extents[3][4][2]); */

    /* // Assign values to the elements */
    /* int values = 0; */
    /* for(index i = 0; i != 3; ++i) */
    /*     for(index j = 0; j != 4; ++j) */
    /*         for(index k = 0; k != 2; ++k) */
    /*             A[i][j][k] = values++; */

    /* // Verify values */
    /* int verify = 0; */
    /* for(index i = 0; i != 3; ++i) */
    /*     for(index j = 0; j != 4; ++j) */
    /*         for(index k = 0; k != 2; ++k) */
    /*             assert(A[i][j][k] == verify++); */
    /* /\* return 0; *\/ */

    /* ---------------------------------------------------------
       boost! no eigenvalue/eigenvector
       --------------------------------------------------------- */
    /* // Create a 3D array that is 3 x 4 x 2 */
    /* typedef boost::multi_array<double, 2> array_type; */
    /* typedef array_type::index index; */
    /* array_type A(boost::extents[3][3]); */
    /* array_type B(boost::extents[3][3]); */

    /* // x row */
    /* A[0][0] = Ixx; */
    /* A[0][1] = Ixy; */
    /* A[0][2] = Ixz; */

    /* // y row */
    /* A[1][0] = Ixy; */
    /* A[1][1] = Iyy; */
    /* A[1][2] = Iyz; */

    /* // y row */
    /* A[2][0] = Ixz; */
    /* A[2][1] = Iyz; */
    /* A[2][2] = Izz; */

    /* /\* // Verify values *\/ */
    /* /\* int verify = 0; *\/ */
    /* /\* for(index i = 0; i != 3; i++) { *\/ */
    /* /\*     for(index j = 0; j != 3; j++) { *\/ */
    /* /\*         /\\* assert(A[i][j] == verify++); *\\/ *\/ */
    /* /\*         // reassign by lowest_value scaling *\/ */
    /* /\*         printf("A: %f\n",A[i][j]); *\/ */
    /* /\*     } *\/ */
    /* /\* } *\/ */

    /* // Verify values A */
    /* for(index i = 0; i != 3; i++) { */
    /*     for(index j = 0; j != 3; j++) { */
    /*         // reassign by lowest_value scaling */
    /*         A[i][j] = A[i][j] / lowest_value; */
    /*         printf("A/d: %f\n",A[i][j]); */
    /*     } */
    /* } */
    /* ---------------------------------------------------------
       boost! no eigenvalue/eigenvector - END
       --------------------------------------------------------- */

    /* gsl_matrix_view m = gsl_matrix_view_array (InertiaTensor, 3, 3); */
    /* gsl_vector *eval = gsl_vector_alloc (3); */
    /* gsl_matrix *evec = gsl_matrix_alloc (3, 3); */

    /* gsl_eigen_symmv_workspace * w = gsl_eigen_symmv_alloc (3); */
    /* gsl_eigen_symmv (&m.matrix, eval, evec, w); */
    /* gsl_eigen_symmv_free (w); */
    /* gsl_eigen_symmv_sort (eval,evec,GSL_EIGEN_SORT_ABS_ASC); */
    /* { */
    /*     int i; */

    /*     for (i = 0; i < 4; i++) */
    /*     { */
    /*         double eval_i */
    /*             = gsl_vector_get (eval, i); */
    /*         gsl_vector_view evec_i */
    /*             = gsl_matrix_column (evec, i); */

    /*         printf ("eigenvalue = %g\n", eval_i); */
    /*         printf ("eigenvector = \n"); */
    /*         gsl_vector_fprintf (stdout, */
    /*                             &evec_i.vector, "%g"); */
    /*     } */
    /* } */

    /* gsl_vector_free (eval); */
    /* gsl_matrix_free (evec); */

    /* return 0; */

    // eigen continued
    A = A / lowest_value;
    std::cout << "Here is the matrix A:\n" << A << std::endl;

    VectorXcd eivals = A.eigenvalues();
    std::cout << "The eigenvalues of the 3x3 matrix of ones are:" << std::endl << eivals << std::endl;

    /* VectorXcd eivecs = A.eigenvectors(); */
    /* std::cout << "The eigenvalues of the 3x3 matrix of ones are:" << std::endl << eivecs << std::endl; */

    EigenSolver<MatrixXd> es(A);
    std::cout << "The first eigenvector of the 3x3 matrix of ones is:"
              << std::endl << es.eigenvectors().col(0) << std::endl;

    std::cout << "The second eigenvector of the 3x3 matrix of ones is:"
              << std::endl << es.eigenvectors().col(1) << std::endl;

    std::cout << "The third eigenvector of the 3x3 matrix of ones is:"
              << std::endl << es.eigenvectors().col(2) << std::endl;

    /* Matrix3d m = Matrix3d::(3,3); */
    /* cout << "m =" << endl << m << endl; */

}
#endif // INERTIA

/* class Map { */
/*     /\* int num_atoms,max_contacts_per_atom; *\/ */
/* public: */

/*     /\* std::vector< std::vector< std::vector<int>>> map; *\/ */
/*     std::vector< std::vector<int> > map; */
/*     /\* std::vector *\/ */
/*     std::vector<Contact> contacts; */

/*     // Constructor */
/*     /\* Map(); *\/ */

/*     /\* declarations: *\/ */
/*     /\* void extend_map(int); *\/ */
/*     void build_map_default(Chain *chain, int cid); // use chain_ref */

/* }; */
/* /\* inline Map::extend_map(int ca){ *\/ */
/* /\*     for (int i=0; i<ca; i++) { *\/ */
/* /\*         map.push_back(i); // ADD the ROW. *\/ */
/* /\*     } *\/ */
/* /\* } *\/ */
/* inline void Map::build_map_default(Chain *chain,int cid){ */
/*     // default: 8A */

/*     /\* int count; *\/ */
/*     double dist; */
/*     Contact x; */
/*         /\* ,y; *\/ */

/*     for (int i=0; i<chain[cid].num_atoms_ca; i++) { */
/*         /\* printf("%d\n",i); *\/ */

/*         for(int j=i; j<chain[cid].num_atoms_ca; j++) { */

/*             if ( i >= j-2 && i <= j+2) { */
/*                 // exclude j = i +/- 2 */
/*                 /\* printf("%d\t no --> %d j: %d %d\n",i,j,j-2,j+2); *\/ */
/*             } else { */

/*                 /\* printf("%d\n",j); *\/ */
/*                 dist = distance(chain[cid].pos[i],chain[cid].pos[j]); */
/*                 /\* printf("%8.4f\n",dist); *\/ */

/*                 if ( dist <= 8.0 ) { */

/*                     x.cresid = j; */
/*                     x.distance = dist; */
/*                     x.index = i + chain[cid].index; */
/*                     x.cindex = j + chain[cid].index; */

/*                     contacts.insert(contacts.end(),x); */
/*                 } /\* END 8.0 *\/ */
/*             } // else/if */

/*         } // j */

/*         /\* map.push_back(contacts); *\/ */
/*         /\* contacts.clear(); *\/ */
/*     } // i */

/*     /\* std::vector<int>::iterator it; *\/ */
/*     /\* for (it=contacts.begin(); it<contacts.end(); it++) { *\/ */
/*     /\*     // std::cout << *it << ' ' << std::endl; *\/ */
/*     /\*     printf("%d ",*it); *\/ */
/*     /\* } *\/ */

/*     /\* std::vector<int>::iterator iter; // declare a general iterator *\/ */
/*     for ( int k = contacts.begin(); k < contacts.end(); k++){ */
/*         printf("%d\n",contacts[k].index); */
/*     } */

/* } */

/* inline Map::Map() { */
/*     // constructor */

/*     for (int i=0; i<num_atoms; i++ ) { */
/*         /\* printf("%d\n",i); *\/ */
/*         for (int j=0; j<max_contacts_per_atom; j++){ */
/*             /\* if ( mol->contacts[i][j].cresid != 0 || mol->contacts[i][j].distance != 0.0) { *\/ */
/*             /\* printf("%d %f\n",mol->contacts[i][j].cresid,mol->contacts[i][j].distance); *\/ */
/*             contacts[i][j].cresid = -1; */
/*             contacts[i][j].distance = 0.0; */
/*         } */
/*     } */
/* } */
/* inline void Map::get_contacts() { */

/* } */
/* inline void Map::set_values (int natoms, int max_contacts) { */
/*     num_atoms = natoms; */
/*     max_contacts_per_atom = max_contacts; */

/* } */

/* ---------------------------------------------------------
   Molecule preparations:
   --------------------------------------------------------- */
void build_contact_map(Molecule *mol);
void build_latlon_contact_map(Chain *chain, int num_chains, int start, int stop,\
                              Contact (*map)[MAX_CONTACTS],int sel);
void build_contact_map_inter_monomer(Chain *chain,int cid, Contact (*map)[MAX_CONTACTS], int sel);
void build_contact_map_inter(Chain *chain,int cid1,int cid2,Contact (*map)[MAX_CONTACTS]);

void compare_contacts(Chain *chref,Chain *chain,int cid, int sel);
void evaluate_original_contacts_now(Chain *chref,Chain *chain,int cid, int sel);

void print_contact_map(Contact (*map)[MAX_CONTACTS]);
void fprintf_count_of_contacts_2x_persist_from_reference( Molecule *m1, Molecule *m2);
void fprintf_bond_vector_angle_with_tension_vector( Molecule *m1, Molecule *m2 );
void fprintf_contacts_for_one_monomer_face(Chain *chain, int num_chains, int sel);
void fprintf_contact_map_gsop(Contact (*map)[MAX_CONTACTS]);
void fprintf_all_indices(Chain *chref,int num_chains);

void tubulin_monomers_angles( Molecule *m1, Molecule *m2, int r1, int r2, FILE *fp);
void get_tubulin_centroid( Molecule *m1, Vector *com);
void get_centroid_from_array( Vector arr_centroid[],int i, int j, Vector *centroid);
Vector get_centroid_from_centroids(Chain *chain,int start_chain,int stop_chain);
void get_centroid_movement(Chain *chain_ref,Chain *chain_later,int num_chains);
void get_centroid_xyzposition_movement(Chain *chain_later,int chains_to_use,FILE *fp_xyzpos);
/* void get_centroid_xyzposition_movement(chain_later,chains_to_use,fp_xyzpos); */

/* void get_latangle_from35centroids(Chain *chain_ref,Chain *chain_later,int m,int sel); */
void get_latangle_from35centroids(Chain *chain_ref,Chain *chain_later,int m,FILE *fp,int sel);
std::vector<double> get_3angles_dimerbyalpha(Chain *chain_ref,Chain *chain_later,int cid);

/* void get_latangle_from35centroids(Chain *chain_ref,Chain *chain_later,int m,int sel,FILE *fp_ang); */

void get_midpoint_reference_indices(Chain *chain,int i,int *index_ref_i,int *index_ref_c,int sel);
/* get_midpoint_reference_indices(chain_ref,i,&index_ref_i,&index_ref_c,1); */

void compute_nearest_neighbors(Chain *chain, int num_chains);
void determine_latlon_neighbors(Chain *chains, int num_chains, Vector axis);
void get_monomer_monomer_info( Molecule *m1, Molecule *m2,\
                               Vector ref_norm, Vector m1centroid,\
                               Vector m2centroid );

/* void load_dcd_coords_to_chain(dcdhandle *dcd,Chain *chain,int num_chains); */
/* void load_chain_coords_to_timestep(Chain *chain,int num_chains,\ */
/*                                    const molfile_timestep_t *ts); */

void print_chain_coords(Chain *chain,int num_chains);

// overloaded.
void assign_contact_eh(Contact (*map)[MAX_CONTACTS],float eh,int low_resid,\
                       int high_resid,int low_cresid,int high_cresid);
void assign_contact_eh(Contact (*map)[MAX_CONTACTS],float eh);

/* void dcd_to_chain(ddc,chain_0); */
/* void dcd_to_chain(dcdhandle *v,Chain *chain); */

/* void load_chain_coords_to_timestep(dcdhandle *dcd,Chain *chain,int num_chains,\ */
/*                                    void *v, const molfile_timestep_t *ts,int natoms); */
/* load_chain_coords_to_timestep(dcdw,chain_later,chains_to_use,vw,&timestep_w,natoms_w); */

/* static int read_next_timestep(void *v, int natoms, molfile_timestep_t *ts) { */
/* static int write_timestep(void *v, const molfile_timestep_t *ts) { */

#endif


text:
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <assert.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_sf_result.h>
#include <fftw3.h>

#include "cosmocalc.h"

double get_k0_fftlog(int N, double mu, double q, double r0, double L, double k0guess);
void um_fftlog(int m, double mu, double q, double k0, double r0, double L, double *realpart, double *imagpart);

void compute_discrete_spherical_fft(double *data, int N, double r0, double L, double q, double *result, double *k0)
{
  double mu,rn,kn,tmp[2],um[2],k0guess;
  int i,m;
  static int initFlag = 1;
  static double *rin;
  static fftw_complex *out;
  static fftw_plan pf,pb;

  if(initFlag)
    {
      initFlag = 0;

      rin = (double*)fftw_malloc(sizeof(double)*N);
      out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*(N/2+1));
      pf = fftw_plan_dft_r2c_1d(N,rin,out,FFTW_MEASURE);
      pb = fftw_plan_dft_c2r_1d(N,out,rin,FFTW_MEASURE);
    }

  for(i=0;i<N;++i)
    {
      rn = r0*exp(L/((double) N)*((double) (i-N/2)));
      rin[i] = data[i]*4.0*M_PI*pow(rn,1.5)*sqrt(M_PI/2.0);
    }
  fftw_execute(pf);

  q = 0.0;
  mu = 0.5;
  k0guess = *k0;
  *k0 = get_k0_fftlog(N,mu,q,r0,L,k0guess);
  for(i=0;i<N/2+1;++i)
    {
      if(i >= N/2)
	m = i-N;
      else
	m = i;

      um_fftlog(m,mu,q,*k0,r0,L,&(um[0]),&(um[1]));
      if(m == N/2 || m == -N/2)
	um[1] = 0.0;

      tmp[0] = out[i][0]*um[0] - out[i][1]*um[1];
      tmp[1] = out[i][1]*um[0] + out[i][0]*um[1];

      out[i][0] = tmp[0];
      out[i][1] = tmp[1];
    }

  fftw_execute(pb);

  for(i=0;i<N;++i)
    {
      kn = (*k0)*exp(L/((double) N)*((double) (N/2-i)));
      result[N-i-1] = rin[i]/((double) N)/pow(kn,1.5)*exp(L/((double) N)*3.0/2.0);
    }

  /* not needed
     fftw_destroy_plan(pf);
     fftw_destroy_plan(pb);

     fftw_free(rin);
     fftw_free(out);
  */
}

double get_k0_fftlog(int N, double mu, double q, double r0, double L, double k0guess)
{
  double zr,zi,arg;
  gsl_sf_result lnr_top,lnr_bot;
  gsl_sf_result arg_top,arg_bot;
  double twox[2];
  double grat[2];
  double tmp[2];
  double lnkr;
  int intval;

  //do top part of gamma function ratio
  zr = (mu + 1.0 + q)/2.0;
  zi = M_PI*((double) N)/L/2.0;
  gsl_sf_lngamma_complex_e(zr,zi,&lnr_top,&arg_top);

  //do bottom
  zr = (mu + 1.0 - q)/2.0;
  zi = -1.0*M_PI*((double) N)/L/2.0;
  gsl_sf_lngamma_complex_e(zr,zi,&lnr_bot,&arg_bot);

  //get real and imag parts of 2^(q+i*2*pi*m/L)
  zr = pow(2.0,q);
  twox[0] = zr*cos(M_PI*((double) N)/L*log(2.0));
  twox[1] = zr*sin(M_PI*((double) N)/L*log(2.0));

  //now do computation of real and imag parts of result
  zr = exp(lnr_top.val-lnr_bot.val);
  grat[0] = zr*cos(arg_top.val-arg_bot.val);
  grat[1] = zr*sin(arg_top.val-arg_bot.val);

  //do product
  tmp[0] = grat[0]*twox[0] - grat[1]*twox[1];
  tmp[1] = grat[0]*twox[1] + grat[1]*twox[0];

  //find k0
  lnkr = log(k0guess*r0);
  arg = atan2(tmp[1],tmp[0])/M_PI*L/((double) N);
  intval = (int) ((lnkr-arg)/(L/((double) N)));

  return exp(arg + ((double) intval)*L/((double) N))/r0;
}

void um_fftlog(int m, double mu, double q, double k0, double r0, double L, double *realpart, double *imagpart)
{
  double zr,zi;
  gsl_sf_result lnr_top,lnr_bot;
  gsl_sf_result arg_top,arg_bot;
  double kr[2];
  double twox[2];
  double grat[2];
  double tmp[2];

  //do top part of gamma function ratio
  zr = (mu + 1.0 + q)/2.0;
  zi = M_PI*((double) m)/L;
  gsl_sf_lngamma_complex_e(zr,zi,&lnr_top,&arg_top);

  //do bottom
  zr = (mu + 1.0 - q)/2.0;
  zi = -1.0*M_PI*((double) m)/L;
  gsl_sf_lngamma_complex_e(zr,zi,&lnr_bot,&arg_bot);

  //get real and imag parts of (k0*r0)^(i*2*pi*m/L)
  zr = -1.0*log(k0*r0)*2.0*M_PI*((double) m)/L;
  kr[0] = cos(zr);
  kr[1] = sin(zr);

  //get real and imag parts of 2^(q+i*2*pi*m/L)
  zr = pow(2.0,q);
  twox[0] = zr*cos(2.0*M_PI*m/L*log(2.0));
  twox[1] = zr*sin(2.0*M_PI*m/L*log(2.0));

  //now do computation of real and imag parts of result
  zr = exp(lnr_top.val-lnr_bot.val);
  grat[0] = zr*cos(arg_top.val-arg_bot.val);
  grat[1] = zr*sin(arg_top.val-arg_bot.val);

  tmp[0] = kr[0]*twox[0] - kr[1]*twox[1];
  tmp[1] = kr[0]*twox[1] + kr[1]*twox[0];

  *realpart = tmp[0]*grat[0] - tmp[1]*grat[1];
  *imagpart = tmp[0]*grat[1] + tmp[1]*grat[0];
}

ext: c

max_stars_repo_name: beckermr/cosmocalc

max_issues_repo_name: beckermr/cosmocalc

max_forks_repo_name: beckermr/cosmocalc


text:
#ifndef KGS_NULLSPACESVD_H
#define KGS_NULLSPACESVD_H

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <string>

#include "math/SVD.h"
#include "math/Nullspace.h"

/**
 * An implementation of Nullspace backed by SVD
 */
class NullspaceSVD: public Nullspace {
 public:
  /** Will construct a nullspace of `matrix` using the SVD decomposition */
  NullspaceSVD(SVD* svd);

  /** Update the Nullspace (and underlying SVD) to reflect an updated state of the matrix */
  void updateFromMatrix() override;

  /** Return the SVD of the nullspace as columns of a matrix */
  SVD *getSVD() const;

  static void setSingularValueTolerance(double val);

private:
  SVD* m_svd;                  ///< SVD underlying this nullspace

  /// These values have to be chosen according to the numerical analysis
  static double SINGVAL_TOL;
//  static constexpr double RIGID_TOL =   1.0e-10; //0.0000000001; //depends on molecule, but 10^-10 seems a good fit!

  friend class Configuration;
};

#endif //KGS_nullptrSPACE_H


text:
#include <math.h>
#include <stdlib.h>
#include <assert.h>
#include <cconfigspace.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_cdf.h>

#define NUM_AREAS 4
#define NUM_SAMPLES 10000

void test_create_roulette_distribution() {
	ccs_distribution_t      distrib = NULL;
	ccs_result_t            err = CCS_SUCCESS;
	int32_t                 refcount;
	ccs_object_type_t       otype;
	ccs_distribution_type_t dtype;
	ccs_numeric_type_t      data_type;
	ccs_interval_t          interval;
	const size_t            num_areas = NUM_AREAS;
	ccs_float_t             areas[NUM_AREAS];
	size_t                  num_areas_ret;
	ccs_float_t             areas_ret[NUM_AREAS];
	const ccs_float_t       epsilon = 1e-15;

	for(size_t i = 0; i < num_areas; i++) {
		areas[i] = (double)(i+1);
	}

	err = ccs_create_roulette_distribution(
		num_areas,
		areas,
		&distrib);
	assert( err == CCS_SUCCESS );

	err = ccs_object_get_type(distrib, &otype);
	assert( err == CCS_SUCCESS );
	assert( otype == CCS_DISTRIBUTION );

	err = ccs_distribution_get_type(distrib, &dtype);
	assert( err == CCS_SUCCESS );
	assert( dtype == CCS_ROULETTE );

	err = ccs_distribution_get_data_types(distrib, &data_type);
	assert( err == CCS_SUCCESS );
	assert( data_type == CCS_NUM_INTEGER );

        err = ccs_distribution_get_bounds(distrib, &interval);
	assert( err == CCS_SUCCESS );
	assert( interval.type == CCS_NUM_INTEGER );
	assert( interval.lower.i == 0 );
	assert( interval.lower_included == CCS_TRUE );
	assert( interval.upper.i == 4 );
	assert( interval.upper_included == CCS_FALSE );

	err = ccs_roulette_distribution_get_num_areas(distrib, &num_areas_ret);
	assert( err == CCS_SUCCESS );
	assert( num_areas_ret == num_areas );

	err = ccs_roulette_distribution_get_areas(distrib, num_areas, areas_ret, &num_areas_ret);
	assert( err == CCS_SUCCESS );
	assert( num_areas_ret == num_areas );

	ccs_float_t inv_sum = 2.0 / (num_areas * (num_areas + 1));
	for (size_t i = 0; i < num_areas; i++) {
		assert( areas_ret[i] <= areas[i] * inv_sum + epsilon &&
		        areas_ret[i] >= areas[i] * inv_sum - epsilon );
	}

	err = ccs_object_get_refcount(distrib, &refcount);
	assert( err == CCS_SUCCESS );
	assert( refcount == 1 );

	err = ccs_release_object(distrib);
	assert( err == CCS_SUCCESS );
}

void test_create_roulette_distribution_errors() {
	ccs_distribution_t      distrib = NULL;
	ccs_result_t            err = CCS_SUCCESS;
	const size_t            num_areas = NUM_AREAS;
	ccs_float_t             areas[NUM_AREAS];

	for(size_t i = 0; i < num_areas; i++) {
		areas[i] = (double)(i+1);
	}

	err = ccs_create_roulette_distribution(
		0,
		areas,
		&distrib);
	assert( err == -CCS_INVALID_VALUE );

	err = ccs_create_roulette_distribution(
		SIZE_MAX,
		areas,
		&distrib);
	assert( err == -CCS_INVALID_VALUE );

	err = ccs_create_roulette_distribution(
		num_areas,
		NULL,
		&distrib);
	assert( err == -CCS_INVALID_VALUE );

	err = ccs_create_roulette_distribution(
		num_areas,
		areas,
		NULL);
	assert( err == -CCS_INVALID_VALUE );

	areas[1] = -2;
	err = ccs_create_roulette_distribution(
		num_areas,
		areas,
		&distrib);
	assert( err == -CCS_INVALID_VALUE );

}

void test_roulette_distribution() {
	ccs_distribution_t distrib = NULL;
	ccs_rng_t          rng = NULL;
	ccs_result_t       err = CCS_SUCCESS;
	const size_t       num_samples = NUM_SAMPLES;
	ccs_numeric_t      samples[NUM_SAMPLES];
	const size_t       num_areas = NUM_AREAS;
	ccs_float_t        areas[NUM_AREAS];
	int                counts[NUM_AREAS];

	for(size_t i = 0; i < num_areas; i++) {
		areas[i] = (double)(i+1);
	        counts[i] = 0;
	}

	err = ccs_rng_create(&rng);
	assert( err == CCS_SUCCESS );
	err = ccs_create_roulette_distribution(
		num_areas,
		areas,
		&distrib);
	assert( err == CCS_SUCCESS );

	err = ccs_distribution_samples(distrib, rng, num_samples, samples);
	assert( err == CCS_SUCCESS );

	ccs_float_t sum = 0.0;
	ccs_float_t inv_sum = 0.0;
	for(size_t i = 0; i < num_areas; i++) {
		sum += areas[i];
	}
	inv_sum = 1.0 / sum;
	for(size_t i = 0; i < num_samples; i++) {
		assert( samples[i].i >=0 && samples[i].i < 4 );
		counts[samples[i].i]++;
	}
	for(size_t i = 0; i < num_areas; i++) {
		ccs_float_t target = num_samples * areas[i] * inv_sum;
		assert( counts[i] >= target * 0.95 && counts[i] <= target * 1.05 );
	}

	err = ccs_release_object(distrib);
	assert( err == CCS_SUCCESS );
	err = ccs_release_object(rng);
	assert( err == CCS_SUCCESS );
}

void test_roulette_distribution_zero() {
	ccs_distribution_t distrib = NULL;
	ccs_rng_t          rng = NULL;
	ccs_result_t       err = CCS_SUCCESS;
	const size_t       num_samples = NUM_SAMPLES;
	ccs_numeric_t      samples[NUM_SAMPLES];
	const size_t       num_areas = NUM_AREAS;
	ccs_float_t        areas[NUM_AREAS];
	int                counts[NUM_AREAS];

	for(size_t i = 0; i < num_areas; i++) {
		areas[i] = (double)(i+1);
	        counts[i] = 0;
	}
	areas[1] = 0.0;

	err = ccs_rng_create(&rng);
	assert( err == CCS_SUCCESS );
	err = ccs_create_roulette_distribution(
		num_areas,
		areas,
		&distrib);
	assert( err == CCS_SUCCESS );

	err = ccs_distribution_samples(distrib, rng, num_samples, samples);
	assert( err == CCS_SUCCESS );

	ccs_float_t sum = 0.0;
	ccs_float_t inv_sum = 0.0;
	for(size_t i = 0; i < num_areas; i++) {
		sum += areas[i];
	}
	inv_sum = 1.0 / sum;
	for(size_t i = 0; i < num_samples; i++) {
		assert( samples[i].i >=0 && samples[i].i < 4 );
		counts[samples[i].i]++;
	}
	for(size_t i = 0; i < num_areas; i++) {
		ccs_float_t target = num_samples * areas[i] * inv_sum;
		assert( counts[i] >= target * 0.95 && counts[i] <= target * 1.05 );
	}

	err = ccs_release_object(distrib);
	assert( err == CCS_SUCCESS );
	err = ccs_release_object(rng);
	assert( err == CCS_SUCCESS );
}

void test_roulette_distribution_strided_samples() {
	ccs_distribution_t distrib1 = NULL;
	ccs_distribution_t distrib2 = NULL;
	ccs_rng_t          rng = NULL;
	ccs_result_t       err = CCS_SUCCESS;
	const size_t       num_samples = NUM_SAMPLES;
	ccs_numeric_t      samples[NUM_SAMPLES*2];
	const size_t       num_areas = NUM_AREAS;
	ccs_float_t        areas1[NUM_AREAS];
	ccs_float_t        areas2[NUM_AREAS];
	int                counts1[NUM_AREAS];
	int                counts2[NUM_AREAS];

	for(size_t i = 0; i < num_areas; i++) {
		areas1[i] = (double)(i+1);
	        counts1[i] = 0;
		areas2[i] = (double)(i+2);
	        counts2[i] = 0;
	}

	err = ccs_rng_create(&rng);
	assert( err == CCS_SUCCESS );
	err = ccs_create_roulette_distribution(
		num_areas,
		areas1,
		&distrib1);
	assert( err == CCS_SUCCESS );

	err = ccs_create_roulette_distribution(
		num_areas,
		areas2,
		&distrib2);
	assert( err == CCS_SUCCESS );

	err = ccs_distribution_strided_samples(distrib1, rng, num_samples, 2, samples);
	assert( err == CCS_SUCCESS );

	err = ccs_distribution_strided_samples(distrib2, rng, num_samples, 2, &(samples[0]) + 1);
	assert( err == CCS_SUCCESS );

	ccs_float_t sum = 0.0;
	ccs_float_t inv_sum = 0.0;
	for(size_t i = 0; i < num_areas; i++) {
		sum += areas1[i];
	}
	inv_sum = 1.0 / sum;
	for(size_t i = 0; i < num_samples; i++) {
		assert( samples[2*i].i >=0 && samples[2*i].i < 4 );
		counts1[samples[2*i].i]++;
	}
	for(size_t i = 0; i < num_areas; i++) {
		ccs_float_t target = num_samples * areas1[i] * inv_sum;
		assert( counts1[i] >= target * 0.95 && counts1[i] <= target * 1.05 );
	}

	sum = 0.0;
	inv_sum = 0.0;
	for(size_t i = 0; i < num_areas; i++) {
		sum += areas2[i];
	}
	inv_sum = 1.0 / sum;
	for(size_t i = 0; i < num_samples; i++) {
		assert( samples[2*i+1].i >=0 && samples[2*i+1].i < 4 );
		counts2[samples[2*i+1].i]++;
	}
	for(size_t i = 0; i < num_areas; i++) {
		ccs_float_t target = num_samples * areas2[i] * inv_sum;
		assert( counts2[i] >= target * 0.95 && counts2[i] <= target * 1.05 );
	}

	err = ccs_release_object(distrib1);
	assert( err == CCS_SUCCESS );
	err = ccs_release_object(distrib2);
	assert( err == CCS_SUCCESS );
	err = ccs_release_object(rng);
	assert( err == CCS_SUCCESS );
}

void test_roulette_distribution_soa_samples() {
	ccs_distribution_t distrib = NULL;
	ccs_rng_t          rng = NULL;
	ccs_result_t       err = CCS_SUCCESS;
	const size_t       num_samples = NUM_SAMPLES;
	ccs_numeric_t      samples[NUM_SAMPLES];
	const size_t       num_areas = NUM_AREAS;
	ccs_float_t        areas[NUM_AREAS];
	int                counts[NUM_AREAS];
	ccs_numeric_t     *p_samples;

	for(size_t i = 0; i < num_areas; i++) {
		areas[i] = (double)(i+1);
	        counts[i] = 0;
	}

	err = ccs_rng_create(&rng);
	assert( err == CCS_SUCCESS );
	err = ccs_create_roulette_distribution(
		num_areas,
		areas,
		&distrib);
	assert( err == CCS_SUCCESS );

	p_samples = &(samples[0]);
	err = ccs_distribution_soa_samples(distrib, rng, num_samples, &p_samples);
	assert( err == CCS_SUCCESS );

	ccs_float_t sum = 0.0;
	ccs_float_t inv_sum = 0.0;
	for(size_t i = 0; i < num_areas; i++) {
		sum += areas[i];
	}
	inv_sum = 1.0 / sum;
	for(size_t i = 0; i < num_samples; i++) {
		assert( samples[i].i >=0 && samples[i].i < 4 );
		counts[samples[i].i]++;
	}
	for(size_t i = 0; i < num_areas; i++) {
		ccs_float_t target = num_samples * areas[i] * inv_sum;
		assert( counts[i] >= target * 0.95 && counts[i] <= target * 1.05 );
	}

	err = ccs_release_object(distrib);
	assert( err == CCS_SUCCESS );
	err = ccs_release_object(rng);
	assert( err == CCS_SUCCESS );
}

int main() {
	ccs_init();
	test_create_roulette_distribution();
	test_create_roulette_distribution_errors();
	test_roulette_distribution();
	test_roulette_distribution_zero();
	test_roulette_distribution_strided_samples();
	test_roulette_distribution_soa_samples();
	ccs_fini();
	return 0;
}

max_stars_repo_name: deephyper/CCS

max_issues_repo_name: deephyper/CCS

max_forks_repo_name: deephyper/CCS


text:
/*
* -----------------------------------------------------------------
*  Thermochemistry Library --- thrm_lib.c
*  Version: 1.6180
*  Date: Mar 15, 2010
* -----------------------------------------------------------------
*  Programmer: Americo Barbosa da Cunha Junior
*              americo.cunhajr@gmail.com
* -----------------------------------------------------------------
*  Copyright (c) 2010 by Americo Barbosa da Cunha Junior
*
*  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.
*
*  A copy of the GNU General Public License is available in
*  LICENSE.txt or http://www.gnu.org/licenses/.
* -----------------------------------------------------------------
*  This is the implementation file of a library
*  with thermochemistry routines.
* -----------------------------------------------------------------
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>

#include "../include/cklib.h"
#include "../include/ckthrm.h"
#include "../include/thrm_lib.h"

/*
*------------------------------------------------------------
*   thrm_title
*
*   This function prints the program title on the screen
*
*   last update: Oct 20, 2009
* -----------------------------------------------------------------
*/

void thrm_title()
{
    time_t curtime = time (NULL);
    struct tm *loctime = localtime (&curtime);

    printf("\n");
    printf("=============================================");
    printf("\n Themochemistry Direct Simulation\n");
    printf("\n Author:");
    printf("\n Americo Barbosa da Cunha Junior --- PUC-Rio");
    printf("\n americo.cunhajr@gmail.com\n");
    printf("\n Date:");
    printf("\n %s", asctime(loctime) );
    printf("=============================================");
    printf("\n");

    return;
}

/*
*------------------------------------------------------------
*   thrm_mech
*
*   This function prints on the screen informations about the
*   chemical kinetics mechanism
*
*   Input:
*   n_e  - # of elements
*   n_s  - # of species
*   n_r  - # of reactions
*   mech - mechanism name
*
*   last update: Jul 22, 2009
*------------------------------------------------------------
*/

void thrm_mech(int n_e,int n_s,int n_r,char *mech)
{
    printf("\n");
    printf("\n Thermochemistry Mechanism: %s", mech);
    printf("\n # of elements  = %d", n_e);
    printf("\n # of species   = %d", n_s);
    printf("\n # of reactions = %d", n_r);
    printf("\n");
    printf("\n");

    return;
}
/*------------------------------------------------------------*/

/*
* -----------------------------------------------------------------
*   thrm_input
*
*   This function receives the pmsr parameters
*   and returns GSL_SUCCESS if there is no error.
*
*   Input:
*   Np       - # of particles
*   k        - # of discrete time steps
*   delta_t  - time step
*   tau_res  - residence time
*   tau_mix  - mixture time
*   tau_pair - pairwize time
*
*   Output:
*   success or error
*
*   last update: Feb 5, 2009
* -----------------------------------------------------------------
*/

int thrm_input(int *steps,double *t0,double *delta_t,
                            double *Tmax,double *Tmin,double *tol)
{
    printf("\n Input simulation parameters:\n");

    printf("\n # of integration steps:");
    scanf("%d", steps);
    printf("\n %d\n", *steps);
    if( *steps <= 0 )
	return GSL_EINVAL;

    printf("\n initial time:");
    scanf("%lf", t0);
    printf("\n %+.1e\n", *t0);
    if( *t0 < 0.0 )
	return GSL_EINVAL;

    printf("\n time step:");
    scanf("%lf", delta_t);
    printf("\n %+.1e\n", *delta_t);
    if( *delta_t <= 0.0 )
	return GSL_EINVAL;

    printf("\n temperature upper bound:");
    scanf("%lf", Tmax);
    printf("\n %+.1e\n", *Tmax);
    if( *Tmax < 0.0 )
	return GSL_EINVAL;

    printf("\n temperature lower bound:");
    scanf("%lf", Tmin);
    printf("\n %+.1e\n", *Tmin);
    if( *Tmin < 0.0 || *Tmin > *Tmax )
	return GSL_EINVAL;

    printf("\n tolerance:");
    scanf("%lf", tol);
    printf("\n %+.1e\n", *tol);
    if( *tol <= 0.0 )
	return GSL_EINVAL;

    return GSL_SUCCESS;
}
/*----------------------------------------------------------------*/

/*
*------------------------------------------------------------
*   thrm_alloc
*
*   This function allocates a struct thermochemistry.
*
*   Output:
*   thrm - pointer to a struct thermochemistry
*
*   last update: May 10, 2009
*------------------------------------------------------------
*/

thrm_wrk* thrm_alloc()
{
    thrm_wrk *thrm = NULL;

    /* memory allocation for ck_wrk struct */
    thrm = (thrm_wrk *) malloc(sizeof(thrm_wrk));
    if ( thrm == NULL )
        return NULL;

    /* setting ck_wrk elements equal NULL and 0.0 */
    thrm->iwk    = NULL;
    thrm->rwk    = NULL;
    thrm->cwk    = NULL;
    thrm->leniwk = 0;
    thrm->lenrwk = 0;
    thrm->lencwk = 0;
    thrm->n_e    = 0;
    thrm->n_s    = 0;
    thrm->n_r    = 0;

    return thrm;
}
/*------------------------------------------------------------*/

/*
*------------------------------------------------------------
*   thrm_free
*
*   This function frees the memory used by CK work vectors
*
*   Input:
*   thrm - pointer to a struct chemkin
*
*   Output:
*   void
*
*   last update: Jul 21, 2009
*------------------------------------------------------------
*/

void thrm_free(void **thrm_bl)
{
    thrm_wrk* thrm = NULL;

    /* checking if thrm_bl is NULL */
    if ( *thrm_bl == NULL )
        return;

    thrm = (thrm_wrk *) (*thrm_bl);

    /* releasing ck work vectors */
    if( thrm->iwk != NULL )
    {
        free(thrm->iwk);
        thrm->iwk = NULL;
    }

    if( thrm->rwk != NULL )
    {
        free(thrm->rwk);
        thrm->rwk = NULL;
    }

    if( thrm->cwk != NULL )
    {
        free(thrm->cwk);
        thrm->cwk = NULL;
    }

    /* releasing thrm struct */
    free(*thrm_bl);
    *thrm_bl = NULL;

    return;
}
/*------------------------------------------------------------*/

/*
*------------------------------------------------------------
*   thrm_init
*
*   This function initiates chemkin work vectors and the
*   mechanism parameters.
*
*   Input:
*   thrm - pointer to a struct thermochemistry
*
*   Output:
*   success or error
*
*   last update: Feb 19, 2010
*------------------------------------------------------------
*/

int thrm_init(thrm_wrk *thrm)
{
    if ( thrm == NULL )
        return GSL_EINVAL;

    /* calling ck_wrk_len subroutine */
    ck_wrk_len(&(thrm->leniwk),&(thrm->lenrwk),&(thrm->lencwk));

    /* memory allocation for ck work vectors */
    thrm->iwk =    (int *) malloc(   (thrm->leniwk)*sizeof(int));
    thrm->rwk = (double *) malloc(   (thrm->lenrwk)*sizeof(double));
    thrm->cwk =   (char *) malloc(16*(thrm->lencwk)*sizeof(char));
    if ( thrm->iwk == NULL || thrm->rwk == NULL || thrm->cwk == NULL )
    {
        thrm_free((void **)&thrm);
        return GSL_ENOMEM;
    }

    /* calling ck_wrk_init subroutine */
    ck_wrk_init(&(thrm->leniwk),&(thrm->lenrwk),&(thrm->lencwk),
        thrm->iwk,thrm->rwk,thrm->cwk);

    /* calling ck_wrk_indx subroutine */
    ck_wrk_indx(thrm->iwk,thrm->rwk,&(thrm->n_e),&(thrm->n_s),&(thrm->n_r));

    return GSL_SUCCESS;
}
/*------------------------------------------------------------*/

/*
*------------------------------------------------------------
*   thrm_composition
*
*   This function sets composition vector with values.
*
*   Input:
*   thrm - pointer to a struct thermochemistry
*
*   Output:
*   phi - composition vector
*   success or error
*
*   last update: Mar 2, 2010
*------------------------------------------------------------
*/

void thrm_composition(thrm_wrk *thrm,gsl_vector *phi)
{
    /* calling ck_wrk_comp subroutine */
    ck_composition(&(thrm->n_s),thrm->iwk,thrm->rwk,thrm->cwk,phi->data);

    return;
}
/*------------------------------------------------------------*/

/*
*------------------------------------------------------------
*   thrm_page_heading0
*
*   This function prints the properties names on the screen.
*
*   Input:
*   thrm - pointer to a struct thermochemistry
*
*   last update: Nov 6, 2009
*------------------------------------------------------------
*/

void thrm_page_heading0(FILE *file)
{
    fprintf(file," t(s)      ");
    fprintf(file," T(K)      ");
    fprintf(file," h         ");
    fprintf(file," p         ");
    fprintf(file," CO2       ");
    fprintf(file," O2        ");
    fprintf(file," CO        ");
    fprintf(file," O         ");

    return;
}
/*------------------------------------------------------------*/

/*
*------------------------------------------------------------
*   thrm_h2T
*
*   This function computes the mixture temperature for a
*   given enthalpy, pressure and species mass fraction.
*
*   Input:
*   ck  - pointer to a struct chemkin
*   phi - composition vector
*
*   Output:
*   T   - mixture temperature
*
*   last update: Mar 2, 2010
*------------------------------------------------------------
*/

double thrm_h2T(thrm_wrk *thrm,double *phi,double Tmax,double Tmin,double tol)
{
    double T = 0.0;

    ck_h2t_bisec(&(thrm->n_s),phi,&Tmax,&Tmin,&tol,thrm->iwk,thrm->rwk,&T);
    /*ck_h2t_newt(&(thrm->n_s),phi,&tol,thrm->iwk,thrm->rwk,&T,&T);*/

    return T;
}
/*------------------------------------------------------------*/

/*
* -----------------------------------------------------------------
*   thrm_temp_meanvar
*
*   This function computes the temperature mean and variance
*   using the algorithm of Knuth/Welford.
*
*   Ref.:
*   D. E. Knuth (1998)
*   The Art of Computer Programming,
*   volume 2: Seminumerical Algorithms
*   3rd ed., Addison-Wesley.
*
*   B. P. Welford (1962).
*   "Note on a method for calculating corrected sums
*   of squares and products". Technometrics 4(3):419–420
*
*   Input:
*   Np   - # of particles
*   ps   - system of particle
*   thrm - pointer to a struct thermochemistry
*   Tmax - temperature upper bound
*   Tmin - temperature lower bound
*   tol  - tolerance
*   mean - temperature mean
*   var  - temperature variance
*
*   Output:
*   success or error
*
*   last update: Mar 4, 2010
* -----------------------------------------------------------------
*/

void thrm_temp_meanvar(int Np,gsl_vector **ps,thrm_wrk *thrm,
        double Tmax,double Tmin,double tol,double *mean,double *var)
{
    int i;
    double T     = 0.0;
    double S     = 0.0;
    double M     = 0.0;
    double delta = 0.0;

    for ( i = 0; i < Np; i++ )
    {
        T     = thrm_h2T(thrm,ps[i]->data,Tmax,Tmin,tol);
	delta = T - M;
	M    += delta/(double) (i+1);
        S    += delta*(T - M);
    }

    /* mean */
    *mean  = M;

    /* variance */
    *var   = S/(double) (i-1);

    return;
}
/*----------------------------------------------------------------*/

/*
*------------------------------------------------------------
*   thrm_rrsum
*
*   This function computes the reaction rates sum given the
*   reaction rate vector.
*
*   Input:
*   phi - composition vector
*
*   Output:
*   sum - reaction rates sum
*
*   last update: Sep 28, 2009
*------------------------------------------------------------
*/

double thrm_rrsum(gsl_vector *phi)
{
    unsigned int i;
    double sum = 0.0;

    for (i = 2; i < phi->size; i++)
        sum += phi->data[i];

    return sum;
}
/*------------------------------------------------------------*/

/*
*------------------------------------------------------------
*   thrm_mfsum
*
*   This function computes the mass fraction sum given a
*   composition vector.
*
*   Input:
*   phi - composition vector
*
*   Output:
*   sum - mass fraction sum
*
*   last update: Sep 24, 2009
*------------------------------------------------------------
*/

double thrm_mfsum(gsl_vector *phi)
{
    unsigned int i;
    double sum = 0.0;

    for (i = 2; i < phi->size; i++)
        sum += phi->data[i];

    return sum;
}
/*------------------------------------------------------------*/

/*
*------------------------------------------------------------
*   thrm_mfsign
*
*   This function checks if the mass fraction
*   components are all positive in a composition vector.
*
*   Input:
*   phi  - composition vector
*
*   Output:
*   flag - success or fail
*
*   last update: Sep 28, 2009
*------------------------------------------------------------
*/

int thrm_mfsign(gsl_vector *phi)
{
    unsigned int i;
    int flag = GSL_SUCCESS;

    for (i = 2; i < phi->size; i++)
        if ( phi->data[i] < 0.0 )
            flag = GSL_EFAILED;

    return flag;
}
/*------------------------------------------------------------*/

/*
*------------------------------------------------------------
*   thrm_eqs
*
*   This function defines the set of governing diferential
*   equations, one for enthalpy, one for pressure and ns
*   for species mass fractions.
*
*   Input:
*   t         - instant time
*   phi       - composition vector
*   Sphi      - reaction rate vector
*   thrm_data - pointer to external data
*
*   Output:
*   success or error
*
*   last update: Feb 19, 2010
*------------------------------------------------------------
*/

int thrm_eqs(double t,N_Vector phi,N_Vector Sphi,void *thrm_data)
{
    int i, Neq;
    double tol = 1.0e-4;
    double Tmax = 5000.0, Tmin = 100.0;
    double p, T, rho;
    double *phi_data  = NULL;
    double *Sphi_data = NULL;
    double *wdot      = NULL;
    double *wt        = NULL;
    thrm_wrk *thrm    = NULL;

    /* initiating ck workspace */
    thrm = (thrm_wrk *)thrm_data;
    if ( thrm == NULL )
        return GSL_EINVAL;

    /* obtaining phi and Rphi length and components */
    Neq       = NV_LENGTH_S(phi);
    phi_data  = NV_DATA_S(phi);
    Sphi_data = NV_DATA_S(Sphi);

    /* memory allocation for species molar mass and reaction rate vectors */
    wt   = (double *) malloc((Neq-2)*sizeof(double));
    wdot = (double *) malloc((Neq-2)*sizeof(double));
    if( wt == NULL || wdot == NULL )
    {
        free(wt);
        free(wdot);
        return GSL_ENOMEM;
    }

    /* pressure (dynes/cm^2) */
    p = phi_data[1];

    /* mixture temperature (K) */
    T = thrm_h2T(thrm,phi_data,Tmax,Tmin,tol);
    if( T < tol )
        return GSL_EFAILED;

    /* mixture density (g/cm^3)*/
    ckrhoy(&p,&T,&(phi_data[2]),thrm->iwk,thrm->rwk,&rho);

    /* species molar weight (g/mol) */
    ckwt(thrm->iwk,thrm->rwk,wt);

    /* species molar production rates (mol/(cm^3.s)) */
    ckwyp(&p,&T,&(phi_data[2]),thrm->iwk,thrm->rwk,wdot);

    /* dh/dt = 0 (isoenthalpic) */
    Sphi_data[0] = 0.0;

    /* dp/dt = 0 (isobaric) */
    Sphi_data[1] = 0.0;

    /* dY[i]/dt = wdot[i-2].wt[i-2]/rho */
    for ( i = 2; i < Neq; i++ )
        Sphi_data[i] = wdot[i-2]*wt[i-2]/rho;

    /* releasing allocated memory */
    free(wdot);
    free(wt);
    phi_data  = NULL;
    Sphi_data = NULL;
    wdot      = NULL;
    wt        = NULL;
    thrm      = NULL;

    return GSL_SUCCESS;
}
/*------------------------------------------------------------*/

/*
*------------------------------------------------------------
*   thrm_eqs2
*
*   This function defines the set of governing diferential
*   equations, one for temperature and ns
*   for species mass fractions.
*
*   Input:
*   t         - instant time
*   phi       - state vector
*   Sphi      - reaction rate vector
*   thrm_data - pointer to external data
*
*   Output:
*   success or error
*
*   last update: Dec 28, 2010
*------------------------------------------------------------
*/

int thrm_eqs2(double t,N_Vector phi,N_Vector Sphi,void *thrm_data)
{
    int i, Neq;
    double T, p, rho, cpb, sum = 0.0;
    double *phi_data  = NULL;
    double *Sphi_data = NULL;
    double *wdot      = NULL;
    double *wt        = NULL;
    double *h         = NULL;
    thrm_wrk *thrm    = NULL;

    /* initiating ck workspace */
    thrm = (thrm_wrk *)thrm_data;
    if ( thrm == NULL )
        return GSL_EINVAL;

    /* obtaining phi and Rphi length and components */
    Neq       = NV_LENGTH_S(phi);
    phi_data  = NV_DATA_S(phi);
    Sphi_data = NV_DATA_S(Sphi);

    /* memory allocation for species molar mass and reaction rate vectors */
    wt   = (double *) malloc((Neq-2)*sizeof(double));
    wdot = (double *) malloc((Neq-2)*sizeof(double));
    h    = (double *) malloc((Neq-2)*sizeof(double));
    if( wt == NULL || wdot == NULL || h == NULL)
    {
        free(wt);
        free(wdot);
	free(h);
        return GSL_ENOMEM;
    }

    /* mixture temperature (K) */
    T = phi_data[0];

    /* pressure (dynes/cm^2) */
    p = phi_data[1];

    /* mixture density (g/cm^3)*/
    ckrhoy(&p,&T,&(phi_data[2]),thrm->iwk,thrm->rwk,&rho);

    /* mean specific heat at constant pressure (ergs/g.K) */
    ckcpbs(&T,&(phi_data[2]),thrm->iwk,thrm->rwk,&cpb);

    /* species molar weight (g/mol) */
    ckwt(thrm->iwk,thrm->rwk,wt);

    /* species molar production rates (mol/(cm^3.s)) */
    ckwyp(&p,&T,&(phi_data[2]),thrm->iwk,thrm->rwk,wdot);

    /* species specific enthalpy (ergs/g) */
    ckhms(&T,thrm->iwk,thrm->rwk,h);

    /* dY[i+2]/dt = wdot[i].wt[i]/rho */
    for ( i = 0; i < thrm->n_s; i++ )
    {
        Sphi_data[i+2] = wdot[i]*wt[i]/rho;
	sum += h[i]*wdot[i]*wt[i];
    }

    /* dT/dt =  -1/(rho.cpb) . sum_i=1^ns h[i].wdot[i].wt[i] */
    Sphi_data[0] = -sum/(rho*cpb);

    /* dp/dt = 0 (isobaric) */
    Sphi_data[1] = 0.0;

    /* releasing allocated memory */
    free(wdot);
    free(wt);
    free(h);
    phi_data  = NULL;
    Sphi_data = NULL;
    wdot      = NULL;
    wt        = NULL;
    h         = NULL;
    thrm      = NULL;

    return GSL_SUCCESS;
}
/*------------------------------------------------------------*/

max_stars_repo_name: americocunhajr/CRFlowLib

max_issues_repo_name: americocunhajr/CRFlowLib

max_forks_repo_name: americocunhajr/CRFlowLib


text:
// randomizer.h
// (C) 2013-2020 Nicholas G Davies

#ifndef RANDOMIZER_H
#define RANDOMIZER_H

#include <vector>
#include <limits>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <cmath>
#include <stdexcept>
#include <string>

class Randomizer
{
public:
    Randomizer(unsigned long int seed = 0);
    Randomizer(const Randomizer& R);
    ~Randomizer();

    void Reset();

    double Uniform(double min = 0.0, double max = 1.0);
    double RoundedUniform(double min = 0.0, double max = 1.0, double shoulder = 0.01);
    double Normal(double mean = 0.0, double sd = 1.0);
    double Normal(double mean, double sd, double clamp);
    double Cauchy(double x0 = 0.0, double gamma = 1.0);
    double LogNormal(double zeta = 0.0, double sd = 1.0);
    double Exponential(double rate = 1.0);
    double Gamma(double shape, double scale);
    double Beta(double alpha, double beta);
    unsigned int Discrete(unsigned int size);
    int Discrete(int min, int max);
    void Multinomial(unsigned int N, std::vector<double>& p, std::vector<unsigned int>& n);
    bool Bernoulli(double p);
    unsigned int Binomial(unsigned int n, double p);
    unsigned int BetaBinomial(unsigned int n, double p, double a_plus_b);
    int Poisson(double mean);
    int Geometric(double p);
    int Round(double x);

    template <typename T>
    void Shuffle(std::vector<T>& vec);

    inline gsl_rng* GSL_RNG() const { return r; }
    inline unsigned long int Seed() const { return seed; }

private:
    unsigned long int seed;
    gsl_rng* r;
};

template <typename T>
void Randomizer::Shuffle(std::vector<T>& vec)
{
    gsl_ran_shuffle(r, &vec[0], vec.size(), sizeof(vec[0]));
}

#endif

max_stars_repo_name: nicholasdavies/covid-tiers


text:
// Copyright (c) 2005 - 2010 Marc de Kamps
// All rights reserved.
//
// Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
//
//    * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
//    * Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation
//      and/or other materials provided with the distribution.
//    * Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software
//      without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY
// DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
// USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
//      If you use this software in work leading to a scientific publication, you should cite
//      the 'currently valid reference', which can be found at http://miind.sourceforge.net
#ifndef _CODE_LIBS_NUMTOOLSLIB_NODE_EXSTATEDVINTEGRATOR_INCLUDE_GUARD
#define _CODE_LIBS_NUMTOOLSLIB_NODE_EXSTATEDVINTEGRATOR_INCLUDE_GUARD

#include <limits>
#include <vector>
#include <gsl/gsl_odeiv.h>
#include <gsl/gsl_errno.h>
#include "AbstractDVIntegratorCode.h"
#include "BasicDefinitions.h"
#include "DVIntegratorException.h"
#include "DVIntegratorStateParameter.h"

using std::vector;
using std::istream;
using std::ostream;
using std::string;

using std::numeric_limits;

namespace NumtoolsLib
{

	//! This is a wrapper class for the GSL numerical ODE solvers. Examples of its usage can be
	//! found in DynamicLib, for example in the Wilson-Cowan Algorithm. Unlike DVIntegrator it operates on a
	//! state vector that is maintained by another class.
	//!
	//! The ExStateDVIntegrator is typically member of a class which requires the solution of a system
	//! of ODEs. In the source file of that class a function must be defined whose return type is int
	//! and whose arguments are given in the typedef for Function. A pointer to this function, (i.e the
	//! function name must be given in the constructor. If the GSL
	//! solver requires a Jacobian, then Derivative must also be defined, but for simple solvers,
	//! such as 4-th order Runge-Kutta this is not required and the derivative argument can be left out.
	//! Often, the ODE system requires parameters that need to be updated during the evolution of the system.
	//! This parameter can be a single number, a vector or anything else. The type of the parameter must
	//! be provided using the template argument. The class of the parameter must have a default constructor.
	template <class ParameterObject=double>
	class ExStateDVIntegrator : public AbstractDVIntegrator<ParameterObject>
	{
	public:

		//! Constructor
		ExStateDVIntegrator
		(
			Number,						//<! maximum number of integrations
			double*,					//<! initial state
			Number,						//<! state size
			TimeStep,					//<! initial time step
			Time,						//<! initial time
			const Precision&,			//<! absolute and relative precision
			Function,					//<! dydt
			Derivative deravative = 0,  //<! Jacobian, if available
			const gsl_odeiv_step_type*	p_step_algorithm = gsl_odeiv_step_rkf45 //<! gsl integrator object
		);

		//! copy constructor must be defined, because the state of the integrator
		//! must be handled correctly
		ExStateDVIntegrator(const ExStateDVIntegrator&);

		//! Ditto for copy operator
		ExStateDVIntegrator& operator=(const ExStateDVIntegrator&);

		//! virtual destructor
		virtual ~ExStateDVIntegrator();

		//! Set a new state
		bool Reconfigure
		(
			const DVIntegratorStateParameter<ParameterObject>&
		);

		//! Current time
		Time CurrentTime() const;

		//! Only use if for efficiency reasons, acces to internal state is necessary
		const double* BeginState() const;

		//! Ditto
		const double* EndState  () const;

		//! streaming tag
		virtual string Tag () const;

		//! Direct access to the parameter object for performance reasons
		ParameterObject& Parameter();

	private:

	}; // end of NodeIntegrator

} // end of namespace

#endif // include guard

max_stars_repo_name: dekamps/miind

max_issues_repo_name: dekamps/miind

max_forks_repo_name: dekamps/miind


text:
// --------------------------------------------------------------------------
//                   OpenMS -- Open-Source Mass Spectrometry
// --------------------------------------------------------------------------
// Copyright The OpenMS Team -- Eberhard Karls University Tuebingen,
// ETH Zurich, and Freie Universitaet Berlin 2002-2013.
//
// This software is released under a three-clause BSD license:
//  * Redistributions of source code must retain the above copyright
//    notice, this list of conditions and the following disclaimer.
//  * Redistributions in binary form must reproduce the above copyright
//    notice, this list of conditions and the following disclaimer in the
//    documentation and/or other materials provided with the distribution.
//  * Neither the name of any author or any participating institution
//    may be used to endorse or promote products derived from this software
//    without specific prior written permission.
// For a full list of authors, refer to the file AUTHORS.
// --------------------------------------------------------------------------
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL ANY OF THE AUTHORS OR THE CONTRIBUTING
// INSTITUTIONS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
// OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
// WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
// OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
// ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// --------------------------------------------------------------------------
// $Maintainer: Lars Nilse $
// $Authors: Steffen Sass, Holger Plattfaut, Bastian Blank $
// --------------------------------------------------------------------------

#ifndef OPENMS_FILTERING_DATAREDUCTION_SILACFILTERING_H
#define OPENMS_FILTERING_DATAREDUCTION_SILACFILTERING_H

#include <OpenMS/KERNEL/StandardTypes.h>
#include <OpenMS/KERNEL/MSExperiment.h>
#include <OpenMS/CONCEPT/ProgressLogger.h>
#include <OpenMS/DATASTRUCTURES/DRange.h>
#include <OpenMS/TRANSFORMATIONS/FEATUREFINDER/PeakWidthEstimator.h>

#include <gsl/gsl_interp.h>
#include <gsl/gsl_spline.h>
#include <list>
#include <map>
#include <vector>

namespace OpenMS
{
  class SILACFilter;

  /**
    @brief Filtering for SILAC data.

    This filtering can be used to extract SILAC features from an MS experiment.
    Several SILACFilters can be added to the filtering to search for specific SILAC patterns.

    @see SILACFilter
  */
  class OPENMS_DLLAPI SILACFiltering :
    public ProgressLogger
  {
public:
    typedef std::vector<SILACFilter> Filters;

    /**
     * @brief holds all filters used in the filtering
     */
    Filters filters_;

    /**
     * @brief Wrapper class for spectrum interpolation
     */
    class OPENMS_DLLAPI SpectrumInterpolation
    {
private:
      gsl_interp_accel * current_;
      gsl_spline * spline_;

public:
      SpectrumInterpolation(const MSSpectrum<> &, const SILACFiltering &);
      ~SpectrumInterpolation();

      DoubleReal operator()(DoubleReal mz) const
      {
        return gsl_spline_eval(spline_, mz, current_);
      }

    };

private:
    /**
     * @brief minimal intensity of SILAC features
     */
    DoubleReal intensity_cutoff_;

    /**
     * @brief raw data
     */
    MSExperiment<Peak1D> & exp_;

    /**
     * @brief picked data
     */
    MSExperiment<Peak1D> picked_exp_;

    /**
     * @brief picked data seeds
     */
    MSExperiment<Peak1D> picked_exp_seeds_;

    /**
     * Filename base for debugging output
     */
    const String debug_filebase_;

    /**
     * @brief pick data seeds
     */
    void pickSeeds_();

    /**
     * @brief apply filtering to picked data seeds
     */
    void filterSeeds_();

public:

    /// peak-width equation
    const PeakWidthEstimator::Result peak_width;

    /**
     * @brief detailed constructor
     * @param exp raw data
     * @param intensity_cutoff minimal intensity of SILAC features
     */
    SILACFiltering(MSExperiment<Peak1D> & exp, const PeakWidthEstimator::Result &, const DoubleReal intensity_cutoff, const String debug_filebase_ = "");

    /**
     * @brief adds a new filter to the filtering
     * @param filter filter to add
     */
    void addFilter(SILACFilter & filter);

    /**
     * @brief starts the filtering based on the added filters
     */
    void filterDataPoints();

    /**
     * @brief structure for blacklist
     * @param range m/z and RT interval to be blacklisted
     * @param charge charge of the generating filter
     * @param mass_separations mass separations of the generating filter
     * @param relative_peak_position m/z position of the blacklisted area relative to the mono-isotopic peak of the unlabelled peptide
     */
    struct BlacklistEntry
    {
      DRange<2> range;
      Int charge;
      std::vector<DoubleReal> mass_separations;
      DoubleReal relative_peak_position;
    };

    /**
     * @brief holds the range that is blacklisted for other filters and the filter that generated the blacklist entry
     */
    std::multimap<DoubleReal, BlacklistEntry> blacklist;
  };
}

#endif /* SILACFILTERING_H_ */

max_stars_repo_name: open-ms/all-svn-branches

max_issues_repo_name: kreinert/OpenMS

max_forks_repo_name: kreinert/OpenMS


text:
/*
** Implementation of LISA algorithm
** for statistical inference of fMRI images
**
** generic plug-in.
** A file containing a list of all 3D permutation images must be supplied.
**
** G.Lohmann, April 2017
*/

#include <viaio/Vlib.h>
#include <viaio/file.h>
#include <viaio/mu.h>
#include <viaio/option.h>
#include <viaio/os.h>
#include <viaio/VImage.h>
#include <via/via.h>

#include <gsl/gsl_cdf.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_histogram.h>

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

#ifdef _OPENMP
#include <omp.h>
#endif /*_OPENMP*/

#define ABS(x) ((x) > 0 ? (x) : -(x))

extern void VIsolatedVoxels(VImage src,float threshold);
extern void VHistogram(gsl_histogram *histogram,VString filename);
extern void VCheckImage(VImage src);
extern void FDR(VImage src,VImage dest,gsl_histogram *nullhist,gsl_histogram *realhist,double alpha);
extern void ImageStats(VImage src,double *,double *,double *hmin,double *hmax);
extern void VBilateralFilter(VImage src,VImage dest,int radius,double var1,double var2,int);
extern double VImageVar(VImage src);
extern void  VImageCount(VImage src);
extern void  VGetHistRange(VImage src,double *hmin,double *hmax);
extern float VGetMode(VImage src);
extern void  VZScale(VImage src,float mode,float stddev);
extern void HistoUpdate(VImage,gsl_histogram *);

/* make sure all input images are in float and have the same number of pixels */
void CheckImageTypes(VImage zmap,VImage *permimages,int numperm)
{
  size_t npixels = VImageNPixels(zmap);
  int i;
  for (i=0; i<numperm; i++) {
    if (npixels != VImageNPixels(permimages[i]))
      VError(" inconsistent number of pixels in permutation image %d",i);
    if (VPixelRepn(permimages[i]) != VFloatRepn) VError(" perm image %d is not in float repn");
  }
}

int main (int argc, char *argv[])
{
  static VString  filename = "";
  static VFloat   alpha = 0.05;
  static VShort   radius = 2;
  static VFloat   rvar = 2.0;
  static VFloat   svar = 2.0;
  static VShort   numiter = 2;
  static VBoolean centering = FALSE;
  static VBoolean cleanup = TRUE;
  static VShort   nproc = 0;
  static VOptionDescRec  options[] = {
    {"permutations",VStringRepn,1,(VPointer) &filename,VRequiredOpt,NULL,"List of all permutation images"},
    {"alpha",VFloatRepn,1,(VPointer) &alpha,VOptionalOpt,NULL,"FDR significance level"},
    {"radius",VShortRepn,1,(VPointer) &radius,VOptionalOpt,NULL,"Bilateral parameter (radius in voxels)"},
    {"rvar",VFloatRepn,1,(VPointer) &rvar,VOptionalOpt,NULL,"Bilateral parameter (radiometric)"},
    {"svar",VFloatRepn,1,(VPointer) &svar,VOptionalOpt,NULL,"Bilateral parameter (spatial)"},
    {"filteriterations",VShortRepn,1,(VPointer) &numiter,VOptionalOpt,NULL,"Bilateral parameter (number of iterations)"},
    {"cleanup",VBooleanRepn,1,(VPointer) &cleanup,VOptionalOpt,NULL,"Whether to remove isolated voxels"},
    {"j",VShortRepn,1,(VPointer) &nproc,VOptionalOpt,NULL,"Number of processors to use, '0' to use all"},
  };
  FILE *out_file=NULL,*in_file=NULL;
  VString in_filename=NULL;
  char *prg_name=GetLipsiaName("vlisa_generic");
  fprintf (stderr, "%s\n", prg_name);
  gsl_set_error_handler_off ();

  /* Parse command line arguments and identify files: */
  VParseFilterCmdZ (VNumber (options), options, argc, argv,&in_file,&out_file,&in_filename);

  /* Read the input zmap file: */
  VAttrList list = VReadAttrListZ(in_file,in_filename,0L,TRUE,FALSE);
  VImage zmap1 = VReadImage(list);
  if (zmap1 == NULL) VError(" no input zmap image found");
  if (VPixelRepn(zmap1) != VFloatRepn) VError(" input pixel repn must be float");
  VAttrList geolist = VGetGeoInfo(list);

  /* Read permutation file containing a list of 3D images */
  VAttrList listperm = VReadAttrList(filename,0L,FALSE,FALSE);
  int numperm  = VAttrListNumImages(listperm);
  VImage *zmap = VAttrListGetImages(listperm,numperm);
  CheckImageTypes(zmap1,zmap,numperm);
  fprintf(stderr," Number of permutation images: %d\n",(int)numperm);

  /* estimate null variance to adjust radiometric parameter, use first 30 permutations */
  double zvar=0,hmin=0,hmax=0;
  float stddev=1.0;
  int nperm=0;
  if (numperm > 0) {
    int tstperm = 30;
    if (tstperm > numperm) tstperm = numperm;
    double varsum=0,nx=0;
    for (nperm = 0; nperm < tstperm; nperm++) {
      zvar = VImageVar(zmap[nperm]);
      varsum += zvar;
      nx++;
    }
    double meanvar = varsum/nx;
    stddev = sqrt(meanvar);
  }

  /* get non-permuted hotspot map */
  float mode=0;
  if (centering) mode = VGetMode(zmap1);
  if (numperm > 0) VZScale(zmap1,mode,stddev);

  VImage dst1 = VCreateImageLike(zmap1);
  VBilateralFilter(zmap1,dst1,(int)radius,(double)rvar,(double)svar,(int)numiter);

  /* ini histograms */
  VGetHistRange(dst1,&hmin,&hmax);
  size_t nbins = 20000;
  gsl_histogram *hist0 = gsl_histogram_alloc (nbins);
  gsl_histogram_set_ranges_uniform (hist0,hmin,hmax);
  gsl_histogram *histz = gsl_histogram_alloc (nbins);
  gsl_histogram_set_ranges_uniform (histz,hmin,hmax);
  HistoUpdate(dst1,histz);

  /* omp-stuff */
#ifdef _OPENMP
  int num_procs=omp_get_num_procs();
  if (nproc > 0 && nproc < num_procs) num_procs = nproc;
  fprintf(stderr," using %d cores\n",(int)num_procs);
  omp_set_num_threads(num_procs);
#endif /* _OPENMP */

  /* do random permutations */
#pragma omp parallel for shared(zmap) schedule(dynamic)
  for (nperm = 0; nperm < numperm; nperm++) {
    if (nperm%20 == 0) fprintf(stderr," perm  %4d  of  %d\r",nperm,(int)numperm);

    float mode=0;
    if (centering) mode = VGetMode(zmap[nperm]);
    VZScale(zmap[nperm],mode,stddev);

    VImage dst = VCreateImageLike (zmap1);
    VBilateralFilter(zmap[nperm],dst,(int)radius,(double)rvar,(double)svar,(int)numiter);

    #pragma omp critical
    {
      HistoUpdate(dst,hist0);
    }
    VDestroyImage(dst);
  }

  /* apply fdr */
  VImage fdrimage = VCopyImage (dst1,NULL,VAllBands);
  if (numperm > 0) {
    FDR(dst1,fdrimage,hist0,histz,(double)alpha);

    if (cleanup && alpha < 1.0) {
      VIsolatedVoxels(fdrimage,(float)(1.0-alpha));
    }
  }

  /* output */
  VAttrList out_list = VCreateAttrList ();
  VHistory(VNumber(options),options,prg_name,&list,&out_list);
  VSetGeoInfo(geolist,out_list);
  VAppendAttr (out_list,"image",NULL,VImageRepn,fdrimage);
  if (! VWriteFile (out_file, out_list)) exit (1);
  fprintf (stderr, "\n%s: done.\n", argv[0]);
  exit(0);
}

max_stars_repo_name: zrajna/lipsia

max_issues_repo_name: zrajna/lipsia

max_forks_repo_name: zrajna/lipsia


text:
#include <stdio.h>
#include <cblas.h>

int
main ()
{
  int i;
  double x[] = { 1, 1, 1 };

  cblas_dscal(3, 4.323, x, 1);

  for (i = 0; i < 3; ++i)
    printf ("%f\n", x[i]);

  return 0;
}

max_stars_repo_name: gkiar/fuzzy-python

max_issues_repo_name: gkiar/fuzzy-python

max_forks_repo_name: gkiar/fuzzy-python


text:
#ifndef _RNG_INTERNAL_H
#define _RNG_INTERNAL_H

#include <gsl/gsl_rng.h>

struct _ccs_rng_data_s {
	const gsl_rng_type *rng_type;
	gsl_rng            *rng;
};

typedef struct _ccs_rng_data_s _ccs_rng_data_t;

struct _ccs_rng_ops_s {
	_ccs_object_ops_t obj_ops;
};
typedef struct _ccs_rng_ops_s _ccs_rng_ops_t;

struct _ccs_rng_s {
	_ccs_object_internal_t obj;
	_ccs_rng_data_t *data;
};

#endif //_RNG_INTERNAL_H