text:
/***************************************************************************
                          misc.h  -  description
                             -------------------
    copyright            : (C) 2005 by MOUSSA
    email                : mmousa@liris.cnrs.fr
 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   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.                                   *
 *                                                                         *
 ***************************************************************************/

#ifndef MISC_H
#define MISC_H

#include "config.h"

#include <gsl/gsl_math.h>
#include <gsl/gsl_monte.h>
#include <gsl/gsl_monte_plain.h>
#include <gsl/gsl_sf.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>

//////////////////////////////////////////////////////////////////////////////
template<class T> void Add_vectors(T *vec1, T *vec2, T *result, int size);
template<class T> void Add_vectors(T *vec1, T *vec2, T *result, int size){
  for(int i=0; i<size; i++){
    result[i] = vec1[i]+vec2[i];
  }
}
//////////////////////////////////////////////////////////////////////////////
template<class T> void Mult_vector_const(T *vec, T c, T *result, int size);
template<class T> void Mult_vector_const(T *vec, T c, T *result, int size){
  for(int i=0; i<size; i++){
    result[i] = vec[i]*c;
  }
}
//////////////////////////////////////////////////////////////////////////////
template<class T>void assign_value(T *vec, T c, int size);
template<class T>void assign_value(T *vec, T c, int size){
  for(int i = 0; i<size; i++) vec[i] = c;
}
//////////////////////////////////////////////////////////////////////////////
template<class T>void copy_vectors(T *src, T *dest, int size);
template<class T>void copy_vectors(T *src, T *dest, int size){
  for(int i = 0; i<size; i++){
    dest[i] = src[i];
  }
}
//////////////////////////////////////////////////////////////////////////////
template<class T> int is_overlaping(T *vec1, T *vec2, int size);
template<class T> int is_overlaping(T *vec1, T *vec2, int size){
  for(int i = 0; i<size; i++){
    if((fabs(vec1[i])>1e-6)&&(fabs(vec2[i])>1e-6)) return true;
  }
  return false;
}
//////////////////////////////////////////////////////////////////////////////
template<class T> void print_vector(T *vec, int size);
template<class T> void print_vector(T *vec, int size){
  for(int i = 0; i<size; i++){
    if(vec[i]>0.0) std::cout<<vec[i]<<std::endl;
  }
}
//////////////////////////////////////////////////////////////////////////////
template<class T> bool AllOnes(T *vec, int size);
template<class T> bool AllOnes(T *vec, int size){
  for(int i = 0; i < size; i++)
    if(fabs(vec[i]) < 1e-6) return false;
  return true;
}
//////////////////////////////////////////////////////////////////////////////
int gettoken(char[30], FILE *);
bool FileName(char[30], char *);
void theta_phi_bounds(Triangle_3 tr, double *mintheta, double *maxtheta, double *minphi, double *maxphi);
void SphericalHarmonic(int, int, double, double, gsl_complex *);
#endif

max_stars_repo_name: kyeonghopark/jgt-code

max_issues_repo_name: kyeonghopark/jgt-code

max_forks_repo_name: kyeonghopark/jgt-code


text:
#include <gsl/gsl_multiset.h>
#include <gsl/gsl_errno.h>
#include <stdio.h>

int mgsl_multiset_fwrite(const char *filename, const gsl_multiset *p)
{
  FILE *fp;
  if((fp = fopen(filename, "w")) == NULL) return GSL_EFAILED;
  if(gsl_multiset_fwrite(fp, p) != GSL_SUCCESS) return GSL_EFAILED;
  fclose(fp);
  return GSL_SUCCESS;
}

int mgsl_multiset_fread(const char *filename, gsl_multiset *p)
{
  FILE *fp;
  if((fp = fopen(filename, "r")) == NULL) return GSL_EFAILED;
  if(gsl_multiset_fread(fp, p) != GSL_SUCCESS) return GSL_EFAILED;
  fclose(fp);
  return GSL_SUCCESS;
}

int mgsl_multiset_fprintf(const char *filename, const gsl_multiset *p, const char *format)
{
  FILE *fp;
  if((fp = fopen(filename, "w")) == NULL) return GSL_EFAILED;
  if(gsl_multiset_fprintf(fp, p, format) != GSL_SUCCESS) return GSL_EFAILED;
  fclose(fp);
  return GSL_SUCCESS;
}

int mgsl_multiset_fscanf(const char *filename, gsl_multiset *p)
{
  FILE *fp;
  if((fp = fopen(filename, "r")) == NULL) return GSL_EFAILED;
  if(gsl_multiset_fscanf(fp, p) != GSL_SUCCESS) return GSL_EFAILED;
  fclose(fp);
  return GSL_SUCCESS;
}

ext: c

max_stars_repo_name: frithnanth/raku-Math-Libgsl-Multiset

max_issues_repo_name: frithnanth/raku-Math-Libgsl-Multiset

max_forks_repo_name: frithnanth/raku-Math-Libgsl-Multiset


text:
//#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_roots.h>

#define NSPLINELD	1001
#define NSPLINEEZ	1001
#define NSPLINERR	10001
#define NSPLINEDF	10001
gsl_interp_accel *accLD, *accRed, *accEz, *accDF;
gsl_spline *splineLD, *splineRed, *splineEz, *splineDF;

#define MPC_IN_M		3.08567758e22
#define GPC_IN_M		(1e3 * MPC_IN_M)
#define C_LIGHT			2.99792458e8
#define HUBBLE0			(1e5 * 0.71 / MPC_IN_M)
#define DH				(C_LIGHT / HUBBLE0 / GPC_IN_M)
#define DH3				(DH * DH * DH)

void load_splines()
{
	long int i;

	// Luminosity distribution

	accLD = gsl_interp_accel_alloc();
	splineLD = gsl_spline_alloc(gsl_interp_cspline, NSPLINELD);

	double z1[NSPLINELD], ld[NSPLINELD];
	FILE *fp1 = fopen("support_data/splines_LuminosityDistance.txt","r");
	for ( i=0; i<NSPLINELD; i++ )
	{
		fscanf(fp1, "%lf\t%lf\n", &z1[i], &ld[i]);
	}
	fclose(fp1);

	gsl_spline_init(splineLD, z1, ld, NSPLINELD);

	// Redshift distribution

	accRed = gsl_interp_accel_alloc();
	splineRed = gsl_spline_alloc(gsl_interp_cspline, NSPLINERR);

	double z2[NSPLINERR], red[NSPLINERR];
	FILE *fp2 = fopen("support_data/splines_RedshiftRate.txt","r");
	for ( i=0; i<NSPLINERR; i++ )
	{
		fscanf(fp2, "%lf\t%lf\n", &z2[i], &red[i]);
	}
	fclose(fp2);

	gsl_spline_init(splineRed, z2, red, NSPLINERR);

	// E(z) cosmology

	accEz = gsl_interp_accel_alloc();
	splineEz = gsl_spline_alloc(gsl_interp_cspline, NSPLINEEZ);

	double z3[NSPLINEEZ], Ez[NSPLINEEZ], temp1, temp2;
	FILE *fp3 = fopen("support_data/splines_Ez.txt","r");
	for ( i=0; i<NSPLINEEZ; i++ )
	{
		fscanf(fp3, "%lf %lf %lf %lf\n", &z3[i], &temp1, &temp2, &Ez[i]);
	}
	fclose(fp3);

	gsl_spline_init(splineEz, z3, Ez, NSPLINEEZ);

	// detection fraction (w.r.t z)

	accDF = gsl_interp_accel_alloc();
	splineDF = gsl_spline_alloc(gsl_interp_cspline, NSPLINEDF);

	double z4[NSPLINEDF], df[NSPLINEDF], temp3;
	FILE *fp4;
	if (runargs.method == NEURALNET)
	{
		fp4 = fopen("support_data/splines_detection_fraction_z.txt","r");
		for ( i=0; i<NSPLINEDF; i++ )
		{
			fscanf(fp4, "%lf %lf %lf\n", &z4[i], &df[i], &temp3);
		}
		fclose(fp4);
	}
	else if (runargs.method == RANDOMFOREST)
	{
		fp4 = fopen("support_data/splines_detection_fraction_z_RF.txt","r");
		for ( i=0; i<NSPLINEDF; i++ )
		{
			fscanf(fp4, "%lf %lf\n", &z4[i], &df[i]);
		}
		fclose(fp4);
	}
	else if (runargs.method == ADABOOST)
	{
		fp4 = fopen("support_data/splines_detection_fraction_z_AB.txt","r");
		for ( i=0; i<NSPLINEDF; i++ )
		{
			fscanf(fp4, "%lf %lf\n", &z4[i], &df[i]);
		}
		fclose(fp4);
	}
	else if (runargs.method == FLUXTHRESH)
	{
		fp4 = fopen("support_data/splines_detection_fraction_z_flux.txt","r");
		for ( i=0; i<NSPLINEDF; i++ )
		{
			fscanf(fp4, "%lf %lf\n", &z4[i], &df[i]);
		}
		fclose(fp4);
	}

	gsl_spline_init(splineDF, z4, df, NSPLINEDF);
}

void unload_splines()
{
	gsl_spline_free(splineLD);
	gsl_interp_accel_free(accLD);

	gsl_spline_free(splineRed);
	gsl_interp_accel_free(accRed);

	gsl_spline_free(splineEz);
	gsl_interp_accel_free(accEz);

	gsl_spline_free(splineDF);
	gsl_interp_accel_free(accDF);
}

static inline double MAX(double v1, double v2)
{
	return (((v1) > (v2)) ? (v1) : (v2));
}

double log_sum_exp(double x, double y)
{
	//convert to natural log
	x /= M_LOG10E;
	y /= M_LOG10E;

	double mx = MAX(x,y);
	double val = log(exp(x - mx) + exp(y - mx)) + mx;

	return val/M_LN10;
}

double log_subtract_exp(double x, double y)
{
	// x >= y must be true!
	if (y>x)
	{
		fprintf(stderr, "%lf < %lf !\n", x, y);
		exit(-1);
	}
	//convert to natural log
	x /= M_LOG10E;
	y /= M_LOG10E;
	double val = x + log1p(-exp(y-x));
	return val/M_LN10;
}

double lumdist_logpowint(double log_lum_min, double log_lum_max, double pwr)
{
	double val;
	if ( pwr > -1.0 )
	{
		val = -pwr * log_lum_star_global - log10(pwr+1.0) + log_subtract_exp((pwr+1.0)*log_lum_max, (pwr+1.0)*log_lum_min);
	}
	else
	{
		val = -pwr * log_lum_star_global - log10(-(pwr+1.0)) + log_subtract_exp((pwr+1.0)*log_lum_min, (pwr+1.0)*log_lum_max);
	}
	//fprintf(stderr, "%lf -> %lf : %lf\n", log_lum_min, log_lum_max, val);
	return val;
}

double lumdist_logpowinv(double prob, double log_lum_min, double log_lum_max, double pwr)
{
	double temp = lumdist_logpowint(log_lum_min, log_lum_max, pwr) + log10(prob);
	temp += pwr * log_lum_star_global;

	if ( pwr > -1.0 )
	{
		temp += log10(pwr+1.0);
		temp = log_sum_exp((pwr+1.0)*log_lum_min, temp);
	}
	else
	{
		temp += log10(-(pwr+1.0));
		temp = log_subtract_exp((pwr+1.0)*log_lum_min, temp);
	}

	return temp/(pwr+1.0);
}

double luminosity_distr_fast(double prob, double log_lum_min, double log_lum_max)
{
	//fprintf(stderr, "%lf, %lf, %lf, %lf, %lf, %lf\n", prob, log_lum_min, log_lum_max, log_lum_star_global, alpha1_global, beta1_global);

	double log_lum_int_lower = lumdist_logpowint(log_lum_min, log_lum_star_global, alpha1_global);
	double log_lum_int_upper = lumdist_logpowint(log_lum_star_global, log_lum_max, beta1_global);
	double log_lum_int_total = log_sum_exp(log_lum_int_lower, log_lum_int_upper);
	double lower_frac = pow(10.0, log_lum_int_lower - log_lum_int_total), upper_frac = 1.0 - lower_frac;

	//fprintf(stderr, "%lf, %lf, %lf, %lf, %lf\n", log_lum_int_lower, log_lum_int_upper, log_lum_int_total, lower_frac, upper_frac);

	if ( prob <= lower_frac )
	{
		return lumdist_logpowinv(prob/lower_frac, log_lum_min, log_lum_star_global, alpha1_global);
	}
	else
	{
		return lumdist_logpowinv((prob - lower_frac)/upper_frac, log_lum_star_global, log_lum_max, beta1_global);
	}
}

double lum2flux_integrand(double E, void *params)
{
	double *vals = (double *)params;
	double alpha = vals[0];
	double beta = vals[1];
	double E0 = vals[2];
	double Ebreak = vals[3];

	if ( E < Ebreak )
	{
		return pow(E/100.0,alpha) * exp(-E/E0) * E;
	}
	else
	{
		return pow((alpha-beta)*E0/100.0,alpha-beta) * exp(beta-alpha) * pow(E/100.0,beta) * E;
	}
}

double lum2flux_integral_numeric(double alpha, double beta, double Epeak, double Emin, double Emax)
{
	double E0 = Epeak / (2.0+alpha);
	double Ebreak = (alpha-beta)*E0;

	//gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000);

	double result, error;
	unsigned long int neval;
	double vals[4] = {alpha, beta, E0, Ebreak};
	gsl_function F;
  	F.function = &lum2flux_integrand;
  	F.params = &vals[0];

  	gsl_set_error_handler_off();
  	//gsl_integration_qag(&F, Emin, Emax, 1e-7, 1e-6, 1000, 1, w, &result, &error);
  	gsl_integration_qng(&F, Emin, Emax, 1e-7, 1e-6, &result, &error, &neval);
  	//fprintf(stderr, "%d evaluations\t", neval);

  	//gsl_integration_workspace_free (w);

  	return result;
}

double Redshift_distribution_unnormalized(double z, double n1, double n2, double z1)
{
	if (z <= z1)
	{
		return pow(1.0 + z, n1);
	} else {
		return pow(1.0 + z1, n1 - n2) * pow(1.0 + z, n2);
	}
}

double Redshift_distribution_normalized(double z, double n0, double n1, double n2, double z1)
{
	return n0 * Redshift_distribution_unnormalized(z, n1, n2, z1);
}

double Redshift_rescaled(double z, void *params)
{
	double *pars = (double *) params;
	double n0 = pars[0];
	double n1 = pars[1];
	double n2 = pars[2];
	double z1 = pars[3];

	double Rprime = Redshift_distribution_normalized(z, n0, n1, n2, z1) / (1.0 + z);

	Rprime *= DH3 * gsl_spline_eval(splineEz, z, accEz);

	return Rprime;
}

double Redshift_rejection_sampler(long int *seed, double n0, double n1, double n2, double z1)
{
	double pars[4] = {n0, n1, n2, z1};

	double Rzmax = Redshift_rescaled(z1, (void *) &pars[0]);

	double z = 0.0, p, x;
	for (;;)
	{
		z = ran2d(seed) * 10.0;

		p = Redshift_rescaled(z, (void *) &pars[0]) / Rzmax;

		x = ran2d(seed);

		if (x < p) break;
	}

	return z;
}

double GRBNumberIntegral(double n0, double n1, double n2, double z1)
{
	double pars[4] = {n0, n1, n2, z1};

	double result, error;
	unsigned long int neval;
	gsl_function F;
  	F.function = &Redshift_rescaled;
  	F.params = (void *) &pars[0];

  	gsl_set_error_handler_off();
  	gsl_integration_qng(&F, ZMIN, ZMAX, 1e-7, 1e-6, &result, &error, &neval);

  	result *= 4.0 * M_PI;

  	return result;
}

double GRBRate(double z, double n0, double n1, double n2, double z1)
{
	double Rprime = Redshift_distribution_normalized(z, n0, n1, n2, z1) / (1.0 + z);

	Rprime *= DH3 * gsl_spline_eval(splineEz, z, accEz);

	Rprime *= 4.0 * M_PI * runargs.tobs * gsl_spline_eval(splineDF, z, accDF) / 6.0;

	return Rprime;
}

double GRBRateFunc(double z, void *params)
{
	double *pars = (double *) params;
	double n0 = pars[0];
	double n1 = pars[1];
	double n2 = pars[2];
	double z1 = pars[3];

	return GRBRate(z, n0, n1, n2, z1);
}

double GRBRateIntegral(double n0, double n1, double n2, double z1)
{
	double pars[4] = {n0, n1, n2, z1};

	double result, error;
	unsigned long int neval;
	gsl_function F;
  	F.function = &GRBRateFunc;
  	F.params = (void *) &pars[0];

  	gsl_set_error_handler_off();
  	gsl_integration_qng(&F, ZMIN, ZMAX, 1e-7, 1e-6, &result, &error, &neval);

  	return result;
}


text:
/****************************************************************************
** Copyright (c) 2021, Fougue Ltd. <http://www.fougue.pro>
** All rights reserved.
** See license at https://github.com/fougue/mayo/blob/master/LICENSE.txt
****************************************************************************/

#pragma once

#include <gsl/span>

namespace Mayo {

template<typename T> using Span = gsl::span<T>;

} // namespace Mayo

ext: h

max_stars_repo_name: Unionfab/mayo

max_issues_repo_name: Unionfab/mayo

max_forks_repo_name: Unionfab/mayo


text:
/*
 * C version of Diffusive Nested Sampling (DNest4) by Brendon J. Brewer
 *
 * Yan-Rong Li, liyanrong@mail.ihep.ac.cn
 * Jun 30, 2016
 *
 */

#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <math.h>
#include <time.h>
#include <string.h>
#include <gsl/gsl_rng.h>

#include "dnestvars.h"

/*! \file dnestpostprocess.c
 *  \brief post process the sample generated by dnest.
 */

int cmp_sample(const void *pa, const void *pb);
typedef struct
{
  double logl, tiebreaker;
  int id;
}SampleType;

/*
 * This function calculates log(exp(x1) - exp(x2)).
 */
double logdiffexp(double x1, double x2)  // x1 is larger
{
  double biggest = x1;
  double xx1 = x1 - biggest, xx2 = x2 - biggest;
  return log(exp(xx1) - exp(xx2)) + biggest;
}
/*
 * This functions calculates log(exp(x1)+...+exp(xn)).
 */
double logsumexp(double *x, int n)
{
  int j;
  double sum, max;

  max = x[0];
  for(j = 0; j < n; j++)
  {
    max = fmax(max, x[j]);
  }
  sum = 0.0;
  for(j=0; j< n; j++)
  {
    sum += exp( x[j] - max);
  }
  return log(sum) + max;
}

void postprocess(double temperature)
{
  printf("# Starts postprocess.\n");
  FILE *fp, *fp_sample;

  double **levels_orig, **sample_info, *logl;
  int *sandwhich;
  double *psample;
  int i, j;
  int num_levels, num_samples;
  char buf[BUF_MAX_LENGTH];
  int moreSample = 1;

  // read number of levels and samples
  fp = fopen(options.sampler_state_file, "r");
  if(fp == NULL)
  {
    fprintf(stderr, "# Error: Cannot open file %s.\n", options.sampler_state_file);
    exit(0);
  }
  fscanf(fp, "%d %d\n", &num_levels, &num_samples);
  fclose(fp);

  // allocate memory for levels
  levels_orig = malloc(num_levels * sizeof(double *));
  for(i=0; i< num_levels; i++)
  {
    levels_orig[i] = malloc(3 * sizeof(double));
  }

  // allocate memory for sample_info
  sample_info = malloc(num_samples * sizeof(double *));
  for(i=0; i< num_samples; i++)
  {
    sample_info[i] = malloc(3 * sizeof(double));
  }

  // allocate memory for samples
  logl = (void *)malloc(num_samples * sizeof(double));
  sandwhich = malloc(num_samples * sizeof(int));
  psample = (double *)malloc(dnest_size_of_modeltype);

  // read levels
  fp = fopen(options.levels_file, "r");
  if(fp == NULL)
  {
    fprintf(stderr, "# Error: Cannot open file %s.\n", options.levels_file);
    exit(0);
  }
  fgets(buf, BUF_MAX_LENGTH, fp);
  for(i=0; i < num_levels; i++)
  {
    if(feof(fp) != 0)
    {
      fprintf(stderr, "# Error: file %s ends at %d.\n", options.levels_file, i);
      exit(0);
    }

    fgets(buf, BUF_MAX_LENGTH, fp);

    if(sscanf(buf, "%lf %lf %lf", &levels_orig[i][0], &levels_orig[i][1], &levels_orig[i][2]) < 3)
    {
      fprintf(stderr, "# Error: Cannot read file %s.\n", options.levels_file);
      exit(0);
    }
    buf[0]='\0';  // clear up buf
  }
  fclose(fp);

  // read sample_info
  if(dnest_flag_sample_info == 0) //no need to recalculate
  {
    fp = fopen(options.sample_info_file, "r");
    if(fp == NULL)
    {
      fprintf(stderr, "# Error: Cannot open file %s.\n", options.sample_info_file);
      exit(0);
    }
    fgets(buf, BUF_MAX_LENGTH, fp);
    for(i=0; i < num_samples; i++)
    {
      if(feof(fp) != 0)
      {
        fprintf(stderr, "# Error: file %s ends at %d.\n", options.sample_info_file, i);
        exit(0);
      }
      fgets(buf, BUF_MAX_LENGTH, fp);
      if(sscanf(buf, "%lf %lf %lf", &sample_info[i][0], &sample_info[i][1], &sample_info[i][2]) < 3)
      {
        fprintf(stderr, "# Error: Cannot read file %s.\n", options.sample_info_file);
        exit(0);
      }
      buf[0]='\0';  // clear buf

      /* reset level assignment for levels larger than the maximum level numbers */
      if(sample_info[i][0] > num_levels -1)
        sample_info[i][0] = num_levels - 1;
    }
    fclose(fp);
  }
  else   //sample_info file doest not exist, need to recalculate.
  {
    fp = fopen(options.sample_info_file, "w");
    if(fp == NULL)
    {
      fprintf(stderr, "# Error: Cannot open file %s.\n", options.sample_info_file);
      exit(0);
    }
    printf("# Dnest starts to recalculate the sample info.\n");
    fprintf(fp, "# level assignment, log likelihood, tiebreaker, ID.\n");

    //read sample
    fp_sample = fopen(options.sample_file, "r");
    if(fp_sample == NULL)
    {
      fprintf(stderr, "# Error: Cannot open file %s.\n", options.sample_file);
      exit(0);
    }
    fgets(buf, BUF_MAX_LENGTH, fp_sample);

    for(i=0; i < num_samples; i++)
    {
      read_particle(fp_sample, (void *)psample);

      sample_info[i][1] = log_likelihoods_cal_initial((void *)psample, dnest_arg);
      sample_info[i][2] = dnest_rand();

      for(j=0; j<num_levels; j++)
      {
        if(sample_info[i][1] < levels_orig[j][1])
          break;
      }

      /*j=num_levels-1;  // find out the highest allowed level
      while( (sample_info[i][1] < levels_orig[j][1]) && (j>=0) )
      {
        j--;
      }*/

      sample_info[i][0] = (double)dnest_rand_int(j); // randomly assign a level [0, j-1]

      fprintf(fp, "%d %e %f %d\n", (int)sample_info[i][0], sample_info[i][1], sample_info[i][2], 1);
    }
    fclose(fp);
    fclose(fp_sample);
  }
  //tempering with a temperature
  for(i=0; i<num_samples; i++)
    logl[i] = sample_info[i][1] / temperature;

  // finding sandwhiching levels for each samples
  for(i=0; i<num_samples; i++)
  {
    sandwhich[i] = (int)sample_info[i][0];

    for(j=sandwhich[i]; j < num_levels; j++)
    {
      if( sample_info[i][1] > levels_orig[j][1] )
        sandwhich[i] = j;
    }
    //printf("%f %d\n", logl[i], sandwhich[i]);
  }

  double *logx_samples, *logp_samples, *logP_samples;
  double logx_min, logx_max, Umin, U;
  int num_samples_thisLevel;
  double *logx_samples_thisLevel;
  SampleType *logl_samples_thisLevel;

  double left, right;

  logx_samples = malloc(num_samples * sizeof(double));
  logp_samples = malloc(num_samples * sizeof(double));
  logP_samples = malloc(num_samples * sizeof(double));

  logx_samples_thisLevel = malloc(num_samples * sizeof(double));
  logl_samples_thisLevel = malloc(num_samples * sizeof(SampleType));

  for(i=0; i<num_levels; i++)
  {
    logx_max = levels_orig[i][0];
    if(i == num_levels - 1)
      logx_min = -1.0E300;
    else
      logx_min = levels_orig[i+1][0];

    Umin = exp( logx_min - logx_max);

    // finding the samples sandwhiched by this levels
    num_samples_thisLevel = 0;
    for(j=0; j<num_samples; j++)
      if( sandwhich[j] == i )
      {
        logl_samples_thisLevel[num_samples_thisLevel].logl = sample_info[j][1]; // logl
        logl_samples_thisLevel[num_samples_thisLevel].tiebreaker = sample_info[j][2]; // tiebreaker
        logl_samples_thisLevel[num_samples_thisLevel].id = j; // id

        num_samples_thisLevel++;
      }

    //printf("%d\n", num_samples_thisLevel);

    for(j=0; j<num_samples_thisLevel; j++)
    {
      U = Umin + (1.0 - Umin) * ( 1.0/(1.0 + num_samples_thisLevel)
           + ( 1.0 - 2.0/(1.0 + num_samples_thisLevel) ) * (num_samples_thisLevel-1 - j)/(num_samples_thisLevel - 1.0) );
      logx_samples_thisLevel[j] = logx_max + log(U);
    }

    qsort(logl_samples_thisLevel, num_samples_thisLevel, sizeof(SampleType), cmp_sample);

    //printf("%f %f %d %f\n", logl_samples_thisLevel[0].logl, logl_samples_thisLevel[0].tiebreaker, logl_samples_thisLevel[0].id, logx_samples_thisLevel[0]);
    //printf("%f %f %d %f\n", logl_samples_thisLevel[1].logl, logl_samples_thisLevel[1].tiebreaker, logl_samples_thisLevel[1].id, logx_samples_thisLevel[1]);

    for(j = 0; j<num_samples_thisLevel; j++)
    {
      if(j != num_samples_thisLevel - 1)
        left = logx_samples_thisLevel[j+1];
      else if (i == num_levels - 1)
        left = -1.0E300;
      else
        left = levels_orig[i+1][0];

      if( j!= 0)
        right = logx_samples_thisLevel[j-1];
      else
        right = levels_orig[i][0];

      //printf("%e %e %e\n", right, left, logdiffexp(right, left));

      logx_samples[logl_samples_thisLevel[j].id] = logx_samples_thisLevel[j];
      logp_samples[logl_samples_thisLevel[j].id] = log(0.5)  + logdiffexp(right, left);
    }
  }

  double sum, max, logz_estimates, H_estimates, ESS;

  sum = logsumexp(logp_samples, num_samples);
  for(j = 0; j < num_samples; j++)
  {
    logp_samples[j] -= sum;
    //logP_samples[j] = logp_samples[j] + sample_info[j][1];
    logP_samples[j] = logp_samples[j] + logl[j];
  }

  logz_estimates = logsumexp(logP_samples, num_samples);

  H_estimates = -logz_estimates;
  ESS = 0.0;
  for(j=0; j<num_samples; j++)
  {
    logP_samples[j] -= logz_estimates;
    //H_estimates += exp(logP_samples[j]) * sample_info[j][1];
    H_estimates += exp(logP_samples[j]) * logl[j];
    ESS += -logP_samples[j]*exp(logP_samples[j]);
  }
  ESS = exp(ESS);

  printf("log(Z) = %f\n", logz_estimates);
  printf("H = %f\n", H_estimates);
  printf("Effective sample size = %f\n", ESS);

  post_logz = logz_estimates;

  // resample to uniform weight

  int num_ps = moreSample*ESS;
  void *posterior_sample;
  double *posterior_sample_info;
  int *posterior_sample_idx;
  int which;

  const gsl_rng_type * dnest_post_gsl_T;
  gsl_rng * dnest_post_gsl_r;

  dnest_post_gsl_T = (gsl_rng_type *) gsl_rng_default;
  dnest_post_gsl_r = gsl_rng_alloc (dnest_post_gsl_T);
#ifndef Debug
  gsl_rng_set(dnest_post_gsl_r, time(NULL));
#else
  gsl_rng_set(dnest_post_gsl_r, 8888);
  printf("# debugging, random seed %d\n", 8888);
#endif

  posterior_sample = malloc(num_ps * dnest_size_of_modeltype);
  posterior_sample_info = malloc(num_ps * sizeof(double));
  posterior_sample_idx = malloc(num_ps * sizeof(int)); // flag for which particle to save

  max = logP_samples[0];
  for(j=0; j<num_samples; j++)
    max = fmax(logP_samples[j], max);
  for(j=0; j<num_samples; j++)
    logP_samples[j] -= max;

  for(j=0; j<num_ps; j++)
  {
    while(true)
    {
      which =  gsl_rng_uniform_int(dnest_post_gsl_r, num_samples);
      if(log(gsl_rng_uniform(dnest_post_gsl_r)) < logP_samples[which])
      {
        posterior_sample_info[j] = logl[which];
        posterior_sample_idx[j]  = which;  // sample this particle
        break;
      }
    }
  }

  // read sample and pick out selected particles

  fp_sample = fopen(options.sample_file, "r");
  if(fp_sample == NULL)
  {
    fprintf(stderr, "# Error: Cannot open file %s.\n", options.sample_file);
    exit(0);
  }
  fgets(buf, BUF_MAX_LENGTH, fp_sample);
  for(i=0; i < num_samples; i++)
  {

    read_particle(fp_sample, (void *)psample);

    for(j=0; j < num_ps; j++)
    {
      if(posterior_sample_idx[j] == i)
      {
        memcpy(posterior_sample+j*dnest_size_of_modeltype, (void *)psample, dnest_size_of_modeltype);
      }
    }
    //printf("%f %f %f\n", sample[i].params[0], sample[i].params[1], sample[i].params[2]);
  }
  fclose(fp_sample);

  //save posterior sample
  fp = fopen(options.posterior_sample_file, "w");
  if(fp == NULL)
  {
    fprintf(stderr, "# Error: Cannot open file %s.\n", options.posterior_sample_file);
    exit(0);
  }
  fprintf(fp, "# %d\n", num_ps);

  for(i=0; i<num_ps; i++)
  {
    print_particle(fp, posterior_sample + i*dnest_size_of_modeltype, dnest_arg);
  }
  fclose(fp);

  //save posterior sample information
  fp = fopen(options.posterior_sample_info_file, "w");
  if(fp == NULL)
  {
    fprintf(stderr, "# Error: Cannot open file %s.\n", options.posterior_sample_info_file);
    exit(0);
  }
  fprintf(fp, "# %d\n", num_ps);
  for(i=0; i<num_ps; i++)
  {
    fprintf(fp, "%e\n", posterior_sample_info[i]);
  }
  fclose(fp);

  for(i=0; i<num_levels; i++)
   free(levels_orig[i]);
  free(levels_orig);

  for(i=0; i<num_samples; i++)
    free(sample_info[i]);
  free(sample_info);
  free(logl);

  free(logx_samples);
  free(logx_samples_thisLevel);
  free(logp_samples);
  free(logP_samples);
  free(logl_samples_thisLevel);
  free(sandwhich);
  free(psample);
  free(posterior_sample);
  free(posterior_sample_info);
  free(posterior_sample_idx);

  gsl_rng_free(dnest_post_gsl_r);

  printf("# Ends dnest postprocess.\n");
}

int cmp_sample(const void *pa, const void *pb)
{
  SampleType *a = (SampleType *)pa;
  SampleType *b = (SampleType *)pb;

  // in acesending order
  if(a->logl > b->logl)
    return true;
  if( a->logl == b->logl && a->tiebreaker > b->tiebreaker)
    return true;

  return false;
}

max_stars_repo_name: LiyrAstroph/PyCALI

max_issues_repo_name: LiyrAstroph/pyCALI

max_forks_repo_name: LiyrAstroph/pyCALI


text:
//===------------------------------------------------------------*- C++ -*-===//
//
// This file is distributed under MIT License. See LICENSE.txt for details.
//
//===----------------------------------------------------------------------===//
//
// Copyright (c) 2017 University of Kaiserslautern.
//

#pragma once

#include <nlopt.h>

namespace gosat {

class OptConfig {
public:
    OptConfig();

    OptConfig(nlopt_algorithm global_alg, nlopt_algorithm local_alg);

    virtual ~OptConfig() = default;

    int MaxEvalCount;
    int MaxLocalEvalCount;
    double RelTolerance;
    double Bound;
    double StepSize;
    unsigned InitialPopulation;
};

class NLoptOptimizer {
public:
    NLoptOptimizer();

    NLoptOptimizer(nlopt_algorithm global_alg,
                   nlopt_algorithm local_alg = NLOPT_LN_BOBYQA);

    virtual ~NLoptOptimizer() = default;

    int optimize
            (nlopt_func func, unsigned dim, double* x,
             double* min) const noexcept;

    double eval
            (nlopt_func func, unsigned dim, const double* x) const noexcept;

    bool existsRoundingError
            (nlopt_func func,
             unsigned int dim,
             const double* x,
             const double* min) const noexcept;

    void fixRoundingErrorNearZero
            (nlopt_func const func,
             unsigned dim,
             double* x,
             double* min) const noexcept;

    int refineResult(nlopt_func func, unsigned dim, double* x, double* min);

    static bool isSupportedGlobalOptAlg(nlopt_algorithm opt_alg) noexcept;

    static bool isSupportedLocalOptAlg(nlopt_algorithm local_opt_alg) noexcept;

    static bool isRequireLocalOptAlg(nlopt_algorithm opt_alg) noexcept;

    static bool isRequirePopulation(nlopt_algorithm opt_alg) noexcept;

private:
    const nlopt_algorithm m_global_opt_alg;
    const nlopt_algorithm m_local_opt_alg;
public:
    OptConfig Config;
};
}

max_stars_repo_name: abenkhadra/gosat

max_issues_repo_name: abenkhadra/gosat

max_forks_repo_name: abenkhadra/gosat


text:
// ${_warning_in_the_generated_file_not_to_edit}

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_errno.h>
#include <math.h>

#include "neqsys.h"

// Global counters
int NFEV = 0;
int NJEV = 0;
int NFJEV = 0;

int
func (const gsl_vector * x, void * params, gsl_vector * f)
{
  /*
    Best is to name all parameters k[0] ... k[P]
   */
  const double *k = (double *) params;
  const double * const y = (double *) x->data;

  /*
    Define variables for common subexpressions
   */
% for cse_token, cse_expr in func_cse_defs:
  const double ${cse_token} = ${cse_expr};
% endfor

  /*
    Calculate residuals
   */
% for i, expr in enumerate(func_new_code):
  gsl_vector_set(f, ${i}, ${expr});
% endfor

  NFEV++;
  return GSL_SUCCESS;
}

int
jac (const gsl_vector * x, void *params, gsl_matrix * J)
{
  const double *k = (double *) params;
  const double * const y = (double *) x->data;

  /*
    Define variables for common subexpressions
   */
% for cse_token, cse_expr in jac_cse_defs:
  const double ${cse_token} = ${cse_expr};
% endfor

  /*
    Populate the NY times NY Jacobian matrix
   */
% for i, expr in enumerate(jac_new_code):
  gsl_matrix_set (J, ${i // NX}, ${i % NX}, ${expr});
% endfor

  NJEV++;
  return GSL_SUCCESS;
}

int
fdf (const gsl_vector * x, void * params, gsl_vector * f, gsl_matrix * J)
{
  const double * const k = (double *) params;
  const double * const y = (double *) x->data;

  /*
    Define variables for common subexpressions
   */
% for cse_token, cse_expr in fj_cse_defs:
  const double ${cse_token} = ${cse_expr};
% endfor

  /*
    Calculate residuals
   */
% for i, expr in enumerate(fj_func_new_code):
  gsl_vector_set(f, ${i}, ${expr});
% endfor

  /*
    Populate the NY times NY Jacobian matrix
   */
% for i, expr in enumerate(fj_jac_new_code):
  gsl_matrix_set (J, ${i // NX}, ${i % NX}, ${expr});
% endfor

  NFJEV++;
  return GSL_SUCCESS;
}

max_stars_repo_name: bjodah/symneqsys

max_issues_repo_name: bjodah/symneqsys

max_forks_repo_name: bjodah/symneqsys


text:
#pragma once

#include <gsl.h>

#include "types/MT3_MathChar.h"
#include "types/MT3_Radical.h"
#include "types/MT3_Accent.h"
#include "types/MT3_VCenter.h"
#include "types/MT3_Overline.h"
#include "types/MT3_Underline.h"
#include "types/MT3_GenFraction.h"
#include "types/MT3_LeftRight.h"
#include "types/MT3_Script.h"
#include "types/MT3_BigOp.h"
#include "types/MT3_SubBox.h"
#include "types/MT3_MList.h"
#include "types/MT3_Kind.h"
#include "types/MT3_MPen.h"
#include "types/MT3_MSpace.h"
#include "types/MT3_Style.h"
#include "types/MT3_Choice.h"


text:
#ifndef FORCE_H
#define FORCE_H

#include <functional>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_errno.h>
#include <vector>
#include "jspec2/electron_beam.h"
#include "jspec2/constants.h"

using std::vector;
enum class ForceFormula {PARKHOMCHUK, NONMAG_DERBENEV, NONMAG_MESHKOV, NONMAG_NUM1D, NONMAG_NUM3D, MESHKOV, DSM};

class FrictionForceSolver{
protected:
    double time_cooler;
    double mag_field = 0;
    static constexpr double k_f = -4*k_pi*k_c*k_c*k_ke*k_ke*k_e*k_e*k_e/(k_me*1e6);
    static constexpr double k_wp = 4*k_pi*k_c*k_c*k_e*k_ke/(k_me*1e6);
    static constexpr double k_rho_min = k_e*k_ke*k_c*k_c/(k_me*1e6);
    virtual void init(const ElectronBeam& ebeam){};
    virtual void fin(){};
    double max_impact_factor(double v_dlt, int charge_number,double density_e);
public:
    void set_time_cooler(double t){time_cooler = t;}
    void set_mag_field(double x){mag_field = x;}
    double t_cooler() const {return time_cooler;}
    virtual void friction_force(int charge_number, int ion_number,
            const vector<double>& v_tr, const vector<double>& v_l, const vector<double>& density,
            const ElectronBeam& ebeam, vector<double>& force_tr, vector<double>& force_long) = 0;
};

class ForcePark: public FrictionForceSolver {
private:
    static constexpr double k_f = -4*k_c*k_c*k_ke*k_ke*k_e*k_e*k_e/(k_me*1e6);
    double t_eff = 0; //Effective temperature.
    double v_eff = 0; //Effective velocity.
    void rho_lamor_dlt2_eff_e(double v2_eff_e, double mag_field, const vector<double>& v_rms_l, const vector<double>& v_rms_t, ElectronBeam::Temperature tpr,
                          int ion_number, vector<double>& dlt2_eff_e, vector<double>& rho_lamor);
    constexpr double dlt(ElectronBeam::Temperature tpr,  double v2, const vector<double>& dlt2_eff_e, int i);
    constexpr double lc(ElectronBeam::Temperature tpr, double rho_max, double rho_min, const vector<double>& rho_lamor, int i);
public:
    void set_t_eff(double x){t_eff = x; v_eff = sqrt(t_eff*k_c*k_c/(k_me*1e6));}
    void set_v_eff(double v){v_eff = v; t_eff = v_eff*v_eff*k_me*1e6/(k_c*k_c);}
    virtual void friction_force(int charge_number, int ion_number,
            const vector<double>& v_tr, const vector<double>& v_l, const vector<double>& density,
            const ElectronBeam& ebeam, vector<double>& force_tr, vector<double>& force_long) override;
};

class ForceNonMag: public FrictionForceSolver {
protected:
    bool smooth_rho_max = false;
    virtual double f_const(int charge_number) const {return charge_number*charge_number*k_f;}
    static constexpr double rho_min_const(int charge_number) {return charge_number*k_rho_min;}
    static constexpr double rho_max_1(int charge_number, double density_e){return pow(3*charge_number/density_e, 1.0/3);}
    double rho_max_2(double v) const {return v*time_cooler;}
    double rho_max(int charge_number, double v2, double ve2, double ne) const;
    virtual void force(double v, double v_tr, double v_l, double v2, double ve_tr, double ve_l, double ve2,
                       double f_const,double rho_min_const, int charge_number, double ne,
                       double& force_tr, double& force_l) = 0;

public:
    void set_smooth_rho_max(bool b){smooth_rho_max = b;}
    virtual void friction_force(int charge_number, int ion_number,
            const vector<double>& v_tr, const vector<double>& v_l, const vector<double>& density,
            const ElectronBeam& ebeam, vector<double>& force_tr, vector<double>& force_long) override;
};

class ForceNonMagDerbenev: public ForceNonMag {
private:
    constexpr double rho_max_ve_tr(double ve_tr, double ne);
    void force(double v, double v_tr, double v_l, double v2, double ve_tr, double ve_l, double ve2,
                               double f_const, double rho_min_const, int charge_number, double ne,
                               double& force_tr, double& force_l) override;
public:

};

class ForceNonMagMeshkov: public ForceNonMag {
private:
    void force(double v, double v_tr, double v_l, double v2, double ve_tr, double ve_l, double ve2,
                               double f_const, double rho_min_const, int charge_number, double ne,
                               double& force_tr, double& force_l) override;
public:

};

class ForceNonMagNumeric1D: public ForceNonMag {
private:
    double f_const(int charge_number) const override {return 0.28209479177387814*charge_number*charge_number*k_f;} //Coef - 1/sqrt(4*pi)
    static constexpr double k_f = 2*sqrt(2*k_pi)*k_pi*k_c*k_c*k_ke*k_ke*k_e*k_e*k_e/(k_me*1e6);

    size_t limit = 100;
    double espabs = 1e-6;
    double esprel = 1e-6;
    struct P{
        double v_tr;
        double v_l;
        double ve_tr;
        double ve_l;
        int flag;   //0: calculate B_tr; else: calculate B_l;
    };

    #ifdef _OPENMP
    vector<gsl_integration_workspace*> gw;
    vector<P> p;
    #else
    gsl_integration_workspace *gw = nullptr;
    P p;
    #endif // _OPENMP
    double b(double q, void* params);
    void force(double v, double v_tr, double v_l, double v2, double ve_tr, double ve_l, double ve2,
                               double f_const, double rho_min_const, int charge_number, double ne,
                               double& force_tr, double& force_l) override;
public:
    void set_espabs(double x){espabs = x;}
    void set_esprel(double x){esprel = x;}
    ForceNonMagNumeric1D(int n=100);
    ~ForceNonMagNumeric1D();
};

//
class ForceNonMagNumeric3D: public ForceNonMag {
private:
    size_t limit = 100;
    double espabs = 1e-5;
    double esprel = 1e-3;
    struct P{
        double v_tr;
        double v_l;
        double ve_tr;
        double ve_l;
        double vtr;
        double vl;
        double rho_max;
        int charge_number;
        int flag;   //0: calculate B_tr; else: calculate B_l;
    };

    bool use_gsl = false;
    bool use_mean_rho_min = false;

    #ifdef _OPENMP
    vector<gsl_integration_workspace*> giw;
    vector<gsl_integration_workspace*> gmw;
    vector<gsl_integration_workspace*> gow;
    vector<P> p;
    #else
    gsl_integration_workspace *giw = nullptr;
    gsl_integration_workspace *gmw = nullptr;
    gsl_integration_workspace *gow = nullptr;
    P p;
    #endif // _OPENMP

    #ifdef _OPENMP
    static double mean_rho_min;
    static double mean_lc;
    #pragma omp threadprivate(mean_rho_min, mean_lc)
    static bool first_run;
    static vector<vector<double>> exp_vtr;
    static vector<double> hlf_v2tr;
    static vector<double> hlf_v2l;
    static vector<vector<double>> vtr_cos;
    static vector<double> vl;
    static vector<double> vtr;
    static vector<vector<double>> v2tr_sin2;
    #pragma omp threadprivate(exp_vtr, hlf_v2tr, hlf_v2l, vtr_cos, vl, vtr, v2tr_sin2, first_run)
    #else
    double mean_rho_min = 0;
    double mean_lc = 0;
    vector<vector<double>> exp_vtr;
    vector<double> hlf_v2tr;
    vector<double> hlf_v2l;
    vector<vector<double>> vtr_cos;
    vector<double> vl;
    vector<double> vtr;
    vector<vector<double>> v2tr_sin2;
    bool first_run = true;
    #endif // _OPENMP
    bool const_tmpr = true;
    int n_tr = 20;
    int n_l = 10;
    int n_phi = 10;
    double d;
    double f_inv_norm;

    void pre_int(double sgm_vtr, double sgm_vl);
    void calc_exp_vtr(double sgm_vtr, double sgm_vl);

    void init(const ElectronBeam& ebeam) override;
    double inner_integrand(double phi, void* params);
    double middle_integrand(double vl, void* params);
    double outter_integrand(double vtr, void* params);
    double inner_norm_integrand(double vl, void* params);
    double outter_norm_integrand(double vtr, void* params);
    void force_grid(double v, double v_tr, double v_l, double v2, double ve_tr, double ve_l, double ve2,
                               double f_const, double rho_min_const, int charge_number, double ne,
                               double& force_tr, double& force_l);
    void force_gsl(double v, double v_tr, double v_l, double v2, double ve_tr, double ve_l, double ve2,
                               double f_const, double rho_min_const, int charge_number, double ne,
                               double& force_tr, double& force_l);
    void force(double v, double v_tr, double v_l, double v2, double ve_tr, double ve_l, double ve2,
                               double f_const, double rho_min_const, int charge_number, double ne,
                               double& force_tr, double& force_l) override;
public:
    void set_espabs(double x){espabs = x;}
    void set_esprel(double x){esprel = x;}
    void set_gsl(bool b) {use_gsl = b;}
    void set_mean_rho_min(bool b) {use_mean_rho_min = b;}
    void set_grid(int ntr, int nl, int nphi){n_tr = ntr; n_l = nl; n_phi = nphi; first_run = true;}
    ForceNonMagNumeric3D(int n=100);
    ~ForceNonMagNumeric3D();
};

//gsl function wrapper for member functions in class.
template< typename F >
    class gsl_function_pp : public gsl_function {
    public:
    gsl_function_pp(const F& func) : _func(func) {
    function = &gsl_function_pp::invoke;
    params=this;
    }
    private:
    const F& _func;
    static double invoke(double x, void *params) {
    return static_cast<gsl_function_pp*>(params)->_func(x);
    }
};

class ForceMeshkov: public FrictionForceSolver {
protected:
    //A fudge factor to smooth the friction force shape. Using the
        // "Classical" default definition here from the BETACOOL documentation
    double k = 2;
    void force(double ve_tr, double ve_l, double ve2_tr, double ve2_l, double v_tr, double v_l, double v2,
               double rho_min_const, int charge_number,  double density, double f_const, double& force_tr,double& force_l);
public:
    void set_smooth_factor(double x){k = x;}
    virtual void friction_force(int charge_number, int ion_number,
            const vector<double>& v_tr, const vector<double>& v_l, const vector<double>& density,
            const ElectronBeam& ebeam, vector<double>& force_tr, vector<double>& force_long) override;
};

class ForceDSM: public FrictionForceSolver {    //Derbenev-Skrinsky-Meshkov formula for magnetized cooling.
    bool mag_only = false;

     //A fudge factor to smooth the friction force shape. Using the
    // "Classical" default definition here from the BETACOOL documentation
    double k = 2;
    int n_a = 100;
    int n_ve = 100;
    bool const_tpr = true;
    #ifdef _OPENMP
    static bool first_run;
    static vector<double> a;
    static vector<double> cos_a;
    static vector<double> tan_a;
    static vector<double> t2;
    static vector<double> ve;
    static vector<double> exp_ve2;
    #pragma omp threadprivate(a, cos_a, tan_a, t2, ve, exp_ve2, first_run)
    #else
    bool first_run = true;
    vector<double> a;
    vector<double> cos_a;
    vector<double> tan_a;
    vector<double> t2;
    vector<double> ve;
    vector<double> exp_ve2;
    #endif // _OPENMP

    int n_tr = 20;
    int n_l = 10;
    int n_phi = 10;
    #ifdef _OPENMP
    static bool first_run_fa;
    static double f_inv_norm;
    static vector<vector<double>> exp_vtr;
    static vector<double> hlf_v2tr;
    static vector<double> hlf_v2l;
    static vector<vector<double>> vtr_cos;
    static vector<double> vl;
    static vector<double> vtr;
    static vector<vector<double>> v2tr_sin2;
    #pragma omp threadprivate(exp_vtr, hlf_v2tr, hlf_v2l, vtr_cos, vl, vtr, v2tr_sin2, first_run_fa, f_inv_norm)
    #else
    bool first_run_fa = true;
    double f_inv_norm;
    vector<vector<double>> exp_vtr;
    vector<double> hlf_v2tr;
    vector<double> hlf_v2l;
    vector<vector<double>> vtr_cos;
    vector<double> vl;
    vector<double> vtr;
    vector<vector<double>> v2tr_sin2;
    #endif // _OPENMP

protected:
    void init(const ElectronBeam& ebeam) override;
    void pre_int(double sgm_vtr, double sgm_vl);
    void calc_exp_vtr(double sgm_vtr, double sgm_vl);
    void calc_alpha();
    void calc_ve();
    void calc_exp_ve2(double ve2_l, vector<double>& t2);
    void force(double ve_tr, double ve_l, double ve2_tr, double ve2_l, double v_tr, double v_l, double v2,
               double rho_min_const, int charge_number,  double density, double f_const, double& force_tr,double& force_l);
public:
    void set_smooth_factor(double x){k = x;}
    void set_steps(int i) {n_a = i; n_ve = i; first_run = true;}
    void set_grid(int ntr, int nl, int nphi){n_tr = ntr; n_l = nl; n_phi = nphi; first_run_fa = true;}
    void set_mag_only(bool b){mag_only = b;}
    virtual void friction_force(int charge_number, int ion_number,
            const vector<double>& v_tr, const vector<double>& v_l, const vector<double>& density,
            const ElectronBeam& ebeam, vector<double>& force_tr, vector<double>& force_long) override;
};
#endif // FORCE_H


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

#pragma once

#include <unordered_map>
#include <mutex>
#include <gsl/gsl>

#include <typeindex>

#include "arcana/containers/ticketed_collection.h"
#include "arcana/string.h"

#if NDEBUG
#define FIRE_OBJECT_TRACE(...)
#else
#define OBJECT_TRACE_ENABLED
#define FIRE_OBJECT_TRACE(channelName, instance, evt) \
    ::mira::object_trace::fire(channelName, instance, evt)
#endif

namespace mira
{
    /*
        Simple static analysis of data for exposing
    */
    class object_trace final
    {
        using channel_callback = std::function<void(std::intptr_t, const void*)>;
        using global_callback = std::function<void(const char*, std::intptr_t, const void*)>;

        using channel_ticketed_collection = ticketed_collection<channel_callback, std::recursive_mutex>;
        using global_ticketed_collection = ticketed_collection<global_callback, std::recursive_mutex>;
    public:
        using channel_ticket = channel_ticketed_collection::ticket;
        using global_ticket = global_ticketed_collection::ticket;

        object_trace() = delete;

        /*
        Adds a listener to type T on the specified channel.

        The channel is useful to split data into logical streams.
        */
        template<typename T>
        static auto listen(const char* channel, std::function<void(std::intptr_t, T&)> callback)
        {
            return add_listener(channel, std::type_index{ typeid(std::decay_t<T>) }, [callback = std::move(callback)](std::intptr_t instance, const void* data)
            {
                callback(instance, const_cast<T&>(*reinterpret_cast<const T*>(data)));
            });
        }

        /*
            Adds a listener to type T to all channels.
        */
        template<typename T>
        static auto listen(std::function<void(const char*, std::intptr_t, T&)> callback)
        {
            return add_listener(std::type_index{ typeid(std::decay_t<T>) },
                [callback = std::move(callback)](const char* channel, std::intptr_t instance, const void* data)
                {
                    callback(channel, instance, const_cast<T&>(*reinterpret_cast<const T*>(data)));
                });
        }

        /*
            Fires an event in the specified channel.
        */
        template<typename T>
        static void fire(const char* channel, const void* instance, const T& evt)
        {
            fire_event(channel, (std::intptr_t)instance, std::type_index{ typeid(std::decay_t<T>) }, &evt);
        }
    private:
        using channels = std::map<std::string, channel_ticketed_collection, string_compare>;

        struct listeners
        {
            channels ChannelListeners;
            global_ticketed_collection GlobalListeners;
        };

        static channel_ticket add_listener(const char* channel, const std::type_index& type, channel_callback callback);
        static global_ticket add_listener(const std::type_index& type, global_callback callback);

        static void fire_event(const char* channel, std::intptr_t instance, const std::type_index& type, const void* data);

        static std::recursive_mutex m_mutex;
        static std::unordered_map<std::type_index, listeners> m_typemap;
    };
}

max_stars_repo_name: syntheticmagus/mageslam

max_issues_repo_name: syntheticmagus/mageslam

max_forks_repo_name: syntheticmagus/mageslam


text:
#pragma once

// View structures
// STATUS: prototype, NSC

#include <vector>

#include <gsl.h>

namespace duck {
/* span<T> represents a reference to a segment of a T array.
 * Equivalent to a (T* base, int len).
 */
using gsl::span;

/* string_view: Const reference (pointer) to a sequence of char.
 * Does not own the data, only points to it.
 * The char sequence may not be null terminated.
 *
 * std::string_view exists in C++17, but not before sadly.
 */
using string_view = gsl::cstring_span;

using gsl::to_string; // Create new std::string from string_view

bool is_prefix_of (string_view prefix, string_view str);
bool is_prefix_of (char prefix, string_view str);

// Split at separator. Does not remove empty parts, does not trim whitespace.
std::vector<string_view> split (char separator, string_view text);
} // namespace duck

max_stars_repo_name: lereldarion/duck

max_issues_repo_name: lereldarion/duck

max_forks_repo_name: lereldarion/duck


text:
/* gsl_histogram2d_maxval.c
 * Copyright (C) 2000  Simone Piccardi
 *
 * This library 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 library; if not, write to the
 * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
 * Boston, MA 02111-1307, USA.
 */
/***************************************************************
 *
 * File gsl_histogram2d_maxval.c:
 * Routine to find maximum and minumum content of a 2D hisogram.
 * Need GSL library and header.
 * Contains the routines:
 * gsl_histogram2d_max_val find max content values
 * gsl_histogram2d_min_val find min content values
 * gsl_histogram2d_bin_max find coordinates of max contents bin
 * gsl_histogram2d_bin_min find coordinates of min contents bin
 *
 * Author: S. Piccardi
 * Jan. 2000
 *
 ***************************************************************/
#include <config.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_histogram2d.h>
/*
 * Return the maximum contents value of a 2D histogram
 */
double
gsl_histogram2d_max_val (const gsl_histogram2d * h)
{
  const size_t nx = h->nx;
  const size_t ny = h->ny;
  size_t i;
  double max = h->bin[0 * ny + 0];

  for (i = 0; i < nx * ny; i++)
    {
      if (h->bin[i] > max)
	{
	  max = h->bin[i];
	}
    }
  return max;
}

/*
 * Find the bin index for maximum value of a 2D histogram
 */

void
gsl_histogram2d_max_bin (const gsl_histogram2d * h, size_t * imax_out, size_t * jmax_out)
{
  const size_t nx = h->nx;
  const size_t ny = h->ny;
  size_t imax = 0, jmax = 0, i, j;
  double max = h->bin[0 * ny + 0];

  for (i = 0; i < nx; i++)
    {
      for (j = 0; j < ny; j++)
	{
	  double x = h->bin[i * ny + j];

	  if (x > max)
	    {
	      max = x;
	      imax = i;
	      jmax = j;
	    }
	}
    }

  *imax_out = imax;
  *jmax_out = jmax;
}

/*
 * Return the minimum contents value of a 2D histogram
 */

double
gsl_histogram2d_min_val (const gsl_histogram2d * h)
{
  const size_t nx = h->nx;
  const size_t ny = h->ny;
  size_t i;
  double min = h->bin[0 * ny + 0];

  for (i = 0; i < nx * ny; i++)
    {
      if (h->bin[i] < min)
	{
	  min = h->bin[i];
	}
    }

  return min;
}
/*
 * Find the bin index for minimum value of a 2D histogram
 */
void
gsl_histogram2d_min_bin (const gsl_histogram2d * h, size_t * imin_out, size_t * jmin_out)
{
  const size_t nx = h->nx;
  const size_t ny = h->ny;
  size_t imin = 0, jmin = 0, i, j;
  double min = h->bin[0 * ny + 0];

  for (i = 0; i < nx; i++)
    {
      for (j = 0; j < ny; j++)
	{
	  double x = h->bin[i * ny + j];

	  if (x < min)
	    {
	      min = x;
	      imin = i;
	      jmin = j;
	    }
	}
    }

  *imin_out = imin;
  *jmin_out = jmin;
}


text:
#include <jni.h>
#include <assert.h>
#include <lapacke.h>

extern void dtrmv_(char *uplo, char *trans, char *diag, int *n, double *A,
                   int *lda, double *x, int *incx);

#define jniRowMajor 101
#define jniColMajor 102

#define jniNoTrans 111
#define jniTrans   112
#define jniConjTrans    113

#define jniUpper   121
#define jniLower   122

#define jniNonUnit 131
#define jniUnit    132

#define jniLeft    141
#define jniRight   142

JNIEXPORT jint Java_JAMAJni_EigenvalueDecomposition_dgeev
(JNIEnv *env, jclass klass, jint matrix_layout, jchar jobvl, jchar jobvr,
 jint n, jdoubleArray a, jint lda, jdoubleArray wr, jdoubleArray wi,
 jdoubleArray vl, jint ldvl, jdoubleArray vr, jint ldvr){

    double *aElems, *wrElems, *wiElems, *vlElems, *vrElems;
    int info;

    aElems = (*env)-> GetDoubleArrayElements (env, a, NULL);
    wrElems = (*env)-> GetDoubleArrayElements (env, wr, NULL);
    wiElems = (*env)-> GetDoubleArrayElements (env, wi, NULL);
    vlElems = (*env)-> GetDoubleArrayElements (env, vl, NULL);
    vrElems = (*env)-> GetDoubleArrayElements (env, vr, NULL);

    assert(aElems && wrElems && wiElems && vlElems && vrElems);

    info = LAPACKE_dgeev((int) matrix_layout, (char) jobvl, (char) jobvr, (lapack_int) n, aElems, (lapack_int) lda, wrElems, wiElems, vlElems, (lapack_int) ldvl, vrElems, (lapack_int) ldvr);

    (*env)-> ReleaseDoubleArrayElements (env, a, aElems, 0);
    (*env)-> ReleaseDoubleArrayElements (env, vl, vlElems, 0);
    (*env)-> ReleaseDoubleArrayElements (env, vr, vrElems, 0);
    (*env)-> ReleaseDoubleArrayElements (env, wr, wrElems, 0);
    (*env)-> ReleaseDoubleArrayElements (env, wi, wiElems, 0);

    return info;

}

JNIEXPORT void JNICALL Java_JAMAJni_EigenvalueDecomposition_dsyev
(JNIEnv *env, jclass obj, jint layout, jchar jobz, jchar uplo, jint n,
 jdoubleArray ja, jint lda, jdoubleArray jw)
{
    jdouble *a = (*env)-> GetDoubleArrayElements (env, ja, NULL);
    jdouble *w = (*env)-> GetDoubleArrayElements (env, jw, NULL);

    assert(a && w);

    LAPACKE_dsyev((int) layout, (char) jobz, (char) uplo, (lapack_int) n, a, (lapack_int) lda, w);

    (*env)-> ReleaseDoubleArrayElements (env, ja, a, 0);
    (*env)-> ReleaseDoubleArrayElements (env, jw, w, 0);
}

JNIEXPORT void JNICALL Java_JAMAJni_EigenvalueDecomposition_dtrsm
(JNIEnv *env, jclass obj, jint layout, jchar jobz, jchar uplo, jint n,
 jdoubleArray ja, jint lda, jdoubleArray jw)
{
    jdouble *a = (*env)-> GetDoubleArrayElements (env, ja, NULL);
    jdouble *w = (*env)-> GetDoubleArrayElements (env, jw, NULL);

    assert(a && w);

    LAPACKE_dsyev((int) layout, (char) jobz, (char) uplo, (lapack_int) n, a, (lapack_int) lda, w);

    (*env)-> ReleaseDoubleArrayElements (env, ja, a, 0);
    (*env)-> ReleaseDoubleArrayElements (env, jw, w, 0);
}


text:
#ifndef NonLinearLSQ_h
#define NonLinearLSQ_h
/**
 * @file
 * $Revision: 4928 $
 * $Date: 2012-12-21 10:44:12 -0700 (Fri, 21 Dec 2012) $
 * $Id: NonLinearLSQ.h 4928 2012-12-21 17:44:12Z slambright@GS.DOI.NET $
 *
 *   Unless noted otherwise, the portions of Isis written by the USGS are
 *   public domain. See individual third-party library and package descriptions
 *   for intellectual property information, user agreements, and related
 *   information.
 *
 *   Although Isis has been used by the USGS, no warranty, expressed or
 *   implied, is made by the USGS as to the accuracy and functioning of such
 *   software and related material nor shall the fact of distribution
 *   constitute any such warranty, and no responsibility is assumed by the
 *   USGS in connection therewith.
 *
 *   For additional information, launch
 *   $ISISROOT/doc//documents/Disclaimers/Disclaimers.html
 *   in a browser or see the Privacy & Disclaimers page on the Isis website,
 *   http://isis.astrogeology.usgs.gov, and the USGS privacy and disclaimers on
 *   http://www.usgs.gov/privacy.html.
 */

#include <iostream>
#include <sstream>
#include <string>
#include <vector>
#include <algorithm>

#include <QString>

#include <tnt_array1d.h>
#include <tnt_array1d_utils.h>
#include <tnt_array2d.h>
#include <tnt_array2d_utils.h>

#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_multifit_nlin.h>

#include "IException.h"

namespace Isis {

  /**
   * @brief NonLinearLSQ Computes a fit using a Levenberg-Marquardt algorithm
   *
   * This virtual base class uses the GSL toolkit to apply the
   * Levenberg-Marquardt algorithm to fit data to a non-linear equation using
   * least squares.
   *
   * @ingroup Utility
   * @author 2007-11-15 Kris Becker
   *
   * @internal
   */
  class NonLinearLSQ {
    public:
      typedef TNT::Array1D<double> NLVector;
      typedef TNT::Array2D<double> NLMatrix;

      //  Constructors and Destructor
      NonLinearLSQ() : _fitParms(), _uncert(), _nIters(0), _maxIters(50),
                       _status(0), _userTerminated(false), _userMessage() {}

      /** Destructor */
      virtual ~NonLinearLSQ() { }

      virtual int nSize() const = 0;
      virtual int nParms() const = 0;

      /**
       * @brief Sets the maximum number of iterations
       *
       * @param m User provides the maximum number iterations
       */
      void setMaxIters(int m) { _maxIters = m; }

      /**
       * @brief Maximum number iterations for valid solution
       *
       *
       * @return int Maximum resolutions
       */
      int maxIters() const { return (_maxIters); }

      virtual NLVector guess()  = 0;
      virtual NLVector f_x(const NLVector &x) = 0;
      virtual NLMatrix df_x(const NLVector &x) = 0;

      virtual double absErr() const { return (1.0E-4); }
      virtual double relErr() const { return (1.0E-4); }

      int curvefit();
      /** Return status of last fit processing */
      inline int status() const { return (_status); }
      /** Determine success from last fit processing */
      inline bool success() const { return (_status == GSL_SUCCESS); }
      /** Check status for success of the given condition  */
      inline bool success(int status) const { return (status == GSL_SUCCESS); }
      /** Return error message pertaining to last fit procesing */
      inline QString statusstr() const {
        return (QString(gsl_strerror(_status)));
      }
      /** Return error message given status condition */
      inline QString statusstr(int status) const {
        return (QString(gsl_strerror(status)));
      }

      /** Default interation test simply returns input status */
      virtual int checkIteration(const int Iter, const NLVector &fitcoefs,
                                 const NLVector &uncerts, double cplxconj,
                                 int Istatus) {
        return (Istatus);
      }

      /** Return coefficients from last fit processing */
      inline NLVector coefs() const { return (_fitParms); }
      /** Return uncertainties from last fit processing */
      inline NLVector uncert() const { return (_uncert); }
      /** Return number of iterations from last fit processing */
      inline int nIterations() const { return (_nIters); }

    protected:
      void Terminate(const QString &message = "");
      void Abort(const QString &reason = "");

      bool doContinue() const { return (!_userTerminated); }

    private:
      NLVector    _fitParms;
      NLVector    _uncert;
      int         _nIters;
      int         _maxIters;
      int         _status;
      bool        _userTerminated;
      QString _userMessage;

      struct _nlsqPointer {
          NonLinearLSQ *nlsq;
      };
      static int f(const gsl_vector *x, void *params, gsl_vector *fx);
      static int df(const gsl_vector *x, void *params, gsl_matrix *J);
      static int fdf(const gsl_vector *x, void *params, gsl_vector *fx,
                     gsl_matrix *J);

      NLVector gslToNlsq(const gsl_vector *v) const;
      NLMatrix gslToNlsq(const gsl_matrix *m) const;
      gsl_vector *NlsqTogsl(const NLVector &v, gsl_vector *gv = 0) const;
      gsl_matrix *NlsqTogsl(const NLMatrix &m, gsl_matrix *gm = 0) const;
      NLVector getUncertainty(const gsl_matrix *m) const;
  };

}     // namespace Isis
#endif


text:
/******************************************************************************

  (c) 2017 - 2019 Scientific Computation Research Center,
      Rensselaer Polytechnic Institute. All rights reserved.

  This work is open source software, licensed under the terms of the
  BSD license as described in the LICENSE file in the top-level directory.

*******************************************************************************/
#ifndef MSI_PETSC_VERSION_H_
#define MSI_PETSC_VERSION_H_
#include <petscversion.h>
#include <petsc.h>
#include <utility>

template <typename ... Args>
auto PetscSetPCBackend(Args&&... args) -> decltype(
#if PETSC_VERSION_GT(3,9,5)
  PCFactorSetMatSolverType(std::forward<Args>(args)...)
#else
  PCFactorSetMatSolverPackage(std::forward<Args>(args)...)
#endif
  )
{
#if PETSC_VERSION_GT(3,9,5)
  return PCFactorSetMatSolverType(std::forward<Args>(args)...);
#else
  return PCFactorSetMatSolverPackage(std::forward<Args>(args)...);
#endif
}
#endif

max_stars_repo_name: SCOREC/msi

max_issues_repo_name: SCOREC/msi

max_forks_repo_name: SCOREC/msi


text:
/*
 * BRAINS
 * (B)LR (R)everberation-mapping (A)nalysis (I)n AGNs with (N)ested (S)ampling
 * Yan-Rong Li, liyanrong@ihep.ac.cn
 * Thu, Aug 4, 2016
 */

#ifdef SpecAstro

#include <math.h>
#include <stddef.h>
#include <string.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <fftw3.h>

#include "brains.h"

/*
 * calculate SA phase and 2D TF
 */
void calculate_sa_transfun2d_from_blrmodel(const void *pm, double *transv, double *trans2d, int n_vel, int flag_save)
{
  int i;
  double *sa_model, *rm_model;

  /* share all parameters */
  if(parset.flag_sa_par_mutual == 0)
  {
    gen_cloud_sample(pm, 5, 0);
    transfun_2d_cal_with_sample(transv, trans2d, n_vel);
    calculate_sa_with_sample(pm);
  }
  /* share mass and inclination */
  else if(parset.flag_sa_par_mutual == 1)
  {
    rm_model = (double *)pm;
    sa_model = rm_model + num_params_blr;
    sa_model[idx_sa_par_mutual[0]] = rm_model[idx_rm_par_mutual[0]]; //mbh
    sa_model[idx_sa_par_mutual[1]] = rm_model[idx_rm_par_mutual[1]]; //inc

    gen_cloud_sample(pm, 2, flag_save);
    transfun_2d_cal_with_sample(transv, trans2d, n_vel);

    gen_sa_cloud_sample((void *)sa_model, 3, flag_save);
    calculate_sa_with_sample(pm);
  }
  /* share inclination, opening angle and all dynamical parameters */
  else
  {
    rm_model = (double *)pm;
    sa_model = rm_model + num_params_blr;

    sa_model[idx_sa_par_mutual[1]] = rm_model[idx_rm_par_mutual[1]]; //inc
    sa_model[idx_sa_par_mutual[2]] = rm_model[idx_rm_par_mutual[2]]; //opening angle
    /* dynamical parameters from mass on */
    for(i=idx_sa_par_mutual[0]; i<num_params_sa_blr_model; i++)
    {
      sa_model[i] = rm_model[i]; /* note the BLR models are identical */
    }

    gen_cloud_sample(pm, 2, flag_save);
    transfun_2d_cal_with_sample(transv, trans2d, n_vel);

    gen_sa_cloud_sample((void *)sa_model, 3, flag_save);
    calculate_sa_with_sample(pm);
  }
}

/*
 * calculate SA phase and 1D TF
 */
void calculate_sa_transfun_from_blrmodel(const void *pm, int flag_save)
{
  double *sa_model, *rm_model;

  if(parset.flag_sa_par_mutual == 0)
  {
    gen_cloud_sample(pm, 4, flag_save);
    transfun_1d_cal_with_sample();
    calculate_sa_with_sample(pm);
  }
  else
  {
    rm_model = (double *)pm;
    sa_model = rm_model + num_params_blr;
    sa_model[idx_sa_par_mutual[1]] = rm_model[idx_rm_par_mutual[1]]; // inc
    sa_model[idx_sa_par_mutual[2]] = rm_model[idx_rm_par_mutual[2]]; // opening angle

    gen_cloud_sample(pm, 1, flag_save);
    transfun_1d_cal_with_sample();

    gen_sa_cloud_sample((void *)sa_model, 3, flag_save);
    calculate_sa_with_sample(pm);
  }
}

/*
 * calculate SA phase and line profile from cloud sample
 */
void calculate_sa_with_sample(const void *pm)
{
  int i, j, k, idV;
  double V, V_offset, dV, y, z, alpha, beta, flux_norm, phase, *phase_norm, *alpha_cent, *beta_cent;
  double DA, PA, FA, CO, cos_PA, sin_PA;
  double *pmodel = (double *)pm;

  phase_norm = workspace_phase;
  alpha_cent = phase_norm + n_vel_sa_data;
  beta_cent = alpha_cent + n_vel_sa_data;

  /* angular size distance */
  DA = exp(pmodel[num_params_blr + num_params_sa_blr_model]);
  /* position angle */
  PA = pmodel[num_params_blr + num_params_sa_blr_model + 1]/180.0 * PI;
  /* line flux scaling factor */
  FA = exp(pmodel[num_params_blr + num_params_sa_blr_model+2]);
  /*
   * line center offset: V+dV = (w-(w0+dw0))/(w0+dw0) ==> dV = - dw0/w0
   * equivalent to redshift offset:  dz = -(1+z) dw0/w0
   */
  CO = -pmodel[num_params_blr + num_params_sa_blr_model+3]/parset.sa_linecenter * C_Unit;

  /* North through East (E of N), 180 deg ambiguity */
  cos_PA = cos(PA);
  sin_PA = sin(PA);

  dV = vel_sa_data[1] - vel_sa_data[0];

  for(i=0; i<n_vel_sa_data; i++)
  {
    for(k=0; k<n_base_sa_data; k++)
    {
      phase_sa[k*n_vel_sa_data + i] = 0.0;
    }
    phase_norm[i] = 0.0;
    alpha_cent[i] = 0.0;
    beta_cent[i] = 0.0;
  }

  for(i=0; i<parset.n_cloud_per_task; i++)
  {
    y = clouds_alpha[i];
    z = clouds_beta[i];

    //alpha = y * cos_PA + z * sin_PA;
    //beta = -y * sin_PA + z * cos_PA;
    alpha = y;
    beta = z;

    for(j=0; j< parset.n_vel_per_cloud; j++)
    {
      V = clouds_vel[i*parset.n_vel_per_cloud + j] + CO;
      V_offset = V + bin_offset * dV; /* bin type: center or left edge */
      if(V_offset<vel_sa_data[0] || V_offset >= vel_sa_data[n_vel_sa_data-1]+dV)
        continue;
      idV = (V_offset - vel_sa_data[0])/dV;

      phase_norm[idV] += clouds_weight[i];
      alpha_cent[idV] += alpha * clouds_weight[i];
      beta_cent[idV] += beta * clouds_weight[i];
    }
  }

  /* normalize line spectrum */
  flux_norm = 0.0;
  for(j=0; j<n_vel_sa_data; j++)
  {
    flux_norm += phase_norm[j];
  }
  flux_norm /= (sa_flux_norm * n_vel_sa_data);
  for(j=0; j<n_vel_sa_data; j++)
  {
    Fline_sa[j] = FA * phase_norm[j]/(flux_norm + EPS);
  }

  for(j=0; j<n_vel_sa_data; j++)
  {
    y = alpha_cent[j];
    z = beta_cent[j];

    alpha_cent[j] = y * cos_PA + z * sin_PA;
    beta_cent[j] = -y * sin_PA + z * cos_PA;

    alpha_cent[j] = (alpha_cent[j]/(phase_norm[j]+EPS)) / DA;
    beta_cent[j] =  (beta_cent[j]/(phase_norm[j]+EPS)) / DA;
  }

  /* phi = -2*pi * f_line * B/lambda * X/DA */
  for(k=0; k<n_base_sa_data; k++)
  {
    for(j=0; j<n_vel_sa_data; j++)
    {
      phase = base_sa_data[k*2] * alpha_cent[j] + base_sa_data[k*2 + 1] * beta_cent[j];
      phase_sa[k*n_vel_sa_data + j] = -Fline_sa[j]/(1.0+Fline_sa[j]) * phase;
    }
  }

  /* smooth */
  sa_smooth_run(vel_sa_data, Fline_sa, n_vel_sa_data, phase_sa, n_base_sa_data);

  return;
}

/*
 * calculate SA phase and line profile from cloud sample at simulated mesh.
 */
void calculate_sa_sim_with_sample(const void *pm, double *vel_sa, int n_vel_sa, double *base_sa, int n_base_sa,
                                  double *p_sa, double *F_sa)
{
  int i, j, k, idV;
  double V, V_offset, dV, y, z, alpha, beta, flux_norm, phase, *phase_norm, *alpha_cent, *beta_cent;
  double DA, PA, FA, CO, cos_PA, sin_PA;
  double *pmodel = (double *)pm;

  phase_norm = workspace_phase;
  alpha_cent = phase_norm + n_vel_sa;
  beta_cent = alpha_cent + n_vel_sa;

  /* angular size distance */
  DA = exp(pmodel[num_params_blr + num_params_sa_blr_model]);
  /* position angle */
  PA = pmodel[num_params_blr + num_params_sa_blr_model + 1]/180.0 * PI;
  /* line flux scaling factor */
  FA = exp(pmodel[num_params_blr + num_params_sa_blr_model+2]);
  /*
   * line center offset: V+dV = (w-(w0+dw0))/(w0+dw0) ==> dV = - dw0/w0
   * equivalent to redshift offset:  dz = -(1+z) dw0/w0
   */
  CO = -pmodel[num_params_blr + num_params_sa_blr_model+3]/parset.sa_linecenter * C_Unit;

  /* North through East (E of N), 180 deg ambiguity */
  cos_PA = cos(PA);
  sin_PA = sin(PA);

  dV = vel_sa[1] - vel_sa[0];

  for(i=0; i<n_vel_sa; i++)
  {
    for(k=0; k<n_base_sa; k++)
    {
      p_sa[k*n_vel_sa + i] = 0.0;
    }
    phase_norm[i] = 0.0;
    alpha_cent[i] = 0.0;
    beta_cent[i] = 0.0;
  }

  for(i=0; i<parset.n_cloud_per_task; i++)
  {
    y = clouds_alpha[i];
    z = clouds_beta[i];

    //alpha = y * cos_PA + z * sin_PA;
    //beta = -y * sin_PA + z * cos_PA;
    alpha = y;
    beta = z;

    for(j=0; j< parset.n_vel_per_cloud; j++)
    {
      V = clouds_vel[i*parset.n_vel_per_cloud + j] + CO;
      V_offset = V + bin_offset * dV; /* bin type: center or left edge */
      if(V_offset<vel_sa[0] || V_offset >= vel_sa[n_vel_sa-1]+dV)
        continue;
      idV = (V_offset - vel_sa[0])/dV;

      phase_norm[idV] += clouds_weight[i];
      alpha_cent[idV] += alpha * clouds_weight[i];
      beta_cent[idV] += beta * clouds_weight[i];
    }
  }

  /* normalize line spectrum */
  flux_norm = 0.0;
  for(j=0; j<n_vel_sa; j++)
  {
    flux_norm += phase_norm[j];
  }
  flux_norm /= (sa_flux_norm * n_vel_sa);
  for(j=0; j<n_vel_sa; j++)
  {
    F_sa[j] = FA * phase_norm[j]/(flux_norm + EPS);
  }

  for(j=0; j<n_vel_sa; j++)
  {
    y = alpha_cent[j];
    z = beta_cent[j];

    alpha_cent[j] = y * cos_PA + z * sin_PA;
    beta_cent[j] = -y * sin_PA + z * cos_PA;

    alpha_cent[j] = (alpha_cent[j]/(phase_norm[j]+EPS)) / DA;
    beta_cent[j] =  (beta_cent[j]/(phase_norm[j]+EPS)) / DA;
  }

  /* phi = -2*pi * f_line * B/lambda * X/DA */
  for(k=0; k<n_base_sa; k++)
  {
    for(j=0; j<n_vel_sa; j++)
    {
      phase = base_sa[k*2] * alpha_cent[j] + base_sa[k*2 + 1] * beta_cent[j];
      p_sa[k*n_vel_sa + j] = -F_sa[j]/(1.0+F_sa[j]) * phase;
    }
  }

  /* smooth */
  sa_smooth_run(vel_sa, F_sa, n_vel_sa, p_sa, n_base_sa);

  return;
}

/*
 * calculate SA phase and line profile from a BLR model.
 */
void calculate_sa_from_blrmodel(const void *pm, int flag_save)
{
  gen_sa_cloud_sample(pm, 3, 0);
  calculate_sa_with_sample(pm);
  return;
}

void set_sa_blr_model()
{
  switch(parset.flag_sa_blrmodel)
  {
    case 0:
      num_params_sa_blr_model = num_params_MyBLRmodel2d;
      gen_sa_cloud_sample = gen_cloud_sample_mymodel;
      break;

    case 1:
      num_params_sa_blr_model = sizeof(BLRmodel1)/sizeof(double);
      gen_sa_cloud_sample = gen_cloud_sample_model1;
      break;

    case 2:
      num_params_sa_blr_model = sizeof(BLRmodel2)/sizeof(double);
      gen_sa_cloud_sample = gen_cloud_sample_model2;
      break;

    case 3:
      num_params_sa_blr_model = sizeof(BLRmodel3)/sizeof(double);
      gen_sa_cloud_sample = gen_cloud_sample_model3;
      break;

    case 4:
      num_params_sa_blr_model = sizeof(BLRmodel4)/sizeof(double);
      gen_sa_cloud_sample = gen_cloud_sample_model4;
      break;

    case 5:
      num_params_sa_blr_model = sizeof(BLRmodel5)/sizeof(double);
      gen_sa_cloud_sample = gen_cloud_sample_model5;
      break;

    case 6:
      num_params_sa_blr_model = sizeof(BLRmodel6)/sizeof(double);
      gen_sa_cloud_sample = gen_cloud_sample_model6;
      break;

    case 7:
      num_params_sa_blr_model = sizeof(BLRmodel7)/sizeof(double);
      gen_sa_cloud_sample = gen_cloud_sample_model7;
      break;

    case 8:
      num_params_sa_blr_model = sizeof(BLRmodel8)/sizeof(double);
      gen_sa_cloud_sample = gen_cloud_sample_model8;
      break;

    case 9:
      num_params_sa_blr_model = sizeof(BLRmodel9)/sizeof(double);
      gen_sa_cloud_sample = gen_cloud_sample_model9;
      break;

    default:
      num_params_sa_blr_model = sizeof(BLRmodel1)/sizeof(double);
      gen_sa_cloud_sample = gen_cloud_sample_model1;
      break;
  }

  return;
}

/*
 * set the index of mutual parameters (BH mass and inclination) in RM and SA BLR models.
 */
void set_idx_par_mutual()
{
  switch(parset.flag_blrmodel)
  {
    case 0:
      idx_rm_par_mutual[0] = offsetof(MyBLRmodel, mbh);
      idx_rm_par_mutual[1] = offsetof(MyBLRmodel, inc);
      idx_rm_par_mutual[2] = offsetof(MyBLRmodel, opn);
      break;

    case 1:
      idx_rm_par_mutual[0] = offsetof(BLRmodel1, mbh);
      idx_rm_par_mutual[1] = offsetof(BLRmodel1, inc);
      idx_rm_par_mutual[2] = offsetof(BLRmodel1, opn);
      break;

    case 2:
      idx_rm_par_mutual[0] = offsetof(BLRmodel2, mbh);
      idx_rm_par_mutual[1] = offsetof(BLRmodel2, inc);
      idx_rm_par_mutual[2] = offsetof(BLRmodel2, opn);
      break;

    case 3:
      idx_rm_par_mutual[0] = offsetof(BLRmodel3, mbh);
      idx_rm_par_mutual[1] = offsetof(BLRmodel3, inc);
      idx_rm_par_mutual[2] = offsetof(BLRmodel3, opn);
      break;

    case 4:
      idx_rm_par_mutual[0] = offsetof(BLRmodel4, mbh);
      idx_rm_par_mutual[1] = offsetof(BLRmodel4, inc);
      idx_rm_par_mutual[2] = offsetof(BLRmodel4, opn);
      break;

    case 5:
      idx_rm_par_mutual[0] = offsetof(BLRmodel5, mbh);
      idx_rm_par_mutual[1] = offsetof(BLRmodel5, inc);
      idx_rm_par_mutual[2] = offsetof(BLRmodel5, opn);
      break;

    case 6:
      idx_rm_par_mutual[0] = offsetof(BLRmodel6, mbh);
      idx_rm_par_mutual[1] = offsetof(BLRmodel6, inc);
      idx_rm_par_mutual[2] = offsetof(BLRmodel6, opn);
      break;

    case 7:
      idx_rm_par_mutual[0] = offsetof(BLRmodel7, mbh);
      idx_rm_par_mutual[1] = offsetof(BLRmodel7, inc);
      idx_rm_par_mutual[2] = offsetof(BLRmodel7, opn);
      break;

    case 8:
      idx_rm_par_mutual[0] = offsetof(BLRmodel8, mbh);
      idx_rm_par_mutual[1] = offsetof(BLRmodel8, inc);
      idx_rm_par_mutual[2] = offsetof(BLRmodel8, theta_min);
      break;

    case 9:
      idx_rm_par_mutual[0] = offsetof(BLRmodel9, mbh);
      idx_rm_par_mutual[1] = offsetof(BLRmodel9, inc);
      idx_rm_par_mutual[2] = offsetof(BLRmodel9, opn);
      break;
  }

  switch(parset.flag_sa_blrmodel)
  {
    case 0:
      idx_sa_par_mutual[0] = offsetof(MyBLRmodel, mbh);
      idx_sa_par_mutual[1] = offsetof(MyBLRmodel, inc);
      idx_sa_par_mutual[2] = offsetof(MyBLRmodel, opn);
      break;

    case 1:
      idx_sa_par_mutual[0] = offsetof(SABLRmodel1, mbh);
      idx_sa_par_mutual[1] = offsetof(SABLRmodel1, inc);
      idx_sa_par_mutual[2] = offsetof(SABLRmodel1, opn);
      break;

    case 2:
      idx_sa_par_mutual[0] = offsetof(SABLRmodel2, mbh);
      idx_sa_par_mutual[1] = offsetof(SABLRmodel2, inc);
      idx_sa_par_mutual[2] = offsetof(SABLRmodel2, opn);
      break;

    case 3:
      idx_sa_par_mutual[0] = offsetof(SABLRmodel3, mbh);
      idx_sa_par_mutual[1] = offsetof(SABLRmodel3, inc);
      idx_sa_par_mutual[2] = offsetof(SABLRmodel3, opn);
      break;

    case 4:
      idx_sa_par_mutual[0] = offsetof(SABLRmodel4, mbh);
      idx_sa_par_mutual[1] = offsetof(SABLRmodel4, inc);
      idx_sa_par_mutual[2] = offsetof(SABLRmodel4, opn);
      break;

    case 5:
      idx_sa_par_mutual[0] = offsetof(SABLRmodel5, mbh);
      idx_sa_par_mutual[1] = offsetof(SABLRmodel5, inc);
      idx_sa_par_mutual[2] = offsetof(SABLRmodel5, opn);
      break;

    case 6:
      idx_sa_par_mutual[0] = offsetof(SABLRmodel6, mbh);
      idx_sa_par_mutual[1] = offsetof(SABLRmodel6, inc);
      idx_sa_par_mutual[2] = offsetof(SABLRmodel6, opn);
      break;

    case 7:
      idx_sa_par_mutual[0] = offsetof(SABLRmodel7, mbh);
      idx_sa_par_mutual[1] = offsetof(SABLRmodel7, inc);
      idx_sa_par_mutual[2] = offsetof(SABLRmodel7, opn);
      break;

    case 8:
      idx_sa_par_mutual[0] = offsetof(SABLRmodel8, mbh);
      idx_sa_par_mutual[1] = offsetof(SABLRmodel8, inc);
      idx_sa_par_mutual[2] = offsetof(SABLRmodel8, theta_min);
      break;

    case 9:
      idx_sa_par_mutual[0] = offsetof(SABLRmodel9, mbh);
      idx_sa_par_mutual[1] = offsetof(SABLRmodel9, inc);
      idx_sa_par_mutual[2] = offsetof(SABLRmodel9, opn);
      break;
  }

  idx_sa_par_mutual[0] /= sizeof(double);
  idx_sa_par_mutual[1] /= sizeof(double);
  idx_sa_par_mutual[2] /= sizeof(double);
  idx_rm_par_mutual[0] /= sizeof(double);
  idx_rm_par_mutual[1] /= sizeof(double);
  idx_rm_par_mutual[2] /= sizeof(double);

  return;
}

/*
 * setup  SA BLR model parameter range.
 */
// model 1
void set_sa_blr_range_model1()
{
  int i;

  i = 0;
  //mu
  sa_blr_range_model[i][0] = log(0.1);
  sa_blr_range_model[i++][1] = log(rcloud_max_set*0.5);
  //beta
  sa_blr_range_model[i][0] = 0.001;
  sa_blr_range_model[i++][1] = 2.0;
  //F
  sa_blr_range_model[i][0] = 0.001;
  sa_blr_range_model[i++][1] = 0.999;
  //inc
  sa_blr_range_model[i][0] = 0.0;  // in cosine
  sa_blr_range_model[i++][1] = 1.0;
  //opn
  sa_blr_range_model[i][0] = 0.0;  // in rad
  sa_blr_range_model[i++][1] = 90.0;
  //k
  sa_blr_range_model[i][0] = -0.5;
  sa_blr_range_model[i++][1] = 0.5;
  //mbh
  sa_blr_range_model[i][0] = log(mass_range[0]);
  sa_blr_range_model[i++][1] = log(mass_range[1]);
  //lambda
  sa_blr_range_model[i][0] = 0.0;
  sa_blr_range_model[i++][1] = 1.5;
  //q
  sa_blr_range_model[i][0] = 0.0;
  sa_blr_range_model[i++][1] = 1.0;

  //sa_blr_range_model[2][1] = fmin(sa_blr_range_model[2][1], log(rcloud_max_set));
  return;
}

// model 2
void set_sa_blr_range_model2()
{
  int i;

  i = 0;
  //mu
  sa_blr_range_model[i][0] = log(0.1);
  sa_blr_range_model[i++][1] = log(rcloud_max_set*0.5);
  //beta
  sa_blr_range_model[i][0] = 0.001;
  sa_blr_range_model[i++][1] = 2.0;
  //F
  sa_blr_range_model[i][0] = 0.001;
  sa_blr_range_model[i++][1] = 0.999;
  //inc
  sa_blr_range_model[i][0] = 0.0;  // in cosine
  sa_blr_range_model[i++][1] = 1.0;
  //opn
  sa_blr_range_model[i][0] = 0.0;  // in rad
  sa_blr_range_model[i++][1] = 90.0;
  //k
  sa_blr_range_model[i][0] = -0.5;
  sa_blr_range_model[i++][1] = 0.5;
  //mbh
  sa_blr_range_model[i][0] = log(mass_range[0]);
  sa_blr_range_model[i++][1] = log(mass_range[1]);
  //sigr
  sa_blr_range_model[i][0] = 0.0;
  sa_blr_range_model[i++][1] = 1.0;
  //sigtheta
  sa_blr_range_model[i][0] = 0.0;
  sa_blr_range_model[i++][1] = 1.0;

  return;
}

// model 3
void set_sa_blr_range_model3()
{
  int i;

  i = 0;
  //Rin
  sa_blr_range_model[i][0] = log(0.1);
  sa_blr_range_model[i++][1] = log(rcloud_max_set*0.5);
  //F
  sa_blr_range_model[i][0] = log(1.0);
  sa_blr_range_model[i++][1] = log(1.0e2);
  //alpha
  sa_blr_range_model[i][0] = -3.0;
  sa_blr_range_model[i++][1] = 3.0;
  //inc
  sa_blr_range_model[i][0] = 0.0;  // in cosine
  sa_blr_range_model[i++][1] = 1.0;
  //opn
  sa_blr_range_model[i][0] = 0.0;   // in rad
  sa_blr_range_model[i++][1] = 90.0;
  //k
  sa_blr_range_model[i][0] = -0.5;
  sa_blr_range_model[i++][1] = 0.5;
  //mbh
  sa_blr_range_model[i][0] = log(mass_range[0]);
  sa_blr_range_model[i++][1] = log(mass_range[1]);
  //xi
  sa_blr_range_model[i][0] = 0.0;
  sa_blr_range_model[i++][1] = 1.0;
  //q
  sa_blr_range_model[i][0] = 0.0;
  sa_blr_range_model[i++][1] = 1.0;

  //rcloud_max_set = fmax(rcloud_max_set, exp(sa_blr_range_model[3][1] + sa_blr_range_model[4][1]));
  return;
}

// model 4
void set_sa_blr_range_model4()
{
  int i;

  i = 0;
  //Rin
  sa_blr_range_model[i][0] = log(0.1);
  sa_blr_range_model[i++][1] = log(rcloud_max_set*0.5);
  //F
  sa_blr_range_model[i][0] = log(1.0);
  sa_blr_range_model[i++][1] = log(1.0e2);
  //alpha
  sa_blr_range_model[i][0] = -3.0;
  sa_blr_range_model[i++][1] = 3.0;
  //inc
  sa_blr_range_model[i][0] = 0.0; // in cosine
  sa_blr_range_model[i++][1] = 1.0;
  //opn
  sa_blr_range_model[i][0] = 0.0;  // in rad
  sa_blr_range_model[i++][1] = 90.0;
  //k
  sa_blr_range_model[i][0] = -0.5;
  sa_blr_range_model[i++][1] = 0.5;
  //mbh
  sa_blr_range_model[i][0] = log(mass_range[0]);
  sa_blr_range_model[i++][1] = log(mass_range[1]);
  //xi
  sa_blr_range_model[i][0] = 0.0;
  sa_blr_range_model[i++][1] = sqrt(2.0)/2.0;
  //q
  sa_blr_range_model[i][0] = 0.0;
  sa_blr_range_model[i++][1] = 1.0;

  //rcloud_max_set = fmax(rcloud_max_set, exp(sa_blr_range_model[3][1] + sa_blr_range_model[4][1]));
  return;
}

// model 5
void set_sa_blr_range_model5()
{
  int i;

  i = 0;
  //mu
  sa_blr_range_model[i][0] = log(0.1);
  sa_blr_range_model[i++][1] = log(rcloud_max_set*0.5);
  //Fin
  sa_blr_range_model[i][0] = 0.0;
  sa_blr_range_model[i++][1] = 1.0;
  //Fout
  sa_blr_range_model[i][0] = log(1.0);
  sa_blr_range_model[i++][1] = log(10.0);
  //alpha
  sa_blr_range_model[i][0] = 1.0;
  sa_blr_range_model[i++][1] = 3.0;
  //inc
  sa_blr_range_model[i][0] = 0.0; // in cosine
  sa_blr_range_model[i++][1] = 1.0;
  //opn
  sa_blr_range_model[i][0] = 0.0;  // in degree
  sa_blr_range_model[i++][1] = 90.0;
  //k
  sa_blr_range_model[i][0] = -0.5;
  sa_blr_range_model[i++][1] = 0.5;
  //beta
  sa_blr_range_model[i][0] = 1.0;
  sa_blr_range_model[i++][1] = 5.0;
  //xi
  sa_blr_range_model[i][0] = 0.0;
  sa_blr_range_model[i++][1] = 1.0;
  //mbh
  sa_blr_range_model[i][0] = log(mass_range[0]);
  sa_blr_range_model[i++][1] = log(mass_range[1]);
  //fellip
  sa_blr_range_model[i][0] = 0.0;
  sa_blr_range_model[i++][1] = 1.0;
  //fflow
  sa_blr_range_model[i][0] = 0.0;
  sa_blr_range_model[i++][1] = 1.0;
  //sigr_circ
  sa_blr_range_model[i][0] = log(0.001);
  sa_blr_range_model[i++][1] = log(0.1);
  //sigthe_circ
  sa_blr_range_model[i][0] = log(0.001);
  sa_blr_range_model[i++][1] = log(1.0);
  //sigr_rad
  sa_blr_range_model[i][0] = log(0.001);
  sa_blr_range_model[i++][1] = log(0.1);
  //sigthe_rad
  sa_blr_range_model[i][0] = log(0.001);
  sa_blr_range_model[i++][1] = log(1.0);
  //theta_rot
  sa_blr_range_model[i][0] = 0.0;
  sa_blr_range_model[i++][1] = 90.0;
  //sig_turb
  sa_blr_range_model[i][0] = log(0.001);
  sa_blr_range_model[i++][1] = log(0.1);

  //rcloud_max_set = fmax(rcloud_max_set, exp(sa_blr_range_model[2][1] + sa_blr_range_model[4][1]));

  return;
}

// model 6
void set_sa_blr_range_model6()
{
  int i;

  i = 0;
  //mu
  sa_blr_range_model[i][0] = log(0.1);
  sa_blr_range_model[i++][1] = log(rcloud_max_set*0.5);
  //beta
  sa_blr_range_model[i][0] = 0.001;
  sa_blr_range_model[i++][1] = 2.0;
  //F
  sa_blr_range_model[i][0] = 0.0;
  sa_blr_range_model[i++][1] = 1.0;
  //inc
  sa_blr_range_model[i][0] = 0.0; //in cosine
  sa_blr_range_model[i++][1] = 1.0;
  //opn
  sa_blr_range_model[i][0] = 0.0;  // in rad
  sa_blr_range_model[i++][1] = 90.0;
  //k
  sa_blr_range_model[i][0] = -0.5;
  sa_blr_range_model[i++][1] = 0.5;
  //gamma
  sa_blr_range_model[i][0] = 1.0;
  sa_blr_range_model[i++][1] = 5.0;
  //xi
  sa_blr_range_model[i][0] = 0.0;
  sa_blr_range_model[i++][1] = 1.0;
  //mbh
  sa_blr_range_model[i][0] = log(mass_range[0]);
  sa_blr_range_model[i++][1] = log(mass_range[1]);
  //fellip
  sa_blr_range_model[i][0] = 0.0;
  sa_blr_range_model[i++][1] = 1.0;
  //fflow
  sa_blr_range_model[i][0] = 0.0;
  sa_blr_range_model[i++][1] = 1.0;
  //sigr_circ
  sa_blr_range_model[i][0] = log(0.001);
  sa_blr_range_model[i++][1] = log(0.1);
  //sigthe_circ
  sa_blr_range_model[i][0] = log(0.001);
  sa_blr_range_model[i++][1] = log(1.0);
  //sigr_rad
  sa_blr_range_model[i][0] = log(0.001);
  sa_blr_range_model[i++][1] = log(0.1);
  //sigthe_rad
  sa_blr_range_model[i][0] = log(0.001);
  sa_blr_range_model[i++][1] = log(1.0);
  //theta_rot
  sa_blr_range_model[i][0] = 0.0;
  sa_blr_range_model[i++][1] = 90.0;
  //sig_turb
  sa_blr_range_model[i][0] = log(0.0001);
  sa_blr_range_model[i++][1] = log(0.1);

  //rcloud_max_set = 44.85;
  //sa_blr_range_model[2][1] = fmin(sa_blr_range_model[2][1], log(rcloud_max_set));
  return;
}

// model 7
void set_sa_blr_range_model7()
{
  int i;

  i = 0;
  //mu
  sa_blr_range_model[i][0] = log(0.1);
  sa_blr_range_model[i++][1] = log(rcloud_max_set*0.5);
  //beta
  sa_blr_range_model[i][0] = 0.001;
  sa_blr_range_model[i++][1] = 2.0;
  //F
  sa_blr_range_model[i][0] = 0.0;
  sa_blr_range_model[i++][1] = 1.0;
  //inc
  sa_blr_range_model[i][0] = 0.0;  // in cosine
  sa_blr_range_model[i++][1] = 1.0;
  //opn
  sa_blr_range_model[i][0] = 0.0;  // in rad
  sa_blr_range_model[i++][1] = 90.0;
  //k
  sa_blr_range_model[i][0] = -0.5;
  sa_blr_range_model[i++][1] = 0.5;
  //gamma
  sa_blr_range_model[i][0] = 1.0;
  sa_blr_range_model[i++][1] = 5.0;
  //xi
  sa_blr_range_model[i][0] = 0.0;
  sa_blr_range_model[i++][1] = 1.0;

  //fsh
  sa_blr_range_model[i][0] = 0.0;
  sa_blr_range_model[i++][1] = 1.0;
  //mu_un
  sa_blr_range_model[i][0] = log(0.1);
  sa_blr_range_model[i++][1] = log(rcloud_max_set*0.5);
  //beta_un
  sa_blr_range_model[i][0] = 0.001;
  sa_blr_range_model[i++][1] = 2.0;
  //F_un
  sa_blr_range_model[i][0] = 0.0;
  sa_blr_range_model[i++][1] = 1.0;
  //opn_un
  sa_blr_range_model[i][0] = 0.0;
  sa_blr_range_model[i++][1] = 90.0;

  //mbh
  sa_blr_range_model[i][0] = log(mass_range[0]);
  sa_blr_range_model[i++][1] = log(mass_range[1]);
  //fellip
  sa_blr_range_model[i][0] = 0.0;
  sa_blr_range_model[i++][1] = 1.0;
  //fflow
  sa_blr_range_model[i][0] = 0.0;
  sa_blr_range_model[i++][1] = 1.0;
  //sigr_circ
  sa_blr_range_model[i][0] = log(0.001);
  sa_blr_range_model[i++][1] = log(0.1);
  //sigthe_circ
  sa_blr_range_model[i][0] = log(0.001);
  sa_blr_range_model[i++][1] = log(1.0);
  //sigr_rad
  sa_blr_range_model[i][0] = log(0.001);
  sa_blr_range_model[i++][1] = log(0.1);
  //sigthe_rad
  sa_blr_range_model[i][0] = log(0.001);
  sa_blr_range_model[i++][1] = log(1.0);
  //theta_rot
  sa_blr_range_model[i][0] = 0.0;
  sa_blr_range_model[i++][1] = 90.0;

  //fellip_un
  sa_blr_range_model[i][0] = 0.0;
  sa_blr_range_model[i++][1] = 1.0;
  //fflow_un
  sa_blr_range_model[i][0] = 0.0;
  sa_blr_range_model[i++][1] = 1.0;

  //sig_turb
  sa_blr_range_model[i][0] = log(0.001);
  sa_blr_range_model[i++][1] = log(0.1);

  return;
}

// model 8
void set_sa_blr_range_model8()
{
  int i;

  i = 0;
  //theta_min
  sa_blr_range_model[i][0] = 20.0;     // in degree
  sa_blr_range_model[i++][1] = 90.0;
  //dtheta_max
  sa_blr_range_model[i][0] = 0.0;      // in degree
  sa_blr_range_model[i++][1] = 90.0;
  //r_min
  sa_blr_range_model[i][0] = log(0.1);
  sa_blr_range_model[i++][1] = log(rcloud_max_set*0.25);
  //fr_max
  sa_blr_range_model[i][0] = 1.0;  //
  sa_blr_range_model[i++][1] = 10.0;
  //gamma
  sa_blr_range_model[i][0] = 0.0;  //
  sa_blr_range_model[i++][1] = 3.0;
  //alpha
  sa_blr_range_model[i][0] = 0.0;
  sa_blr_range_model[i++][1] = 3.0;
  //lambda
  sa_blr_range_model[i][0] = -3.0;
  sa_blr_range_model[i++][1] = 0.0;
  //k
  sa_blr_range_model[i][0] = -0.5;
  sa_blr_range_model[i++][1] = 0.5;
  //xi
  sa_blr_range_model[i][0] = 0.0;
  sa_blr_range_model[i++][1] = 1.0;
  //Rv
  sa_blr_range_model[i][0] = log(10.0);
  sa_blr_range_model[i++][1] = log(50.0);
  //Rblr
  sa_blr_range_model[i][0] = log(10.0);
  sa_blr_range_model[i++][1] = log(rcloud_max_set*0.5);
  //inc
  sa_blr_range_model[i][0] = 0.0;   // cos(inc)
  sa_blr_range_model[i++][1] = 1.0;
  //mbh
  sa_blr_range_model[i][0] = log(mass_range[0]);
  sa_blr_range_model[i++][1] = log(mass_range[1]);

  return;
}

// model 9
void set_sa_blr_range_model9()
{
  int i;

  i = 0;
  //mu
  sa_blr_range_model[i][0] = log(0.1);
  sa_blr_range_model[i++][1] = log(rcloud_max_set*0.5);
  //beta
  sa_blr_range_model[i][0] = 0.001;
  sa_blr_range_model[i++][1] = 2.0;
  //F
  sa_blr_range_model[i][0] = 0.001;
  sa_blr_range_model[i++][1] = 0.999;
  //inc
  sa_blr_range_model[i][0] = 0.0;  // in cosine
  sa_blr_range_model[i++][1] = 1.0;
  //opn
  sa_blr_range_model[i][0] = 0.0;  // in degree
  sa_blr_range_model[i++][1] = 90.0;
  //mbh
  sa_blr_range_model[i][0] = log(mass_range[0]);
  sa_blr_range_model[i++][1] = log(mass_range[1]);

  return;
}

/*!
 * This function copes with parameter fixing.\n
 * Only fix BLR model parameters.
 */
void set_par_fix_sa_blrmodel()
{
  int i;
  char *pstr;

  npar_fix = 0;

  if(thistask == roottask)
  {
    pstr = parset.sa_str_par_fix_val;
    // set the default value if not provided.
    for(i=strlen(parset.sa_str_par_fix); i<num_params_sa_blr_model; i++)
      parset.sa_str_par_fix[i] = '0';

    for(i=0; i<num_params_sa_blr_model; i++)
    {
      if(parset.sa_str_par_fix[i] == '0')
      {
        par_fix[num_params_blr + i] = 0;
        par_fix_val[num_params_blr + i] = -DBL_MAX;  /* set to be the smallest value */
      }
      else if(parset.sa_str_par_fix[i] == '1')
      {
        if(pstr == NULL)
        {
          printf("# %d-th SA BLR parameter value is not provided (counting from 0).\n", i);
          exit(0);
        }
        par_fix[num_params_blr + i] = 1;
        sscanf(pstr, "%lf", &par_fix_val[num_params_blr + i]);
        npar_fix++;
        printf("# %d-th SA BLR parameter fixed, value= %f.\n", i, par_fix_val[num_params_blr + i]);
        pstr = strchr(pstr, ':'); /* values are separated by ":" */
        if(pstr!=NULL)
        {
          pstr++;
        }
      }
      else   // default value
      {
        par_fix[num_params_blr + i] = 0;
        par_fix_val[num_params_blr + i] = -DBL_MAX;
      }
    }
  }

  MPI_Bcast(par_fix + num_params_blr, num_params_sa_blr_model, MPI_INT, roottask, MPI_COMM_WORLD);
  MPI_Bcast(par_fix_val + num_params_blr, num_params_sa_blr_model, MPI_DOUBLE, roottask, MPI_COMM_WORLD);

  return;
}

int nd_fft_sa, nd_fft_sa_cal, npad_sa;

fftw_complex *sa_data_fft, *sa_resp_fft, *sa_conv_fft;
fftw_plan sa_pdata, sa_presp, sa_pback;
double *sa_real_data, *sa_real_resp, *sa_real_conv;

/*!
 * This function initiates workspace for FFT.
 */
void sa_smooth_init(int n_v_sa, const double *v_sa, double sigV)
{
  int i;
  double dV;

  npad_sa = fmin(n_v_sa * 0.1, 20);
  nd_fft_sa = n_v_sa+npad_sa*2;
  nd_fft_sa_cal = nd_fft_sa/2+1;

  sa_data_fft = (fftw_complex *) fftw_malloc((nd_fft_sa_cal) * sizeof(fftw_complex));
  sa_resp_fft = (fftw_complex *) fftw_malloc((nd_fft_sa_cal) * sizeof(fftw_complex));
  sa_conv_fft = (fftw_complex *) fftw_malloc((nd_fft_sa_cal) * sizeof(fftw_complex));

  sa_real_data = (double *)fftw_malloc(nd_fft_sa * sizeof(double));
  //real_resp = (double *)fftw_malloc(nd_fft_sa * sizeof(double));
  sa_real_conv = (double *)fftw_malloc(nd_fft_sa * sizeof(double));

  sa_pdata = fftw_plan_dft_r2c_1d(nd_fft_sa, sa_real_data, sa_data_fft, FFTW_PATIENT);
  //presp = fftw_plan_dft_r2c_1d(nd_fft_sa, sa_real_resp, sa_resp_fft, FFTW_PATIENT);
  sa_pback = fftw_plan_dft_c2r_1d(nd_fft_sa, sa_conv_fft, sa_real_conv, FFTW_PATIENT);

  dV = v_sa[1] - v_sa[0];
  for(i=0; i<nd_fft_sa_cal; i++)
  {
    sa_resp_fft[i][0] = exp(-2.0 * PI*PI * sigV/dV*sigV/dV * i*i*1.0/nd_fft_sa/nd_fft_sa)/nd_fft_sa;
    sa_resp_fft[i][1] = 0.0;
  }
  return;
}

/*!
 * This function finalizes FFT.
 */
void sa_smooth_end()
{
  fftw_destroy_plan(sa_pdata);
  fftw_destroy_plan(sa_pback);

  fftw_free(sa_data_fft);
  fftw_free(sa_resp_fft);
  fftw_free(sa_conv_fft);

  fftw_free(sa_real_data);
  //fftw_free(sa_real_resp);
  fftw_free(sa_real_conv);
  return;
}

void sa_smooth_run(double *v_sa, double *F_sa, int n_v_sa, double *p_sa, int n_base_sa)
{
  int i, j;
  /* first profile */
  memcpy(sa_real_data+npad_sa, F_sa, n_v_sa*sizeof(double));
  for(i=0; i<npad_sa; i++)
    sa_real_data[i] = sa_real_data[nd_fft_sa-1-i] = 0.0;

  /* FFT of line */
  fftw_execute(sa_pdata);
  /* complex multiply and inverse FFT
   * note that for FFT of real data, FFTW outputs n/2+1 complex numbers.
   * similarly, for complex to real transform, FFTW needs input of n/2+1 complex numbers.
   */
  for(i=0; i<nd_fft_sa_cal; i++)
  {
    sa_conv_fft[i][0] = sa_data_fft[i][0]*sa_resp_fft[i][0] - sa_data_fft[i][1]*sa_resp_fft[i][1];
    sa_conv_fft[i][1] = sa_data_fft[i][0]*sa_resp_fft[i][1] + sa_data_fft[i][1]*sa_resp_fft[i][0];
  }
  fftw_execute(sa_pback);
  memcpy(F_sa, sa_real_conv+npad_sa, n_v_sa*sizeof(double));

  /* then phase */
  for(j=0; j<n_base_sa; j++)
  {
    memcpy(sa_real_data+npad_sa, &p_sa[j*n_v_sa], n_v_sa*sizeof(double));
    for(i=0; i<npad_sa; i++)
      sa_real_data[i] = sa_real_data[nd_fft_sa-1-i] = 0.0;

    /* FFT of line */
    fftw_execute(sa_pdata);
    /* complex multiply and inverse FFT
     * note that for FFT of real data, FFTW outputs n/2+1 complex numbers.
     * similarly, for complex to real transform, FFTW needs input of n/2+1 complex numbers.
     */
    for(i=0; i<nd_fft_sa_cal; i++)
    {
      sa_conv_fft[i][0] = sa_data_fft[i][0]*sa_resp_fft[i][0] - sa_data_fft[i][1]*sa_resp_fft[i][1];
      sa_conv_fft[i][1] = sa_data_fft[i][0]*sa_resp_fft[i][1] + sa_data_fft[i][1]*sa_resp_fft[i][0];
    }
    fftw_execute(sa_pback);
    memcpy(&p_sa[j*n_v_sa], sa_real_conv+npad_sa, n_v_sa*sizeof(double));
  }
  return;
}
#endif

max_stars_repo_name: LiyrAstroph/BRAINS

max_issues_repo_name: LiyrAstroph/BRAINS

max_forks_repo_name: LiyrAstroph/BRAINS


text:
/*#include <fftw.h>*/
#include <fftw3.h>
#include <stdio.h>
#include <stdlib.h>
#define N 16
main()
{
     int i;
     fftw_complex in[N], out[N];
     float fact;
     fftw_plan p,p2;
printf(" stuff in data\n");
     for (i=0;i<N;i++) {
/*      in[i].re=(i*i); in[i].im=1; */
     	in[i][0]=(i*i); in[i][1]=1;
     }
printf(" create plans\n");
     p  = (fftw_plan)fftw_plan_dft_1d(N,in,out,FFTW_FORWARD,FFTW_ESTIMATE);
     p2 = (fftw_plan)fftw_plan_dft_1d(N,out,in,FFTW_BACKWARD,FFTW_ESTIMATE);
printf(" do it\n");
     fftw_execute(p);
     for (i=0;i<N;i++) {
/*	printf("%12.4f %12.4f\n",out[i].re,out[i].im); */
  	printf("%12.4f %12.4f\n",out[i][0],out[i][1]);
     }
printf(" \n");
printf(" undo it\n");
     fftw_execute(p2);
     fact=1.0/(float)N;
     for (i=0;i<N;i++) {
/*	printf("%10.2f %10.2f\n",in[i].re*fact,in[i].im*fact); */
	printf("%10.2f %10.2f\n",in[i][0]*fact,in[i][1]*fact);
     }
printf(" clean up\n");
     fftw_destroy_plan(p);
     fftw_destroy_plan(p2);
}

max_stars_repo_name: timkphd/examples

max_issues_repo_name: timkphd/examples

max_forks_repo_name: timkphd/examples


text:
#include <viaio/Vlib.h>
#include <viaio/VImage.h>
#include <viaio/mu.h>
#include <viaio/option.h>

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_cblas.h>
#include <gsl/gsl_blas.h>

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

int DNorm(double *data,size_t n)
{
  size_t j;
  double s1=0,s2=0,mean,sd,u,tiny=1.0e-6;
  double nx = (double)n;

  for (j=0; j<n; j++) {
    u = data[j];
    s1 += u;
    s2 += u*u;
  }
  mean = s1/nx;
  sd = sqrt((s2 - nx * mean * mean) / (nx - 1.0));
  if (sd < tiny) return -1;
  for (j=0; j<n; j++) {
    u = data[j];
    data[j] = (u-mean)/sd;
  }
  return 1;
}

double DCorr(double *data1,double *data2,size_t n)
{
  size_t j;
  double sum=0,nx=(double)n;
  for (j=0; j<n; j++) {
    sum += data1[j]*data2[j];
  }
  return sum/nx;
}

void DMN(VImage map,float *A,size_t nvox)
{
  size_t i,j,k;
  int b,r,c;

  fprintf(stderr," DMN...\n");

  VAttrList list = VReadAttrList("dmn.v",0L,TRUE,FALSE);
  VImage dmn = VReadImage(list);
  if (dmn == NULL) VError(" no DMN found");
  VImage dst = VCreateImageLike(dmn);
  VFillImage(dst,VAllBands,0);

  /*
  size_t i0=0;
  for (i=0; i<nvox; i++) {
    b = VPixel(map,0,0,i,VInteger);
    r = VPixel(map,0,1,i,VInteger);
    c = VPixel(map,0,2,i,VInteger);
    if (c==31 && r==25 && b==33) i0 = i;
    if (c==20 && r==18 && b==32) i0 = i;
  }
  i = i0;

  for (j=0; j<nvox; j++) {
    b = VPixel(map,0,0,j,VInteger);
    r = VPixel(map,0,1,j,VInteger);
    c = VPixel(map,0,2,j,VInteger);
    if (j < i) k = j+i*(i+1)/2;
    else k=i+j*(j+1)/2;
    VPixel(dst,b,r,c,VFloat) = A[k];
  }
  */

  double *data1 = (double *)VCalloc(nvox,sizeof(double));
  double *data2 = (double *)VCalloc(nvox,sizeof(double));

  for (j=0; j<nvox; j++) {
    b = VPixel(map,0,0,j,VInteger);
    r = VPixel(map,0,1,j,VInteger);
    c = VPixel(map,0,2,j,VInteger);
    data1[j] = VGetPixel(dmn,b,r,c);
  }
  DNorm(data1,nvox);

  for (i=0; i<nvox; i++) {
    if (i%100 == 0) fprintf(stderr," i: %7lu  of %lu\r",i,nvox);
    k=0;
    for (j=0; j<nvox; j++) {
      b = VPixel(map,0,0,j,VInteger);
      r = VPixel(map,0,1,j,VInteger);
      c = VPixel(map,0,2,j,VInteger);

      if (j < i) k = j+i*(i+1)/2;
      else k=i+j*(j+1)/2;
      data2[j] = A[k];
    }
    if (DNorm(data2,nvox) < 0) continue;

    b = VPixel(map,0,0,i,VInteger);
    r = VPixel(map,0,1,i,VInteger);
    c = VPixel(map,0,2,i,VInteger);
    VPixel(dst,b,r,c,VFloat) = DCorr(data1,data2,nvox);
  }
  VFree(data1);
  VFree(data2);

  VAttrList out_list = VCreateAttrList();
  VAttrList geolist = VGetGeoInfo(list);
  VSetGeoInfo(geolist,out_list);
  VAppendAttr(out_list,"image",NULL,VImageRepn,dst);
  FILE *fp = fopen("dst.v","w");
  VWriteFile (fp, out_list);
  fclose(fp);
  fprintf(stderr,"\n DMN, ok.\n");
}

max_stars_repo_name: zrajna/lipsia

max_issues_repo_name: zrajna/lipsia

max_forks_repo_name: zrajna/lipsia