text:
/*
    Code for the article 'Landscape heterogeneity buffers biodiversity of meta-food-webs under global change through rescue and drainage effects'

    Copyright (C) 2016 Christian Guill & Florian D. Schneider
    Copyright (C) 2020 Remo Ryser

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.

    last edit: 01.03.2020 by RR

 */

/*
// ************** Local foodweb code extended into space ***************************

Remo Ryser, Johanna Häussler, Markus Stark

*/

#include <stdio.h>
#include <math.h>
#include <string.h>
#include <fstream>      // for input/output from or to files
//#include <stdlib.h>
#include <iostream>		// for input/output on terminal
#include <sstream>
#include <vector>
#include <cstdlib>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_rng.h>				// random number generator
#include <gsl/gsl_randist.h>			// random number distributions
#include <gsl/gsl_blas.h>				// linear algebra routines
#include <gsl/gsl_linalg.h>				// linear algebra
#include <gsl/gsl_sort_vector.h>		// vector operations
#include <gsl/gsl_odeiv.h>              // ODE solver

#include <cvode/cvode.h>             /* prototypes for CVODE fcts., consts. */
#include <nvector/nvector_serial.h>  /* serial N_Vector types, fcts., macros */
#include <cvode/cvode_dense.h>       /* prototype for CVDense */

double *web_calc(gsl_rng *r);
static void dynamics(double B[], double Bdot[], void *params, double t);
static void pdef_structure(gsl_rng *r,gsl_matrix *Ap, gsl_vector *mass, gsl_matrix *Up, gsl_vector *DmaxPlant);
static void output(gsl_matrix *SW, gsl_vector *mass, gsl_vector *con, double iniB[], double meanB[], double meanB_tot[], double CV[], double CV_tot[], double params[], int paramslength, double Biomass_Temp[], double RGG_params[], double net_growth_Temp[]);
static void getInputSpp();
static void getInputRicker();
static void getInputMass(gsl_vector *mass, gsl_matrix *Up,gsl_vector *DmaxPlant);
static void getInputLandscape();
static void getInputCoords(gsl_vector *coords_in_v);
static void getInputSW(gsl_matrix *SW_input);
static void getInputNutrient(gsl_vector *NS_input);

double get_random_parameter(double mean, double sigma, double low_cutoff, double high_cutoff ,gsl_rng *r);
static void mean_body_mass(double meanB[], gsl_vector *mass, double mean_body_masses[]);

static void set_parameters(gsl_rng *r, gsl_matrix *Ap, gsl_vector *mass, double params[], gsl_matrix *Up);
static void set_spatial_parameters(gsl_rng *r, gsl_matrix *SW, gsl_vector *Loc, gsl_vector *NS, double params[]);
static void initialise_biomass_densities(double B[], double iniB[], double params[], gsl_rng *r);
static void solve_ode(double B[], double meanB[],  double meanB_tot[], double CV[],  double CV_tot[], double params[], gsl_matrix *CovM, double Biomass_Temp[], double net_growth_Temp[]);
static void Extinction(double B[], int Num);
static void Extinct_Species(double B[], int Num, int Pat);
static int Cvode_rates(realtype t, N_Vector y, N_Vector ydot, void *params);

static void net_rates(gsl_vector *TBvec, gsl_vector *Ivec, gsl_vector *Dvec, gsl_vector *Xvec, gsl_vector *Emvec, gsl_vector *N_in, gsl_vector *N_out, gsl_vector *N_up, void *params, double t);

static void Prepare_timeseries_file(FILE *timeseries);
static void Write_timeseries_to_file(double B[], double t, FILE *timeseries);
static void Prepare_netrates_file(FILE *netrates);
static void Write_netrates_to_file(int p, gsl_vector *netvec, gsl_vector *Emvec_p, FILE *netrates);
static void Prepare_disprates_file(FILE *disprates);
static void Write_disprates_to_file(int i, int j, gsl_vector *Imvec_i, FILE *disprates);
static void Write_params_to_file(double params[], int paramslength);
static void Write_landscape_to_file(double params[]);
static void Write_SW_matrix_to_file(gsl_matrix *SW);

static void Generate_RGG_structure(gsl_rng *ran, gsl_vector *mass, gsl_vector *Loc, gsl_matrix *SW, gsl_matrix *SW_input, gsl_vector *con, double params[], double RGG_params[], gsl_vector *DmaxPlant);

static void calc_patch_distances(gsl_matrix *P_Dist, gsl_vector *Loc, double RGG_params[]);
static void calc_spatial_connectance(gsl_matrix *SW, gsl_vector *con, double RGG_params[]);

static void show_matrix(gsl_matrix *A, int Num, int N2);
static void show_vector(gsl_vector *A, int Num);
double calc_sd(gsl_vector *con, int START, int P, int div);
double calc_mean(gsl_vector *con, int START, int P, int div);

max_stars_repo_name: RemoRyser/Metafoodweb

max_issues_repo_name: RemoRyser/Metafoodweb

max_forks_repo_name: RemoRyser/Metafoodweb


text:
#include "readArray_MPI.h"
#include "../include/paralleltt.h"

#include <mpi.h>
#include <stdio.h>
#include <stdlib.h>
#include <lapacke.h>
#include <time.h>
#include <math.h>

void* p_hilbert_init(int d, int* n){
    p_hilbert* parameters = (p_hilbert*) malloc(sizeof(p_hilbert));
    parameters->d = d;
    parameters->n = (int*) malloc(d*sizeof(int));
    for (int ii = 0; ii < d; ++ii){
        parameters->n[ii] = n[ii];
    }
    parameters->n_stream = (int*) calloc(d, sizeof(int));
    parameters->n_prod = (long*) malloc(d * sizeof(long));
    return (void*) parameters;
}

void p_hilbert_free(void* parameters){
    p_hilbert* p_casted = (p_hilbert*) parameters;
    free(p_casted->n); p_casted->n = NULL;
    free(p_casted->n_stream); p_casted->n_stream = NULL;
    free(p_casted->n_prod); p_casted->n_prod = NULL;
    free(p_casted);
}

void f_hilbert(double* restrict X, int* ind1, int* ind2, const void* parameters)
{
    p_hilbert* p_casted = (p_hilbert*) parameters;
    int d = p_casted->d;
    int* n = p_casted->n;
    int* n_stream = p_casted->n_stream;
    long* n_prod = p_casted->n_prod;

    long N_stream = 1;
    for (int ii = 0; ii < d; ++ii){
        n_stream[ii] = ind2[ii] - ind1[ii];
        N_stream = N_stream*n_stream[ii];
    }

    n_prod[0] = 1;
    for (int ii = 1; ii < d; ++ii){
        n_prod[ii] = (long) n_prod[ii-1] * n_stream[ii-1];
    }

    int bias = 1;
    for (int jj = 0; jj < d; ++jj){
        bias = bias + ind1[jj];
    }

    for (long ii = 0; ii < N_stream; ++ii){
        int xinv = bias;
		long tmp = ii;
		for (int jj = d-1; jj >= 0; jj--){
			xinv += tmp/n_prod[jj];
			tmp = tmp % n_prod[jj];
		}

        X[ii] = (double) 1/xinv;
    }
}

void* p_gaussian_bumps_init(int d, int* n, int M, double gamma, double* region, double* centers)
{
    p_gaussian_bumps* parameters = (p_gaussian_bumps*) malloc(sizeof(p_gaussian_bumps));
    parameters->d = d;
    parameters->n = (int*) calloc(d, sizeof(int));
    for (int ii = 0; ii < d; ++ii){
        parameters->n[ii] = n[ii];
    }

    parameters->M = M;
    parameters->gamma = gamma;

    parameters->region = (double*) calloc(2*d, sizeof(double));
    for (int ii = 0; ii < 2*d; ++ii){
        parameters->region[ii] = region[ii];
    }
    parameters->centers = (double*) calloc(d*M, sizeof(double));
    for (int ii = 0; ii < d*M; ++ii){
        parameters->centers[ii] = centers[ii];
    }

//    parameters->ten_ind_ii = (int*) calloc(d, sizeof(int));
    parameters->x_ii = (double*) calloc(d, sizeof(double));
    parameters->n_stream = (int*) calloc(d, sizeof(double));
    parameters->n_prod = (long*) calloc(d, sizeof(long));

    return (void*) parameters;
}

void* unit_random_p_gaussian_bumps_init(int d, int* n, int M, double gamma, int seed)
{
    p_gaussian_bumps* parameters = (p_gaussian_bumps*) malloc(sizeof(p_gaussian_bumps));
    parameters->d = d;
    parameters->n = (int*) calloc(d, sizeof(int));
    for (int ii = 0; ii < d; ++ii){
        parameters->n[ii] = n[ii];
    }

    parameters->M = M;
    parameters->gamma = gamma;

    parameters->region = (double*) calloc(2*d, sizeof(double));
    for (int ii = 0; ii < d; ++ii){
        parameters->region[2*ii] = -1.0;
        parameters->region[2*ii+1] = 1.0;
    }

    parameters->centers = (double*) calloc(d*M, sizeof(double));
    srand(seed);
    int r1 = rand()%4096, r2 = rand()%4096, r3 = rand()%4096, r4 = rand()%4096;
    int iseed[4] = {r1, r2, r3, r4+(r4%2 == 0?1:0)};
    LAPACKE_dlarnv(2, iseed, d*M, parameters->centers);

    parameters->x_ii = (double*) calloc(d, sizeof(double));
    parameters->n_stream = (int*) calloc(d, sizeof(double));
    parameters->n_prod = (long*) calloc(d, sizeof(long));

    return (void*) parameters;
}

void p_gaussian_bumps_free(void* parameters)
{
    p_gaussian_bumps* p_casted = (p_gaussian_bumps*) parameters;

    free(p_casted->n);          p_casted->n = NULL;
    free(p_casted->region);     p_casted->region = NULL;
    free(p_casted->centers);    p_casted->centers = NULL;
    free(p_casted->x_ii);       p_casted->x_ii = NULL;
    free(p_casted->n_stream);   p_casted->n_stream = NULL;
    free(p_casted->n_prod);     p_casted->n_prod = NULL;
    free(p_casted);
}

void f_gaussian_bumps(double* restrict X, int* ind1, int* ind2, const void* parameters)
{
    p_gaussian_bumps* p_casted = (p_gaussian_bumps*) parameters;
    int d = p_casted->d;
    int* n = p_casted->n;
    int M = p_casted->M;
    double gamma = p_casted->gamma;
    double* region = p_casted->region;
    double* centers = p_casted->centers;
    double* x_ii = p_casted->x_ii;
    int* n_stream = p_casted->n_stream;
    long* n_prod = p_casted->n_prod;

    long N_stream = 1;
    for (int ii = 0; ii < d; ++ii){
        n_stream[ii] = ind2[ii] - ind1[ii];
        N_stream = N_stream*n_stream[ii];
    }

    n_prod[0] = 1;
    for (int ii = 1; ii < d; ++ii){
        n_prod[ii] = (long) n_prod[ii-1] * n_stream[ii-1];
    }

    for (long ii = 0; ii < N_stream; ++ii){
        long tmp = ii;
        for (int jj = d-1; jj >= 0; jj--){
            int ind_jj = tmp/n_prod[jj];
            x_ii[jj] = region[2*jj] + (ind_jj + ind1[jj]) * (region[2*jj + 1] - region[2*jj]) / (n[jj] - 1);
			tmp = tmp % n_prod[jj];
		}

        X[ii] = 0;
        for (int kk = 0; kk < M; ++kk){
            double exponent = 0;
            for (int jj = 0; jj < d; ++jj){
                exponent += (x_ii[jj] - centers[kk*d + jj]) * (x_ii[jj] - centers[kk*d + jj]);
            }
            X[ii] += exp(-gamma * exponent);
        }
    }
}

void* p_arithmetic_init(int d, int* n)
{
    p_arithmetic* parameters = (p_arithmetic*) malloc(sizeof(p_arithmetic));
    parameters->d = d;
    parameters->n = (int*) malloc(d*sizeof(int));
    for (int ii = 0; ii < d; ++ii){
        parameters->n[ii] = n[ii];
    }
    return (void*) parameters;
}

void p_arithmetic_free(void* parameters)
{
    p_arithmetic* p_casted = (p_arithmetic*) parameters;
    free(p_casted->n); p_casted->n = NULL;
    free(p_casted);
}

void f_arithmetic(double* restrict X, int* ind1, int* ind2, const void* parameters)
{
    p_arithmetic* p_casted = (p_arithmetic*) parameters;
    int d = p_casted->d;
    int* n = p_casted->n;

    int* n_stream = (int*) calloc(d, sizeof(int));
    long N_stream = 1;
    for (int ii = 0; ii < d; ++ii){
        n_stream[ii] = ind2[ii] - ind1[ii];
        N_stream = N_stream*n_stream[ii];
    }

    int* ind_ii = (int*) calloc(d, sizeof(int));
    for (int ii = 0; ii < N_stream; ++ii){
        to_tensor_ind(ind_ii, ii, n_stream, d);
        for (int jj = 0; jj < d; ++jj){
            ind_ii[jj] = ind_ii[jj] + ind1[jj];
        }
        X[ii] = 1 + to_vec_ind(ind_ii, n, d);
//        X[ii] = X[ii] + (0.000000001 / X[ii]);
    }

    free(n_stream);
    free(ind_ii);
}

void* p_tt_init(tensor_train* tt)
{
    p_tt* parameters = (p_tt*) malloc(sizeof(p_tt));
    parameters->tt = tt;
    return (void*) parameters;
}

// Doesn't actually free the tensor train. You gotta do that yourself
void p_tt_free(void* parameters)
{
    p_tt* p_casted = (p_tt*) parameters;
    p_casted->tt = NULL;
    free(p_casted);
}

void f_tt(double* restrict X, int* ind1, int* ind2, const void* parameters)
{
    p_tt* p_casted = (p_tt*) parameters;
    tensor_train* tt = p_casted->tt;
    int d = tt->d;
    int* n = tt->n;
    int* r = tt->r;

    matrix_tt** train_submats = (matrix_tt**) calloc(d, sizeof(matrix_tt*));
    for (int ii = 0; ii < d; ++ii){
        matrix_tt* train_mat_ii = matrix_tt_wrap(r[ii]*n[ii], r[ii+1], tt->trains[ii]);
        train_submats[ii] = submatrix(train_mat_ii, ind1[ii]*r[ii], ind2[ii]*r[ii], 0, r[ii+1]);
        free(train_mat_ii);
    }

    matrix_tt* mult = submatrix_copy(train_submats[d-1]);
    matrix_tt_reshape(r[d-1], ind2[d-1] - ind1[d-1], mult);

    for (int jj = d-2; jj >= 1; --jj){
        matrix_tt* new_mult = matrix_tt_init(train_submats[jj]->m, mult->n);
        matrix_tt_dgemm(train_submats[jj], mult, new_mult, 1.0, 0.0);
        matrix_tt_reshape(r[jj], (new_mult->n) * (new_mult->m) / r[jj], new_mult);
        matrix_tt_free(mult); mult = new_mult;
    }

    long X_m = 1;
    for (int ii = 1; ii < d; ++ii){
        X_m = X_m*(ind2[ii] - ind1[ii]);
    }

    matrix_tt* X_mat = matrix_tt_wrap(ind2[0] - ind1[0], X_m, X);
    matrix_tt_dgemm(train_submats[0], mult, X_mat, 1.0, 0.0);

    for (int ii = 0; ii < d; ++ii){
        free(train_submats[ii]); train_submats[ii] = NULL;
    }

    free(train_submats); train_submats = NULL;
    matrix_tt_free(mult);   mult = NULL;
    free(X_mat);         X_mat = NULL;
}

double tt_error(tensor_train* tt, MPI_tensor* ten)
{
    MPI_Comm comm = ten->comm;
    tt_broadcast(comm, tt);

    void* parameters_tt = p_tt_init(tt);
    MPI_tensor* ten_tt = MPI_tensor_init(ten->d, ten->n, ten->nps, ten->comm, &f_tt, parameters_tt);
    int rank = ten_tt->rank;
    int* schedule_rank = ten_tt->schedule[rank];

    double true_norm_squared = 0;
    double diff_norm_squared = 0;

    flattening_info* fi = flattening_info_init(ten, 0, 1, 0);
    for (int ii = 0; ii < ten->n_schedule; ++ii){
        int block = schedule_rank[ii];

        if (block != -1){
            stream(ten, block);
            stream(ten_tt, block);

            flattening_info_update(fi, ten, block);
            long N = fi->t_N;

            for (int jj = 0; jj < N; ++jj){
                ten_tt->X[jj] = ten_tt->X[jj] - ten->X[jj];
            }

            matrix_tt* mat_true = matrix_tt_wrap(N, 1, ten->X);
            double tmp = frobenius_norm(mat_true);
            true_norm_squared = true_norm_squared + tmp*tmp;

            matrix_tt* mat_diff = matrix_tt_wrap(N, 1, ten_tt->X);
            tmp = frobenius_norm(mat_diff);
            diff_norm_squared = diff_norm_squared + tmp*tmp;
//            printf("r%d ii%d true_norm_squared = %e, diff_norm_squared = %e\n", rank, ii, true_norm_squared, diff_norm_squared);

            free(mat_true); mat_true = NULL;
            free(mat_diff); mat_diff = NULL;
        }
    }
    int head = 0;

    double true_norm_reduced;
    double diff_norm_reduced;
    MPI_Reduce(&true_norm_squared, &true_norm_reduced, 1, MPI_DOUBLE, MPI_SUM, head, comm);
    true_norm_reduced = sqrt(true_norm_reduced);

    MPI_Reduce(&diff_norm_squared, &diff_norm_reduced, 1, MPI_DOUBLE, MPI_SUM, head, comm);
    diff_norm_reduced = sqrt(diff_norm_reduced);

    double rel_err = diff_norm_reduced/true_norm_reduced;

    p_tt_free(parameters_tt); parameters_tt = NULL;
    flattening_info_free(fi); fi = NULL;
    MPI_tensor_free(ten_tt); ten_tt = NULL;

    return rel_err;
}

ext: c


text:
#include <stdlib.h>
#include <stdio.h>
#include <stdint.h>
#include <dc1394/dc1394.h>
#include <math.h>
#include <fitsio.h>
#include <gsl/gsl_statistics.h>
#include <xpa.h>
#include <time.h>
#include <sys/time.h>
#include <string.h>
#include <inttypes.h>

#define NXPA 10

typedef struct _Box {
    double  x;
    double  y;
    double  fwhm;
    double  cenx;
    double  ceny;
    double  counts;
    double  background;
    double  noise;
    double  sigmaxythresh;
    double  snr;
    double  sigmaxy;
    double  sigmafwhmthresh;
    double  sigmafwhm;
    double  strehl;
    double  max;
    int     r;
} Box;

typedef struct _Back {
    double x;
    double y;
    double background;
    double sigma;
    int r;
    int width;
} Back;

Box box[3];
Back back;
long nelements, naxes[2], fpixel;
int boxsize = 17;

/* MASSDIMM values */
/*  20101228:  determined DIMM pixel scale to be 1.78"/pixel using HR1866
    and HR8895 as astrometric double stars
*/
double pixel_scale = 1.78;
double d = 0.080;
double r = 0.170;

double lambda = 0.6e-6;

double stardist(int i, int j) {
    return( sqrt( (box[i].cenx-box[j].cenx)*(box[i].cenx-box[j].cenx) +
		  (box[i].ceny-box[j].ceny)*(box[i].ceny-box[j].ceny) ) );
}

/* this routine uses a. tokovinin's modified equation given in
   2002, PASP, 114, 1156
 */
double seeing(double var) {
    double b, K, seeing;

    b = r/d;

    /* pixel scale in "/pixel, convert to rad */
    var = var * pow(pixel_scale/206265.0, 2);

    K = 0.364 * (1.0 - 0.532 * pow(b, -1/3) - 0.024 * pow(b, -7/3));

    seeing = 206265.0 * 0.98 * pow(d/lambda, 0.2) * pow(var/K, 0.6);

    return seeing;
}

/* this routine uses the classic DIMM equation */
double old_seeing(double var) {
    double r0;

    /* pixel scale in "/pixel, convert to rad */
    var = var * pow(pixel_scale/206265.0, 2);

    r0 = pow(2.0 * (lambda * lambda) * (
		 ( 0.1790 * pow(d, -1.0/3.0) -
		   0.0968 * pow(r, -1.0/3.0) ) / var
		 ), 0.6);

    // return 206265.0*0.98*lambda/r0;
    return r0;
}

/* measure the background in an annulus around the spot pattern */
int background(unsigned char *image, int imwidth, int imheight) {
    int i, j, backpix;
    int low_y, up_y, low_x, up_x;
    double dist, sum, sumsq;

    backpix = 0;
    sum = 0.0;
    sumsq = 0.0;

    low_y = back.y - back.r - back.width;
    up_y  = back.y + back.r + back.width;
    low_x = back.x - back.r - back.width;
    up_x  = back.x + back.r + back.width;
    if (low_y < 0) {
	low_y = 0;
    }
    if (up_y >= imheight) {
	up_y = imheight;
    }
    if (low_x < 0) {
	low_x = 0;
    }
    if (up_x >= imwidth) {
	up_x = imwidth;
    }

    for (i=low_y; i<up_y; i++) {
	for (j=low_x; j<up_x; j++) {
	    dist = sqrt(pow(back.x-j, 2) + pow(back.y-i, 2));
	    if (dist >= back.r && dist <= back.r+back.width) {
		sum += image[i*naxes[0]+j];
		backpix++;
	    }
	}
    }

    back.background = sum/backpix;
    //back.background = 0.0;

    for (i=low_y; i<up_y; i++) {
	for (j=low_x; j<up_x; j++) {
	    dist = sqrt(pow(back.x-j, 2) + pow(back.y-i, 2));
	    if (dist >= back.r && dist <= back.r+back.width) {
		sumsq += (image[i*naxes[0]+j]-back.background)*
		    (image[i*naxes[0]+j]-back.background);
	    }
	}
    }

    back.sigma = sqrt(sumsq/backpix);
    //back.sigma = 1.0;

    return 1;
}

/* measure centroid using center-of-mass algorithm */
int centroid(unsigned char *image, int imwidth, int imheight, int num) {
    int i, j;
    double  max    = 0.0;
    double  sum    = 0.0;
    double  sumx   = 0.0;
    double  sumxx  = 0.0;
    double  sumy   = 0.0;
    double  sumyy  = 0.0;
    double  val = 0.0;
    double  gain = 0.9;
    double  rmom;
    double  dist, dx;
    double  nsigma = 3.0;
    int low_y, up_y, low_x, up_x;
    int sourcepix = 0;

    dx = pixel_scale/206265.0;

    low_y = box[num].y - box[num].r;
    up_y  = box[num].y + box[num].r;
    low_x = box[num].x - box[num].r;
    up_x  = box[num].x + box[num].r;
    if (low_y < 0) {
	low_y = 0;
    }
    if (up_y >= imheight) {
	up_y = imheight;
    }
    if (low_x < 0) {
	low_x = 0;
    }
    if (up_x >= imwidth) {
	up_x = imwidth;
    }

    for (i=low_y; i<up_y; i++) {
	for (j=low_x; j<up_x; j++) {
	    dist = sqrt(pow(box[num].x-j, 2) + pow(box[num].y-i, 2));
	    if (dist <= box[num].r) {
		val = image[i*naxes[0]+j] - back.background;
		if (val >= nsigma*back.sigma) {
		    if (val > max) {
			max = val;
		    }
		    sum   += val;
		    sumx  += val*j;
		    sumxx += val*j*j;
		    sumy  += val*i;
		    sumyy += val*i*i;
		    sourcepix++;
		}
	    }
	}
    }

    if ( sum <= 0.0 ) {
	box[num].sigmaxy = -1.0;
	box[num].sigmafwhm = -1.0;
	box[num].snr = 0.0;
	box[num].fwhm = -1.0;
	box[num].counts = 1.0;
	box[num].strehl = 0.0;
    } else {
	box[num].max = max;
	box[num].strehl = (max/sum)*(4.0/3.14159)*pow( lambda/(d*dx), 2);
	rmom = ( sumxx - sumx * sumx / sum + sumyy - sumy * sumy / sum ) / sum;

	if ( rmom <= 0 ) {
	    box[num].fwhm = -1.0;
    } else {
	    box[num].fwhm = sqrt(rmom)  * 2.354 / sqrt(2.0);
    }

    box[num].counts = sum;
    box[num].cenx   = sumx / sum;
    box[num].ceny   = sumy / sum;
    box[num].noise  = back.sigma*sourcepix;
    box[num].x += gain*(box[num].cenx - box[num].x);
    box[num].y += gain*(box[num].ceny - box[num].y);
    box[num].snr = sum/sqrt(box[num].noise*box[num].noise + sum);
    box[num].sigmaxy = 1.0 / box[num].snr / sqrt(6.0);
    box[num].sigmafwhm = back.sigma * pow(sourcepix,1.5) / 10.
	/ box[num].fwhm / box[num].counts
	* 2.354 * 2.354 / 2.0;
  }

  return 1;

}

int grab_frame(dc1394camera_t *c, unsigned char *buf, int nbytes) {
    dc1394video_frame_t *frame=NULL;
    dc1394error_t err;

    err = dc1394_capture_dequeue(c, DC1394_CAPTURE_POLICY_WAIT, &frame);
    if (err != DC1394_SUCCESS) {
	dc1394_log_error("Unable to capture.");
	dc1394_capture_stop(c);
	dc1394_camera_free(c);
	exit(1);
    }
    memcpy(buf, frame->image, nbytes);
    dc1394_capture_enqueue(c, frame);
    return 1;
}

int add_gaussian(unsigned char *buffer, float cenx, float ceny, float a, float sigma) {
    float gauss, rsq;
    int i, j, low_x, up_x, low_y, up_y, size;
    double xoff, yoff;

    xoff = 2.0*drand48() - 1.0;
    yoff = 2.0*drand48() - 1.0;

    cenx += 0.5*xoff;
    ceny += 0.5*yoff;

    size = 30;
    low_x = (int)(cenx-size);
    up_x = (int)(cenx+size);
    low_y = (int)(ceny-size);
    up_y = (int)(ceny+size);

    for (i=low_y; i<up_y; i++) {
	for (j=low_x; j<up_x; j++) {
	    rsq = (cenx - j)*(cenx - j) + (ceny - i)*(ceny - i);
	    gauss = a*expf(-1.0*rsq/(sigma*sigma));
	    if (gauss > 255)
		gauss = 255;
	    buffer[i*naxes[0]+j] += (unsigned char)gauss;
	}
    }

    return 1;
}

int main(int argc, char *argv[]) {
    dc1394_t * dc;
    dc1394camera_t * camera;
    dc1394camera_list_t * list;
    dc1394error_t err;
    dc1394video_mode_t mode;
    unsigned int max_height, max_width, winleft, wintop;
    uint64_t total_bytes = 0;
    unsigned char *buffer, *buffer2, *average;
    char *names[NXPA];
    char *messages[NXPA];
    char filename[256], xpastr[256];
    char *froot;
    fitsfile *fptr;
    int i, j, f, fstatus, status, nimages, nboxes, test, xsize, ysize;
    int nbad = 0, nbad_l = 0;
    FILE *init, *out, *cenfile;
    float xx = 0.0, yy = 0.0, xsum = 0.0, ysum = 0.0;
    float elapsed_time, fps, avemax;
    double *dist, *sig, *dist_l, *sig_l, *weight, *weight_l;
    double mean, var, var_l, avesig, airmass, exptime;
    double r0, seeing_short, seeing_long, seeing_ave;
    struct timeval start_time, end_time;
    struct tm ut;
    time_t start_sec, end_sec;
    suseconds_t start_usec, end_usec;

    srand48((unsigned)time(NULL));

    XPA xpa;
    xpa = XPAOpen(NULL);

    stderr = freopen("measure_seeing.log", "w", stderr);

    if (argc <= 3) {
	printf("Usage: measure_seeing <n> <airmass> <exptime>\n");
	exit(-1);
    }

    nimages = atoi(argv[1]);
    airmass = atof(argv[2]);
    exptime = atof(argv[3]);
    fstatus = 0;
    status = 0;
    xsize = 320;
    ysize = 240;
    naxes[0] = xsize;
    naxes[1] = ysize;
    fpixel = 1;

    fps = 1.0/exptime;
    if (fps > 330.0) {
	fps = 330.0;
    }

    nelements = naxes[0]*naxes[1];

    mode = DC1394_VIDEO_MODE_FORMAT7_1;
    dc = dc1394_new();
    if (!dc)
	return -1;
    err = dc1394_camera_enumerate(dc, &list);
    DC1394_ERR_RTN(err, "Failed to enumerate cameras.");

    if (list->num == 0) {
	dc1394_log_error("No cameras found.");
	return -1;
    }

    camera = dc1394_camera_new(dc, list->ids[0].guid);
    if (!camera) {
	dc1394_log_error("Failed to initialize camera with guid %"PRIx64".",
			 list->ids[0].guid);
	return -1;
    }
    dc1394_camera_free_list(list);

    printf("Using camera with GUID %"PRIx64"\n", camera->guid);

    // need to use legacy firewire400 mode for now.  800 not quite reliable.
    dc1394_video_set_iso_speed(camera, DC1394_ISO_SPEED_400);

    // configure camera for format7
    err = dc1394_video_set_mode(camera, mode);
    DC1394_ERR_CLN_RTN(err, dc1394_camera_free(camera), "Can't choose video mode.");
    err = dc1394_format7_get_max_image_size(camera, mode, &max_width, &max_height);
    DC1394_ERR_CLN_RTN(err, dc1394_camera_free(camera), "cannot get max image size.");

    winleft = 0;
    wintop = 0;

    err = dc1394_format7_set_roi(camera,
				 mode,
				 DC1394_COLOR_CODING_MONO8,
				 DC1394_USE_MAX_AVAIL,
				 winleft, wintop, // left, top
				 xsize, ysize);
    DC1394_ERR_CLN_RTN(err, dc1394_camera_free(camera), "Can't set ROI.");

    // set the frame rate to absolute value in frames/sec
    err = dc1394_feature_set_mode(camera, DC1394_FEATURE_FRAME_RATE,
				  DC1394_FEATURE_MODE_MANUAL);
    DC1394_ERR_CLN_RTN(err, dc1394_camera_free (camera), "cannot set framerate to manual");
    err = dc1394_feature_set_absolute_control(camera,
					      DC1394_FEATURE_FRAME_RATE, DC1394_TRUE);
    DC1394_ERR_CLN_RTN(err, dc1394_camera_free (camera),
		       "cannot set framerate to absolute mode");
    err = dc1394_feature_set_absolute_value(camera, DC1394_FEATURE_FRAME_RATE, fps);
    DC1394_ERR_CLN_RTN(err,dc1394_camera_free (camera),"cannot set framerate");
    printf("I: framerate is %f fps\n", fps);

    // set the shutter speed to absolute value in seconds
    err = dc1394_feature_set_mode(camera, DC1394_FEATURE_SHUTTER, DC1394_FEATURE_MODE_MANUAL);
    DC1394_ERR_CLN_RTN(err, dc1394_camera_free(camera), "cannot set shutter to manual");
    err = dc1394_feature_set_absolute_control(camera, DC1394_FEATURE_SHUTTER, DC1394_TRUE);
    DC1394_ERR_CLN_RTN(err, dc1394_camera_free(camera), "cannot set shutter to absolute mode");
    err = dc1394_feature_set_absolute_value(camera, DC1394_FEATURE_SHUTTER, exptime);
    DC1394_ERR_CLN_RTN(err, dc1394_camera_free(camera), "cannot set shutter");
    printf("I: exptime is %f s\n", exptime);

    // set gain manually.  use relative value here in range 48 to 730.
    err = dc1394_feature_set_mode(camera, DC1394_FEATURE_GAIN, DC1394_FEATURE_MODE_MANUAL);
    DC1394_ERR_CLN_RTN(err, dc1394_camera_free(camera), "cannot set gain to manual");
    err = dc1394_feature_set_value(camera, DC1394_FEATURE_GAIN, 400);
    DC1394_ERR_CLN_RTN(err, dc1394_camera_free(camera), "cannot set gain");

    // set brightness manually.  use relative value in range 0 to 1023.
    err = dc1394_feature_set_mode(camera, DC1394_FEATURE_BRIGHTNESS, DC1394_FEATURE_MODE_MANUAL);
    DC1394_ERR_CLN_RTN(err, dc1394_camera_free(camera), "cannot set brightness to manual");
    err = dc1394_feature_set_value(camera, DC1394_FEATURE_BRIGHTNESS, 100);
    DC1394_ERR_CLN_RTN(err, dc1394_camera_free(camera), "cannot set brightness");

    err = dc1394_format7_get_total_bytes(camera, DC1394_VIDEO_MODE_FORMAT7_1, &total_bytes);
    DC1394_ERR_CLN_RTN(err, dc1394_camera_free(camera), "Can't get total bytes.");

    err = dc1394_capture_setup(camera, 16, DC1394_CAPTURE_FLAGS_DEFAULT);
    DC1394_ERR_CLN_RTN(err, dc1394_camera_free(camera), "Error capturing.");

    // start the camera up
    err = dc1394_video_set_transmission(camera, DC1394_ON);
    if (err != DC1394_SUCCESS) {
	dc1394_log_error("Unable to start camera iso transmission.");
	dc1394_capture_stop(camera);
	dc1394_camera_free(camera);
	return -1;
    }

    printf("Camera successfully initialized.\n");

    out = fopen("seeing.dat", "a");
    cenfile = fopen("centroids.dat", "w");
    init = fopen("init_cen_all", "r");

    i = 0;
    while (fscanf(init, "%f %f\n", &xx, &yy) != EOF) {
	box[i].x = xx;
	box[i].cenx = xx;
	box[i].y = yy;
	box[i].ceny = yy;
	box[i].r = boxsize/2.0;
	i++;
    }
    nboxes = i;
    fclose(init);

    back.r = 60;
    back.width = 10;

    /* allocate the buffers */
    if (!(dist = calloc(nimages, sizeof(double)))) {
	printf("Couldn't allocate dist array.\n");
	exit(-1);
    }
    if (!(sig = calloc(nimages, sizeof(double)))) {
	printf("Couldn't allocate sig array.\n");
	exit(-1);
    }
    if (!(weight = calloc(nimages, sizeof(double)))) {
	printf("Couldn't allocate sig array.\n");
	exit(-1);
    }
    if (!(weight_l = calloc(nimages, sizeof(double)))) {
	printf("Couldn't allocate sig array.\n");
	exit(-1);
    }
    if (!(dist_l = calloc(nimages, sizeof(double)))) {
	printf("Couldn't allocate dist_l array.\n");
	exit(-1);
    }
    if (!(sig_l = calloc(nimages, sizeof(double)))) {
	printf("Couldn't allocate sig_l array.\n");
	exit(-1);
    }
    if (!(buffer = malloc(nelements*sizeof(char)))) {
	printf("Couldn't Allocate Image Buffer\n");
	exit(-1);
    }
    if (!(buffer2 = malloc(nelements*sizeof(char)))) {
	printf("Couldn't Allocate 2nd Image Buffer\n");
	exit(-1);
    }
    if (!(average = malloc(nelements*sizeof(char)))) {
	printf("Couldn't Allocate Average Image Buffer\n");
	exit(-1);
    }

    froot = "seeing.fits";
    avemax = 0.0;

    /* get initial frame */
    grab_frame(camera, buffer2, nelements*sizeof(char));
    grab_frame(camera, buffer2, nelements*sizeof(char));
    grab_frame(camera, buffer2, nelements*sizeof(char));
    grab_frame(camera, buffer2, nelements*sizeof(char));
    // add_gaussian(buffer2, 195.0, 130.0, 140.0, 1.5);
    // add_gaussian(buffer2, 140.0, 115.0, 140.0, 1.5);

    gettimeofday(&start_time, NULL);

    for (f=0; f<nimages; f++) {
	grab_frame(camera, buffer, nelements*sizeof(char));
	// add_gaussian(buffer, 195.0, 130.0, 140.0, 1.5);
	// add_gaussian(buffer, 140.0, 115.0, 140.0, 1.5);

	// find center of star images and calculate background
	xsum = 0.0;
	ysum = 0.0;
	for (i=0; i<nboxes; i++) {
	    xsum += box[i].cenx;
	    ysum += box[i].ceny;
	}
	back.x = xsum/nboxes;
	back.y = ysum/nboxes;
	background(buffer, naxes[0], naxes[1]);

	for (i=0; i<nboxes; i++) {
	    box[i].r = boxsize/2.0;
	    centroid(buffer, naxes[0], naxes[1], i);
	}

	/* fix both boxes to the same Y */
	ysum = (box[0].y + box[1].y)/2.0;
	box[0].y = ysum;
	box[1].y = ysum;

	if (box[0].snr < box[1].snr) {
	    weight[f] = (box[0].snr/box[0].counts)*(box[0].snr/box[0].counts);
	} else {
	    weight[f] = (box[1].snr/box[1].counts)*(box[1].snr/box[1].counts);
	}

	if (weight[f] == 0.0) {
	    nbad++;
	}

	dist[f] = stardist(0, 1);
	if (dist[f] > 55.0 || dist[f] < 15.0) {
	    weight[f] = 0.0;
	}

	if (box[0].fwhm > 0.0 && box[1].fwhm > 0.0) {
	    for (i=0; i<nboxes; i++) {
		fprintf(cenfile,
			"%6.2f %6.2f %5.2f %9.4f %8.4f %7.4f %7.4f %8.4f %5.3f\t ",
			box[i].cenx,
			box[i].ceny,
			box[i].fwhm,
			box[i].counts,
			back.background,
			box[i].noise,
			box[i].sigmaxy,
			box[i].snr,
			box[i].strehl
		    );
	    }
	    sig[f] = sqrt(box[0].sigmaxy*box[0].sigmaxy + box[1].sigmaxy*box[1].sigmaxy);
	    fprintf(cenfile, "%6.2f %6.2f %.1e\n", dist[f], sig[f], exptime);
	    if (box[0].max > box[1].max) {
		avemax += box[0].max;
	    } else {
		avemax += box[1].max;
	    }
	}

	/* now average two exposures */
	for (j=0; j<nelements; j++) {
	    test = (int)buffer[j] + (int)buffer2[j];
	    if (test <= 254) {
		average[j] = buffer[j]+buffer2[j];
	    } else {
		average[j] = 254;
	    }
	}
	memcpy(buffer2, buffer, nelements*sizeof(char));

	xsum = 0.0;
	ysum = 0.0;
	for (i=0; i<nboxes; i++) {
	    xsum += box[i].cenx;
	    ysum += box[i].ceny;
	}
	back.x = xsum/nboxes;
	back.y = ysum/nboxes;
	background(average, naxes[0], naxes[1]);

	for (i=0; i<nboxes; i++) {
	    box[i].r = boxsize/2.0;
	    //centroid(average, i);
	    //box[i].r = boxsize/3.0;
	    //centroid(average, i);
	    //box[i].r = boxsize/4.0;
	    centroid(average, naxes[0], naxes[1], i);
	}

	if (box[0].snr < box[1].snr) {
	    weight_l[f] = (box[0].snr/box[0].counts)*(box[0].snr/box[0].counts);
	} else {
	    weight_l[f] = (box[1].snr/box[1].counts)*(box[1].snr/box[1].counts);
	}

	if (weight_l[f] == 0.0) {
	    nbad_l++;
	}

	dist_l[f] = stardist(0, 1);
	sig_l[f] = sqrt(box[0].sigmaxy*box[0].sigmaxy + box[1].sigmaxy*box[1].sigmaxy);

	if (f % 160 == 0) {
	    sprintf(filename, "!%s", froot);
	    fits_create_file(&fptr, "!seeing.fits", &fstatus);
	    fits_create_img(fptr, BYTE_IMG, 2, naxes, &fstatus);
	    fits_write_img(fptr, TBYTE, fpixel, nelements, buffer, &fstatus);
	    fits_close_file(fptr, &fstatus);
	    fits_report_error(stdout, fstatus);

	    status = XPASet(xpa, "timDIMM", "array [xdim=320,ydim=240,bitpix=8]", "ack=false",
			    (char *)buffer, nelements, names, messages, NXPA);
	    sprintf(xpastr, "image; box %f %f %d %d 0.0",
		    box[0].x, box[0].y, boxsize, boxsize);
	    status = XPASet(xpa, "timDIMM", "regions", "ack=false",
			    xpastr, strlen(xpastr), names, messages, NXPA);
	    sprintf(xpastr, "image; box %f %f %d %d 0.0",
		    box[1].x, box[1].y, boxsize, boxsize);
	    status = XPASet(xpa, "timDIMM", "regions", "ack=false",
			    xpastr, strlen(xpastr), names, messages, NXPA);
	}

    }

    gettimeofday(&end_time, NULL);
    printf("End capture.\n");

    /*-----------------------------------------------------------------------
     *  stop data transmission
     *-----------------------------------------------------------------------*/
    err = dc1394_video_set_transmission(camera, DC1394_OFF);
    DC1394_ERR_RTN(err,"couldn't stop the camera?");

    start_sec = start_time.tv_sec;
    start_usec = start_time.tv_usec;
    end_sec = end_time.tv_sec;
    end_usec = end_time.tv_usec;

    elapsed_time = (float)((end_sec + 1.0e-6*end_usec) - (start_sec + 1.0e-6*start_usec));
    fps = nimages/elapsed_time;
    printf("Elapsed time = %g seconds.\n", elapsed_time);
    printf("Framerate = %g fps.\n", fps);

    avemax /= nimages;
    printf("Avemax is %.3f\n", avemax);

    /* analyze short exposure */
    printf("\t SHORT EXPOSURE\n");
    mean = gsl_stats_wmean(weight, 1, dist, 1, nimages);
    avesig = gsl_stats_wmean(weight, 1, sig, 1, nimages);
    printf("mean = %f, avesig = %f\n", mean, avesig);

    printf("\n");

    var = gsl_stats_wvariance_m(weight, 1, dist, 1, nimages, mean);
    var = var - avesig*avesig;
    if (var < 0.0) {
	var = 0.0;
    }
    seeing_short = seeing(var)/pow(airmass,0.6);
    r0 = old_seeing(var);
    printf("sigma = %f, seeing = %f\n", sqrt(var), seeing_short);

    /* analyze long exposure */
    printf("\t LONG EXPOSURE\n");
    mean = gsl_stats_wmean(weight_l, 1, dist_l, 1, nimages);
    avesig = gsl_stats_wmean(weight_l, 1, sig_l, 1, nimages);
    printf("mean_l = %f, avesig_l = %f\n", mean, avesig);

    printf("\n");

    var_l = gsl_stats_wvariance_m(weight_l, 1, dist_l, 1, nimages, mean);
    var_l = var_l - avesig*avesig;
    if (var_l < 0.0) {
	var_l = 0.0;
    }
    seeing_long = seeing(var_l)/pow(airmass,0.6);
    printf("sigma_l = %f, seeing_l = %f\n", sqrt(var_l), seeing_long);

    printf("Bad samples:  %d for short, %d for long.\n", nbad, nbad_l);

    /* reduce exposure when too bright */
    if (avemax > 100.0) {
      if (exptime > 1.0e-3) {
	exptime -= 1.0e-3;
      } else {
	exptime /= 2.0;
      }
      printf("\033[0;33mStar too bright, reducing exposure time to %.1e seconds.\033[0;39m\n", exptime);
    }

    /* increase if too faint, but max out at 8 ms */
    if (avemax < 25.0) {
      if (exptime >= 1.0e-3) {
	exptime += 1.0e-3;
      } else {
	exptime *= 2.0;
      }
      if (exptime > 8.0e-3) {
	exptime = 8.0e-3;
	printf("\033[0;33mStar too faint, but exposure time max'ed out at %.1e seconds.\033[0;39m\n",
	       exptime);
      } else {
	printf("\033[0;33mStar too faint, exposure time increased to %.1e seconds.\033[0;39m\n",
	       exptime);
      }
    }

    init = fopen("exptime", "w");
    fprintf(init, "%.2e\n", exptime);
    fclose(init);

    if (nbad < 30) {
	seeing_ave = pow(seeing_short, 1.75)*pow(seeing_long,-0.75);
	printf("\033[0;32mAirmass corrected seeing = %4.2f\"\033[0;39m\n\n", seeing_short);
	printf("\033[0;32mFried Parameter, R0 = %.2f cm\033[0;39m\n\n", 100*r0);

	gmtime_r(&end_sec, &ut);
	fprintf(out, "%d-%02d-%02d %02d:%02d:%02d %f %f %f %f %f %f\n",
		ut.tm_year+1900,
		ut.tm_mon+1,
		ut.tm_mday,
		ut.tm_hour,
		ut.tm_min,
		ut.tm_sec,
		var,
		var_l,
		seeing_short,
		seeing_long,
		seeing_ave,
		airmass);

	init = fopen("init_cen_all", "w");
	for (i=0; i<nboxes; i++) {
	    fprintf(init, "%f %f\n", box[i].cenx, box[i].ceny);
	}
	fclose(init);
	init = fopen("seeing.out", "w");
	fprintf(init, "%.2f\n", seeing_short);
	fclose(init);
	init = fopen("r0.out", "w");
	fprintf(init, "%.1f\n", 100*r0);
	fclose(init);
    } else {
	printf("\n\n\033[0;31mABORTING measurement!  Lost at least one box.\033[0;39m\n\n");
    }

    fclose(cenfile);
    fclose(out);

    return (status);
}

max_stars_repo_name: marissakotze/timDIMM

max_issues_repo_name: marissakotze/timDIMM

max_forks_repo_name: marissakotze/timDIMM


text:
#pragma once

#include "Utils.h"

#include <algorithm>
#include <array>
#include <atomic>
#include <deque>
#include <functional>
#include <initializer_list>
#include <limits>
#include <memory>
#include <mutex>
#include <set>
#include <vector>

#include <gsl/gsl>

namespace indexes::utils {
template <typename epoch_t, typename ReclaimedType,
          int ReclamationThreshold = 1000, typename Enable = void>
class EpochManager;

// Provides support for epoch based reclaimation.
template <typename epoch_t, typename ReclaimedType, int ReclamationThreshold>
class EpochManager<
    epoch_t, ReclaimedType, ReclamationThreshold,
    typename std::enable_if_t<
        (ReclamationThreshold > 0) && std::is_integral_v<epoch_t> &&
        !std::is_same_v<epoch_t, bool> && !std::is_same_v<epoch_t, char>>> {
  // Calling thread must be registered by calling `RegisterThread` before using
  // epoch based safe memory reclamation.

  using ReclaimedPtrType = ReclaimedType *;

public:
  // It is undefined behaviour to use epoch manager, without registering the
  // thread. Returns false, if active registered thread count is more than
  // `MAX_THREADS`. If returned false, thread is not registered.

  // enter_epoch guarantees that all shared objects, accessed by the calling
  // thread, after enter_epoch is called are safe.
  // Recursive calls are allowed
  inline void enter_epoch() {
    auto &pd = m_local_epoch[ThreadRegistry::ThreadID()];

    if (pd.nested_level == 0) {
      pd.epoch.store(now(), std::memory_order_release);
      std::atomic_thread_fence(std::memory_order_acquire);
    }
    pd.nested_level++;
  }

  // exit_epoch marks quiescent state of the calling thread, when nesting
  // becomes 0.
  // Enables reclaimation of objects retired before calling thread's epoch.
  inline void exit_epoch() {
    auto &pd = m_local_epoch[ThreadRegistry::ThreadID()];

    pd.nested_level--;
    if (pd.nested_level == 0)
      pd.epoch.store(QUIESCENT_STATE, std::memory_order_release);
  }

  // Switch to new epoch and return old epoch
  inline epoch_t switch_epoch() { return m_global_epoch++; }

  // Returns current threads epoch
  inline epoch_t my_epoch() {
    return m_local_epoch[ThreadRegistry::ThreadID()].epoch.load(
        std::memory_order_relaxed);
  }

  // Return current epoch
  inline epoch_t now() const { return m_global_epoch; }

  // retire objects and start a new epoch.
  // Objects will be reclaimed at a suitable and safe epoch.
  // When reclaimed `reclaimer` will be called for each reclaimed object.
  void
  retire_in_new_epoch(std::function<void(ReclaimedPtrType object)> reclaimer,
                      gsl::span<ReclaimedPtrType> objects) {
    retire(reclaimer, objects, switch_epoch());
  }

  inline void
  retire_in_new_epoch(std::function<void(ReclaimedPtrType object)> reclaimer,
                      ReclaimedPtrType object) {
    retire_in_new_epoch(reclaimer, {&object, 1});
  }

  // retire objects in current (without starting new epoch) epoch.
  // Objects will be reclaimed at a suitable and safe epoch.
  // When reclaimed `reclaimer` will be called for each reclaimed object.
  void retire_in_current_epoch(
      std::function<void(ReclaimedPtrType object)> reclaimer,
      gsl::span<ReclaimedPtrType> objects) {
    retire(reclaimer, objects, now());
  }

  inline void retire_in_current_epoch(
      std::function<void(ReclaimedPtrType object)> reclaimer,
      ReclaimedPtrType object) {
    retire_in_current_epoch(reclaimer, {&object, 1});
  }

  // Reclaim objects which are safe to reclaim.
  // Returns # objects still not reclaimed (but retired).
  // A long running thread using an epoch could prevent
  // reclaimation of objects visible to that thread.
  size_t do_reclaim() {
    return reclaim_in_retire_list(m_retire_list[ThreadRegistry::ThreadID()],
                                  get_min_used_epoch());
  }

  void reclaim_all() {
    epoch_t min_used_epoch = get_min_used_epoch();

    for (auto &retire_list : m_retire_list) {
      reclaim_in_retire_list(retire_list, min_used_epoch);
    }
  }

  // Getter/Setter for reclaimation threshold.
  // # objects to collect before triggering reclaimation.
  // A low value means possibly reduced peak memory consumption, but would less
  // performance, because of frequent calls to `DoReclaim`.
  void set_reclamation_threshold(int threshold) {
    if (threshold > 0)
      m_reclaimation_threshold = threshold;
  }

  int get_reclamation_threshold(int threshold) {
    return m_reclaimation_threshold;
  }

  EpochManager()
      : m_reclaimation_threshold{ReclamationThreshold}, m_global_epoch{0},
        m_local_epoch(ThreadRegistry::MAX_THREADS),
        m_retire_list(ThreadRegistry::MAX_THREADS) {}

  EpochManager(const EpochManager &) = delete;
  EpochManager(EpochManager &&) = delete;

private:
  class Retiree {
  public:
    Retiree(ReclaimedPtrType object,
            std::function<void(ReclaimedPtrType object)> reclaimer,
            epoch_t retired_epoch)
        : m_object(object), m_retired_epoch(retired_epoch),
          m_reclaimer(reclaimer) {}

    // Can reclaim this object after the given `epoch`
    bool can_reclaim(epoch_t epoch) const { return epoch > m_retired_epoch; }

    void reclaim() { return m_reclaimer(m_object); }

  private:
    // Object to reclaim safely.
    ReclaimedPtrType m_object;
    // Epoch on which the object was retired.
    epoch_t m_retired_epoch;
    // Callback to call during reclaimation.
    std::function<void(ReclaimedPtrType object)> m_reclaimer;
  };

  static size_t reclaim_in_retire_list(std::deque<Retiree> &retire_list,
                                       epoch_t min_used_epoch) {
    auto begin = std::begin(retire_list);
    auto reclaim_upto = begin;

    for (auto &retiree : retire_list) {
      if (!retiree.can_reclaim(min_used_epoch))
        break;

      retiree.reclaim();
      reclaim_upto++;
    }

    size_t num_reclaimed = reclaim_upto - begin;

    if (num_reclaimed) {
      retire_list.erase(begin, reclaim_upto);
      retire_list.shrink_to_fit();
    }

    return retire_list.size();
  }

  epoch_t get_min_used_epoch() {
    return std::min_element(std::begin(m_local_epoch),
                            std::begin(m_local_epoch) +
                                ThreadRegistry::MaxThreadID() + 1,
                            [](const auto &a, const auto &b) {
                              return a.epoch.load() < b.epoch.load();
                            })
        ->epoch;
  }

  void retire(std::function<void(ReclaimedPtrType object)> reclaimer,
              gsl::span<ReclaimedPtrType> objects, epoch_t retired_epoch) {
    auto &retire_list = m_retire_list[ThreadRegistry::ThreadID()];

    for (auto object : objects)
      retire_list.emplace_back(object, reclaimer, retired_epoch);

    if (static_cast<int>(retire_list.size()) >= m_reclaimation_threshold)
      do_reclaim();
  }

  // Reclaimation threshold. Default 1000 objects
  std::atomic<int> m_reclaimation_threshold;

  // Global epoch
  std::atomic<epoch_t> m_global_epoch;

  // Quiescent state
  static constexpr epoch_t QUIESCENT_STATE =
      std::numeric_limits<epoch_t>::max();

  struct alignas(128) PrivateData {
    std::atomic<epoch_t> epoch = QUIESCENT_STATE;
    int nested_level = 0;
  };

  // Thread local epoch, denotes the epoch of each thread.
  // Accessed using `slot` by each thread.
  std::vector<PrivateData> m_local_epoch;

  // Thread local retire list. Accessed using `slot` by each thread.
  std::vector<std::deque<Retiree>> m_retire_list;
};
} // namespace indexes::utils

ext: h


text:
/*
 * Kaczmarz down-up sweep on diagonally banded matrix. The essential loop is:
 *
 * for each row i
 *     x = x + w*(b(i) - A(i,:)*x)*A(i,:)'
 * end
 *
 * The matrix is given in band storage format, where each row (stored
 * contiguously in memory) of the array R stores a diagonal of the matrix with
 * offset idx(j), such that
 *
 * 	   A(i,i+idx(j)) = R(i,j)
 *
 * use (from MATLAB):
 * 		y = sweepR_mex(R,idx,x,b,w,dir)
 *
 * 		R			- matrix of diagonals of matrix A
 *		idx			- offsets of diagonals
 *		x			- initial guess
 * 		b 			- right hand side (source)
 *		w 			- relaxation parameter (0 <= w <= 2)
 *		ns 		    - number of RHS
 * 		n_threads	- OPTIONAL argument to control the number of execution
 * 		threads solving CARP blocks in parallel. The number of threads can also
 * 		be defined via an environment variable (OMP_NUM_THREADS), but this
 * 		optional argument takes precedence. The default number of threads is
 * 		one. Take care if using more than one MATLAB worker per node: each
 * 		MATLAB worker will use OMP_NUM_THREADS, so if there are four workers on
 * 		a node, there will be 4 x OMP_NUM_THREADS parallel CARP sweeps.
 *  compile :
 *	   mex -largeArrayDims sweep_MT_mex.c -DDEFINEUNIX -lmwblas CFLAGS="\$CFLAGS -fopenmp" LDFLAGS="\$LDFLAGS -fopenmp"
 *
 * Author: Mathias Louboutin from Art Petrenko sweepR_mex.c
 *         Seismic Laboratory for Imaging and Modeling
 *         Department of Earth, Ocean, and Atmosperic Sciences
 *         The University of British Columbia
 *
 * Date: March, 2015

 * You may use this code only under the conditions and terms of the
 * license contained in the file LICENSE provided with this source
 * code. If you do not agree to these terms you may not use this
 * software.
*/

#include <stdlib.h> /* for getenv */
#include <stddef.h> /* for size_t type */
#include <string.h> /* for memcpy */
#include <pthread.h> /* for threading */
#include <omp.h>

/* The following section allows this file to compile on Mac OS X 10.8.5. Pass
 * the flag -DARCH_MACI64 to the compiler to activate it. */
#ifdef ARCH_MACI64
#include <mach/error.h>
typedef wchar_t char16_t;
#else /* not ARCH_MACI64 */
#include <error.h>
#endif /* ARCH_MACI64 */

#include <math.h>
#include <mex.h>
#include <matrix.h>
#include <blas.h>

struct copy_init_guess_data_t
{
	double *copy_src_real, *copy_dst_real;
	double *copy_src_imag, *copy_dst_imag;
	long n_to_copy;
};

struct sweep_data_t
{
	long start_row, end_row, ncol, ny, nx, haloWidth, main_diagonal_offset;
	double *Rr, *Ri, *yr, *yi, *br, *bi;
	long *idx;
	double w;
	int dir, n_threads, ns;
};

struct average_data_t
{
	double *copy_src_real, *copy_dst_real;
	double *copy_src_imag, *copy_dst_imag;
	double *halo_1_real, *halo_2_real;
	double *halo_1_imag, *halo_2_imag;
	double *halo_dst_real, *halo_dst_imag;
	long n_to_copy, n_in_halo;
};

struct thread_data_t
{
	struct copy_init_guess_data_t copy_init_guess_data;
	struct sweep_data_t sweep_data;
	struct average_data_t average_data;
	pthread_barrier_t *barrier;
};

void *do_sweep(void *thread_args_void)
{
    struct sweep_data_t *thread_args;
	/* Variables contained in thread_args_void struct */
	long start_row, end_row, ncol, ny, nx;
	/* Rr and Ri are pointers to short fat ncol-by-N matrices */
	double *Rr, *Ri, *yr, *yi, *br, *bi;
	long *idx;
	double w;
	int n_threads,ns;
	/* Temporary storage variables */
	double cr = 0, ci = 0;
	long offset, main_diagonal_offset;

	/* Assign local pointers to data locations in shared memory */
	thread_args = (struct sweep_data_t *) thread_args_void;
	start_row 	= thread_args->start_row;
	end_row 	= thread_args->end_row;
	ncol 		= thread_args->ncol;
	ny 			= thread_args->ny;
	nx			= thread_args->nx;
	main_diagonal_offset = thread_args->main_diagonal_offset;
	Rr 			= thread_args->Rr;
	Ri 			= thread_args->Ri;
	idx 		= thread_args->idx;
	yr 			= thread_args->yr;
	yi 			= thread_args->yi;
	br 			= thread_args->br;
	bi 			= thread_args->bi;
	w 			= thread_args->w;
	n_threads 	= thread_args->n_threads;
	ns 			= thread_args->ns;
	offset = (start_row == 0 ? 0 : - main_diagonal_offset);

	long i;
	int s;
	long j;
	long k;
	long indj;
	long toto;

	#pragma omp parallel for schedule(static,1) private(k,indj,i,j,cr,ci) num_threads(n_threads)
	for(s=0 ; s<ns; s++)
	{
		/* Kaczmarz sweep on one row block */
		for(i = start_row ; i<end_row ;i++ )
		{

			if (0 <= i + main_diagonal_offset && i + main_diagonal_offset < nx){
				cr = br[i + main_diagonal_offset+s*nx];
				ci = bi[i + main_diagonal_offset+nx*s];
			}
			else{
				error(1,0,"Discovery of whether the iterate vector is haloed failed.");
			}
			/* First loop over non-zero row elements calculates inner product
			 * of matrix row and CARP iterate */
			long diff = i - start_row + offset;
			long icol=i*ncol;

			for(j=0 ; j<ncol;j++)
			{
				/* i + idx[j] is the column index for the full Helmholtz matrix.
				 * k is the index into the vector representing the CARP iterate
				 * of the given block. */
				k = diff + idx[j];
				indj=icol + j;
				if(0<=k && k<ny)
				{
					cr -= Rr[indj]*yr[k+s*nx] - Ri[indj]*yi[k+s*nx];
					ci -= Rr[indj]*yi[k+s*nx] + Ri[indj]*yr[k+s*nx];

				}
			}

			cr*=w;
			ci*=w;
			/* Second loop over non-zero row elements updates Karkzmarz iterate */

			for(j=0 ; j<ncol;j++)
			{
				k = diff + idx[j];
				indj=icol + j;
				if(0<=k && k<ny)
				{
					yr[k+s*nx] +=   cr*Rr[indj] + ci*Ri[indj];
					yi[k+s*nx] +=  -cr*Ri[indj] + ci*Rr[indj];

				}
			}

		}
		cr = 0;
		ci = 0;
		for(i = end_row-1 ; i>start_row-1 ;i-- )
		{

			if (0 <= i + main_diagonal_offset && i + main_diagonal_offset < nx){
				cr = br[i + main_diagonal_offset+s*nx];
				ci = bi[i + main_diagonal_offset+nx*s];
			}
			else{
				error(1,0,"Discovery of whether the iterate vector is haloed failed.");
			}
			/* First loop over non-zero row elements calculates inner product
			 * of matrix row and CARP iterate */
			long diff = i - start_row + offset;
			long icol=i*ncol;

			for(j=0 ; j<ncol;j++)
			{
				/* i + idx[j] is the column index for the full Helmholtz matrix.
				 * k is the index into the vector representing the CARP iterate
				 * of the given block. */
				k = diff + idx[j];
				indj=icol + j;
				if(0<=k && k<ny)
				{
					cr -= Rr[indj]*yr[k+s*nx] - Ri[indj]*yi[k+s*nx];
					ci -= Rr[indj]*yi[k+s*nx] + Ri[indj]*yr[k+s*nx];

				}
			}

			cr*=w;
			ci*=w;
			/* Second loop over non-zero row elements updates Karkzmarz iterate */

			for(j=0 ; j<ncol;j++)
			{
				k = diff + idx[j];
				indj=icol + j;
				if(0<=k && k<ny)
				{
					yr[k+s*nx] +=   cr*Rr[indj] + ci*Ri[indj];
					yi[k+s*nx] +=  -cr*Ri[indj] + ci*Rr[indj];

				}
			}

		}
	}

	return NULL;
}

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
	/* structs to hold all arguments to each thread in one variable */
	struct sweep_data_t **thread_args_sweep = NULL;
	struct copy_init_guess_data_t **copy_init_guess_data = NULL;

	mwSize ncol, nx;
	ptrdiff_t ncolBlas, idxIncBlas = 1, maxIdxLoc = 0;

	double *Rr,*Ri,*idxd = NULL,*xr,*xi,*br,*bi,*yr,*yi;
	long *idx = NULL;
    double w = 0;
	int ns;
	mwSize n_threads = 1;
	char *n_threads_str = NULL;
	mwSize N=1, numGridPointsPerBlock, main_diagonal_offset;

	/* Flags that are set if memory is allocated within the MEX file */
	int Ri_alloc=0, xi_alloc=0, bi_alloc=0;
	mwSize *seg_bounds_hi, *seg_bounds_mid, *seg_bounds_row, *seg_bounds_lo;

	/* Read input arguments; initialize complex part to zero if input is real. */
	ns = lrint(mxGetScalar(prhs[5]));
	N = mxGetN(prhs[0]);
    ncol = mxGetM(prhs[0]);
	ncolBlas = (ptrdiff_t)ncol;
	Rr  = mxGetPr(prhs[0]);
    if(mxIsComplex(prhs[0])){
        Ri  = mxGetPi(prhs[0]);
    }
    else{
        Ri = mxCalloc(N*ncol,sizeof(double));
		Ri_alloc = 1;
    }
	idxd = mxGetPr(prhs[1]);
    nx  = mxGetM(prhs[2]);
	xr  = mxGetPr(prhs[2]);
    if(mxIsComplex(prhs[2])){
        xi  = mxGetPi(prhs[2]);
    }
    else{
        xi = mxCalloc(nx*ns,sizeof(double));
		xi_alloc = 1;
    }
	br  = mxGetPr(prhs[3]);
    if(mxIsComplex(prhs[3])){
        bi  = mxGetPi(prhs[3]);
    }
    else{
        bi = mxCalloc(nx*ns,sizeof(double));
		bi_alloc = 1;
    }
    if (mxGetM(prhs[3]) != nx){
    mexErrMsgIdAndTxt("SLIM_release_apps:tools:algorithms:ThreeDFreqModeling:sweepR_mex:NumElements",
				"The number of elements in the iterate and right hand side vectors must be equal.");
	}

	/* Allocate the final output vector */
	plhs[0]  = mxCreateDoubleMatrix(nx,ns, mxCOMPLEX);
	yr       = mxGetPr(plhs[0]);
    yi       = mxGetPi(plhs[0]);

	/* Check to make sure memory was allocated correctly */
	if (Rr==NULL || Ri==NULL || idxd==NULL || xr==NULL || xi==NULL || br==NULL || bi==NULL ||
	    yr==NULL || yi==NULL){
		mexErrMsgIdAndTxt("SLIM_release_apps:tools:algorithms:ThreeDFreqModeling:sweepR_mex:OutOfMemory",
				"Could not allocate memory for main computational variables.");
	}

	if ((idx = (long *) mxCalloc(ncol, sizeof(long))) == NULL){
		mexErrMsgIdAndTxt("SLIM_release_apps:tools:algorithms:ThreeDFreqModeling:sweepR_mex:OutOfMemory",
				"Could not allocate memory for main computational variables.");
	}
	mwSize i;
	for (i=0; i < ncol; i++){
		idx[i] = lrint(idxd[i]);
	}

	/* The default value for the number of threads can be overridden by an
	 * environment variable. */
	n_threads_str = getenv("OMP_NUM_THREADS");
	if (n_threads_str == NULL){
		n_threads = 1;
	}
	else{
		n_threads = strtol(n_threads_str, NULL, 10);
		if(n_threads < 1){
			n_threads = 1;
		}
	}
	/* The environment variable can in turn be overridden by an optional
	 * argument to the mexFunction. */
	if (nrhs >= 7){
		if(1 <= lrint(mxGetScalar(prhs[6]))){
			n_threads = lrint(mxGetScalar(prhs[6]));
		}
	}

	/* printf("Using %d threads \n",n_threads);
	/* Partition the iterate vector into blocks. Note that the below
	 * partitioning scheme is slighlty different from that in pCARPCG.m in this
	 * directory. The partitioning scheme of pCARPCG corresponds to
	 * distributing a three dimensional array with dimensions given by n
	 * according to Matlab's codistributor1d.defaultPartition(n(3)), and then
	 * vectorizing it. The partition scheme of the present file instead uses
	 * Matlab's codistributor1d.defaultPartion(prod(n)). In other words,
	 * pCARPCG divides the iterate into blocks along the slow dimension,
	 * whereas sweepR_mex.c does not take dimensionality into account, only the
	 * total number of gridpoints. This is done to avoid needing an extra input
	 * parameter with the the system dimensions. The seg_bounds_hi, _lo and
	 * _mid arrays contain indices into non-haloed vectors, while the
	 * seg_bounds_row array contains indices to the rows of the system matrix.
	 *
		 * yr_seg[i_thread-1]        overlap      yr_seg[i_thread]
	 * ------------------------|-----|-----|-------------------------------
	 *         .----------------^     |     ^-------------------.
	 *    seg_bounds_lo[i_thread], seg_bounds_mid[i_thread], seg_bounds_hi[i_thread]
	 */
	numGridPointsPerBlock = N;
	seg_bounds_hi  = (mwSize *)mxCalloc(2,sizeof(mwSize));
	seg_bounds_mid = (mwSize *)mxCalloc(2,sizeof(mwSize));
	seg_bounds_lo  = (mwSize *)mxCalloc(2,sizeof(mwSize));
	seg_bounds_row = (mwSize *)mxCalloc(2,sizeof(mwSize));
	if (N == nx){
		main_diagonal_offset = 0;
	}
	else{
		/* The vector is haloed. We are only able to correctly process matrices
		 * with a non-zero main diagonal and symmetric off-main diagonal offsets. */
		if (ncol % 2 != 1){
			mexErrMsgIdAndTxt("SLIM_release_apps:tools:algorithms:ThreeDFreqModeling:sweepR_mex:EvenNumberOfDiags",
					"Input iterate vector appears to be haloed but there is an even number of non-zero diagonals in the system matrix.");
		}
		main_diagonal_offset = idx[ncol/2];
		mwSize i;
		for (i = 1; i <= ncol/2; i++){
			if (idx[ncol/2 + i] - main_diagonal_offset != -(idx[ncol/2 - i] - main_diagonal_offset)){
				mexErrMsgIdAndTxt("SLIM_release_apps:tools:algorithms:ThreeDFreqModeling:sweepR_mex:DiagsNotSymmetric",
						"Input iterate vector appears to be haloed but the pattern of non-zero diagonals in the system matrix is not symmetric.");
			}
		}
	}

	seg_bounds_hi[0]  = 0;
	seg_bounds_mid[0] = 0;
	seg_bounds_row[0] = 0;
	seg_bounds_lo[0]  = 0;
	seg_bounds_lo[1]  = nx;
	seg_bounds_hi[1]  = nx;
	seg_bounds_mid[1] = nx;
	seg_bounds_row[1] = N;

	thread_args_sweep = (struct sweep_data_t **)mxCalloc(1, sizeof(struct sweep_data_t *));

	/* Set thread arguments */
	thread_args_sweep[0] = (struct sweep_data_t *)mxCalloc(1,sizeof(struct sweep_data_t));
	thread_args_sweep[0]->start_row 	= seg_bounds_row[0];
	thread_args_sweep[0]->end_row   	= seg_bounds_row[1];
	thread_args_sweep[0]->ncol 	 	= ncol;
	thread_args_sweep[0]->ny 	 	= seg_bounds_hi[1]-seg_bounds_lo[0];
	thread_args_sweep[0]->nx			= nx;
	thread_args_sweep[0]->main_diagonal_offset = main_diagonal_offset;
	thread_args_sweep[0]->Rr  		= Rr;
	thread_args_sweep[0]->Ri  		= Ri;
	thread_args_sweep[0]->idx 		= idx;
	thread_args_sweep[0]->yr	 	= yr;
	thread_args_sweep[0]->yi	 	= yi;
	thread_args_sweep[0]->br 		= br;
	thread_args_sweep[0]->bi 		= bi;
	thread_args_sweep[0]->w 		 = w;
	thread_args_sweep[0]->n_threads = n_threads;
	thread_args_sweep[0]->ns = ns;
	/* Set the initial guess directly in the output array too */
	memcpy((void *)yr, (void *)xr, sizeof(double)*nx*ns);
	memcpy((void *)yi, (void *)xi, sizeof(double)*nx*ns);
	do_sweep((void *)thread_args_sweep[0]);

	/* Free memory if it was allocated within the MEX file. */
	if (Ri_alloc){
		mxFree(Ri);
	}
	if (xi_alloc){
		mxFree(xi);
	}
	if (bi_alloc){
		mxFree(bi);
	}

	mxFree(idx);
	mxFree(thread_args_sweep);

	/* Don't think I need pthread_exit() here, because pthread_join is called above */
	return;
}

max_stars_repo_name: liaman/SLIM-release-apps-public


text:
#ifndef GBPINTERPOLATE_AWAKE
#define GBPINTERPOLATE_AWAKE
#include <gbpLib.h>
#include <gsl/gsl_interp.h>
#include <gsl/gsl_integration.h>

typedef struct interp_struct interp_info;
struct interp_struct {
    gsl_interp *           interp;
    gsl_interp_accel *     accel;
    size_t                 n;
    const gsl_interp_type *T;
    double *               x;
    double *               y;
};

void   free_interpolate(void **interp, void *params);
void   init_interpolate(double *x, double *y, size_t n, const gsl_interp_type *T, interp_info **interp);
void   ADaPS_store_interp(ADaPS **list, void *data, const char *name, ...);
double interpolate(interp_info *interp, double x);
double interpolate_derivative(interp_info *interp, double x);
double interpolate_integral(interp_info *interp, double x_lo, double x_hi);
double interpolate_maximum_function(double x, void *params);
void   interpolate_maximum(interp_info *interp,
                           double       x_lo_in,
                           double       x_guess_in,
                           double       x_hi_in,
                           double       threshold,
                           double *     x_maxima,
                           double *     y_maxima);
double interpolate_minimum_function(double x, void *params);
void   interpolate_minimum(interp_info *interp,
                           double       x_lo_in,
                           double       x_guess_in,
                           double       x_hi_in,
                           double       threshold,
                           double *     x_minima,
                           double *     y_minima);

#endif

max_stars_repo_name: gbpoole/gbpCode

max_issues_repo_name: gbpoole/gbpCode

max_forks_repo_name: gbpoole/gbpCode


text:
/**
 * @file      LoadFromBST.h
 *
 * @brief     Definition of function to load voxel volume from BST files
 *
 * @author    Stefan Reinhold
 * @copyright Copyright (C) 2018 Stefan Reinhold  -- All Rights Reserved.
 *            You may use, distribute and modify this code under the terms of
 *            the AFL 3.0 license; see LICENSE for full license details.
 */

#pragma once

#include "VolumeSize.h"
#include "VoxelSize.h"

#include <Eigen/Core>
#include <ImageStack/ImageStack.h>
#include <ImageStack/ImageStackLoaderBST.h>
#include <ImageStack/ResolutionDecorator.h>
#include <gsl/gsl>

namespace CortidQCT {

namespace IO {

/**
 * @brief Loads a BST file and calles the given functional with the loaded
 * data.
 *
 * The functional is called with a pointer to the raw data as the first
 * argument, the 3d size of the loaded data as the second argument a with the
 * 3d resolution (i.e. voxel size) as the third argument.
 *
 * @tparam F functional to be called by the function
 * @param filename Path to the file to be loaded
 * @param f functional that is beeing called
 * @throw std::runtime_error if the loading failed
 */
template <class F> void loadFromBST(std::string const &filename, F &&f) {
  using Decorator = ImageStack::ResolutionDecorator;
  using Img = ImageStack::ImageStack<float, ImageStack::HostStorage, Decorator>;
  using Loader = ImageStack::ImageStackLoaderBST<Img>;

  auto const img = Img(Loader(filename));

  auto const mapping = img.map();
  auto const size = img.size().template cast<std::size_t>().eval();
  auto const resol = img.resolution.template cast<float>().eval();
  VoxelSize voxelSize{resol.x(), resol.y(), resol.z()};
  VolumeSize volumeSize{size.x(), size.y(), size.z()};

  f(mapping.data(), volumeSize, voxelSize);
}

} // namespace IO

} // namespace CortidQCT

max_stars_repo_name: ithron/CortidQCT

max_issues_repo_name: ithron/CortidQCT

max_forks_repo_name: ithron/CortidQCT