text:
/*
# Program to run a Monte Carlo radiation transfer through the 2D
# simulations of GRB jets.
#
# Python code written by D. Lazzati at Oregonstate, C code written by Tyler Parsotan @ Oregon State
# ver 0.1 July 8, 2015
# ver 1.1 July 20, 2015: added record of number of scatterings, included
# 	all terms in weight. Should now give correct light curves.
# ver 1.2 July 21, 2015: added parameter file to keep track of input
# 	params of each simulation

# ver 2.0 July 22, 2015: corrected the problem that arises when there is
# 	no scattering in the time span of one frame. Fixed output arrays dimension.

# ver 2.1 July 25, 2015: fixed bug that did not make the number of
# 	scattering grow with the number of photons.

# ver 3.0 July 28, 2015: using scipy nearest neighbor interpolation to
# 	speed things up. Gained about factor 2

# ver 3.1 July 29, 2015: added radial spread of photon injection points
# ver 3.2 July 31, 2015: added Gamma to the weight of photons!!!

# ver 4.0 Aug 5, 2015: try to speed up by inverting cycle
# ver 4.1 Aug 8, 2015: add spherical test as an option
# ver 4.2 Aug 9, 2015: saving files appending rather than re-writing
# ver 4.3 Aug 11, 2015: corrected error in the calculation of the local temperature
# ver 4.4 Aug 13, 2015: added cylindrical test
# ver 4.5 Aug 18, 2015: fixd various problems pointed by the cylindrical test
# ver 4.6 Aug 21, 2015: corrected mean free path for large radii

# ver 5.0 Aug 25, 2015: corrected problem with high-T electrons and excess scatterings
# ver 5.1 Aug 25, 2015: cleaned-up coding
# ver 5.2 Sept 3, 2015: fixed problem with number of scatterings for multiple injections
 *
 * ver 6.0 Dec 28, 2016: rewrote the code in C, added checkpoint file so if the code is interrupted all the progress wont be lost, made the code only need to be compiled once for a given MC_XXX directory path
                                            so you just need to supply the sub directory of MC_XXX as a command line argument
*/

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <unistd.h>
#include <dirent.h>
#include <math.h>
#include <gsl/gsl_rng.h>
//#include "mclib.h"
#include "mclib_3d.h"
#include <omp.h>

#define THISRUN "Science"
#define FILEPATH "/Users/Tylerparsotan/Documents/Box\ Sync/RIKEN_HYDRO_DATA/JP_HYDRODATA/"
#define FILEROOT "u"
#define MC_PATH "PHOTON_TEST/"
//#define MC_PATH "MC_16OI/Single_Photon_Cy_mc_total/"
#define MCPAR "riken_mc.par"

int main(int argc, char **argv)
{
    //compile each time a macro is changed

	// Define variables
	char flash_prefix[200]="";
	char mc_file[200]="" ;
    char this_run[200]=THISRUN;
    char *cyl="Cylindrical";
    char *sph="Spherical";
    char spect;//type of spectrum
    char restrt;//restart or not
    double fps, theta_jmin, theta_jmax ;//frames per second of sim, min opening angle of jet, max opening angle of jet in radians
    double inj_radius_small, inj_radius_large,  ph_weight_suggest ;//radius at chich photons are injected into sim
    int frm0_small, frm0_large,last_frm, frm2_small, frm2_large, j=0, min_photons, max_photons ;//frame starting from, last frame of sim, frame of last injection

    int num_thread=0, angle_count=0;
    int half_threads=floor((num_thread/2));
    int num_angles=0;
    int dim_switch=0;
    double *thread_theta=NULL; //saves ranges of thetas for each thread to go through
    double delta_theta=0;

    //new muliple threads injecting and propagating photons
    const gsl_rng_type *rng_t;
    gsl_rng **rng;
    gsl_rng_env_setup();
    rng_t = gsl_rng_ranlxs0;

    //want to break up simulation by angle and injection frame & have each thread save data in its own folder
    //have each thread check if its directory is made and if its restarting (delete evrything) or if its continuing with a previous simulation
    //the angle and the injection frames will be the names of mc_dir, therefore read mc.par first in MC_XXX directory

    //make strings of proper directories etc.
	snprintf(flash_prefix,sizeof(flash_prefix),"%s%s",FILEPATH,FILEROOT );
	//snprintf(mc_dir,sizeof(flash_prefix),"%s%s",FILEPATH,MC_PATH);
    //snprintf(mc_file,sizeof(flash_prefix),"%s%s",mc_dir, MCPAR);
    snprintf(mc_file,sizeof(flash_prefix),"%s%s%s",FILEPATH, MC_PATH,MCPAR);

    //printf(">> mc.py:  Reading mc.par\n");

    readMcPar(mc_file, &fps, &theta_jmin, &theta_jmax, &delta_theta, &inj_radius_small,&inj_radius_large, &frm0_small , &frm0_large ,&last_frm ,&frm2_small, &frm2_large, &ph_weight_suggest, &min_photons, &max_photons, &spect, &restrt, &num_thread, &dim_switch); //thetas that comes out is in degrees

    //printf("small: %d large: %d weights: %e, %c\n", min_photons, max_photons, ph_weight_suggest, restrt);
    if (num_thread==0)
    {
        num_thread=omp_get_max_threads(); //if user specifies 0 in mc.par, default to max number of threads possible
    }

    half_threads=floor((num_thread/2));
    //thread_theta=malloc(((half_threads)+1)*sizeof(double) );

    rng = (gsl_rng **) malloc((num_thread ) * sizeof(gsl_rng *));
     rng[0] = gsl_rng_alloc (rng_t); //initalize first random number generator to seed the others with random numbers

    for(j=1;j<num_thread;j++)
    {
        rng[j] = gsl_rng_alloc (rng_t);
        gsl_rng_set(rng[j],gsl_rng_get(rng[0]));
    }

    //divide up angles and frame injections among threads DONT WANT NUMBER OF THREADS TO BE ODD
    //assign ranges to array that hold them
    /*
    delta_theta=(theta_jmax-theta_jmin)/(half_threads);
    *(thread_theta+0)=theta_jmin;
    //printf("%e\n", *(thread_theta+0));

    for (j=1;j<(half_threads+1); j++)
    {
        *(thread_theta+j)=*(thread_theta+(j-1))+delta_theta;
        //printf("%e\n", *(thread_theta+j)*180/M_PI);
    }
    */

     //leave angles in degrees here
    num_angles=(int) (((theta_jmax-theta_jmin)/delta_theta)) ;//*(180/M_PI));
    thread_theta=malloc( num_angles *sizeof(double) );
    *(thread_theta+0)=theta_jmin;//*(180/M_PI);
    printf("%e\n", *(thread_theta+0));

    for (j=1;j<(num_angles); j++)
    {
        *(thread_theta+j)=*(thread_theta+(j-1))+delta_theta;
        printf("%e\n", *(thread_theta+j));
    }

    //printf("half_threads: %d\n", half_threads);
    //start parallel section with half the threads possible and assign each thread to its value of theta by rewriting its private value of theta min and max that was read from mcrat
    omp_set_nested(1); //allow for nested parallelization
    //omp_set_dynamic(0);
    //omp_set_num_threads(num_thread);
    //#pragma omp parallel firstprivate(theta_jmin_thread, theta_jmax_thread, mc_dir, phPtr, framestart, scatt_framestart, num_ph, restrt, time_now, st, dirp, file_count, mc_filename, mc_operation, dt_max)
    //printf("%d\n", omp_get_num_threads() );
    //    for (angle_count=(int) (theta_jmin); angle_count< (int) (theta_jmax+1)  ;angle_count=angle_count+delta_theta )

    #pragma omp parallel for num_threads(num_thread) private(angle_count)
    for (angle_count=0; angle_count< num_angles ;angle_count++ )
    {

        printf("%d\t%lf\n", omp_get_thread_num(), delta_theta );
        double inj_radius;
        int frm2, frm0;
        char mc_filename[200]="";
        char mc_filename_2[200]="";
        char mc_operation[200]="";
        char mc_dir[200]="" ;
        int file_count = 0;
        DIR * dirp;
        struct dirent * entry;
        struct stat st = {0};
        double theta_jmin_thread=0, theta_jmax_thread=0;

        theta_jmin_thread= (*(thread_theta+angle_count))*(M_PI/180);//(*(thread_theta+omp_get_thread_num() ));
        theta_jmax_thread= ((*(thread_theta+angle_count))+delta_theta)*(M_PI/180);//(*(thread_theta+omp_get_thread_num()+1 ));
        printf("Thread %d: %0.1lf, %0.1lf \n %d %d\n", omp_get_thread_num(),  theta_jmin_thread*180/M_PI,  theta_jmax_thread*180/M_PI, frm2_small, frm2_large );

        snprintf(mc_dir,sizeof(flash_prefix),"%s%s%0.1lf-%0.1lf/",FILEPATH,MC_PATH, theta_jmin_thread*180/M_PI, theta_jmax_thread*180/M_PI ); //have to add angle into this

        printf(">> Thread %d in  MCRaT: I am working on path: %s \n",omp_get_thread_num(), mc_dir  );

        if ((theta_jmin_thread >= 0) &&  (theta_jmax_thread <= (2*M_PI/180) )) //if within small angle (0-2 degrees) use _small inj_radius and frm2
        {
            inj_radius=inj_radius_small;
            frm2=frm2_small;
            frm0=frm0_small;
        }
        else
        {
            inj_radius=inj_radius_large;
            frm2=frm2_large;
            frm0=frm0_large;
        }
        printf("Thread %d: %0.1lf, %0.1lf \n %d %e %d\n", omp_get_thread_num(),  theta_jmin_thread*180/M_PI,  theta_jmax_thread*180/M_PI, frm2, inj_radius, frm0 );

        //want to also have another set of threads that each has differing ranges of frame injections therefore do nested parallelism here so each thread can read its own checkpoint file
        //#pragma omp parallel num_threads(2) firstprivate(restrt)
        {
            char flash_file[200]="";
            char log_file[200]="";
            FILE *fPtr=NULL; //pointer to log file for each thread
            double *xPtr=NULL,  *yPtr=NULL,  *rPtr=NULL,  *thetaPtr=NULL,  *velxPtr=NULL,  *velyPtr=NULL,  *densPtr=NULL,  *presPtr=NULL,  *gammaPtr=NULL,  *dens_labPtr=NULL;
            double *phiPtr=NULL, *velzPtr=NULL, *zPtr=NULL;
            double *szxPtr=NULL,*szyPtr=NULL, *tempPtr=NULL; //pointers to hold data from FLASH files
            int num_ph=0, array_num=0, ph_scatt_index=0, max_scatt=0, min_scatt=0,i=0; //number of photons produced in injection algorithm, number of array elleemnts from reading FLASH file, index of photon whch does scattering, generic counter
            double dt_max=0, thescatt=0, accum_time=0;
            double  gamma_infinity=0, time_now=0, time_step=0, avg_scatt=0; //gamma_infinity not used?
            double ph_dens_labPtr=0, ph_vxPtr=0, ph_vyPtr=0, ph_tempPtr=0, ph_vzPtr=0;// *ph_cosanglePtr=NULL ;
            double min_r=0, max_r=0;
            int frame=0, scatt_frame=0, frame_scatt_cnt=0, scatt_framestart=0, framestart=0;
            struct photon *phPtr=NULL; //pointer to array of photons

            //if (omp_get_thread_num()==0)
            //{
                //printf( "A:%d Im Thread: %d with ancestor %d working in %s\n", omp_get_num_threads(),  omp_get_thread_num(), omp_get_ancestor_thread_num(1), mc_dir);
            //}

            if (restrt=='c')
            {
                printf(">> mc.py:  Reading checkpoint\n");
                //#pragma omp critical
                {
                    readCheckpoint(mc_dir, &phPtr, frm0, &framestart, &scatt_framestart, &num_ph, &restrt, &time_now);

                /*
                for (i=0;i<num_ph;i++)
                {
                    printf("%e,%e,%e, %e,%e,%e, %e, %e\n",(phPtr+i)->p0, (phPtr+i)->p1, (phPtr+i)->p2, (phPtr+i)->p3, (phPtr+i)->r0, (phPtr+i)->r1, (phPtr+i)->r2, (phPtr+i)->num_scatt );
                }
                */
                if (restrt=='c')
                {
                    printf(">> Thread %d with ancestor %d: Starting from photons injected at frame: %d out of %d\n", omp_get_thread_num(), omp_get_ancestor_thread_num(1),framestart, frm2);
                    printf(">> Thread %d with ancestor %d: Continuing scattering %d photons from frame: %d\n", omp_get_thread_num(), omp_get_ancestor_thread_num(1),num_ph, scatt_framestart);
                    printf(">> Thread %d with ancestor %d: The time now is: %e\n", omp_get_thread_num(), omp_get_ancestor_thread_num(1),time_now);
                }
                else
                {
                    printf(">> Thread %d with ancestor %d: Continuing simulation by injecting photons at frame: %d out of %d\n", omp_get_thread_num(), omp_get_ancestor_thread_num(1),framestart, frm2); //starting with new photon injection is same as restarting sim
                }
                }
            }
            else if (stat(mc_dir, &st) == -1)
            {
                mkdir(mc_dir, 0777); //make the directory with full permissions

                framestart=frm0; //if restarting then start from parameters given in mc.par file
                scatt_framestart=frm0;
            }
            else
            {
                //remove everything from MC directory to ensure no corruption of data if theres other files there besides the mc.par file
                //for a checkpoint implementation, need to find the latest checkpoint file and read it and not delete the files
                #pragma omp critical
                {
                    printf(">> Thread %d with ancestor %d:  Cleaning directory \n",omp_get_thread_num(), omp_get_ancestor_thread_num(1));
                    dirp = opendir(mc_dir);
                    while ((entry = readdir(dirp)) != NULL)
                    {
                        if (entry->d_type == DT_REG) { /* If the entry is a regular file */
                             file_count++; //count how many files are in dorectory
                        }
                    }
                    printf("File count %d\n", file_count);

                    if (file_count>0)
                    {
                        for (i=0;i<=last_frm;i++)
                        {

                            snprintf(mc_filename,sizeof(mc_filename),"%s%s%d%s", mc_dir,"mcdata_",i,"_P0.dat");
                            if(( access( mc_filename, F_OK ) != -1 ) )
                            {
                                snprintf(mc_operation,sizeof(flash_prefix),"%s%s%s%d%s","exec rm ", mc_dir,"mcdata_",i,"_*.dat"); //prepares string to remove *.dat in mc_dir
                                printf("%s\n",mc_operation);
                                system(mc_operation);

                                //snprintf(mc_operation,sizeof(flash_prefix),"%s%s%s%d%s","exec rm ", mc_dir,"mcdata_",i,"_*");
                                //system(mc_operation);
                            }
                        }
                        //snprintf(mc_operation,sizeof(flash_prefix),"%s%s%s","exec rm ", mc_dir,"mcdata_PW_*.dat"); //prepares string to remove *.dat in mc_dir
                        //system(mc_operation);

                        snprintf(mc_operation,sizeof(flash_prefix),"%s%s%s","exec rm ", mc_dir,"mcdata_PW.dat"); //prepares string to remove *.dat in mc_dir
                        system(mc_operation);

                        //snprintf(mc_operation,sizeof(flash_prefix),"%s%s%s","exec rm ", mc_dir,"mc_chkpt_*.dat"); //prepares string to remove *.dat in mc_dir
                        //system(mc_operation);

                        snprintf(mc_operation,sizeof(flash_prefix),"%s%s%s","exec rm ", mc_dir,"mc_output_*.log"); //prepares string to remove *.log in mc_dir
                        system(mc_operation);

                    }
                }
                framestart=frm0; //if restarting then start from parameters given in mc.par file
                scatt_framestart=frm0;

            }

            dt_max=1.0/fps;
            //#pragma omp barrier

            snprintf(log_file,sizeof(log_file),"%s%s",mc_dir,"mc_output.log" );
            printf("%s\n",log_file);
            fPtr=fopen(log_file, "w");

            fprintf(fPtr, "%d Im Thread: %d with ancestor %d Starting on Frame: %d scatt_framestart: %d\n", omp_get_num_threads(),  omp_get_thread_num(), omp_get_ancestor_thread_num(1), framestart, scatt_framestart);
            fflush(fPtr);
            //fclose(fPtr);
            //#pragma omp barrier
            //exit(0);
            //loop over frames
            //for a checkpoint implementation, start from the last saved "frame" value and go to the saved "frm2" value

            //#pragma omp for
            for (frame=framestart;frame<=frm2;frame++)
            {
                 if (restrt=='r')
                 {
                    time_now=frame/fps; //for a checkpoint implmentation, load the saved "time_now" value when reading the ckeckpoint file otherwise calculate it normally
                 }

                //printf(">> mc.py: Working on Frame %d\n", frame);
                fprintf(fPtr,"%d Im Thread: %d with ancestor %d Working on Frame: %d\n", omp_get_num_threads(),  omp_get_thread_num(), omp_get_ancestor_thread_num(1), frame);
                fflush(fPtr);

                if (restrt=='r')
                {

                    //exit(0);
                    //read in FLASH file
                    //for a checkpoint implmentation, dont need to read the file yet
                    if (dim_switch==0)
                    {
                        //put proper number at the end of the flash file
                        modifyFlashName(flash_file, flash_prefix, frame, dim_switch);
                        fprintf(fPtr,">> Im Thread: %d with ancestor %d: Opening FLASH file %s\n",omp_get_thread_num(), omp_get_ancestor_thread_num(1), flash_file);
                        fflush(fPtr);

                        #pragma omp critical
                        {
                            readAndDecimate(flash_file, inj_radius, fps, &xPtr,  &yPtr,  &szxPtr, &szyPtr, &rPtr,\
                                &thetaPtr, &velxPtr,  &velyPtr,  &densPtr,  &presPtr,  &gammaPtr,  &dens_labPtr, &tempPtr, &array_num, 1, min_r, max_r, fPtr);
                        }
                    }
                    else
                    {
                        read_hydro(FILEPATH, frame, inj_radius, &xPtr,  &yPtr, &zPtr,  &szxPtr, &szyPtr, &rPtr,\
                                &thetaPtr, &phiPtr, &velxPtr,  &velyPtr, &velzPtr,  &densPtr,  &presPtr,  &gammaPtr,  &dens_labPtr, &tempPtr, &array_num, 1, min_r, max_r, fps, fPtr);
                    }

                    //check for run type
                    if(strcmp(cyl, this_run)==0)
                    {
                        //printf("In cylindrical prep\n");
                        cylindricalPrep(gammaPtr, velxPtr, velyPtr, densPtr, dens_labPtr, presPtr, tempPtr, array_num);
                    }
                    else if (strcmp(sph, this_run)==0)
                    {
                        sphericalPrep(rPtr, xPtr, yPtr,gammaPtr, velxPtr, velyPtr, densPtr, dens_labPtr, presPtr, tempPtr, array_num );
                    }

                    //determine where to place photons and how many should go in a given place
                    //for a checkpoint implmentation, dont need to inject photons, need to load photons' last saved data
                    fprintf(fPtr,">>  Thread: %d with ancestor %d: Injecting photons\n",omp_get_thread_num(), omp_get_ancestor_thread_num(1));
                    fflush(fPtr);

                    if (dim_switch==0)
                    {
                        photonInjection(&phPtr, &num_ph, inj_radius, ph_weight_suggest, min_photons, max_photons,spect, array_num, fps, theta_jmin_thread, theta_jmax_thread, xPtr, yPtr, szxPtr, szyPtr,rPtr,thetaPtr, tempPtr, velxPtr, velyPtr,rng[omp_get_thread_num()] );
                    }
                    else
                    {
                        photonInjection3D(&phPtr, &num_ph, inj_radius, ph_weight_suggest, min_photons, max_photons,spect, array_num, fps, theta_jmin_thread, theta_jmax_thread, xPtr, yPtr, zPtr, szxPtr, szyPtr,rPtr,thetaPtr, phiPtr, tempPtr, velxPtr, velyPtr, velzPtr, rng[omp_get_thread_num()] );
                    }

                    //printf("%d\n",num_ph); //num_ph is one more photon than i actually have
                    /*
                    for (i=0;i<num_ph;i++)
                        printf("%e,%e,%e \n",(phPtr+i)->r0, (phPtr+i)->r1, (phPtr+i)->r2 );
                    */
                }

                //scatter photons all the way thoughout the jet
                //for a checkpoint implmentation, start from the last saved "scatt_frame" value eh start_frame=frame or start_frame=cont_frame
                if (restrt=='r')
                {
                    scatt_framestart=frame; //have to make sure that once the inner loop is done and the outer loop is incrememnted by one the inner loop starts at that new value and not the one read by readCheckpoint()
                }

                for (scatt_frame=scatt_framestart;scatt_frame<=last_frm;scatt_frame++)
                {
                    fprintf(fPtr,">>\n");
                    fprintf(fPtr,">> Thread %d with ancestor %d : Working on photons injected at frame: %d out of %d\n", omp_get_thread_num(), omp_get_ancestor_thread_num(1),frame, frm2);
                    fprintf(fPtr,">> Thread %d with ancestor %d: %s - Working on frame %d\n",omp_get_thread_num(), omp_get_ancestor_thread_num(1), THISRUN, scatt_frame);
                    fprintf(fPtr,">> Thread %d with ancestor %d: Opening file...\n", omp_get_thread_num(), omp_get_ancestor_thread_num(1));
                    fflush(fPtr);

                    if (dim_switch==0)
                    {
                        //put proper number at the end of the flash file
                        modifyFlashName(flash_file, flash_prefix, scatt_frame, dim_switch);
                        #pragma omp critical
                        {
                            phMinMax(phPtr, num_ph, &min_r, &max_r);
                            readAndDecimate(flash_file, inj_radius, fps, &xPtr,  &yPtr,  &szxPtr, &szyPtr, &rPtr,\
                                &thetaPtr, &velxPtr,  &velyPtr,  &densPtr,  &presPtr,  &gammaPtr,  &dens_labPtr, &tempPtr, &array_num, 0, min_r, max_r, fPtr);
                        }
                    }
                    else
                    {
                        phMinMax(phPtr, num_ph, &min_r, &max_r);

                        read_hydro(FILEPATH, scatt_frame, inj_radius, &xPtr,  &yPtr, &zPtr,  &szxPtr, &szyPtr, &rPtr,\
                                &thetaPtr, &phiPtr, &velxPtr,  &velyPtr, &velzPtr,  &densPtr,  &presPtr,  &gammaPtr,  &dens_labPtr, &tempPtr, &array_num, 0, min_r, max_r, fps, fPtr);
                    }

                    //check for run type
                    if(strcmp(cyl, this_run)==0)
                    {
                        //printf("In cylindrical prep\n");
                        cylindricalPrep(gammaPtr, velxPtr, velyPtr, densPtr, dens_labPtr, presPtr, tempPtr, array_num);
                    }
                    else if (strcmp(sph, this_run)==0)
                    {
                        sphericalPrep(rPtr, xPtr, yPtr,gammaPtr, velxPtr, velyPtr, densPtr, dens_labPtr, presPtr, tempPtr, array_num );
                    }
                        //printf("The result of read and decimate are arrays with %d elements\n", array_num);

                    fprintf(fPtr,">> Thread %d with ancestor %d: propagating and scattering %d photons\n", omp_get_thread_num(), omp_get_ancestor_thread_num(1),num_ph);
                    fflush(fPtr);

                    frame_scatt_cnt=0;
                    while (time_now<((scatt_frame+1)/fps))
                    {
                        //if simulation time is less than the simulation time of the next frame, keep scattering in this frame
                        //go through each photon and find blocks closest to each photon and properties of those blocks to calulate mean free path
                        //and choose the photon with the smallest mfp and calculate the timestep

                        ph_scatt_index=findNearestPropertiesAndMinMFP(phPtr, num_ph, array_num, &time_step, xPtr,  yPtr, zPtr, velxPtr,  velyPtr,  velzPtr, dens_labPtr, tempPtr,\
                            &ph_dens_labPtr, &ph_vxPtr, &ph_vyPtr, &ph_vzPtr, &ph_tempPtr, rng[omp_get_thread_num()], dim_switch);

                        //printf("In main: %e, %d, %e, %e\n", *(ph_num_scatt+ph_scatt_index), ph_scatt_index, time_step, time_now);
                        printf("In main: %e, %d, %e, %e\n",((phPtr+ph_scatt_index)->num_scatt), ph_scatt_index, time_step, time_now);

                         if (time_step<dt_max)
                        {

                            //update number of scatterings and time
                            //(*(ph_num_scatt+ph_scatt_index))+=1;
                            ((phPtr+ph_scatt_index)->num_scatt)+=1;
                            frame_scatt_cnt+=1;
                            time_now+=time_step;

                            updatePhotonPosition(phPtr, num_ph, time_step);

                            //scatter the photon
                            //printf("Passed Parameters: %e, %e, %e\n", (ph_vxPtr), (ph_vyPtr), (ph_tempPtr));

                            photonScatter( (phPtr+ph_scatt_index), (ph_vxPtr), (ph_vyPtr), ph_vzPtr, (ph_tempPtr), rng[omp_get_thread_num()] , dim_switch, fPtr);

                            //if (frame_scatt_cnt%1000 == 0)
                            {
                                fprintf(fPtr,"Scattering Number: %d\n", frame_scatt_cnt);
                                fprintf(fPtr,"The local temp is: %e\n", (ph_tempPtr));
                                fprintf(fPtr,"Average photon energy is: %e\n", averagePhotonEnergy(phPtr, num_ph)); //write function to average over the photons p0 and then do (*3e10/1.6e-9)
                                fflush(fPtr);
                            }

                        }
                        else
                        {
                            time_now+=dt_max;

                            //for each photon update its position based on its momentum

                            updatePhotonPosition(phPtr, num_ph, dt_max);
                        }

                        //printf("In main 2: %e, %d, %e, %e\n", ((phPtr+ph_scatt_index)->num_scatt), ph_scatt_index, time_step, time_now);

                    }

                //get scattering statistics
                phScattStats(phPtr, num_ph, &max_scatt, &min_scatt, &avg_scatt);

                fprintf(fPtr,"The number of scatterings in this frame is: %d\n", frame_scatt_cnt);
                fprintf(fPtr,"The last time step was: %lf.\nThe time now is: %lf\n", time_step,time_now);
                fprintf(fPtr,"The maximum number of scatterings for a photon is: %d\nThe minimum number of scattering for a photon is: %d\n", max_scatt, min_scatt);
                fprintf(fPtr,"The average number of scatterings thus far is: %lf\n", avg_scatt);
                fflush(fPtr);

                printPhotons(phPtr, num_ph,  scatt_frame , frame, mc_dir);
                exit(0);
                //for a checkpoint implmentation,save the checkpoint file here after every 5 frames or something
                //save the photons data, the scattering number data, the scatt_frame value, and the frame value
                //WHAT IF THE PROGRAM STOPS AFTER THE LAST SCATT_FRAME, DURING THE FIRST SCATT_FRAME OF NEW FRAME VARIABLE - save restrt variable as 'r'
                fprintf(fPtr, ">> Thread %d with ancestor %d: Making checkpoint file\n", omp_get_thread_num(), omp_get_ancestor_thread_num(1));
                fflush(fPtr);
                saveCheckpoint(mc_dir, frame, scatt_frame, num_ph, time_now, phPtr, last_frm);

                free(xPtr);free(yPtr);free(szxPtr);free(szyPtr);free(rPtr);free(thetaPtr);free(velxPtr);free(velyPtr);free(densPtr);free(presPtr);
                free(gammaPtr);free(dens_labPtr);free(tempPtr);
                xPtr=NULL; yPtr=NULL;  rPtr=NULL;thetaPtr=NULL;velxPtr=NULL;velyPtr=NULL;densPtr=NULL;presPtr=NULL;gammaPtr=NULL;dens_labPtr=NULL;
                szxPtr=NULL; szyPtr=NULL; tempPtr=NULL;
                }
                restrt='r';//set this to make sure that the next iteration of propogating photons doesnt use the values from the last reading of the checkpoint file
                free(phPtr);
                phPtr=NULL;

            }

        }//end omp parallel inner section

        //merge files from each worker thread within a directory
        //dirFileMerge(mc_dir, frm0, last_frm);

    } //end omp parallel section

    //gsl_rng_free (rand);
    //free rand number generator
    for (j=0;j<num_thread;j++)
    {
        gsl_rng_free(rng[j]);
    }
    free(rng);
    free(thread_theta);

	return 0;

}

max_stars_repo_name: outflows/MCRaT

max_issues_repo_name: outflows/MCRaT

max_forks_repo_name: outflows/MCRaT


text:
#ifndef PROXIMAL_MATRIX_H_GUARD
#define PROXIMAL_MATRIX_H_GUARD

#include <math.h>
#include <float.h>

#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_eigen.h>

void project_sdc(gsl_matrix *X);
void project_spectral_norm_ball(gsl_matrix *X);
void prox_neg_log_det(gsl_matrix *X, const double rho);
void prox_quad(gsl_vector *x, const double rho, gsl_matrix *A, gsl_matrix *b);

#endif

ext: h

max_stars_repo_name: askuyue/proximal

max_issues_repo_name: zhiyuxzh/proximal

max_forks_repo_name: zhiyuxzh/proximal


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

#pragma once

#include <string>

#include <gsl/gsl>

#include "onnx/onnx_pb.h"

#include "core/graph/basic_types.h"

namespace onnxruntime::utils {

// keep these signatures in sync with DECLARE_MAKE_ATTRIBUTE_FNS below
/** Creates an AttributeProto with the specified name and value. */
ONNX_NAMESPACE::AttributeProto MakeAttribute(std::string attr_name, int64_t value);
/** Creates an AttributeProto with the specified name and values. */
ONNX_NAMESPACE::AttributeProto MakeAttribute(std::string attr_name, gsl::span<const int64_t> values);

#define DECLARE_MAKE_ATTRIBUTE_FNS(type)                                           \
  ONNX_NAMESPACE::AttributeProto MakeAttribute(std::string attr_name, type value); \
  ONNX_NAMESPACE::AttributeProto MakeAttribute(std::string attr_name, gsl::span<const type> values)

DECLARE_MAKE_ATTRIBUTE_FNS(float);
DECLARE_MAKE_ATTRIBUTE_FNS(std::string);
DECLARE_MAKE_ATTRIBUTE_FNS(ONNX_NAMESPACE::TensorProto);
#if !defined(DISABLE_SPARSE_TENSORS)
DECLARE_MAKE_ATTRIBUTE_FNS(ONNX_NAMESPACE::SparseTensorProto);
#endif
DECLARE_MAKE_ATTRIBUTE_FNS(ONNX_NAMESPACE::TypeProto);
DECLARE_MAKE_ATTRIBUTE_FNS(ONNX_NAMESPACE::GraphProto);

#undef DECLARE_MAKE_ATTRIBUTE_FNS

// The below overload is made so the compiler does not attempt to resolve
// string literals with the gsl::span overload
inline ONNX_NAMESPACE::AttributeProto MakeAttribute(std::string attr_name, const char* value) {
  return MakeAttribute(std::move(attr_name), std::string{value});
}

/**
 * Sets an attribute in `node_attributes` with key `attribute.name()` and value `attribute`.
 * If an attribute with the same name exists, it will be overwritten.
 * @return Pair of (iterator to attribute, whether attribute was added (true) or updated (false)).
 */
std::pair<NodeAttributes::iterator, bool> SetNodeAttribute(ONNX_NAMESPACE::AttributeProto attribute,
                                                           NodeAttributes& node_attributes);

}  // namespace onnxruntime::utils

max_stars_repo_name: SiriusKY/onnxruntime

max_issues_repo_name: SiriusKY/onnxruntime

max_forks_repo_name: SiriusKY/onnxruntime


text:
#include <stdio.h>
#include <gsl/gsl_statistics.h>

int main(void)
{
  double data[5] = {17.2, 18.1, 16.5, 18.3, 12.6};
  double mean, variance, largest, smallest;

  mean     = gsl_stats_mean(data, 1, 5);
  variance = gsl_stats_variance(data, 1, 5);
  largest  = gsl_stats_max(data, 1, 5);
  smallest = gsl_stats_min(data, 1, 5);

  printf ("The dataset is %g, %g, %g, %g, %g\n",
         data[0], data[1], data[2], data[3], data[4]);

  printf ("The sample mean is %g\n", mean);
  printf ("The estimated variance is %g\n", variance);
  printf ("The largest value is %g\n", largest);
  printf ("The smallest value is %g\n", smallest);
  return 0;
}

ext: c


text:
#include <pygsl/error_helpers.h>
#include <pygsl/function_helpers.h>
#include <pygsl/pygsl_features.h>

#ifdef PyGSL_DERIV_MODULE
#if ( _PYGSL_GSL_HAS_DERIV == 0 )
#error "The deriv module was only introduced by GSL 1.5. You seem to compile against an older verion!"
#endif
#include <gsl/gsl_deriv.h>
#endif /* PyGSL_DERIV_MODULE */

#include <gsl/gsl_diff.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <setjmp.h>
/*
 * callback functions
 * - python function passed by user
 * - actual C callback
 * - GSL wrapper struct
 */

/* Used for traceback */
static PyObject *module = NULL;

typedef struct{
	PyObject * callback;
	PyObject * args;
	jmp_buf  buffer;
}pygsl_diff_args;

static double
diff_callback(double x, void *p)
{
	double value;
	int flag;
	pygsl_diff_args *pargs = NULL;

	pargs = (pygsl_diff_args *) p;

	assert(pargs->callback);
	assert(pargs->args);
	flag = PyGSL_function_wrap_helper(x, &value, NULL, pargs->callback,
					  pargs->args, (char *)__FUNCTION__);

	if(GSL_SUCCESS != flag){
		longjmp(pargs->buffer, flag);
		return gsl_nan();
	}
	return value;
}

/* wrapper function */
typedef int pygsl_deriv_func(const gsl_function *, double, double, double *, double *);
typedef int pygsl_diff_func(const gsl_function *, double, double *, double *);

static PyObject *
PyGSL_diff_generic(PyObject *self, PyObject *args,
#ifndef PyGSL_DIFF_MODULE
		   pygsl_deriv_func func
#else
		   pygsl_diff_func func
#endif
)
{
	PyObject *result=NULL, *myargs=NULL;
	PyObject *cb=NULL;

	pygsl_diff_args pargs = {NULL, NULL};
	/* Changed to compile using Sun's Compiler */
	gsl_function diff_gsl_callback = {NULL, NULL};

	double x, value, abserr;
	int flag;
#ifndef PyGSL_DIFF_MODULE
	double h;
	if(! PyArg_ParseTuple(args, "Odd|O", &cb, &x, &h, &myargs)){
		return NULL;
	}
#else
	if(! PyArg_ParseTuple(args, "Od|O", &cb, &x, &myargs)){
		return NULL;
	}
#endif

	/* Changed to compile using Sun's Compiler */
	diff_gsl_callback.function = diff_callback;
	diff_gsl_callback.params = (void *) &pargs;

	if(! PyCallable_Check(cb)) {
		PyErr_SetString(PyExc_TypeError,
				"The first parameter must be callable");
		return NULL;
	}
	Py_INCREF(cb);               /* Add a reference to new callback */

	pargs.callback = cb;        /* Remember new callback */

	/* Did I get arguments? If so handle them */
	if(NULL == myargs){
		Py_INCREF(Py_None);
		pargs.args= Py_None;
	}else{
		Py_INCREF(myargs);
		pargs.args= myargs;
	}

	if((flag=setjmp(pargs.buffer)) == 0){
		/* Jmp buffer set, call the function */
#ifndef PyGSL_DIFF_MODULE
	     flag = func(&diff_gsl_callback, x, h, &value, &abserr);
#else
	     flag = func(&diff_gsl_callback, x, &value, &abserr);
#endif
	}else{
		DEBUG_MESS(2, "CALLBACK called longjmp! flag =%d", flag);
	}

	/* Arguments no longer used */
	Py_DECREF(pargs.args);

	/* Dispose of callback */
	Py_DECREF(pargs.callback);

	if(flag != GSL_SUCCESS){
		PyGSL_ERROR_FLAG(flag);
		return NULL;
	}

	result = Py_BuildValue("(dd)", value, abserr);
	return result;
}

max_stars_repo_name: juhnowski/FishingRod

max_issues_repo_name: juhnowski/FishingRod

max_forks_repo_name: juhnowski/FishingRod


text:
#include <asm/types.h>
#include <glib.h>
#include <stdio.h>
#include <assert.h>
#include <inttypes.h>
#include <float.h>

#include <gsl/gsl_histogram.h>

#include <blktrace_api.h>
#include <blktrace.h>
#include <plugins.h>
#include <utils.h>
#include <list_plugins.h>
#include <reqsize.h>

#define DECL_ASSIGN_I2C(name,data)				\
	struct i2c_data *name = (struct i2c_data *)data

#define OIO_ALLOC	(8)

#define N_BINS		(21)
#define BINS_SEP	(8)

struct oio_data
{
#define READ	0
#define WRITE	1
	gsl_histogram *op[2];
	__u64 time;
};

struct i2c_data
{
	GTree *is;
	FILE *oio_f;

	__u32 outstanding;
	__u32 maxouts;

	/* oio hist */
	struct oio_data *oio;
	__u32 oio_size;
	__u64 oio_prev_time;
	FILE *oio_hist_f;
};

static void write_outs(struct i2c_data *i2c, struct blk_io_trace *t)
{
	if(i2c->oio_f)
		fprintf(i2c->oio_f, "%f %u\n",
				NANO_ULL_TO_DOUBLE(t->time),
				i2c->outstanding);
}

static gsl_histogram *new_hist()
{
	gsl_histogram *h;

	/* allocate bins and set ranges uniformally */
	h = gsl_histogram_alloc(N_BINS);
	gsl_histogram_set_ranges_uniform(h,0,N_BINS*BINS_SEP);

	/* last bin should hold everything greater
	 * than (N_BINS-1)*BINS_SEP */
	h->range[h->n] = DBL_MAX;

	return h;
}

static void init_oio_data(struct oio_data *oio, int n)
{
	int i;

	for (i = 0; i < n; i++) {
		oio[i].time = 0;
		oio[i].op[READ] = new_hist();
		oio[i].op[WRITE] = new_hist();
	}
}

static gboolean add_to_matrix(__u64 *__unused, struct blk_io_trace *t, struct i2c_data *i2c)
{
	gsl_histogram_increment(i2c->oio[i2c->outstanding].op[IS_WRITE(t)], (double)(t->bytes/BLK_SIZE));

	return FALSE;
}

static void oio_change(struct i2c_data *i2c, struct blk_io_trace *t, int inc)
{
	/* allocate oio space if the one I had is over */
	if(i2c->outstanding + 1 >= i2c->oio_size) {
		i2c->oio = realloc(i2c->oio, (i2c->oio_size + OIO_ALLOC)*sizeof(struct oio_data));
		init_oio_data(i2c->oio + i2c->oio_size, OIO_ALLOC);
		i2c->oio_size += OIO_ALLOC;
	}

	/* increase the time */
	if(i2c->oio_prev_time != UINT64_MAX) {
		i2c->oio[i2c->outstanding].time += t->time - i2c->oio_prev_time;
	}
	i2c->oio_prev_time = t->time;

	if(inc) {
		i2c->outstanding++;
	} else {
		i2c->outstanding--;
	}
	i2c->maxouts = MAX(i2c->maxouts, i2c->outstanding);

	g_tree_foreach(i2c->is, (GTraverseFunc)add_to_matrix, i2c);

	write_outs(i2c, t);
}

static void C(struct blk_io_trace *t, void *data)
{
	DECL_ASSIGN_I2C(i2c,data);

	if(g_tree_lookup(i2c->is, &t->sector)!=NULL) {
		g_tree_remove(i2c->is, &t->sector);

		oio_change(i2c, t, FALSE);
	}
}

static void I(struct blk_io_trace *t, void *data)
{
	DECL_ASSIGN_I2C(i2c,data);

	if(g_tree_lookup(i2c->is, &t->sector)==NULL) {
		DECL_DUP(struct blk_io_trace,new_t,t);
		g_tree_insert(i2c->is,&new_t->sector,new_t);

		oio_change(i2c, t, TRUE);
	}
}

static void add_histogram(gsl_histogram *h1, gsl_histogram *h2)
{
	unsigned int i;

	for (i = 0; i < h1->n && i < h2->n; i++) {
		h1->bin[i] += h2->bin[i];
	}
}

void i2c_add(void *data1, const void *data2)
{
	DECL_ASSIGN_I2C(i2c1,data1);
	DECL_ASSIGN_I2C(i2c2,data2);
	__u32 i;

	i2c1->maxouts = MAX(i2c1->maxouts, i2c2->maxouts);
	i2c1->oio_prev_time = MAX(i2c1->oio_prev_time, i2c2->oio_prev_time);

	if(i2c1->oio_size < i2c2->oio_size) {
		__u32 diff = i2c2->oio_size - i2c1->oio_size;
		i2c1->oio = realloc(i2c1->oio, i2c2->oio_size*sizeof(struct oio_data));
		init_oio_data(i2c1->oio + diff, diff);
		i2c1->oio_size = i2c2->oio_size;
	}

	for (i = 0; i <= i2c2->maxouts; i++) {
		i2c1->oio[i].time += i2c2->oio[i].time;
		add_histogram(i2c1->oio[i].op[READ],i2c2->oio[i].op[READ]);
		add_histogram(i2c1->oio[i].op[WRITE],i2c2->oio[i].op[WRITE]);
	}
}

void i2c_print_results(const void *data)
{
	DECL_ASSIGN_I2C(i2c,data);
	double p;
	double avg = 0;
	__u32 i;
	__u64 tot_time = 0;

	for (i = 0; i <= i2c->maxouts; i++) {
		tot_time += i2c->oio[i].time;
	}

	for (i = 0; i <= i2c->maxouts; i++) {
		p = ((double)i2c->oio[i].time)/((double)tot_time);

		if(i2c->oio_hist_f)
			fprintf(i2c->oio_hist_f, "%u\t%.2lf\n",i,100*p);

		avg += p*i;
	}

	/* print all histograms */
	if(i2c->oio_hist_f) {
		for (i = 1; i <= i2c->maxouts; i++) {
			fprintf(i2c->oio_hist_f, "\nread: %d\n",i);
			gsl_histogram_fprintf(i2c->oio_hist_f,
					i2c->oio[i].op[READ], "%g", "%g");
			fprintf(i2c->oio_hist_f, "\nwrite: %d\n",i);
			gsl_histogram_fprintf(i2c->oio_hist_f,
					i2c->oio[i].op[WRITE], "%g", "%g");
		}
	}

	printf("I2C Max. OIO: %u, Avg: %.2lf\n",i2c->maxouts,avg);
}

void i2c_init(struct plugin *p, struct plugin_set *__un1, struct plug_args *pa)
{
	char filename[FILENAME_MAX];
	struct i2c_data *i2c = p->data = g_new(struct i2c_data,1);

	i2c->is = g_tree_new(comp_int64);
	i2c->outstanding = 0;
	i2c->maxouts = 0;

	i2c->oio_f = NULL;
	if(pa->i2c_oio_f) {
		get_filename(filename, "i2c_oio", pa->i2c_oio_f, pa->end_range);
		i2c->oio_f = fopen(filename,"w");
		if(!i2c->oio_f) perror_exit("Opening I2C detail file");
	}

	i2c->oio_hist_f = NULL;
	if(pa->i2c_oio_hist_f) {
		get_filename(filename, "i2c_oio_hist", pa->i2c_oio_hist_f, pa->end_range);
		i2c->oio_hist_f = fopen(filename,"w");
		if(!i2c->oio_hist_f) perror_exit("Opening I2C detail file");
	}

	i2c->oio = NULL;
	i2c->oio_size = 0;
	i2c->oio_prev_time = UINT64_MAX;
}

void i2c_ops_init(struct plugin_ops *po)
{
	po->add = i2c_add;
	po->print_results = i2c_print_results;

	/* association of event int and function */
	g_tree_insert(po->event_tree,(gpointer)__BLK_TA_COMPLETE,C);
	g_tree_insert(po->event_tree,(gpointer)__BLK_TA_INSERT,I);
}

void i2c_destroy(struct plugin *p)
{
	unsigned int i;

	DECL_ASSIGN_I2C(i2c,p->data);

	if(i2c->oio_f)
		fclose(i2c->oio_f);

	if(i2c->oio_hist_f)
		fclose(i2c->oio_hist_f);

	for (i = 0; i <= i2c->maxouts; i++) {
		gsl_histogram_free(i2c->oio[i].op[READ]);
		gsl_histogram_free(i2c->oio[i].op[WRITE]);
	}
	free(i2c->oio);

	g_free(p->data);
}

max_stars_repo_name: useche/btstats

max_issues_repo_name: useche/btstats

max_forks_repo_name: useche/btstats


text:
/*****************************************************************************/
#define this_is "FDBinary v.3 beta (test version of 21 Feb 2010)"
/*****************************************************************************/

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

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_sf.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_multimin.h>

#include "mxfuns.h"
#include "fd3sep.h"
#include "triorb.h"

/*****************************************************************************/

/* function and macro to kill the program with a short message */
#define FDBErrorString "\nError in FDBianry"
#define DIE(s) {fprintf(stderr,"%s: %s\n",FDBErrorString,s);fdbfailure();};
void fdbfailure(void) { exit ( EXIT_FAILURE ); }

/*****************************************************************************/

/* macros for reading values from keyboard */
#define GETDBL(x) {if(1!=scanf("%lg",x)) DIE("failed reading double");}
#define GETINT(x) {if(1!=scanf("%d", x)) DIE("failed reading int");}
#define GETLNG(x) {if(1!=scanf("%ld",x)) DIE("failed reading long");}
#define GETSTR(x) {if(1!=scanf("%s", x)) DIE("failed reading string");}

/*****************************************************************************/

#define SPEEDOFLIGHT 2.998E5 /* speed of light in km/s */

static char *triorb_strings[] = {
  /*  0 */  "period of AB--C",                    "day", "%.9lg",
  /*  1 */  "time of periast pass in AB--C",      "day", "%.9lg",
  /*  2 */  "eccentricity of AB--C",                "1", "%.9lg",
  /*  3 */  "periast long of AB in AB--C",        "deg", "%.9lg",
  /*  4 */  "a sin(incl) of AB in AB--C",   "light-day", "%.9lg",
  /*  5 */  "a sin(incl) of C in AB--C",    "light-day", "%.9lg",
  /*  6 */  "period A--B",                        "day", "%.9lg",
  /*  7 */  "time of periast pass in A--B",       "day", "%.9lg",
  /*  8 */  "eccentricity of A--B",                 "1", "%.9lg",
  /*  9 */  "periast long of A in A--B",          "deg", "%.9lg",
  /* 10 */  "rv semiamp of A in A--B",           "km/s", "%.9lg",
  /* 11 */  "rv semiamp of B in A--B",           "km/s", "%.9lg",
  /* 12 */  "periast long adv per cycle in A--B", "deg", "%.9lg"
};

static char *triorb_pname   (long j) { return triorb_strings[3*j];   }
static char *triorb_punit   (long j) { return triorb_strings[3*j+1]; }
static char *triorb_pformat (long j) { return triorb_strings[3*j+2]; }

/*****************************************************************************/

static long   K, M, N, Ndft, ksw[3], nfp, opsw[TRIORB_NP];
static double **dftobs, **dftmod, **dftres;
static double rvstep, *otimes, *rvcorr, *sig, **lfm, **rvm;
static double op0[TRIORB_NP], dop0[TRIORB_NP];
static double meritfngsl ( const gsl_vector *v, void *params );
static double meritfn ( double *op );

#define MX_FDBINARY_FORMAT "%15.8E   "
static char *mxfd3fmts=MX_FDBINARY_FORMAT;

/* FDBinary */

int main ( void ) {

  long i, i0, i1, j, k, vc, vlen, nruns, niter;
  double **masterobs, **mod, **res, **obs, z0, z1, stoprat;
  char obsfn[1024], resfn[1024], modfn[1024], rvsfn[1024], logfn[1024];
  char *starcode[] = {"A","B","C"};
  FILE *logfp;

  setbuf ( stdout, NULL );
  MxError( FDBErrorString, stdout, fdbfailure );
  MxFormat( mxfd3fmts );

  printf ( "\n  %s\n", this_is );

  printf ( "\n  SECTION 1: LOADING OBSERVED SPECTRA\n\n" );
  printf ( "  observed spectra master file = " );
  GETSTR ( obsfn ); printf ( "%s\n", obsfn );
  vc=0; vlen=0; masterobs = MxLoad ( obsfn, &vc, &vlen ); M = vc - 1;
  printf ( "  loaded %ld spectra with %ld data pts per spectrum\n", M, vlen );
  z0 = **masterobs; z1 = *(*masterobs+vlen-1);
  rvstep = SPEEDOFLIGHT * ( - 1 + exp ((z1-z0)/(vlen-1)) );
  printf ( "  lnlambda range from %lf to %lf\n", z0, z1 );
  printf ( "  rv step %lg km/s per data point\n", rvstep );
  printf ( "  use data from lnlambda = " );
  GETDBL ( &z0 ); printf ( "%lg\n", z0 );
  printf ( "             to lnlambda = " );
  GETDBL ( &z1 ); printf ( "%lg\n", z1 );
  i0 = 0; while ( *(*masterobs+i0) < z0 ) i0++;
  i1 = vlen-1; while ( z1 < *(*masterobs+i1) ) i1--;
  N = i1 - i0 + 1;
  printf ( "  selected %ld data points per spectrum\n", N );
  obs = MxAlloc ( M+1, N );
  for ( i = 0 ; i < N ; i++ ) for ( j = 0 ; j <= M ; j++ )
    *(*(obs+j)+i) = *(*(masterobs+j)+i0+i);
  MxFree ( masterobs, vc, vlen );
  printf ( "  write selected data to file = " );
  GETSTR ( obsfn ); printf ( "%s\n", obsfn ); MxWrite ( obs, M+1, N, obsfn );

  printf ( "\n  SECTION 2: COMPONENT SPECTRA MODELING SWITCHES (0/1)\n\n" );
  for ( K = i = 0 ; i < 3 ; i++ ) {
    int sw;
    GETINT(&sw);
    sw = sw ? 1 : 0;
    printf ( "  %s = %d\n", starcode[i], sw );
    if ( sw ) ksw[K++] = i;
  }
  printf ( "\n  number of components to be resolved is %ld\n", K );

  /* allocating memory */
  Ndft = 2*(N/2 + 1);
  dftobs = MxAlloc ( M, Ndft ); dft_fwd ( M, N, obs+1, dftobs );
  otimes = *MxAlloc ( 1, M );
  rvcorr = *MxAlloc ( 1, M );
  sig = *MxAlloc ( 1, M );
  res = MxAlloc ( M+1, N ); dftres = MxAlloc ( M, Ndft );
  mod = MxAlloc ( K+1, N ); dftmod = MxAlloc ( K, Ndft );
  rvm = MxAlloc ( K, M );
  lfm = MxAlloc ( K, M );
  for ( i = 0 ; i < N ; i++ ) { *(*res+i) = *(*obs+i); *(*mod+i) = *(*obs+i); }

  printf ( "\n  SECTION 3: DESCRIPTORS TO OBSERVED SPECTRA\n" );
  printf ( "             AND LIGHT-FACTORS ASSIGNED TO COMPONENTS\n\n" );
  printf ( "  %14s%12s%12s", "t_obs", "rv_corr", "noise_rms" );
  for ( k = 0 ; k < K ; k++ ) printf ( "      lf_%s", starcode[ksw[k]] );
  printf ( "\n" );
  printf ( "  %14s%12s%12s", "[day]", "[km/s]", "[1]" );
  for ( k = 0 ; k < K ; k++ ) printf ( "       [1]" );
  printf ( "\n\n" );
  for ( j = 0 ; j < M ; j++ )
  {
    GETDBL(otimes+j); printf ( "  %14.5lf", *(otimes+j) );
    GETDBL(rvcorr+j); printf ( "%12.4lf", *(rvcorr+j) );
    GETDBL(sig+j); printf ( "%12.4lf", *(sig+j) );
    for ( k = 0; k < K ; k++ )
      { GETDBL(*(lfm+k)+j); printf ( "%10.4lf", *(*(lfm+k)+j) ); }
    putchar ( '\n' );
  }

  printf ( "\n  SECTION 4: PARAMETERS OF ORBITS\n" );
  for ( nfp = i = 0 ; i < TRIORB_NP ; i++ ) {
    if ( 0 == i ) printf ( "\n  wide (AB--C) orbit\n\n" );
    if ( 6 == i ) printf ( "\n  tight (A--B) orbit\n\n" );
    printf ( "  %s [%s] = ", triorb_pname(i), triorb_punit(i) );
    GETDBL(op0+i); printf ( triorb_pformat(i), *(op0+i) );
    printf ( " +/- " );
    GETDBL(dop0+i); printf ( triorb_pformat(i), *(dop0+i) );
    if ( *(dop0+i) ) { opsw[nfp++] = i; printf ( " *** free ***" ); }
    if ( 2 == i || 8 == i ) *(op0+i) = *(op0+i) / (1-*(op0+i)); /* ecc */
    printf ( "\n" );
  }
  printf ( "\n  number of free orbital parameters is %ld\n", nfp );

  printf ( "\n  SECTION 5: OPTIMISATION DETAILS\n\n" );
  if ( 0 == nfp ) {
    printf ( "  WARNING: *** no free orbital parameters *** \n" );
    printf ( "  optimisation of orb. parameters will not be performed\n" );
    printf ( "  all input in this section is ignored\n\n" );
  }
  GETLNG( &nruns );
  printf ( "  number of independent optimisation runs = %ld\n", nruns );
  GETLNG( &niter );
  printf ( "  number of allowed iterations per run = %ld\n", niter );
  GETDBL( &stoprat );
  printf ( "  stop when simplex shrinked by factor = %lg\n", stoprat );

  printf ( "\n  SECTION 6: OUTPUT FILES\n\n" );

  printf ( "  write model spectra to file = " );
  GETSTR ( modfn ); printf ( "%s\n", modfn );
  printf ( "  write residuals to file = " );
  GETSTR ( resfn ); printf ( "%s\n", resfn );
  printf ( "  write radial velocities to file = " );
  GETSTR ( rvsfn ); printf ( "%s\n", rvsfn );
  printf ( "  write optimisation log to file = " );
  scanf ( "%s", logfn );
  if ( strlen ( logfn ) ) {
    printf ( "%s\n", logfn );
    if ( NULL == ( logfp = fopen ( logfn, "w" ) ) )
      DIE ( "could not open log file" );
    setbuf ( logfp, NULL );
  } else {
    printf ( "[no log file will be written]\n" );
  }

  printf ( "\n  SECTION 7: OPTIMISATION\n\n" );

  if ( 0 == nfp ) {                                            /* separation */

    double chi2;

    printf ( "  WARNING: *** no free orbital parameters *** \n" );
    chi2 = meritfn ( op0 );
    printf ( "  separation at the starting point:  chi2=%lg  gof=%.2lf\n",
      chi2, gsl_sf_gamma_inc_Q ( N*(M-K)/2.0, chi2/2.0 ) );
    dft_bck ( K, N, dftmod, mod+1 ); MxWrite ( mod, K+1, N, modfn );
    dft_bck ( M, N, dftres, res+1 ); MxWrite ( res, M+1, N, resfn );
    MxWrite( rvm, K, M, rvsfn );

  } else {                                                  /* disentangling */

    const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex;
    gsl_multimin_fminimizer *s = NULL; gsl_multimin_function minex_func;
    gsl_vector *ss, *x;
    size_t iter = 0, irun = 0, i; int status;
    double size0 = 1.0, size, chi2;
    const gsl_rng_type * R; gsl_rng * r;

    gsl_rng_env_setup (); R = gsl_rng_default; r = gsl_rng_alloc (R);

    minex_func.f = &meritfngsl;
    minex_func.n = nfp;
    minex_func.params = (void *) NULL;

    s = gsl_multimin_fminimizer_alloc (T, nfp);
    x  = gsl_vector_alloc (nfp); /* starting point */
    ss = gsl_vector_alloc (nfp); /* initial vertex size vector */

    for ( i = 0 ; i < nfp ; i++ ) {
      gsl_vector_set ( x, i, op0[opsw[i]] );
      gsl_vector_set ( ss, i, dop0[opsw[i]] );
      size0 *= dop0[opsw[i]];
    }

    chi2 = meritfngsl ( x, NULL );
    printf ( "  separation at the starting point:  chi2=%lg  gof=%.2lf\n",
      chi2, gsl_sf_gamma_inc_Q ( N*(M-K)/2.0, chi2/2.0 ) );
    dft_bck ( K, N, dftmod, mod+1 ); MxWrite ( mod, K+1, N, modfn );
    dft_bck ( M, N, dftres, res+1 ); MxWrite ( res, M+1, N, resfn );
    MxWrite( rvm, K, M, rvsfn );

    printf
      ( "  converged disentangling runs reported only if chi2 decreases\n" );

    for ( irun = 0 ; irun < nruns ; irun++ )  {

      for ( i = 0 ; i < nfp ; i++ ) {
        double tmp = 2.0 * ( gsl_rng_uniform (r) - 0.5 );
        gsl_vector_set ( x,  i, op0[opsw[i]] + tmp * dop0[opsw[i]] );
        gsl_vector_set ( ss, i, dop0[opsw[i]] );
      }

      gsl_multimin_fminimizer_set ( s, &minex_func, x, ss );
      status = GSL_FAILURE; size = size0;

      for ( iter = 0 ; iter < niter ; iter++ ) {
        if ( gsl_multimin_fminimizer_iterate (s) ) break;
        size = gsl_multimin_fminimizer_size (s);
        status = gsl_multimin_test_size (size, stoprat*size0);
        if ( GSL_SUCCESS == status ) break;
      }

      if ( (GSL_SUCCESS == status) && (s->fval < chi2) ) {
        printf ( "\n  irun=%d  iter=%d  sxshrnkf=%.2lg  chi2=%lg  ",
          irun+1, iter, size/size0, chi2 = s->fval );
        printf ( "gof=%.2lf\n",
          gsl_sf_gamma_inc_Q ( (N*(M-K)-nfp)/2.0, chi2/2.0 ) );
        for ( i = 0 ; i < nfp ; i++ ) {
          double xi = gsl_vector_get (s->x, i);
          if ( 2 == opsw[i] || 8 == opsw[i] ) xi = fabs ( xi / (1 + xi) );
          printf ( "  %40s = ", triorb_pname(opsw[i]) );
          printf ( triorb_pformat(opsw[i]), xi );
          if ( strcmp ( triorb_punit(opsw[i]), "1" ) )
            printf ( " %s", triorb_punit(opsw[i]) );
          printf ( "\n" );
        }
        dft_bck ( K, N, dftmod, mod+1 ); MxWrite ( mod, K+1, N, modfn );
        dft_bck ( M, N, dftres, res+1 ); MxWrite ( res, M+1, N, resfn );
        MxWrite( rvm, K, M, rvsfn );
      }

      if ( (GSL_SUCCESS == status) && strlen ( logfn ) ) {
        fprintf ( logfp, "%d  %d  %.2lg  %lg  ",
          irun+1, iter, size/size0, s->fval );
        fprintf ( logfp, "%.2lf  ",
          gsl_sf_gamma_inc_Q ( (N*(M-K)-nfp)/2.0, (s->fval)/2.0 ) );
        for ( i = 0 ; i < nfp ; i++ ) {
          double xi = gsl_vector_get (s->x, i);
          if ( 2 == opsw[i] || 8 == opsw[i] ) xi = fabs ( xi / (1 + xi) );
          fprintf ( logfp, "  " );
          fprintf ( logfp, triorb_pformat(i), xi );
        }
        fprintf ( logfp, "\n" );
      }

    }

    gsl_multimin_fminimizer_free (s);
    gsl_vector_free(x);
    gsl_vector_free(ss);

    printf ( "\n  completed %ld optimisation runs\n\n", nruns );

  }

  printf ( "  EXITING REGULARLY\n\n" );

  return EXIT_SUCCESS;

}

/*****************************************************************************/

double meritfngsl ( const gsl_vector *v, void *params )
{

  long i, j;
  double op[TRIORB_NP];

  for ( i = j = 0 ; i < TRIORB_NP ; i++ )
    op[i] = dop0[i] ? gsl_vector_get ( v, j++ ) : op0[i];

  return meritfn ( op );
}

/*****************************************************************************/

double meritfn ( double *opin )
{

  long j, k;
  double op[TRIORB_NP], rv[3];

  op[ 0] = opin[ 0];
  op[ 1] = opin[ 1];
  op[ 2] = fabs ( opin[2] / (1 + opin[2]) );
  op[ 3] = opin[ 3] * (M_PI/180);
  op[ 4] = opin[ 4];
  op[ 5] = opin[ 5];
  op[ 6] = opin[ 6];
  op[ 7] = opin[ 7];
  op[ 8] = fabs ( opin[8] / (1 + opin[8]) );
  op[ 9] = opin[ 9] * (M_PI/180);
  op[10] = opin[10] / rvstep;
  op[11] = opin[11] / rvstep;
  op[12] = opin[12] * (M_PI/180);

  for ( j = 0 ; j < M ; j++ ) {
    triorb_rv ( op, otimes[j], rv );
    for ( k = 0 ; k < K ; k++ )
      *(*(rvm+k)+j) = rv[ksw[k]] + *(rvcorr+j) / rvstep;
  }

  return fd3sep ( K, M, N, dftobs, sig, rvm, lfm, dftmod, dftres );
}

/*****************************************************************************/

max_stars_repo_name: mrawls/FDBinary-tools

max_issues_repo_name: mrawls/FDBinary-tools

max_forks_repo_name: mrawls/FDBinary-tools


text:
#pragma once

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

#include <MathLib/MathLibDll.h>

#include <MathLib/MathLib.h>
#include <MathLib/ThreeTuple.h>

#include <gsl/matrix/gsl_matrix.h>
#include <gsl/vector/gsl_vector.h>

#include <MathLib/Matrix.h>

#define VECTOR_AT(v, i) (*((v->data + (i * v->tda ))))

/*====================================================================================================================================================================*
 | This class will be used to represent matrices of arbitrary sizes (m rows by n columns) that have elements of type double. The underlying data strucutre used by    |
 | this class is gsl's (Gnu Scientific Library) matrix class. This class also makes use of the ATLAS implementation of BLAS for some operations such as matrix-matrix |
 | multiplication. This class is meant to improve performance, not necessarily ease of use.                                                                           |
 *====================================================================================================================================================================*/

class MATHLIB_DECLSPEC Vector : public Matrix{
public:
	/**
		constructor - creates an n row vector that is not initialized to any particular values
	*/
	Vector(int n);

	/**
		default constructor
	*/
	Vector();

	/**
		copy constructor - performs a deep copy of the matrix passed in as a parameter.
	*/
	Vector(const Vector& other);

	/**
		destructor.
	*/
	~Vector();

	/**
		copy operator - performs a deep copy of the Vector passed in as a parameter.
	*/
	Vector& operator=(const Vector &other);

	/**
		copy operator - performs a deep copy of the Vector passed in as a parameter.
	*/
	Vector& operator=(const Matrix &other);

	/**
		this method performs a shallow copy of the Vector that is passed in as a parameter.
	*/
	void shallowCopy(const Matrix& other);

	/**
		this method performs a deep copy of the vector that is passed in as a paramerer.
	*/
	void deepCopy(const Matrix& other);

	/**
		This method sets the current vector to be equal to one of the products: A * b or A'*b.
		The value of transA indicates if A is transposed or not
	*/
	void setToProductOf(const Matrix& A, const Matrix& B, bool transA = false, bool transB = false);

	/**
		This method sets the current vector to be equal to one of the rows of A - shallow column only!
	*/
	void setToRow(const Matrix& A, int row, int start = 0, int howMany = -1);

	/**
		This method sets the current vector to be equal to one of the cols of A - shallow column only!
	*/
	void setToCol(const Matrix& A, int col, int start = 0, int howMany = -1);

	/**
		This method prints the contents of the matrix - testing purpose only.
	*/
	void printVector() const;

	/**
		This method returns a copy of the value of the matrix at (i,j)
	*/
	double get(int i) const;

	/**
		This method sets the value of the matrix at (i,j) to newVal.
	*/
	void set(int i, double newVal);

	/**
		Computes the 2-norm squared for the current vector.
	*/
	inline double normSquared(){
		int r = getRowCount();
		double result = 0;
		for (int i=0;i<r;i++){
			double n = MATRIX_AT(matrix, i, 0) ;
			result += n*n;
		}
		return result;
	}

};

MATHLIB_DECLSPEC
void testVectorClass();

max_stars_repo_name: MontyThibault/centre-of-mass-awareness

max_issues_repo_name: MontyThibault/centre-of-mass-awareness

max_forks_repo_name: MontyThibault/centre-of-mass-awareness


text:
#pragma once

#include "arcana/expected.h"
#include "arcana/functional/inplace_function.h"
#include "arcana/iterators.h"
#include "arcana/type_traits.h"

#include "cancellation.h"

#include <gsl/gsl>
#include <memory>
#include <stdexcept>

#include <atomic>

namespace arcana
{
    template<typename ResultT, typename ErrorT>
    class task_completion_source;

    //
    // A scheduler that will invoke the continuation inline
    // right after the previous task.
    //
    namespace
    {
        constexpr auto inline_scheduler = [](auto&& callable)
        {
            callable();
        };
    }
}

#include "internal/internal_task.h"

namespace arcana
{
    //
    //  Generic task system to run work with continuations on a generic scheduler.
    //
    //  The scheduler on which tasks are queued must satisfy this contract:
    //
    //      struct scheduler
    //      {
    //          template<CallableT>
    //          void queue(CallableT&& callable)
    //          {
    //              callable must be usable like this: callable();
    //          }
    //      };
    //
    template<typename ResultT, typename ErrorT>
    class task
    {
        using payload_t = internal::task_payload_with_return<ResultT, ErrorT>;
        using payload_ptr = std::shared_ptr<payload_t>;

        static_assert(std::is_same<typename as_expected<ResultT, ErrorT>::value_type, ResultT>::value,
            "task can't be of expected<T>");

    public:
        using result_type = ResultT;
        using error_type = ErrorT;

        task() = default;

        task(const task& other) = default;
        task(task&& other) = default;

        task(const task_completion_source<ResultT, ErrorT>& source)
            : m_payload{ source.m_payload }
        {}

        task(task_completion_source<ResultT, ErrorT>&& source)
            : m_payload{ std::move(source.m_payload) }
        {}

        task& operator=(task&& other) = default;
        task& operator=(const task& other) = default;

        bool operator==(const task& other)
        {
            return m_payload == other.m_payload;
        }

        //
        // Executes a callable on this scheduler once this task is finished and
        // returns a task that represents the callable.
        //
        // Calling .then() on the returned task will queue a task to run after the
        // callable is run.
        //
        template<typename SchedulerT, typename CallableT>
        auto then(SchedulerT& scheduler, cancellation& token, CallableT&& callable)
        {
            using traits = internal::callable_traits<CallableT, result_type>;
            using wrapper = internal::input_output_wrapper<result_type, error_type, traits::handles_expected::value>;

            static_assert(error_priority<error_type>::value <= error_priority<typename traits::error_propagation_type>::value,
                "The error of a parent task needs to be convertible to the error of a child task.");

            static_assert(std::is_same_v<typename expected_error_or<typename traits::input_type, error_type>::type, error_type>,
                "Continuation expected input parameter needs to use the same error type as the parent task");

            auto factory{ internal::make_task_factory(
                internal::make_work_payload<typename traits::expected_return_type::value_type, typename traits::error_propagation_type>(
                    [callable = wrapper::wrap_callable(std::forward<CallableT>(callable), token)]
                    (internal::base_task_payload* self) mutable noexcept
                    {
                        return callable(*static_cast<payload_t*>(self)->Result);
                    })
            ) };

            m_payload->create_continuation([&scheduler](auto&& c)
            {
                scheduler(std::forward<decltype(c)>(c));
            }, m_payload, std::move(factory.to_run.m_payload));

            return factory.to_return;
        }

    private:
        explicit task(payload_ptr payload)
            : m_payload{ std::move(payload) }
        {}

        template<typename OtherResultT, typename OtherErrorT>
        friend class task;

        friend class task_completion_source<ResultT, ErrorT>;

        template<typename OtherErrorT, typename OtherResultT>
        friend struct internal::task_factory;

        template<typename SchedulerT, typename CallableT>
        friend auto make_task(SchedulerT& scheduler, cancellation& token, CallableT&& callable)
            -> typename internal::task_factory<
                            typename internal::callable_traits<CallableT, void>::error_propagation_type,
                            typename internal::callable_traits<CallableT, void>::expected_return_type::value_type>::task_t;

        payload_ptr m_payload;
    };
}

namespace arcana
{
    template<typename ResultT, typename ErrorT>
    class task_completion_source
    {
        using payload_t = internal::task_payload_with_return<ResultT, ErrorT>;
        using payload_ptr = std::shared_ptr<payload_t>;

    public:
        using result_type = ResultT;
        using error_type = ErrorT;

        task_completion_source()
            : m_payload{ std::make_shared<payload_t>() }
        {}

        //
        // Completes the task this source represents.
        //
        void complete()
        {
            static_assert(std::is_same<ResultT, void>::value,
                "complete with no arguments can only be used with a void completion source");
            m_payload->complete(basic_expected<void, ErrorT>::make_valid());
        }

        //
        // Completes the task this source represents.
        //
        template<typename ValueT>
        void complete(ValueT&& value)
        {
            m_payload->complete(std::forward<ValueT>(value));
        }

        //
        // Returns whether or not the current source has already been completed.
        //
        bool completed() const
        {
            return m_payload->completed();
        }

        //
        // Converts this task_completion_source to a task object for consumers to use.
        //
        task<ResultT, ErrorT> as_task() const &
        {
            return task<ResultT, ErrorT>{ m_payload };
        }

        task<ResultT, ErrorT> as_task() &&
        {
            return task<ResultT, ErrorT>{ std::move(m_payload) };
        }

    private:
        explicit task_completion_source(std::shared_ptr<payload_t> payload)
            : m_payload{ std::move(payload) }
        {}

        friend class task<ResultT, ErrorT>;
        friend class abstract_task_completion_source;

        template<typename E, typename R>
        friend struct internal::task_factory;

        payload_ptr m_payload;
    };

    //
    // a type erased version of task_completion_source.
    //
    class abstract_task_completion_source
    {
        using payload_t = internal::base_task_payload;
        using payload_ptr = std::shared_ptr<payload_t>;

    public:
        abstract_task_completion_source()
            : m_payload{}
        {}

        template<typename T, typename E>
        explicit abstract_task_completion_source(const task_completion_source<T, E>& other)
            : m_payload{ other.m_payload }
        {}

        template<typename T, typename E>
        explicit abstract_task_completion_source(task_completion_source<T, E>&& other)
            : m_payload{ std::move(other.m_payload) }
        {}

        //
        // Returns whether or not the current source has already been completed.
        //
        bool completed() const
        {
            return m_payload->completed();
        }

        template<typename T, typename E>
        bool operator==(const task_completion_source<T, E>& other)
        {
            return m_payload == other.m_payload;
        }

        template<typename T, typename E>
        task_completion_source<T, E> unsafe_cast()
        {
            return task_completion_source<T, E>{ std::static_pointer_cast<internal::task_payload_with_return<T, E>>(m_payload) };
        }

    private:
        payload_ptr m_payload;
    };

    //
    // creates a task and queues it to run on the given scheduler
    //
    template<typename SchedulerT, typename CallableT>
    inline auto make_task(SchedulerT& scheduler, cancellation& token, CallableT&& callable)
        -> typename internal::task_factory<
                            typename internal::callable_traits<CallableT, void>::error_propagation_type,
                            typename internal::callable_traits<CallableT, void>::expected_return_type::value_type>::task_t
    {
        using traits = internal::callable_traits<CallableT, void>;
        using wrapper = internal::input_output_wrapper<void, typename traits::error_propagation_type, false>;

        auto factory{ internal::make_task_factory(
            internal::make_work_payload<typename traits::expected_return_type::value_type, typename traits::error_propagation_type>(
                [callable = wrapper::wrap_callable(std::forward<CallableT>(callable), token)]
                (internal::base_task_payload*) mutable noexcept
                {
                    return callable(basic_expected<void, typename traits::error_propagation_type>::make_valid());
                })
        ) };

        scheduler([to_run = std::move(factory.to_run)]
        {
            to_run.m_payload->run(nullptr);
        });

        return factory.to_return;
    }

    //
    // creates a completed task from the given result
    //
    template<typename ErrorT, typename ResultT>
    inline task<typename as_expected<ResultT, ErrorT>::value_type, ErrorT> task_from_result(ResultT&& value)
    {
        task_completion_source<typename as_expected<ResultT, ErrorT>::value_type, ErrorT> result;
        result.complete(std::forward<ResultT>(value));
        return std::move(result);
    }

    template<typename ErrorT>
    inline task<void, ErrorT> task_from_result()
    {
        task_completion_source<void, ErrorT> result;
        result.complete();
        return std::move(result);
    }

    template<typename ResultT, typename ErrorT>
    inline task<ResultT, std::error_code> task_from_error(const ErrorT& error)
    {
        task_completion_source<ResultT, std::error_code> result;
        result.complete(make_unexpected(make_error_code(error)));
        return std::move(result);
    }

    template<typename ResultT>
    inline task<ResultT, std::error_code> task_from_error(const std::error_code& error)
    {
        task_completion_source<ResultT, std::error_code> result;
        result.complete(make_unexpected(error));
        return std::move(result);
    }

    template<typename ResultT>
    inline task<ResultT, std::exception_ptr> task_from_error(const std::exception_ptr& error)
    {
        task_completion_source<ResultT, std::exception_ptr> result;
        result.complete(make_unexpected(error));
        return std::move(result);
    }

    template<typename ErrorT>
    inline task<void, ErrorT> when_all(gsl::span<task<void, ErrorT>> tasks)
    {
        if (tasks.empty())
        {
            return task_from_result<ErrorT>();
        }

        struct when_all_data
        {
            std::mutex mutex;
            size_t pendingCount;
            std::error_code error;
        };

        task_completion_source<void, ErrorT> result;
        auto data = std::make_shared<when_all_data>();
        data->pendingCount = tasks.size();

        for (task<void, ErrorT>& task : tasks)
        {
            task.then(inline_scheduler, cancellation::none(), [data, result](const basic_expected<void, ErrorT>& exp) mutable noexcept {
                bool last = false;
                {
                    std::lock_guard<std::mutex> guard{ data->mutex };

                    data->pendingCount -= 1;
                    last = data->pendingCount == 0;

                    if (exp.has_error() && !data->error)
                    {
                        // set the first error, as it might have cascaded
                        data->error = exp.error();
                    }
                }

                if (last) // we were the last task to complete
                {
                    if (data->error)
                    {
                        result.complete(make_unexpected(data->error));
                    }
                    else
                    {
                        result.complete();
                    }
                }
            });
        }

        return std::move(result);
    }

    template<typename T, typename ErrorT>
    inline task<std::vector<T>, ErrorT> when_all(gsl::span<task<T, ErrorT>> tasks)
    {
        if (tasks.empty())
        {
            return task_from_result<ErrorT, std::vector<T>>(std::vector<T>());
        }

        struct when_all_data
        {
            std::mutex mutex;
            size_t pendingCount;
            std::error_code error;
            std::vector<T> results;
        };

        task_completion_source<std::vector<T>, ErrorT> result;
        auto data = std::make_shared<when_all_data>();
        data->pendingCount = tasks.size();
        data->results.resize(tasks.size());

        //using forloop with index to be able to keep proper order of results
        for (auto idx = 0U; idx < data->results.size(); idx++)
        {
            tasks[idx].then(arcana::inline_scheduler, cancellation::none(), [data, result, idx](const basic_expected<T, ErrorT>& exp) mutable noexcept {
                bool last = false;
                {
                    std::lock_guard<std::mutex> guard{ data->mutex };

                    data->pendingCount -= 1;
                    last = data->pendingCount == 0;

                    if (exp.has_error() && !data->error)
                    {
                        // set the first error, as it might have cascaded
                        data->error = exp.error();
                    }
                    if (exp.has_value())
                    {
                        data->results[idx] = exp.value();
                    }
                }

                if (last) // we were the last task to complete
                {
                    if (data->error)
                    {
                        result.complete(make_unexpected(data->error));
                    }
                    else
                    {
                        result.complete(data->results);
                    }
                }
            });
        }

        return std::move(result);
    }

    template<typename ErrorT, typename... ArgTs>
    inline task<std::tuple<typename arcana::void_passthrough<ArgTs>::type...>, ErrorT> when_all(task<ArgTs, ErrorT>... tasks)
    {
        using void_passthrough_tuple = std::tuple<typename arcana::void_passthrough<ArgTs>::type...>;

        struct when_all_data
        {
            std::mutex mutex;
            int pending;
            ErrorT error;
            void_passthrough_tuple results;
        };

        task_completion_source<void_passthrough_tuple, ErrorT> result;

        auto data = std::make_shared<when_all_data>();
        data->pending = std::tuple_size<void_passthrough_tuple>::value;

        std::tuple<task<ArgTs, ErrorT>&...> taskrefs = std::make_tuple(std::ref(tasks)...);

        iterate_tuple(taskrefs, [&](auto& task, auto idx) {
            using task_t = std::remove_reference_t<decltype(task)>;

            task.then(inline_scheduler,
                cancellation::none(),
                [data, result](const basic_expected<typename task_t::result_type, ErrorT>& exp) mutable noexcept {
                bool last = false;
                {
                    std::lock_guard<std::mutex> guard{ data->mutex };

                    data->pending -= 1;
                    last = data->pending == 0;

                    internal::write_expected_to_tuple<decltype(idx)::value, ErrorT>(data->results, exp);

                    if (exp.has_error() && !data->error)
                    {
                        // set the first error, as it might have cascaded
                        data->error = exp.error();
                    }
                }

                if (last) // we were the last task to complete
                {
                    if (data->error)
                    {
                        result.complete(make_unexpected(data->error));
                    }
                    else
                    {
                        result.complete(std::move(data->results));
                    }
                }

                return basic_expected<void, ErrorT>::make_valid();
            });
        });
        return std::move(result);
    }
}


text:
/* randist/exppow.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 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., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

#include <config.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

/* The exponential power probability distribution is

   p(x) dx = (1/(2 a Gamma(1+1/b))) * exp(-|x/a|^b) dx

   for -infty < x < infty. For b = 1 it reduces to the Laplace
   distribution.

   The exponential power distribution is related to the gamma
   distribution by E = a * pow(G(1/b),1/b), where E is an exponential
   power variate and G is a gamma variate.

   We use this relation for b < 1. For b >=1 we use rejection methods
   based on the laplace and gaussian distributions which should be
   faster.

   See P. R. Tadikamalla, "Random Sampling from the Exponential Power
   Distribution", Journal of the American Statistical Association,
   September 1980, Volume 75, Number 371, pages 683-686.

*/

double
gsl_ran_exppow (const gsl_rng * r, const double a, const double b)
{
  if (b < 1)
    {
      double u = gsl_rng_uniform (r) ;
      double v = gsl_ran_gamma (r, 1/b, 1.0) ;
      double z = a * pow(v, 1/b) ;

      if (u > 0.5)
	{
	  return z ;
	}
      else
	{
	  return -z ;
	}
    }
  else if (b == 1)
    {
      /* Laplace distribution */
      return gsl_ran_laplace (r, a) ;
    }
  else if (b < 2)
    {
      /* Use laplace distribution for rejection method */

      double x, y, h, ratio, u ;

      /* Scale factor chosen by upper bound on ratio at b = 2 */

      double s = 1.4489 ;
      do
	{
	  x = gsl_ran_laplace (r, a) ;
	  y = gsl_ran_laplace_pdf (x,a) ;
	  h = gsl_ran_exppow_pdf (x,a,b) ;
	  ratio = h/(s * y) ;
	  u = gsl_rng_uniform (r) ;
	}
      while (u > ratio) ;

      return x ;
    }
  else if (b == 2)
    {
      /* Gaussian distribution */
      return gsl_ran_gaussian (r, a/sqrt(2.0)) ;
    }
  else
    {
      /* Use gaussian for rejection method */

      double x, y, h, ratio, u ;
      const double sigma = a / sqrt(2.0) ;

      /* Scale factor chosen by upper bound on ratio at b = infinity.
	 This could be improved by using a rational function
	 approximation to the bounding curve. */

      double s = 2.4091 ;  /* this is sqrt(pi) e / 2 */

      do
	{
	  x = gsl_ran_gaussian (r, sigma) ;
	  y = gsl_ran_gaussian_pdf (x, sigma) ;
	  h = gsl_ran_exppow_pdf (x, a, b) ;
	  ratio = h/(s * y) ;
	  u = gsl_rng_uniform (r) ;
	}
      while (u > ratio) ;

      return x;
    }
}

double
gsl_ran_exppow_pdf (const double x, const double a, const double b)
{
  double p ;
  double lngamma = gsl_sf_lngamma (1+1/b) ;
  p = (1/(2*a)) * exp(-pow(fabs(x/a),b) - lngamma);
  return p;
}