text:
// mimetic_bayes program
// Alessandro Casalino, University of Trento
// For licence informations, see the Github repository license
//
// Compile with: gcc-8 -O2 -DHAVE_INLINE -DGSL_RANGE_CHECK_OFF Ligo_abc_gsl.c -o Ligo_abc_gsl.o -lgsl -lgslcblas -fopenmp
// Run with: ./Ligo_abc_gsl.o
// Test running time: time ./Ligo_abc_gsl.o

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <sys/time.h>
#include <omp.h>

// PROGRAM VARIABLES

// Variables for parallelization with OpenMP
// Info on OpenMP: https://bisqwit.iki.fi/story/howto/openmp/#Syntax
// Enable OpenMP
int OPENMP = 1;
// Number of threads used by OpenMP (can't be bigger than the number of threads available on your computer)
// If 0, the number of threads available on the computer is used
int OMP_THREAD_NUM = 0;

// Number of Chains (for higher efficiency, consider N_CHAIN as a multiple of OMP_THREAD_NUM)
int N_CHAINS = 8;

// Make the Lag test (can make the computation very long!)
int LAG_TEST = 0;
// Ratio of lag function computed with respect to the total number of points in every chain
// Smaller values make the computation faster
double LAG_TEST_STOP_RATIO = 0.001;

// Write (in the terminal) debug informations and informations that might be considered in order to improve results
int DEBUG = 1;

// Write .csv files
#define CSV 1
// Write covariant matrix in .csv file
int CSV_COV = 1;
// Write correlation matrix (Pearson coefficients) in .csv file
int CSV_CORR = 1;

// NUMERICAL VARIABLES
// Chain evaluation dimension (every step found from proposal distribution is multiplied by this quantity)
// This is also used to define the covariance matrix for the Gaussian proposal
double SPEED[] = {1e-1,5e-16,4e-16};
// Mean value for Gaussian proposal (PROPOSAL=1)
double MEAN_P1[] = {0.,0.,0.};
// Number of points in every chain
int N = 1e6;

// PHYSICAL VARIABLES
// Choose to use the constraint on Ct2 in the prior
int CT2_CONSTR = 1;
// Choose the proposal distribution
int PROPOSAL = 1;

// Note that the constraints are defined in rho(pos)

// PROGRAM

typedef struct chain {
  gsl_vector * mean;
  gsl_matrix * covariance;
  double acceptance_ratio;
} chain_res;

inline double cT2 (gsl_vector * pos) {

  double a=gsl_vector_get(pos,0), b=gsl_vector_get(pos,1), c=gsl_vector_get(pos,2);

  return (2.0-2.0*a)/(2.0-2.0*a+b);

}

inline double cS2 (gsl_vector * pos) {

  double a=gsl_vector_get(pos,0), b=gsl_vector_get(pos,1), c=gsl_vector_get(pos,2);

  return (b-c)*(2.0*a-2.0)/(2.0*a-b-2.0)/(4.0-4.0*a-b+3.0*c);

}

inline double MPl2 (gsl_vector * pos) {

  double a=gsl_vector_get(pos,0), b=gsl_vector_get(pos,1), c=gsl_vector_get(pos,2);

  return 4.0-4.0*a-b+3.0*c;

}

inline double rho(gsl_vector * pos){

  double a=gsl_vector_get(pos,0), b=gsl_vector_get(pos,1), c=gsl_vector_get(pos,2);

  // Value of the sound speed squared
  double cs2=cS2(pos);
  // Value of the tensor speed squared
  double ct2=cT2(pos);
  // Value of the Planck mass rescaling factor
  double Mpl2=MPl2(pos);

  // Value of the argument of the Gaussian in the Likelihood
  double delta=fabs(sqrt(ct2)-1.0);

  // This is the experimental sigma on cT2
  double sigma = 4.5e-16;

  double result = 0.;

  // With these ifs I'm introducing the constraints on cT2, cS2 and on the parameters a and c
  // To eliminate these constraints, put CT2_CONSTR=0 at the beginning of the file
  //
  // OLD constraint: CT2_CONSTR==1 && a>=-1. && a<=1. && b>=0. && b<=1. && ct2>=0.0 && ct2<= 1.0 && cs2>=0.0 && cs2<=1.0
  //
  // NEW (correct) constraint: CT2_CONSTR==1 && a>=-1. && a<=1. && Mpl2>=0. && c<=1. && c>=0 && ct2>=0.0 && ct2<= 1.0 && cs2<=0.0
  //
  if (CT2_CONSTR==1 && a>=-1. && a<=1. && c>=0. && c<=1. && ct2>=0. && ct2<=1. && cs2>=0. && cs2<=1.){
    result=exp(-0.5 * pow(delta/sigma,2.0) );
  }
  else if(CT2_CONSTR==1){
    result=exp(-1e90);
  }
  else{
    result=exp(-0.5 * pow(delta/sigma,2.0) );
  }

  return result;
}

chain_res chain_analysis(double (* chain)[3]){

  chain_res result;

  gsl_vector * mean = gsl_vector_alloc(3);

  int i,k,s;
  for(k=0;k<3;k++){
    double mean_temp = 0.;
    for(i=0;i<N;i++){
      mean_temp = mean_temp + chain[i][k];
    }
    gsl_vector_set(mean, k, mean_temp);
  }
  gsl_vector_scale(mean,1./N);

  gsl_matrix * cm = gsl_matrix_alloc(3,3);

  for(k=0;k<3;k++){
    double mean_row = gsl_vector_get(mean,k);
    for(s=0;s<3;s++){
      double mean_column = gsl_vector_get(mean,s);
      double cm_temp = 0.;
      for(i=0;i<N;i++){
        cm_temp = cm_temp + (chain[i][k]-mean_row)*(chain[i][s]-mean_column);
      }
      gsl_matrix_set(cm, k, s, cm_temp);
    }
  }
  gsl_matrix_scale(cm,1./(N-1.));

  result.mean = gsl_vector_alloc(3);
  gsl_vector_memcpy(result.mean,mean);
  result.covariance = gsl_matrix_alloc(3,3);
  gsl_matrix_memcpy(result.covariance,cm);

  gsl_matrix_free(cm);
  gsl_vector_free(mean);

  return result;

}

inline double Gaussian_proposal(double mean, double sigma, double u1, double u2) {

  // Box Muller transform from flat distribution to Guassian distribution
  return mean + sigma * sqrt(-2.*log(u1)) * cos(2.*M_PI*u2);

}

void C_Delta(double (* chain)[3], chain_res results, double f_stop, int chain_number){

  int id_thread = omp_get_thread_num();

  double mean_X = gsl_vector_get(results.mean,0);
  double mean_Y = gsl_vector_get(results.mean,1);
  double mean_Z = gsl_vector_get(results.mean,2);

  double cm_XX = gsl_matrix_get(results.covariance,0,0);
  double cm_YY = gsl_matrix_get(results.covariance,1,1);
  double cm_ZZ = gsl_matrix_get(results.covariance,2,2);

  int stop = N*f_stop;

  double * C_Delta_X; double * C_Delta_Y; double * C_Delta_Z;
  C_Delta_X = (double *) malloc(sizeof(double) * stop);
  C_Delta_Y = (double *) malloc(sizeof(double) * stop);
  C_Delta_Z = (double *) malloc(sizeof(double) * stop);

  int i,j;
  for(i=0;i<stop;i++){
    for(j=0;j<N-i;j++){
      C_Delta_X[i]=C_Delta_X[i]+(chain[j][0]-mean_X)*(chain[j+i][0]-mean_X);
      C_Delta_Y[i]=C_Delta_Y[i]+(chain[j][1]-mean_Y)*(chain[j+i][1]-mean_Y);
      C_Delta_Z[i]=C_Delta_Z[i]+(chain[j][2]-mean_Z)*(chain[j+i][2]-mean_Z);
    }
    C_Delta_X[i]=C_Delta_X[i]/(N-i)/cm_XX;
    C_Delta_Y[i]=C_Delta_Y[i]/(N-i)/cm_YY;
    C_Delta_Z[i]=C_Delta_Z[i]/(N-i)/cm_ZZ;
  }

  FILE *fp;
  char filename[25];

  sprintf(filename,"Ligo_abc_DeltaC_t%d.csv",chain_number);
  fp = fopen (filename, "w+");

  for(i=0;i<stop;i++){
    fprintf(fp, "%d,%.8e,%.8e,%.8e\n", i, C_Delta_X[i], C_Delta_Y[i], C_Delta_Z[i]);
  }

  fclose(fp);

  free(C_Delta_X);
  free(C_Delta_Y);
  free(C_Delta_Z);

}

chain_res LHSampling(int proposal, int chain_number){

  int k;

  int id_thread=omp_get_thread_num();

  const gsl_rng_type * T;
  gsl_rng * r;

  gsl_rng_env_setup();

  T = gsl_rng_default;
  r = gsl_rng_alloc (T);

  unsigned long seed = time(NULL) + chain_number;

  gsl_rng_set(r, seed);

  gsl_vector * pos_ini = gsl_vector_alloc(3);
  for(k=0;k<3;k++){
    if(k==0){
      double random_number = 2. * gsl_rng_uniform (r) - 1.;
      gsl_vector_set(pos_ini,k,SPEED[k]*random_number);
    }
    else if (k==1){
      double random_number = gsl_rng_uniform (r);
      gsl_vector_set(pos_ini,k,SPEED[k]*random_number);
    }
    else {
      double random_number = gsl_rng_uniform (r);
      gsl_vector_set(pos_ini,k,SPEED[k]*random_number);
    }
  }

  // Here definitions only for the Gaussian proposal distribution
  gsl_vector * pos_mean_p1 = gsl_vector_alloc(3);
  for(k=0;k<3;k++){
    gsl_vector_set(pos_mean_p1, k, MEAN_P1[k]);
  }

  gsl_matrix * S = gsl_matrix_calloc(3,3);
  gsl_matrix * L = gsl_matrix_alloc(3,3);
  gsl_vector * pos_temp2 = gsl_vector_alloc(3);

  gsl_matrix_set(S,0,0,pow(SPEED[0],2.));
  gsl_matrix_set(S,1,1,pow(SPEED[1],2.));
  gsl_matrix_set(S,2,2,pow(SPEED[2],2.));
  gsl_matrix_memcpy(L,S);
  gsl_linalg_cholesky_decomp1(L);
  // End of Gaussian proposal definitions

  gsl_vector * pos = gsl_vector_alloc(3);
  gsl_vector_memcpy(pos,pos_ini);

  #if CSV==1
    FILE *fp;
    char filename[25];
    sprintf(filename,"Ligo_abc_%d.txt",chain_number);
    fp = fopen (filename, "w+");
  #endif

  double (* chain)[3] = malloc(sizeof(* chain) * N);

  int i = 0; int j = 0; int rc = 1;

  while (i<N) {

    double f0 = rho(pos);

    gsl_vector * pos_temp = gsl_vector_calloc(3);

    // NOTE: Introducing the proposal distributions here

    if(proposal==0){

      for(k=0;k<3;k++){

        double random_number = 2. * gsl_rng_uniform (r) - 1.;
        gsl_vector_set(pos_temp, k, SPEED[k]*random_number);

      }

      gsl_vector_add(pos_temp,pos);

    }

    else if(proposal=!0){

      for(k=0;k<3;k++){

        gsl_vector_set(pos_temp2, k, Gaussian_proposal(0., 1., gsl_rng_uniform (r), gsl_rng_uniform (r)));

      }

      gsl_blas_dgemv (CblasNoTrans,
                  1., L, pos_temp2,
                  0., pos_temp);
      gsl_vector_add(pos_temp,pos_mean_p1);
      gsl_vector_add(pos_temp,pos);

    }

    double f1 = rho(pos_temp);

    // NOTE: here we also print data in .csv file (if CSV==1) in Getdist ready format
    // Getdist output with (no commas): repetitions count log(rho) {parameters} {(optional) other outputs - not considered by Getdist}

    if (f1>f0) {

      gsl_vector_memcpy(pos,pos_temp);
      j++;

      #if CSV==1
      fprintf(fp, "%d %.8e %.8e %.8e %.8e\n", rc, -log(rho(pos)), gsl_vector_get(pos,0), gsl_vector_get(pos,1), gsl_vector_get(pos,2));
      #endif

      rc = 1;

    }
    else {

      double ra = gsl_rng_uniform (r);

      if (ra<f1/f0) {
        gsl_vector_memcpy(pos,pos_temp);
        j++;

        #if CSV==1
        fprintf(fp, "%d %.8e %.8e %.8e %.8e\n", rc, -log(rho(pos)), gsl_vector_get(pos,0), gsl_vector_get(pos,1), gsl_vector_get(pos,2));
        #endif

        rc = 1;
      }
      else{
        rc++;
      }

    }

    for(k=0;k<3;k++){
      chain[i][k]=gsl_vector_get(pos,k);
    }

    gsl_vector_free(pos_temp);
    i++;

  }

  chain_res result;

  result = chain_analysis(chain);

  result.acceptance_ratio = (double) j/i;

  if(DEBUG==1 && result.acceptance_ratio<0.05) printf("\n WARNING (chain %d): Acceptance ratio at %e. The chain seems stuck at initial point. \n", chain_number, result.acceptance_ratio);
  if(DEBUG==1 && result.acceptance_ratio>0.8) printf("\n WARNING (chain %d): Acceptance ratio at %e. The chain seems evolving too slowly. \n", chain_number, result.acceptance_ratio);

  // Check of correlation (with lag)
  if(LAG_TEST==1){

    printf(" Computing Lag function..\n");
    if(DEBUG==1) printf(" Note: This computation can take a lot of computational time.\n");
    if(DEBUG==1) printf("       Lower LAG_TEST_STOP_RATIO for faster results.\n");
    C_Delta(chain,result,LAG_TEST_STOP_RATIO, chain_number);

  }

  #if CSV==1
    if(CSV_COV==1){
      FILE *fp2;
      sprintf(filename,"Ligo_abc_%d.covmat",chain_number);
      fp2 = fopen (filename, "w+");
      fprintf(fp2, "%.8e %.8e %.8e\n", gsl_matrix_get(result.covariance,0,0),gsl_matrix_get(result.covariance,0,1),gsl_matrix_get(result.covariance,0,2));
      fprintf(fp2, "%.8e %.8e %.8e\n", gsl_matrix_get(result.covariance,1,0),gsl_matrix_get(result.covariance,1,1),gsl_matrix_get(result.covariance,1,2));
      fprintf(fp2, "%.8e %.8e %.8e\n", gsl_matrix_get(result.covariance,2,0),gsl_matrix_get(result.covariance,2,1),gsl_matrix_get(result.covariance,2,2));
      fclose(fp2);
    }
  #endif

  #if CSV==1
    if(CSV_CORR==1){
      FILE *fp3;
      sprintf(filename,"Ligo_abc_%d.corr",chain_number);
      fp3 = fopen (filename, "w+");
      fprintf(fp3, "%.8e %.8e %.8e\n", 1.,gsl_matrix_get(result.covariance,0,1)/sqrt(fabs(gsl_matrix_get(result.covariance,0,0)))/sqrt(fabs(gsl_matrix_get(result.covariance,1,1))),gsl_matrix_get(result.covariance,0,2)/sqrt(fabs(gsl_matrix_get(result.covariance,0,0)))/sqrt(fabs(gsl_matrix_get(result.covariance,2,2))));
      fprintf(fp3, "%.8e %.8e %.8e\n", gsl_matrix_get(result.covariance,1,0)/sqrt(fabs(gsl_matrix_get(result.covariance,0,0)))/sqrt(fabs(gsl_matrix_get(result.covariance,1,1))),1.,gsl_matrix_get(result.covariance,1,2)/sqrt(fabs(gsl_matrix_get(result.covariance,1,1)))/sqrt(fabs(gsl_matrix_get(result.covariance,2,2))));
      fprintf(fp3, "%.8e %.8e %.8e\n", gsl_matrix_get(result.covariance,2,0)/sqrt(fabs(gsl_matrix_get(result.covariance,0,0)))/sqrt(fabs(gsl_matrix_get(result.covariance,2,2))),gsl_matrix_get(result.covariance,2,1)/sqrt(fabs(gsl_matrix_get(result.covariance,1,1)))/sqrt(fabs(gsl_matrix_get(result.covariance,2,2))),1.);
      fclose(fp3);
    }
  #endif

  // Free all the variables
  free(chain);

  gsl_vector_free(pos_ini);
  gsl_vector_free(pos_mean_p1);
  gsl_vector_free(pos);
  gsl_rng_free(r);

  gsl_vector_free(pos_temp2);
  gsl_matrix_free(S);
  gsl_matrix_free(L);

  #if CSV==1
    fclose(fp);
  #endif

  return result;

}

inline double R(double sigma2_chain, double sigma2_mean){

  int noc = N_CHAINS;

  return sqrt(((N-1.)/N*sigma2_chain + (noc+1.)/N/noc*sigma2_mean)/sigma2_chain);

}

gsl_vector * GR_Test(chain_res * chain_results) {

  int i, k;
  int noc = N_CHAINS;

  gsl_vector * mean = gsl_vector_alloc(3);
  gsl_vector * sigma2_chain = gsl_vector_alloc(3);
  gsl_vector * sigma2_mean = gsl_vector_alloc(3);

  for(k=0;k<3;k++){
    double mean_temp = 0.;
    double sigma2_temp = 0.;
    for(i=0;i<noc;i++){
      mean_temp = mean_temp + gsl_vector_get(chain_results[i].mean,k)/noc;
      sigma2_temp = sigma2_temp + fabs(gsl_matrix_get(chain_results[i].covariance,k,k))/noc;
    }
    gsl_vector_set(mean, k, mean_temp);
    gsl_vector_set(sigma2_chain, k, sigma2_temp);
  }

  for(k=0;k<3;k++){
    double mean_temp = 0.;
    double sigma2_temp = 0.;
    for(i=0;i<noc;i++){
      sigma2_temp = sigma2_temp + pow(gsl_vector_get(mean,k)-gsl_vector_get(chain_results[i].mean,k),2.)*N/(noc-1.);
    }
    gsl_vector_set(sigma2_mean, k, sigma2_temp);
  }

  double R_a = R(gsl_vector_get(sigma2_chain,0), gsl_vector_get(sigma2_mean,0));
  double R_b = R(gsl_vector_get(sigma2_chain,1), gsl_vector_get(sigma2_mean,1));
  double R_c = R(gsl_vector_get(sigma2_chain,2), gsl_vector_get(sigma2_mean,2));

  printf("\n Gelman-Rubin test.\n");
  printf(" \t- R_a : %f\n", R_a);
  printf(" \t- R_b : %f\n", R_b);
  printf(" \t- R_c : %f\n", R_c);

  if(DEBUG==1) printf(" Converge if (approximately) R < 1.2. For better results R < 1.01.\n");

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

  gsl_vector_free(mean);
  gsl_vector_free(sigma2_chain);
  gsl_vector_free(sigma2_mean);

  gsl_vector * result = gsl_vector_alloc(3);
  gsl_vector_set(result,0,R_a);gsl_vector_set(result,1,R_b);gsl_vector_set(result,2,R_c);

  return result;

}

void init_msg(int proposal){

  printf(" Starting the MCMC computation with Metropolis-Hastings algorithm.\n");
  printf(" Properties of the chains:\n");
  if(proposal==0){
    printf(" \t- Proposal distribution: flat.\n");
  }
  if(proposal!=0){
    printf(" \t- Proposal distribution: Gaussian.\n");
  }
  if(OPENMP==0){
    printf(" \t- The chain will be %d points long.\n", N);
  }
  if(OPENMP==1){
    printf(" \t- Each of the %d chains will be %d points long.\n", N_CHAINS, N);
  }
  printf(" ------------------------------------------------------------------------ \n");

}

int main(){

  int proposal=PROPOSAL;

  printf(" ------------------------------------------------------------------------ \n");
  printf("\t       mimetic_bayes program ");
  if(DEBUG==1) printf("(with DEBUG mode on)");
  printf("\n");
  printf(" ------------------------------------------------------------------------ \n");

  if(OPENMP==1){

    int i;

    chain_res * chain_results;
    chain_results = (chain_res *) malloc(sizeof(chain_res) * N_CHAINS);

    if(OMP_THREAD_NUM==0) OMP_THREAD_NUM = omp_get_max_threads();
    if(DEBUG==1) {
      printf(" Info on OpenMP on this computer. \n");
      printf(" \t- Number of processors available = %d\n", omp_get_num_procs ( ) );
      printf(" \t- Number of threads              = %d\n", omp_get_max_threads ( ) );
      printf(" \t- Number of working threads      = %d\n", OMP_THREAD_NUM );
      printf(" Specify the number of working threads in .c file.\n");
      printf(" ------------------------------------------------------------------------ \n");
    }

    printf(" Number of chains requested = %d. \n", N_CHAINS);

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

    if(OMP_THREAD_NUM>omp_get_max_threads()){
      printf(" Fatal Error. Can't set number of threads (OMP_THREAD_NUM = %d) bigger than the available number of threads (%d). \n",OMP_THREAD_NUM, omp_get_max_threads());
      printf(" \t Change OMP_THREAD_NUM and try again. \n");
      printf(" The program will be terminated. \n");
      exit(0);
    }

    omp_set_num_threads(OMP_THREAD_NUM);

    init_msg(proposal);

    printf("\n Doing parallelized workflow. Please wait. \n\n");
    if(DEBUG==1) printf(" NOTE: Some WARNING might be shown. \n");
    if(DEBUG==1) printf("       Consider to stop the program and check the parameters if WARNINGs appear. \n\n");

    #pragma omp parallel for
    for(i=1;i<=N_CHAINS;i++)
    {
      int current_thread = omp_get_thread_num();
      if(DEBUG==1) printf(" Thread %d working on chain %d .\n", current_thread, i);
      chain_results[i-1]=LHSampling(proposal, i);
      if(DEBUG==1) printf(" Thread %d finished working on chain %d. \n", current_thread, i);
    }

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

    for(i=0;i<N_CHAINS;i++){

      printf(" CHAIN %d RESULTS. \n", i+1);

      printf("\n Acceptance ratio: %f \n\n", chain_results[i].acceptance_ratio);
      printf(" Results: Mean value pm sqrt(covariant_ii).\n");
      printf(" \ta: %e pm %e\n", gsl_vector_get(chain_results[i].mean,0), sqrt(fabs(gsl_matrix_get(chain_results[i].covariance,0,0))));
      printf(" \tb: %e pm %e\n", gsl_vector_get(chain_results[i].mean,1), sqrt(fabs(gsl_matrix_get(chain_results[i].covariance,1,1))));
      printf(" \tc: %e pm %e\n", gsl_vector_get(chain_results[i].mean,2), sqrt(fabs(gsl_matrix_get(chain_results[i].covariance,2,2))));

      printf("\n Results: Correlation coefficients.\n");
      printf(" \tab: %e \n", gsl_matrix_get(chain_results[i].covariance,0,1)/sqrt(fabs(gsl_matrix_get(chain_results[i].covariance,0,0)))/sqrt(fabs(gsl_matrix_get(chain_results[i].covariance,1,1))));
      printf(" \tac: %e \n", gsl_matrix_get(chain_results[i].covariance,0,2)/sqrt(fabs(gsl_matrix_get(chain_results[i].covariance,0,0)))/sqrt(fabs(gsl_matrix_get(chain_results[i].covariance,2,2))));
      printf(" \tbc: %e \n", gsl_matrix_get(chain_results[i].covariance,1,2)/sqrt(fabs(gsl_matrix_get(chain_results[i].covariance,1,1)))/sqrt(fabs(gsl_matrix_get(chain_results[i].covariance,2,2))));
      printf(" These are the Pearson coefficients : covariant_ij / sigma_i sigma_j\n\n");

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

    }

    gsl_vector * GR_test_result = gsl_vector_alloc(3);
    GR_test_result = GR_Test(chain_results);
    double R_a = gsl_vector_get(GR_test_result,0);
    double R_b = gsl_vector_get(GR_test_result,1);
    double R_c = gsl_vector_get(GR_test_result,2);

    if(DEBUG==1){

      int k = 0;

      double C00 = sqrt(gsl_matrix_get(chain_results[0].covariance,0,0));
      double C11 = sqrt(gsl_matrix_get(chain_results[0].covariance,1,1));
      double C22 = sqrt(gsl_matrix_get(chain_results[0].covariance,2,2));

      printf("\n INFORMATIONS for better convergence. \n\n");
      printf(" WARNING: These are hints based on the results and thumb rules. \n");
      printf("          Follow these hints at your own risk. \n\n");

      if(R_a>1.2||R_b>1.2||R_c>1.2) {
        printf(" - An R_i>1.2 : Chains are sampling different parameter spaces.\n");
        printf("                There might be some troubles in your SPEED vector definition.\n");
        k++;
      }
      if((R_a<1.2&&R_a>1.01)||(R_b<1.2&&R_b>1.01)||(R_c<1.2&&R_c>1.01)) {
        printf(" - An R_i>1.01 : The results might be good, but the results can be improved.\n");
        printf("                 Try different values for the SPEED vector to achieve better results.\n");
        k++;
      }

      double p_value = 0.5;
      if(fabs(SPEED[0]-C00/4.)/fabs(SPEED[0]+C00/4.)>p_value){
        printf(" - From covariance matrix: Try with a SPEED[0] = %e .\n",C00/4.);
        k++;
      }
      if(fabs(SPEED[1]-C11/4.)/fabs(SPEED[1]+C11/4.)>p_value){
        printf(" - From covariance matrix: Try with a SPEED[1] = %e .\n",C11/4.);
        k++;
      }
      if(fabs(SPEED[2]-C22/4.)/fabs(SPEED[2]+C22/4.)>p_value){
        printf(" - From covariance matrix: Try with a SPEED[2] = %e .\n",C22/4.);
        k++;
      }

      int j=0, r=0;
      for(i=0;i<N_CHAINS;i++){
        if(chain_results[i].acceptance_ratio>0.5) j++;
        if(chain_results[i].acceptance_ratio<0.1) r++;
      }
      if(j>0){
        printf(" - From acceptance rations: Maybe some acceptance ratios are too big.\n");
        printf("                            The configuration space might be explored slowly.\n");
        printf("                            This is not necessarily a problem,\n");
        printf("                            but if the results are not good:\n");
        printf("                            try with different SPEED vector values.\n");
        k++;
      }
      if(r>0){
        printf(" - From acceptance rations: Maybe some acceptance ratios are too small.\n");
        printf("                            The configuration space might be explored very inefficiently.\n");
        printf("                            This is not necessarily a problem,\n");
        printf("                            but if the results are not good:\n");
        printf("                            try with different SPEED vector values.\n");
        k++;
      }

      double mean_a = gsl_vector_get(chain_results[0].mean,0);
      double sigma_a = sqrt(fabs(gsl_matrix_get(chain_results[0].covariance,0,0)));
      if(PROPOSAL==1 && (fabs(mean_a-MEAN_P1[0])/fabs(mean_a+MEAN_P1[0])>3.*sigma_a)) {
        printf(" - A burn-in removal might be needed: Try with MEAN_P1[0] = %e.\n",mean_a);
        printf("   Note that this can improve but also worsen the results.\n");
        printf("   Revert the result and make a burn-in removal instead if the results get worse.\n");
      }
      double mean_b = gsl_vector_get(chain_results[0].mean,1);
      double sigma_b = sqrt(fabs(gsl_matrix_get(chain_results[0].covariance,1,1)));
      if(PROPOSAL==1 && (fabs(mean_b-MEAN_P1[1])/fabs(mean_b+MEAN_P1[1])>3.*sigma_b)) {
        printf(" - A burn-in removal might be needed: Try with MEAN_P1[1] = %e.\n",mean_b);
        printf("   Note that this can improve but also worsen the results.\n");
        printf("   Revert the result and make a burn-in removal instead if the results get worse.\n");
      }
      double mean_c = gsl_vector_get(chain_results[0].mean,2);
      double sigma_c = sqrt(fabs(gsl_matrix_get(chain_results[0].covariance,2,2)));
      if(PROPOSAL==1 && (fabs(mean_a-MEAN_P1[2])/fabs(mean_a+MEAN_P1[2])>3.*sigma_c)) {
        printf(" - A burn-in removal might be needed: Try with MEAN_P1[2] = %e .\n",mean_c);
        printf("   Note that this can improve but also worsen the results.\n");
        printf("   Revert the result and make a burn-in removal instead if the results get worse.\n");
      }

      if(k==0){
        printf(" Congratulations! The convercence seems good. \n");
      }

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

    }

    free(chain_results);

  }
  else{

    init_msg(proposal);

    LHSampling(proposal, 0);
    printf(" ------------------------------------------------------------------------ \n");

  }

  return 0;

}


text:
#ifndef MITRO_GRAPH_H
#define MITRO_GRAPH_H

#include <algorithm>
#include <cstdint>
#include <functional>
#include <memory>
#include <set>
#include <string>
#include <vector>
#include <gsl.h>

class Graph;

struct GraphEdge
{
    double weight = 0;
    int from = -1;
    int to = -1;
    GraphEdge(double weight, int from, int to) : weight(weight), from(from), to(to) {}
    bool operator<(const GraphEdge& other) const { return std::make_tuple(weight, from, to) < std::make_tuple(other.weight, other.from, other.to); }
    int otherEnd(int oneEnd) const
    {
        if (oneEnd != from && oneEnd != to)
            throw "Invalid end";
        return (oneEnd == from) ? to : from;
    }
};

struct EdgePtrCompare
{
    bool operator() (const GraphEdge* lhs, const GraphEdge* rhs)
    {
        return *lhs < *rhs;
    }
};

typedef std::set<gsl::not_null<GraphEdge*>, EdgePtrCompare> EdgePtrSet;
typedef std::vector<GraphEdge*> EdgePtrVector;

template <class ValueType, class WrappedValueType, bool is_const_iterator = false>
class basic_graph_iterator
{
    using vector_in_creator = std::vector<WrappedValueType>;
    using iterator_in_creator = typename std::conditional<is_const_iterator,
        typename vector_in_creator::const_iterator,
        typename vector_in_creator::iterator
    >::type;
public:
    using iterator_category = std::random_access_iterator_tag;
    using difference_type = typename iterator_in_creator::difference_type;
    using value_type = ValueType;
    using pointer = typename std::conditional<is_const_iterator, const value_type*, value_type*>::type;
    using reference = typename std::conditional<is_const_iterator, const value_type&, value_type&>::type;

    basic_graph_iterator(const basic_graph_iterator& other) = default;
    basic_graph_iterator& operator=(const basic_graph_iterator& other) = default;
    explicit basic_graph_iterator(const iterator_in_creator &it) : it(it) {}

    friend bool operator==(const basic_graph_iterator& l, const basic_graph_iterator& r)  { return l.it == r.it; }
    friend bool operator!=(const basic_graph_iterator& l, const basic_graph_iterator& r)  { return l.it != r.it; }
    friend bool operator<(const basic_graph_iterator& l, const basic_graph_iterator& r)   { return l.it < r.it; }

    reference operator*() { return **it; }
    pointer  operator->() { return *it; }

    basic_graph_iterator& operator++()
    {
        ++it;
        return *this;
    }

    basic_graph_iterator operator++(int)
    {
        basic_graph_iterator tmp(*this);
        ++it;
        return tmp;
    }

private:
    iterator_in_creator it;
};

class GraphNode
{
public:
	friend class Graph;
    using edge_iterator  = basic_graph_iterator<GraphEdge, gsl::not_null<GraphEdge*>, false>;
    using const_edge_iterator = basic_graph_iterator<GraphEdge, gsl::not_null<GraphEdge*>, true>;

    GraphNode(double weight = 0.0): weight(weight) {};
	~GraphNode() = default;

	GraphNode &operator=(const GraphNode &other) = delete;
	GraphNode(const GraphNode &other) = default;

    void setWeight(double w) { weight = w; }
    double getWeight() const { return weight; }

	const_edge_iterator begin() const { return const_edge_iterator{ edges.begin() }; }
	const_edge_iterator end() const { return const_edge_iterator{ edges.end() }; }

	edge_iterator begin() { return edge_iterator{ edges.begin() }; }
	edge_iterator end() { return edge_iterator{ edges.end() }; }

    int getEdgeCount() const { return int(edges.size()); }
	const GraphEdge& getEdge(int idx) const { return *edges[idx]; }
private:
    void addEdge(gsl::not_null<GraphEdge*> edge) { edges.push_back(edge); }
    void removeEdge(gsl::not_null<GraphEdge*> edge) { edges.erase(std::remove(edges.begin(), edges.end(), edge), edges.end()); }

    std::vector<gsl::not_null<GraphEdge*>> edges;
	double weight;
};

class Graph
{
public:
    friend class GraphNode;
    friend std::istream &operator>>(std::istream &is, Graph &g);
    friend std::ostream &operator<<(std::ostream &os, const Graph &g);

    using node_iterator = basic_graph_iterator<GraphNode, std::unique_ptr<GraphNode>, false>;
    using const_node_iterator = basic_graph_iterator<GraphNode, std::unique_ptr<GraphNode>, true>;
    using edge_iterator = basic_graph_iterator<GraphEdge, std::unique_ptr<GraphEdge>, false>;
    using const_edge_iterator = basic_graph_iterator<GraphEdge, std::unique_ptr<GraphEdge>, true>;

    explicit Graph(bool dir) : directed(dir) {}
	~Graph() = default;

	Graph &operator=(const Graph &other) = delete;
	Graph(const Graph &other) = delete;

	int addNode(double weight = 0.0);

	void addEdge(int from, int to, double weight = 0.0);
	void removeEdge(int from, int to);
    bool hasEdge(int from, int to);

    const_node_iterator begin() const { return const_node_iterator{nodes.begin()}; }
    const_node_iterator end() const { return const_node_iterator{nodes.end()}; }

    node_iterator begin() { return node_iterator{nodes.begin()}; }
    node_iterator end() { return node_iterator{nodes.end()}; }

    const_edge_iterator edges_begin() const { return const_edge_iterator{edges.begin()}; }
    const_edge_iterator edges_end() const { return const_edge_iterator{edges.end()}; }

    edge_iterator edges_begin() { return edge_iterator{edges.begin()}; }
    edge_iterator edges_end() { return edge_iterator{edges.end()}; }

	void resize(int newNodes);
	GraphNode &getNode(int idx) { return *nodes[idx]; }
	const GraphNode &getNode(int idx) const { return *nodes[idx]; }
	int getNodeCount() const { return int(nodes.size()); }

    bool isDirected() const { return directed; }

    bool hasWeightedNodes() const { return weightedNodes; }
    void setWeightedNodes(bool weighted) { weightedNodes = weighted; }

    bool hasWeightedEdges() const { return weightedEdges; }
    void setWeightedEdges(bool weighted) { weightedEdges = weighted; }

    void clear(bool newDirected);
private:
    bool directed;
    bool weightedNodes = false;
    bool weightedEdges = true;
    using node_vector = std::vector<std::unique_ptr<GraphNode>>;
    using edge_vector = std::vector<std::unique_ptr<GraphEdge>>;
    node_vector nodes;
    edge_vector edges;
};

std::istream &operator>>(std::istream &is, Graph &g);
std::ostream &operator<<(std::ostream &os, const Graph &g);

#endif // MITRO_UTIL_GRAPH_H

ext: h


text:
/* cdf/inverse_normal.c
 *
 * Copyright (C) 2002 Przemyslaw Sliwa and Jason H. Stover.
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307, USA.
 */

/*
 * Computes the inverse normal cumulative distribution function
 * according to the algorithm shown in
 *
 *      Wichura, M.J. (1988).
 *      Algorithm AS 241: The Percentage Points of the Normal Distribution.
 *      Applied Statistics, 37, 477-484.
 */

#include <config.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_cdf.h>

#include "rat_eval.h"

static double
small (double q)
{
  const double a[8] = { 3.387132872796366608, 133.14166789178437745,
    1971.5909503065514427, 13731.693765509461125,
    45921.953931549871457, 67265.770927008700853,
    33430.575583588128105, 2509.0809287301226727
  };

  const double b[8] = { 1.0, 42.313330701600911252,
    687.1870074920579083, 5394.1960214247511077,
    21213.794301586595867, 39307.89580009271061,
    28729.085735721942674, 5226.495278852854561
  };

  double r = 0.180625 - q * q;

  double x = q * rat_eval (a, 8, b, 8, r);

  return x;
}

static double
intermediate (double r)
{
  const double a[] = { 1.42343711074968357734, 4.6303378461565452959,
    5.7694972214606914055, 3.64784832476320460504,
    1.27045825245236838258, 0.24178072517745061177,
    0.0227238449892691845833, 7.7454501427834140764e-4
  };

  const double b[] = { 1.0, 2.05319162663775882187,
    1.6763848301838038494, 0.68976733498510000455,
    0.14810397642748007459, 0.0151986665636164571966,
    5.475938084995344946e-4, 1.05075007164441684324e-9
  };

  double x = rat_eval (a, 8, b, 8, (r - 1.6));

  return x;
}

static double
tail (double r)
{
  const double a[] = { 6.6579046435011037772, 5.4637849111641143699,
    1.7848265399172913358, 0.29656057182850489123,
    0.026532189526576123093, 0.0012426609473880784386,
    2.71155556874348757815e-5, 2.01033439929228813265e-7
  };

  const double b[] = { 1.0, 0.59983220655588793769,
    0.13692988092273580531, 0.0148753612908506148525,
    7.868691311456132591e-4, 1.8463183175100546818e-5,
    1.4215117583164458887e-7, 2.04426310338993978564e-15
  };

  double x = rat_eval (a, 8, b, 8, (r - 5.0));

  return x;
}

double
gsl_cdf_ugaussian_Pinv (const double P)
{
  double r, x, pp;

  double dP = P - 0.5;

  if (P == 1.0)
    {
      return GSL_POSINF;
    }
  else if (P == 0.0)
    {
      return GSL_NEGINF;
    }

  if (fabs (dP) <= 0.425)
    {
      x = small (dP);

      return x;
    }

  pp = (P < 0.5) ? P : 1.0 - P;

  r = sqrt (-log (pp));

  if (r <= 5.0)
    {
      x = intermediate (r);
    }
  else
    {
      x = tail (r);
    }

  if (P < 0.5)
    {
      return -x;
    }
  else
    {
      return x;
    }

}

double
gsl_cdf_ugaussian_Qinv (const double Q)
{
  double r, x, pp;

  double dQ = Q - 0.5;

  if (Q == 1.0)
    {
      return GSL_NEGINF;
    }
  else if (Q == 0.0)
    {
      return GSL_POSINF;
    }

  if (fabs (dQ) <= 0.425)
    {
      x = small (dQ);

      return -x;
    }

  pp = (Q < 0.5) ? Q : 1.0 - Q;

  r = sqrt (-log (pp));

  if (r <= 5.0)
    {
      x = intermediate (r);
    }
  else
    {
      x = tail (r);
    }

  if (Q < 0.5)
    {
      return x;
    }
  else
    {
      return -x;
    }
}

double
gsl_cdf_gaussian_Pinv (const double P, const double sigma)
{
  return sigma * gsl_cdf_ugaussian_Pinv (P);
}

double
gsl_cdf_gaussian_Qinv (const double Q, const double sigma)
{
  return sigma * gsl_cdf_ugaussian_Qinv (Q);
}

ext: c


text:
#pragma once

#include <cassert>
#include <exception>
#include <functional>
#include <set>
#include <map>

#include <nod/nod.hpp>

#include <gsl-lite/gsl-lite.hpp>

#include "mutils/time/profiler.h"

#include "threadpool/threadpool.h"

#include "nodegraph/model/node.h"
#include "nodegraph/model/pin.h"

namespace NodeGraph
{

// A collection of nodes that can be computed
class Graph
{
public:
    Graph();
    virtual ~Graph();
    virtual void Clear();

    // Use this method to create nodes and add them to the m_graph
    template <typename T, typename... Args>
    T* CreateNode(Args&&... args)
    {
        PreModify();

        auto pNode = new T(*this, std::forward<Args>(args)...);
        nodes.insert(pNode);
        m_displayNodes.insert(pNode);
        m_mapIdToNode[pNode->GetId()] = pNode;

        PostModify();
        return pNode;
    }

    void DestroyNode(Node* pNode);

    bool IsType(Node& node, ctti::type_id_t type) const;

    template <class T>
    std::set<T*> Find(ctti::type_id_t type) const
    {
        std::set<T*> found;
        for (auto& pNode : nodes)
        {
            if (IsType(*pNode, type))
            {
                found.insert(static_cast<T*>(pNode));
            }
        }
        return found;
    }

    template <class T>
    std::set<T*> Find(const std::vector<ctti::type_id_t>& nodeTypes) const
    {
        std::set<T*> found;
        for (auto& t : nodeTypes)
        {
            auto f = Find<T>(t);
            found.insert(f.begin(), f.end());
        }
        return found;
    }

    void Visit(Node& node, PinDir dir, ParameterType type, std::function<bool(Node&)> fn);

    // Get the list of pins that could be on the UI
    virtual std::vector<Pin*> GetControlSurface() const;

    virtual void Compute(const std::set<Node*>& nodes, int64_t numTicks);

    const std::set<Node*>& GetNodes() const
    {
        return nodes;
    }

    const std::set<Node*>& GetDisplayNodes() const
    {
        return m_displayNodes;
    }
    void SetDisplayNodes(const std::set<Node*>& nodes)
    {
        m_displayNodes = nodes;
    }

    const std::set<Node*>& GetOutputNodes() const
    {
        return m_outputNodes;
    }
    void SetOutputNodes(const std::set<Node*>& nodes)
    {
        m_outputNodes = nodes;
    }

    void PreModify()
    {
        if (m_modifyTracker == 0)
        {
            sigBeginModify(this);
        }

        m_modifyTracker++;
    }

    void PostModify()
    {
        assert(m_modifyTracker > 0);
        m_modifyTracker--;
        if (m_modifyTracker == 0)
        {
            sigEndModify(this);
        }
    }

    void SetName(const std::string& name);
    std::string Name() const;

    // Called to notify that this graph is about to be destroyed
    void NotifyDestroy(Graph* pGraph);

    const std::map<uint64_t, Node*>& GetNodesById() const
    {
        return m_mapIdToNode;
    }

    // Signals
    nod::signal<void(Graph*)> sigBeginModify;
    nod::signal<void(Graph*)> sigEndModify;
    nod::signal<void(Graph*)> sigDestroy;

protected:
    uint32_t m_modifyTracker = 0;

    std::map<uint64_t, Node*> m_mapIdToNode;

    std::set<Node*> nodes;
    std::set<Node*> m_displayNodes;
    std::set<Node*> m_outputNodes;

    uint64_t currentGeneration = 1;
    std::string m_strName;
}; // Graph

} // namespace NodeGraph

max_stars_repo_name: Rezonality/nodegraph

max_issues_repo_name: cmaughan/nodegraph

max_forks_repo_name: cmaughan/nodegraph


text:
static char help[] = "Assemble a Mat sparsely.\n";

#include <petsc.h>

int main(int argc,char **args) {
  PetscErrorCode ierr;
  Mat    A;
  int    i1[3] = {0, 1, 2},
         j1[3] = {0, 1, 2},
         i2 = 3,
         j2[3] = {1, 2, 3},
         i3 = 1,
         j3 = 3;
  double aA1[9] = { 1.0,  2.0,  3.0,
                    2.0,  1.0, -2.0,
                   -1.0,  1.0,  1.0},
         aA2[3] = { 1.0,  1.0, -1.0},
         aA3 = -3.0;

  PetscInitialize(&argc,&args,NULL,help);

  ierr = MatCreate(PETSC_COMM_WORLD,&A); CHKERRQ(ierr);
  ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,4,4); CHKERRQ(ierr);
  ierr = MatSetFromOptions(A); CHKERRQ(ierr);
  ierr = MatSetUp(A); CHKERRQ(ierr);
  ierr = MatSetValues(A,3,i1,3,j1,aA1,INSERT_VALUES); CHKERRQ(ierr);
  ierr = MatSetValues(A,1,&i2,3,j2,aA2,INSERT_VALUES); CHKERRQ(ierr);
  ierr = MatSetValue(A,i3,j3,aA3,INSERT_VALUES); CHKERRQ(ierr);
  ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
  ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);

  MatDestroy(&A);
  return PetscFinalize();
}


text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gbpLib.h>
#include <gbpMath.h>
#include <gbpCosmo_core.h>
#include <gbpCosmo_NFW_etc.h>
#include <gsl/gsl_sf_expint.h>

// Cole and Lacey (1996) V_c(r) profile
double V_circ_NFW(double r, double M_vir, double z, int mode, cosmo_info **cosmo) {
    double c_vir;
    double R_vir;
    double V2_c_characteristic;
    double g_c;
    double r_o;
    double x;
    double r_val = 0.;

    set_NFW_params(M_vir, z, mode, cosmo, &c_vir, &R_vir);

    r_o = R_vir / c_vir;
    x   = r / r_o;
    if(x > 0.) {
        V2_c_characteristic = G_NEWTON * M_vir / R_vir;
        g_c                 = 1. / (log(1. + c_vir) - c_vir / (1. + c_vir));
        r_val               = sqrt(V2_c_characteristic * g_c * (log(1. + x) - x / (1. + x)) * c_vir / x);
    }
    return (r_val);
}

max_stars_repo_name: gbpoole/gbpCode

max_issues_repo_name: gbpoole/gbpCode

max_forks_repo_name: gbpoole/gbpCode


text:
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_vector.h>
#include "MeshAttributes.h"
#include "mathfunctions.h"
#include "oneTimeStep.h"

extern const double g;

void oneTimeStep2D (double time, double dt)
{

	// store H, Qx, Qy
	gsl_vector *gzeta = gsl_vector_alloc(3*NumEl);
	gsl_vector *gQx = gsl_vector_alloc(3*NumEl);
	gsl_vector *gQy = gsl_vector_alloc(3*NumEl);

	int j = 0;
	for (int i=0; i < NumEl; ++i) {
		for (int k=0; k < 3; ++k){
			gsl_vector_set(gzeta,j,zeta[index(i,k,3)]);
			gsl_vector_set(gQx,j, Qx[index(i,k,3)]);
			gsl_vector_set(gQy,j, Qy[index(i,k,3)]);
			++j;
		}
	}

	/************ print out zeta, Qx and  Qy for debugging ********************/
	#ifdef DEBUG
	printf("Initial values:\n zeta \n");
	gsl_vector_fprintf(stdout,gzeta,"%e");
	printf("\n Qx \n");
	gsl_vector_fprintf(stdout,gQx,"%e");
	printf("\n Qy \n");
	gsl_vector_fprintf(stdout,gQy,"%e");
	printf("\n");
	#endif
	/***************************************************************************/

	// Calculate the wet-dry status of the elements
//	wetDryStatus2D();

	// Intermediate RK step
	compute2DL(time);

	gsl_vector *gRHSZeta = gsl_vector_alloc(3*NumEl);
	gsl_vector *gRHSQx = gsl_vector_alloc(3*NumEl);
	gsl_vector *gRHSQy = gsl_vector_alloc(3*NumEl);

	j = 0;
	for (int i=0; i < NumEl; ++i) {
		for (int k=0; k < 3; ++k){
			gsl_vector_set(gRHSZeta,j,RHSZeta[index(i,k,3)]);
			gsl_vector_set(gRHSQx,j,RHSQx[index(i,k,3)]);
			gsl_vector_set(gRHSQy,j,RHSQy[index(i,k,3)]);
			++j;
		}
	}

	/************ print out RHSZeta, RHSQx and RHSQy for debugging ********************/
	#ifdef DEBUG
	printf("After first computeL:\n RHSZeta \n");
	gsl_vector_fprintf(stdout,gRHSZeta,"%e");
	printf("\n RHSQx \n");
	gsl_vector_fprintf(stdout,gRHSQx,"%e");
	printf("\n RHSQy \n");
	gsl_vector_fprintf(stdout,gRHSQy,"%e");
	printf("\n");
	#endif
	/***************************************************************************/

	// w_1 = w + dt*L;
	gsl_vector *gzeta_1 = gsl_vector_alloc(3*NumEl);
	gsl_vector *gQx_1 = gsl_vector_alloc(3*NumEl);
	gsl_vector *gQy_1 = gsl_vector_alloc(3*NumEl);
	for (int i=0; i<3*NumEl; ++i)
	{
		gsl_vector_set(gzeta_1,i,gsl_vector_get(gzeta,i));
		gsl_vector_set(gQx_1,i,gsl_vector_get(gQx,i));
		gsl_vector_set(gQy_1,i,gsl_vector_get(gQy,i));
	}

	gsl_vector_scale(gRHSZeta,dt);
	gsl_vector_scale(gRHSQx,dt);
	gsl_vector_scale(gRHSQy,dt);
	gsl_vector_add(gzeta_1,gRHSZeta);
	gsl_vector_add(gQx_1, gRHSQx);
	gsl_vector_add(gQy_1, gRHSQy);

	/************ print out zeta, Qx and  Qy for debugging ********************/
	#ifdef DEBUG
	printf("After first computeL:\n Zeta \n");
	gsl_vector_fprintf(stdout,gzeta_1,"%e");
	printf("\n Qx \n");
	gsl_vector_fprintf(stdout,gQx_1,"%e");
	printf("\n Qy \n");
	gsl_vector_fprintf(stdout,gQy_1,"%e");
	printf("\n");
	#endif
	/***************************************************************************/

	j = 0;
	for (int i=0; i<NumEl; ++i){
		for (int k =0; k < 3; ++k){
			zeta[index(i,k,3)] = gsl_vector_get(gzeta_1,j);
			Qx[index(i,k,3)] = gsl_vector_get(gQx_1,j);
			Qy[index(i,k,3)] = gsl_vector_get(gQy_1,j);
			++j;
		}
	}

/*	printf("zeta\n");
	for (int i=0; i < NumEl; ++i)
	{
		for (int k =0; k <3; ++k)
		{
			printf("%1.15f\t",junc->zeta[index(i,k,3)]);
		}
		printf("\n");

	}

	printf("Qx\n");
	for (int i=0; i < NumEl; ++i)
	{
		for (int k =0; k <3; ++k)
		{
			printf("%1.15f\t",junc->Qx[index(i,k,3)]);
		}
		printf("\n");

	}

	printf("Qy\n");
	for (int i=0; i < NumEl; ++i)
	{
		for (int k =0; k <3; ++k)
		{
			printf("%1.15f\t",junc->Qy[index(i,k,3)]);
		}
		printf("\n");

	}
*/
	// Apply minmod slope limiter on w_1
	SlopeLimiter();

	// Apply the positive-depth operator on w_1
//	PDop2D();

	j = 0;
	for (int i=0; i<NumEl; ++i)
	{
		for(int k=0; k <3; ++k)
		{
			gsl_vector_set(gzeta_1,j,zeta[index(i,k,3)]);
			gsl_vector_set(gQx_1,j,Qx[index(i,k,3)]);
			gsl_vector_set(gQy_1,j,Qy[index(i,k,3)]);
			++j;
		}
	}

	/************ print out zeta, Qx and  Qy for debugging ********************/
	#ifdef DEBUG
	printf("After first slope limiting:\n zeta \n");
	gsl_vector_fprintf(stdout,gzeta_1,"%e");
	printf("\n Qx \n");
	gsl_vector_fprintf(stdout,gQx_1,"%e");
	printf("\n Qy \n");
	gsl_vector_fprintf(stdout,gQy_1,"%e");
	printf("\n");
	#endif
	/***************************************************************************/

	// w_1 = w_1 + w;
	gsl_vector_add(gzeta_1,gzeta);
	gsl_vector_add(gQx_1,gQx);
	gsl_vector_add(gQy_1,gQy);

	// Calcualte the wet-dry status of the elements again
//	wetDryStatus2D();

	// Compute w
	compute2DL(time);

	j = 0;
	for (int i=0; i < NumEl; ++i) {
		for (int k=0; k < 3; ++k){
			gsl_vector_set(gRHSZeta,j,RHSZeta[index(i,k,3)]);
			gsl_vector_set(gRHSQx,j,RHSQx[index(i,k,3)]);
			gsl_vector_set(gRHSQy,j,RHSQy[index(i,k,3)]);
			++j;
		}
	}

	/************ print out zeta, Qx and  Qy for debugging ********************/
	#ifdef DEBUG
	printf("After second computeL:\n RHSZeta \n");
	gsl_vector_fprintf(stdout,gRHSZeta,"%e");
	printf("\n RHSQx \n");
	gsl_vector_fprintf(stdout,gRHSQx,"%e");
	printf("\n RHSQy \n");
	gsl_vector_fprintf(stdout,gRHSQy,"%e");
	printf("\n");
	#endif
	/***************************************************************************/

	gsl_vector_scale(gzeta_1,0.5);
	gsl_vector_scale(gQx_1,0.5);
	gsl_vector_scale(gQy_1,0.5);
	gsl_vector_scale(gRHSZeta,dt/2);
	gsl_vector_scale(gRHSQx,dt/2);
	gsl_vector_scale(gRHSQy,dt/2);

	gsl_vector_add(gzeta_1,gRHSZeta);
	gsl_vector_add(gQx_1,gRHSQx);
	gsl_vector_add(gQy_1,gRHSQy);

	j = 0;
	for (int i=0; i<NumEl; ++i){
		for (int k =0; k < 3; ++k){
			zeta[index(i,k,3)] = gsl_vector_get(gzeta_1,j);
			Qx[index(i,k,3)] = gsl_vector_get(gQx_1,j);
			Qy[index(i,k,3)] = gsl_vector_get(gQy_1,j);
			++j;
		}
	}

	// Apply minmod slope limiter to w
	SlopeLimiter();

	// Apply the positive depth operator on w
//	PDop2D();

	gsl_vector_free(gzeta);
	gsl_vector_free(gQx);
	gsl_vector_free(gQy);
	gsl_vector_free(gzeta_1);
	gsl_vector_free(gQx_1);
	gsl_vector_free(gQy_1);
	gsl_vector_free(gRHSZeta);
	gsl_vector_free(gRHSQx);
	gsl_vector_free(gRHSQy);

}

void time_evolution(double FinalTime)
{

	// output the file with x y and z information for the elements
	boundary_conditions();

	double time = 0;
	double dt = 1e-7;
	int Nstep = 0;
	int fileNumber = 1;

	// write out initial conditions
	FILE* Initfile1;
	FILE* Initfile2;
	FILE* Initfile3;
//	FILE* Initfile4;

	char flname1[100];
	char flname2[100];
	char flname3[100];
//	char flname4[100];
	sprintf(flname1, "000_zeta.dat");
	sprintf(flname2, "000_Qx.dat");
	sprintf(flname3, "000_Qy.dat");
//	sprintf(flname4, "000_Qdotn.dat");
	Initfile1 = fopen(flname1, "w");
	Initfile2 = fopen(flname2, "w");
	Initfile3 = fopen(flname3, "w");
//	Initfile4 = fopen(flname4, "w");

/*	fprintf(Initfile1, "H at time 0\n");
	fprintf(Initfile2, "Qx at time 0 \n");
	fprintf(Initfile3, "Qy at time 0 \n");
	fprintf(Initfile4, "Qdotn at time 0 \n");
*/
	int k;
	for (int j=0; j<NumEl; ++j)
	{
		for(k=0; k<3; ++k)
		{
			double zeta_val = zeta[index(j,k,3)];
			double zval = (z[EltoVert[j*3+k]]+z[EltoVert[j*3+((k+1)%3)]])/2;
			double H = zeta_val + zval;
			double Qx_val = Qx[index(j,k,3)];
			double Qy_val = Qy[index(j,k,3)];
		//	double Qdotn = Qx_val*normvecx[index(j,k,3)] + Qy_val*normvecy[index(j,k,3)];
			fprintf(Initfile1, "%1.15f\t\t",zeta_val);
			fprintf(Initfile2, "%1.15f\t\t",Qx_val);
			fprintf(Initfile3, "%1.15f\t\t",Qy_val);
		//	fprintf(Initfile4, "%1.15f\t\t",Qdotn);
		}
		fprintf(Initfile1, "\n");
		fprintf(Initfile2, "\n");
		fprintf(Initfile3, "\n");
//		fprintf(Initfile4, "\n");
	}
	fclose(Initfile1);
	fclose(Initfile2);
	fclose(Initfile3);
//	fclose(Initfile4);

	//wetDryStatus2D();
	while (time<FinalTime)
	{
		time = time + dt;
		printf("dt = %e\n", dt);
		printf("time = %e\n", time);

		oneTimeStep2D(time,dt);

		boundary_conditions();

		/********* Output data every twentienth step *************/
		Nstep = Nstep + 1;
		int NstepinSec = floor(1./dt);
		if (time == FinalTime)
	//	if (Nstep%500 == 0 || time == FinalTime)
//		if ((Nstep%(10*NstepinSec) == 0) || time == FinalTime)
		{

			FILE* file3;
			FILE* file4;
			FILE* file5;
//			FILE* file6;

			char filename1[100];
			char filename2[100];
			char filename3[100];
//			char filename4[100];

/*			sprintf(filename1, "Zeta%d", NumEl);
			sprintf(filename2, "Qx%d", NumEl);
			sprintf(filename3, "Qy%d", NumEl);
*/

			sprintf(filename1, "%3.3d_H.dat", fileNumber);
			sprintf(filename2, "%3.3d_Qx.dat", fileNumber);
			sprintf(filename3, "%3.3d_Qy.dat", fileNumber);

		//	sprintf(filename1, "%3.3d_H.dat", fileNumber);
		//	sprintf(filename2, "%3.3d_Qx.dat", fileNumber);
		//	sprintf(filename3, "%3.3d_Qy.dat", fileNumber);
//			sprintf(filename4, "%3.3d_Q.n.dat", fileNumber);
			file3 = fopen(filename1, "w");
			file4 = fopen(filename2, "w");
			file5 = fopen(filename3, "w");
//			file6 = fopen(filename4, "w");

			//fprintf(file3, "H at time %e\n", time);
			//fprintf(file4, "Qx at time %e \n", time);
			//fprintf(file5, "Qy at time %e \n", time);
//			fprintf(file6, "Qdotn at time %e \n", time);

			int k;
			for (int j=0; j<NumEl; ++j)
			{
				for(k=0; k<3; ++k)
				{
					double zeta_val = zeta[index(j,k,3)];
					double Qx_val = Qx[index(j,k,3)];
					double Qy_val = Qy[index(j,k,3)];
					double zval = (z[EltoVert[j*3+k]]+z[EltoVert[j*3+((k+1)%3)]])/2;
					double hval = zeta_val + zval;

//					double Qdotn = Qx_val*normvecx[index(j,k,3)] + Qy_val*normvecy[index(j,k,3)];
					fprintf(file3, "%1.15f\t\t",zeta_val);
					fprintf(file4, "%1.15f\t\t",Qx_val);
					fprintf(file5, "%1.15f\t\t",Qy_val);
//					fprintf(file6, "%1.15f\t\t",Qdotn);
				}
				fprintf(file3, "\n");
				fprintf(file4, "\n");
				fprintf(file5, "\n");
//				fprintf(file6, "\n");
			}
			fclose(file3);
			fclose(file4);
			fclose(file5);
//			fclose(file6);

			++fileNumber;
		}

		/**************** Set next time step ****************/

		// set next time step
	//	printf("max_lambda = %e\n", max_lambda);
	//	printf("minEdgLength = %e\n", minEdgLength);
		dt = 0.1*minEdgLength/max_lambda;
		if ((time + dt) > FinalTime)
			dt = FinalTime - time;

	}
}

max_stars_repo_name: evalseth/DG-RAIN

max_issues_repo_name: evalseth/DG-RAIN

max_forks_repo_name: evalseth/DG-RAIN


text:
static char help[] = "This code implements the distribution federate of the T-D power flow use-case. A global convergence for the federation is implemented in this use-case. The The distribution federate obtains the upstream distribution substation voltage from transmission federate, solves power flow, and sends the updated injection into transmission. This code does not use HELICS iterative API\n\n";

#include <petsc.h>
#include <ValueFederate.h>
#include <opendss.h>
#include <complex.h>

typedef struct {
  helics_federate vfed;
  helics_publication pub;
  helics_subscription sub;
  helics_time_t currenttime;
  helics_iteration_status currenttimeiter;
  char       pub_topic[PETSC_MAX_PATH_LEN];
} UserData;

helics_status CreateDistributionFederate(helics_federate *distfed,char* pub_topic)
{
  helics_federate_info_t fedinfo;
  helics_status   status;
  const char*    fedinitstring="--federates=1";
  helics_federate vfed;
  char            fedname[PETSC_MAX_PATH_LEN];
  PetscErrorCode ierr;

  /* Create Federate Info object that describes the federate properties */
  fedinfo = helicsFederateInfoCreate();

  ierr = PetscStrcpy(fedname,"Distribution Federate ");CHKERRQ(ierr);
  ierr = PetscStrcat(fedname,pub_topic);CHKERRQ(ierr);
  /* Set Federate name */
  status = helicsFederateInfoSetFederateName(fedinfo,fedname);

  /* Set core type from string */
  status = helicsFederateInfoSetCoreTypeFromString(fedinfo,"zmq");

  /* Federate init string */
  status = helicsFederateInfoSetCoreInitString(fedinfo,fedinitstring);

  /* Set the message interval (timedelta) for federate. Note that
     HELICS minimum message time interval is 1 ns and by default
     it uses a time delta of 1 second. What is provided to the
     setTimedelta routine is a multiplier for the default timedelta.
  */
  /* Set one second message interval */
  status = helicsFederateInfoSetTimeDelta(fedinfo,1.0);

  status = helicsFederateInfoSetLoggingLevel(fedinfo,1);

  status = helicsFederateInfoSetMaxIterations(fedinfo,100);
  /* Create value federate */
  vfed = helicsCreateValueFederate(fedinfo);

  *distfed = vfed;

  return status;
}

int main(int argc,char **argv)
{
  PetscErrorCode ierr;
  const char*    helicsversion;
  UserData       user;
  helics_status   status;
  PetscBool      flg;
  char           netfile[PETSC_MAX_PATH_LEN],dcommand[PETSC_MAX_PATH_LEN];
  char           sub_topic[PETSC_MAX_PATH_LEN];
  double complex *Stotal;
  PetscInt      iter=0;
  PetscScalar   Vm,Va;
  PetscScalar   Pg,Qg,Pgprv,Qgprv;
  PetscInt      pflowT_conv=0,pflowD_conv=0,global_conv=0;
  char          pubstr[PETSC_MAX_PATH_LEN],substr[PETSC_MAX_PATH_LEN];
  PetscInt      strlen;
  PetscReal     tol=1E-6,mis;

  PetscInitialize(&argc,&argv,NULL,help);

  helicsversion = helicsGetVersion();

  printf("D FED: Helics version = %s\n",helicsversion);

  /* Get network data file from command line */
  ierr = PetscOptionsGetString(NULL,NULL,"-netfile",netfile,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr);

  ierr = PetscSNPrintf(dcommand,PETSC_MAX_PATH_LEN-1,"Redirect ");CHKERRQ(ierr);
  ierr = PetscStrcat(dcommand,netfile);CHKERRQ(ierr);

  /* Load the D file */
  ierr = OpenDSSRunCommand(dcommand);CHKERRQ(ierr);

  flg = PETSC_FALSE;
  /* Get the distribution topic for publication stream */
  ierr = PetscOptionsGetString(NULL,NULL,"-dtopic",user.pub_topic,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr);
  if(!flg) {
    SETERRQ(PETSC_COMM_SELF,0,"Need to specify the publication name, option -dtopic <topic_name>.\n This is same as the distribution feeder file name without the extension");
  }

    /* Create distribution federate */
  status = CreateDistributionFederate(&user.vfed,user.pub_topic);
  printf("Created D FEDERATE %s\n",user.pub_topic);

  /* Register the publication */
  user.pub = helicsFederateRegisterGlobalPublication(user.vfed,user.pub_topic,"string","");
  printf("D FEDERATE %s: Publication registered\n",user.pub_topic);

  /* Subscribe to transmission federate's publication */
  ierr = PetscStrcpy(sub_topic,"Trans/");
  ierr = PetscStrcat(sub_topic,user.pub_topic);
  user.sub = helicsFederateRegisterSubscription(user.vfed,sub_topic,"string","");
  printf("D FEDERATE %s: Subscription registered\n",user.pub_topic);

  status = helicsFederateEnterInitializationMode(user.vfed);
  printf("D FEDERATE %s: Entered initialization mode\n",user.pub_topic);

  user.currenttime = 0.0;
  user.currenttimeiter = iterating;

  PetscInt solve;
  /* Solve power flow */
  //  printf("D FEDERATE %s running power flow\n",user.pub_topic);
  ierr = OpenDSSSolutionGetSolve(&solve);CHKERRQ(ierr);

  /* Send power injection to transmission */
  /* Get the net injection at the boundary bus */
  ierr = OpenDSSCircuitGetTotalPower(&Stotal);CHKERRQ(ierr);
  Pgprv = -creal(Stotal[0])/1000.0; /* Conversion to MW */
  Qgprv = -cimag(Stotal[0])/1000.0;

  ierr = PetscSNPrintf(pubstr,PETSC_MAX_PATH_LEN-1,"%18.16f,%18.16f,%d",Pgprv,Qgprv,pflowD_conv);CHKERRQ(ierr);
  status = helicsPublicationPublishString(user.pub,pubstr);
  //  printf("D FEDERATE %s sent Pg = %4.3f, Qg = %4.3f, conv = %d from T FEDERATE\n",user.pub_topic,Pgprv,Qgprv,pflowD_conv);

  status = helicsFederateEnterExecutionMode(user.vfed);
  printf("D FEDERATE %s: Entered execution mode\n",user.pub_topic);

  while(user.currenttimeiter == iterating) {
    iter++;

    status = helicsSubscriptionGetString(user.sub,substr,PETSC_MAX_PATH_LEN-1,&strlen);
    sscanf(substr,"%lf,%lf,%d",&Vm,&Va,&pflowT_conv);
    //    printf("D FEDERATE %s received Vm = %4.3f, Va = %4.3f, conv = %d from T FEDERATE\n",user.pub_topic,Vm,Va,pflowT_conv);

    global_conv = pflowT_conv & pflowD_conv;

    if(global_conv) {
      helicsFederateRequestTimeIterative(user.vfed,user.currenttime,no_iteration,&user.currenttime,&user.currenttimeiter);
    } else {
      /* Set source bus voltage */
      ierr = OpenDSSVsourcesSetPU(Vm);CHKERRQ(ierr);
      ierr = OpenDSSVsourcesSetAngleDeg(Va);CHKERRQ(ierr);

      /* 2. Solve power flow */
      //      printf("D FEDERATE %s running power flow\n",user.pub_topic);
      ierr = OpenDSSSolutionGetSolve(&solve);CHKERRQ(ierr);

      /* Send power injection to transmission */
      /* Get the net injection at the boundary bus */
      ierr = OpenDSSCircuitGetTotalPower(&Stotal);CHKERRQ(ierr);
      Pg = -creal(Stotal[0])/1000.0; /* Conversion to MW */
      Qg = -cimag(Stotal[0])/1000.0;

      mis = PetscSqrtScalar((Pg-Pgprv)/100*(Pg-Pgprv)/100 + (Qg-Qgprv)/100*(Qg-Qgprv)/100); /* Divide by 100 for conversion to pu */
      if(mis < tol) {
	pflowD_conv = 1;
      } else {
	pflowD_conv = 0;
	Pgprv = Pg;
	Qgprv = Qg;
      }
      ierr = PetscSNPrintf(pubstr,PETSC_MAX_PATH_LEN-1,"%18.16f,%18.16f,%d",Pg,Qg,pflowD_conv);CHKERRQ(ierr);
      status = helicsPublicationPublishString(user.pub,pubstr);
      //      printf("D FEDERATE %s sent Pg = %4.3f, Qg = %4.3f, conv = %d to T FEDERATE\n",user.pub_topic,Pg,Qg,pflowD_conv);

      fflush(NULL);

      /*3. Publish Pg, Qg, and convergence status to transmission */
      status = helicsFederateRequestTimeIterative(user.vfed,user.currenttime,force_iteration,&user.currenttime,&user.currenttimeiter);

      printf("Iteration %d: D Federate %s mis = %g,converged = %d\n",iter,user.pub_topic,mis,pflowD_conv);

    }

  }

  status = helicsFederateEnterExecutionModeComplete(user.vfed);
  status = helicsFederateFinalize(user.vfed);

  PetscFinalize();
  return 0;
}

max_stars_repo_name: GMLC-TDC/Use-Cases

max_issues_repo_name: GMLC-TDC/Use-Cases

max_forks_repo_name: GMLC-TDC/Use-Cases


text:
#include <stdio.h>
#include <stdlib.h>
#include <cblas.h>
#include <math.h>
#include <string.h>
#include "../inc/knnring.h"
#include "mpi.h"

// Application Entry Point
/*
 * X is this node's points from which distances will be calculated
 * n is the number of points in X and each incoming-outcoming block
 * d is the number of dimensions of each point
 * k is the number of minimum points to calculate
 */
knnresult distrAllkNN(double * X, int n, int d, int k){

	// Timers
	double start = MPI_Wtime();
	double end;

	// Get processes number and process id
	int p, id;
	MPI_Comm_rank(MPI_COMM_WORLD, &id); // Task ID
	MPI_Comm_size(MPI_COMM_WORLD, &p); // # tasks

	// Find IDs of previous and next node in ring
	int prev, nxt;
	if(id == 0)
		prev = p - 1;
	else
		prev = id - 1;
	if(id == p-1)
		nxt = 0;
	else
		nxt = id + 1;

	// MPI variables used for async communication
    MPI_Request reqSend, reqReceive;
    //MPI_Status status;

    // Y matrix for receiving points
	double* Y = calloc(n*d, sizeof(double));

	// Temporary Y matrix, used for trading in nodes with even ids so that we don't overwrite Y before sending it while receiving
	double* tempY = calloc(n*d, sizeof(double));

	// Start pre-fetching and sending data asynchronously
    MPI_Isend(X, n*d, MPI_DOUBLE, nxt, 2, MPI_COMM_WORLD, &reqSend);
    MPI_Irecv(tempY, n*d, MPI_DOUBLE, prev, 2, MPI_COMM_WORLD, &reqReceive);

	// This is the id that corresponds to the node that initially had the current block of points Y
	// It's used to reconstruct the ids array we receive so we don't have to also send ids
	int blockID = id;

	// Generate ids matrix
	int* ids = calloc(n, sizeof(int));
	for(int i = 0; i < n; i++){
		// The first node gets points [0, n], second [n+1, 2n] etc. Node 0 gets the last n points
		if(id != 0)
			ids[i] = (id - 1) * n + i;
		else
			ids[i] = (p-1) * n + i;
	}

	// First calculate using the original dataset, then move on to sending/receiving blocks from others
	knnresult results = kNN(X, X, n, n, d, k);

	// IDs in the knnresult structure are relative. They start from 0 and go up to n-1. We need to map them to the actual ids generated above
	for(int i = 0; i < n * k; i++)
		results.nidx[i] = ids[results.nidx[i]];

	// Temporary results struct, used for calculations in incoming block
	knnresult newResults;

	// Pointers and temporary arrays used for merging results.
	int ptrResults = 0;
	int ptrNewResults = 0;
	double* tempResDis = calloc(n * k, sizeof(double));
	int* tempResIDs = calloc(n * k, sizeof(int));

	// Variables to hold minimum and maximum reduction
	double mindis, maxdis, tempmindis, tempmaxdis;

	// Trade blocks p-1 times in the ring
	for(int i = 0; i < p-1; i++){

		// Wait for sending/receiving to complete before proceeding
		MPI_Wait(&reqSend, MPI_STATUS_IGNORE);
		MPI_Wait(&reqReceive, MPI_STATUS_IGNORE);

		// Write the received points into Y so we can start pre-fetching the next ones in tempY
		memcpy(Y, tempY, n * d * sizeof(double));

		// Last time we go in data doesn't need to be prefetched
		if(i < p-2){
			// Start pre-fetching
			MPI_Isend(Y, n*d, MPI_DOUBLE, nxt, 2, MPI_COMM_WORLD, &reqSend);
			MPI_Irecv(tempY, n*d, MPI_DOUBLE, prev, 2, MPI_COMM_WORLD, &reqReceive);
		}

		// Reconstruct ids array of received block of points based on the node and the number of blocks already traded
		blockID--;
		if(blockID < 0)
			blockID = p-1;
		for(int j = 0; j < n; j++){
			if(blockID != 0)
				ids[j] = (blockID - 1) * n + j;
			else
				ids[j] = (p-1) * n + j;
		}

		// Run calculations on received points
		newResults = kNN(Y, X, n, n, d, k);

		// Again map ids as done initially
		for(int l = 0; l < n * k; l++)
			newResults.nidx[l] = ids[newResults.nidx[l]];

		// Copy points and ids array of current results to temporary arrays to safely alter the results struct bellow
		memcpy(tempResDis, results.ndist, n * k * sizeof(double));
		memcpy(tempResIDs, results.nidx, n * k * sizeof(int));

		// Merge results and newResults arrays in a merge-sort fashion
		// Iterate for each query point
		for(int r = 0; r < n; r++){
			ptrResults = 0;
			ptrNewResults = 0;

			// Iterate for each neighbor
			for(int j = 0; j < k; j++){

				// If the point in older results is closer than the new one, add it, else add the other one and increment the appropriate array pointer
				if(tempResDis[r*k + ptrResults] < newResults.ndist[r*k + ptrNewResults]){
					results.ndist[r*k + j] = tempResDis[r*k + ptrResults];
					results.nidx[r*k + j] = tempResIDs[r*k + ptrResults];
					ptrResults++;
				}else{
					results.ndist[r*k + j] = newResults.ndist[r*k + ptrNewResults];
					results.nidx[r*k + j] = newResults.nidx[r*k + ptrNewResults];
					ptrNewResults++;
				}

			}

		}

	}

	// Find overall minimum and maximum distances of knn neighbors - excluding initial 0s
	mindis = results.ndist[1];
	maxdis = results.ndist[k-1];

	// Locally reduce minimum and maximum
	for(int i = 0; i < n; i++){
		if(results.ndist[k*i] < mindis)
			mindis = results.ndist[k*i + 1];
		if(results.ndist[k*i + k - 1] > maxdis)
			maxdis = results.ndist[k*i + k - 1];
	}

	MPI_Reduce(&mindis, &tempmindis, 1, MPI_DOUBLE, MPI_MIN, 0, MPI_COMM_WORLD);
	MPI_Reduce(&maxdis, &tempmaxdis, 1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD);

	end = MPI_Wtime();
	printf("Execution time for node %d: %f s\n", id, end-start);

	if(id == 0){

		if(mindis < tempmindis)
			tempmindis = mindis;
		if(maxdis > tempmaxdis)
			tempmaxdis = maxdis;
		printf("Minimum distance is %f, maximum is %f\n", tempmindis, tempmaxdis);
		printf("Note: this message is printed twice because the tester runs two tests: one for column major and one for row major in this order.\n");
	}

	return results;
}

// Application Entry Point
knnresult kNN(double * X, double * Y, int n, int m, int d, int k)
{

	// Calculate distances matrix D - D is row-major and nxm
	double* tempD = calculateD(X, Y, n, m, d, k);

	double* D = calloc(m*n, sizeof(double));

	// Transpose D to mxn
	//cblas_dimatcopy(CblasRowMajor, CblasTrans, n, m, 1.0, D, m, n);

	double* identityMat = calloc(n*n, sizeof(double));
	for(int f=0; f<n; f++)
		identityMat[f*n + f] = 1;

	cblas_dgemm(CblasRowMajor, CblasTrans, CblasNoTrans, m, n, n, 1, tempD, m, identityMat, n, 0, D, n);

	// Free memory
	free(identityMat);
	free(tempD);

	// Create results struct
	knnresult results;
	results.k = k;
	results.m = m;
	results.ndist = calloc(m*k, sizeof(double));
	results.nidx = calloc(m*k, sizeof(int));

	// Create ids array for the X array
	int* ids = calloc(n, sizeof(int));

	// K-Select using quickselect to find k smallest distances for each point of Y
	for(int j = 0; j < m; j++){

		// Re-set ids vector before executing quickselect each time
		for(int i = 0; i < n; i++)
				ids[i] = i;

		// Call quickselect on the row of D to sort it up to its k-th element
		kselect(D + j * n, ids, n, k);

		// Quicksort the results
		quickSort(D + j * n, ids, 0, n-1);

		// Write results (point id and distance)
		for(int l = 0; l < k; l++){

			results.ndist[j * k + l] = D[j * n + l];
			results.nidx[j * k + l] = ids[l];
		}
	}

	// Free memory and return results
	free(ids);
	return results;

}

// Function used to calculate D = sqrt(sum(X.^2,2) -2* X*Y.' + sum(Y.^2,2).');
double* calculateD(double * X, double * Y, int n, int m, int d, int k){

		// Temporary variable
		double temp = 0;

		// Distance matrix for results
		double* D = calloc(n*m, sizeof(double));

		// Euclidean norm vectors
		double* normX = calloc(n, sizeof(double));
		double* normY = calloc(m, sizeof(double));

		// Matrice to store -2*X*Y'
		double* XY = calloc(n*m, sizeof(double));

		// Calculate -2*XY (https://software.intel.com/en-us/mkl-tutorial-c-multiplying-matrices-using-dgemm)
		cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasTrans, n, m, d, -2, X, d, Y, d, 0, XY, m);

		// Calculate sum(X.^2,2), sum(Y.^2,2)
		for(int i = 0; i < n; i++){
			temp = cblas_dnrm2(d, X+i*d, 1);
			normX[i] = temp * temp;
		}
		for(int i = 0; i < m; i++){
			temp = cblas_dnrm2(d, Y+i*d, 1);
			normY[i] = temp * temp;
		}

		// XY = sum(X.^2,2) -2* X*Y.'
	        for (int i=0; i<n; i++)
	    		for(int j=0; j<m; j++)
				XY[i*m+j] += normX[i] + normY[j];

		// D = sqrt(sum(X.^2,2) -2* X*Y.' + sum(Y.^2,2).');
		for(int i = 0; i < n*m; i++)
			D[i] = sqrt(fabs(XY[i]));

		// Free memory
		free(normX);
		free(normY);
		free(XY);

		return D;
}

// A utility function to swap two elements
void swap_double(double *a, double *b){
    double t = *a;
    *a = *b;
    *b = t;
}

// A utility function to swap two elements
void swap_int(int *a, int *b){
    int t = *a;
    *a = *b;
    *b = t;
}

// QuickSort Partition function. low and high are the range of indexes in arr where partition should work
int partition (double arr[], int *ids, int low, int high){

    // Select a pivot and initialize flag to position of smallest element before pivot
    double pivot = arr[high];
    int i = (low - 1);

    // Go through the array examining each element
    for (int j = low; j <= high - 1; j++)
    {
        // If current element is smaller than the pivot, increment i and swap it out with the one currently pointed by i
        if (arr[j] < pivot)
        {
            i++;
            // Swap distances and corresponding point ids
            swap_double(&arr[i], &arr[j]);
            swap_int(&ids[i], &ids[j]);
        }
    }

    // Finally place pivot in its correct position in the array and return the position as the middle point
    swap_double(&arr[i + 1], &arr[high]);
    swap_int(&ids[i + 1], &ids[high]);
    return (i + 1);
}

// Returns the median using the QuickSelect algorithm
double kselect(double arr[], int *ids, int length, int k){
	return quickSelect(arr, ids, length, k);
}

// Returns the idx-th element of arr when arr is sorted. idx is the index (starting from 1) of the point we want to find when the array is sorted.
double quickSelect(double arr[], int *ids, int length, int idx){

    // Check to end recursion
    if (length == 1){
        return arr[0];
    }

    // Select last array element as pivot
    double pivot = arr[length - 1];
    // Get index of pivot after we partition the array
    int pivotIndex = partition(arr, ids, 0, length - 1);

    // Create the higher and lower arrays that occur after partitioning in QuickSort fashion
    int lowerLength = pivotIndex;
    pivotIndex++;
    int higherLength = (length - (lowerLength + 1));
    // At this point pivotIndex, lowerLength and higherLength all start from 1 not 0

    // Variable to store result of following recursive calls
    double result = 0;

    // This means that the point we're looking (median in our case) is in the lower partition
    if (idx <= lowerLength){
        result = quickSelect(arr, ids, lowerLength, idx);
    }
    // This means that idx-th element is our pivot point
    else if(idx == pivotIndex){
        result = pivot;
    }
    // This means that the median is in the higher partition
    else{
        result = quickSelect(arr + pivotIndex, ids + pivotIndex, higherLength, idx - pivotIndex);
    }

    // Return result
    return result;
}

// Simple quicksort implementation that also sorts the ids array according to the way the main array was sorted
void quickSort(double arr[], int *ids, int low, int high)
{
    if (low < high)
    {
        // idx is the partition point
        int idx = partition(arr, ids, low, high);

        // Divide and conquer
        quickSort(arr, ids, low, idx - 1);
        quickSort(arr, ids, idx + 1, high);
    }
}


text:
/*
 * Copyright (c) 2016-2021 lymastee, All rights reserved.
 * Contact: lymastee@hotmail.com
 *
 * This file is part of the gslib project.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 * SOFTWARE.
 */

#pragma once

#ifndef framesys_300a49a6_c103_4a1f_8166_358e823a432b_h
#define framesys_300a49a6_c103_4a1f_8166_358e823a432b_h

#include <gslib/type.h>
#include <ariel/sysop.h>
#include <ariel/config.h>
#include <ariel/rendersys.h>

__ariel_begin__

enum frame_strategy
{
    fs_busy_loop,
    fs_lazy_passive,
};

struct app_config;
struct app_env;
class framesys;
struct frame_event;
class frame_listener;
class rose;

enum frame_event_id
{
    feid_invalid,
    feid_draw,
    feid_timer,
    feid_paint,
    feid_mouse_down,
    feid_mouse_up,
    feid_mouse_move,
    feid_key_down,
    feid_key_up,
    feid_char,
    feid_show,
    feid_create,
    feid_close,
    feid_resize,
    feid_halt,
    feid_resume,
    feid_custom,
};

struct __gs_novtable frame_event abstract
{
    virtual uint get_id() const = 0;
    virtual uint get_event_size() const = 0;
};

struct frame_draw_event:
    public frame_event
{
public:
    virtual uint get_id() const override { return feid_draw; }
    virtual uint get_event_size() const override { return sizeof(*this); }
};

struct frame_close_event:
    public frame_event
{
public:
    virtual uint get_id() const override { return feid_close; }
    virtual uint get_event_size() const override { return sizeof(*this); }
};

struct frame_halt_event:
    public frame_event
{
public:
    virtual uint get_id() const override { return feid_halt; }
    virtual uint get_event_size() const override { return sizeof(*this); }
};

struct frame_resume_event:
    public frame_event
{
public:
    virtual uint get_id() const override { return feid_resume; }
    virtual uint get_event_size() const override { return sizeof(*this); }
};

struct frame_paint_event:
    public frame_event
{
    rect            boundary;

public:
    virtual uint get_id() const override { return feid_paint; }
    virtual uint get_event_size() const override { return sizeof(*this); }
};

struct frame_timer_event:
    public frame_event
{
    uint            timerid;

public:
    virtual uint get_id() const override { return feid_timer; }
    virtual uint get_event_size() const override { return sizeof(*this); }
};

struct frame_show_event:
    public frame_event
{
    bool            show;

public:
    virtual uint get_id() const override { return feid_show; }
    virtual uint get_event_size() const override { return sizeof(*this); }
};

struct frame_resize_event:
    public frame_event
{
    rect            boundary;

public:
    virtual uint get_id() const override { return feid_resize; }
    virtual uint get_event_size() const override { return sizeof(*this); }
};

struct frame_mouse_down_event:
    public frame_event
{
    uint            modifier;
    unikey          key;
    point           position;

public:
    virtual uint get_id() const override { return feid_mouse_down; }
    virtual uint get_event_size() const override { return sizeof(*this); }
};

struct frame_mouse_up_event:
    public frame_event
{
    uint            modifier;
    unikey          key;
    point           position;

public:
    virtual uint get_id() const override { return feid_mouse_up; }
    virtual uint get_event_size() const override { return sizeof(*this); }
};

struct frame_mouse_move_event:
    public frame_event
{
    uint            modifier;
    point           position;

public:
    virtual uint get_id() const override { return feid_mouse_move; }
    virtual uint get_event_size() const override { return sizeof(*this); }
};

struct frame_key_down_event:
    public frame_event
{
    uint            modifier;
    unikey          key;

public:
    virtual uint get_id() const override { return feid_key_down; }
    virtual uint get_event_size() const override { return sizeof(*this); }
};

struct frame_key_up_event:
    public frame_event
{
    uint            modifier;
    unikey          key;

public:
    virtual uint get_id() const override { return feid_key_up; }
    virtual uint get_event_size() const override { return sizeof(*this); }
};

struct frame_char_event:
    public frame_event
{
    uint            modifier;
    uint            charactor;

public:
    virtual uint get_id() const override { return feid_char; }
    virtual uint get_event_size() const override { return sizeof(*this); }
};

struct frame_create_event:
    public frame_event
{
    system_driver*  driver;
    rect            boundary;

public:
    virtual uint get_id() const override { return feid_create; }
    virtual uint get_event_size() const override { return sizeof(*this); }
};

template<frame_event_id _feid>
struct frame_event_table;

#define register_frame_event(feid, fevt) \
    template<> \
    struct frame_event_table<feid> { \
        typedef fevt type; \
    };
register_frame_event(feid_draw, frame_draw_event);
register_frame_event(feid_timer, frame_timer_event);
register_frame_event(feid_paint, frame_paint_event);
register_frame_event(feid_mouse_down, frame_mouse_down_event);
register_frame_event(feid_mouse_up, frame_mouse_up_event);
register_frame_event(feid_mouse_move, frame_mouse_move_event);
register_frame_event(feid_key_down, frame_key_down_event);
register_frame_event(feid_key_up, frame_key_up_event);
register_frame_event(feid_char, frame_char_event);
register_frame_event(feid_show, frame_show_event);
register_frame_event(feid_create, frame_create_event);
register_frame_event(feid_close, frame_close_event);
register_frame_event(feid_resize, frame_resize_event);
register_frame_event(feid_halt, frame_halt_event);
register_frame_event(feid_resume, frame_resume_event);
#undef register_frame_event

class __gs_novtable frame_listener abstract
{
public:
    virtual ~frame_listener() {}
    virtual void on_frame_start() = 0;
    virtual void on_frame_end() = 0;
    virtual bool on_frame_event(const frame_event& event) = 0;
};

class frame_dispatcher:
    public system_notify
{
public:
    virtual void on_show(bool b) override;
    virtual void on_create(system_driver* ptr, const rect& rc) override;
    virtual void on_close() override;
    virtual void on_resize(const rect& rc) override;
    virtual void on_paint(const rect& rc) override;
    virtual void on_halt() override;
    virtual void on_resume() override;
    virtual bool on_mouse_down(uint um, unikey uk, const point& pt) override;
    virtual bool on_mouse_up(uint um, unikey uk, const point& pt) override;
    virtual bool on_mouse_move(uint um, const point& pt) override;
    virtual bool on_key_down(uint um, unikey uk) override;
    virtual bool on_key_up(uint um, unikey uk) override;
    virtual bool on_char(uint um, uint ch) override;
    virtual void on_timer(uint tid) override;

protected:
    frame_listener*     _listener;

public:
    frame_dispatcher() { _listener = 0; }
    void set_listener(frame_listener* listenter) { _listener = listenter; }
    void do_draw();
    void on_frame_start();
    void on_frame_end();
};

struct frame_configs
{
    uint                handles;
    int                 width;
    int                 height;
    /* more to come.. */
};

class framesys:
    public system_driver
{
public:
    friend class frame_dispatcher;
    typedef void (framesys::*fn_empty_frame)();

public:
    void set_frame_strategy(frame_strategy stt);
    void set_notify(system_notify* notify) { _notify = notify; }
    void set_frame_listener(frame_listener* lis) { _dispatcher.set_listener(lis); }
    system_notify* get_notify() const { return _notify; }
    rendersys* get_rendersys() const { return _rendersys; }
    rose* get_rose() const { return _rose; }
    const frame_configs& get_configs() const { return _configs; }
    void initialize(const rect& rc);
    void refresh();

public:
    virtual ~framesys();
    virtual void initialize(const system_context& ctx) override;
    virtual void setup() override;
    virtual void close() override;
    virtual void set_timer(uint tid, int t) override;
    virtual void kill_timer(uint tid) override;
    virtual void update() override;
    virtual void emit(int msgid, void* msg, int size) override;
    virtual void set_ime(point pt, const font& ft) override;
    virtual void set_cursor(cursor_type curty) override;

private:
    framesys();

public:
    static int run();
    static framesys* get_framesys()
    {
        static framesys inst;
        return &inst;
    }

#if defined(WIN32) || defined(_WINDOWS)
    static void set_default_config(app_config& cfg);
    static void on_app_initialized(const app_env& env);
    static void on_app_windowed(HWND);
#endif

protected:
    void idle()
    {
        assert(_empty_frame_proc);
        (this->*_empty_frame_proc)();
    }
    void empty_frame_lazy();
    void empty_frame_busy();

protected:
    frame_strategy      _strategy;
    frame_configs       _configs;
    frame_dispatcher    _dispatcher;
    system_notify*      _notify;
    rendersys*          _rendersys;
    rose*               _rose;
    fn_empty_frame      _empty_frame_proc;
};

__ariel_end__

#endif

max_stars_repo_name: lymastee/gslib

max_issues_repo_name: lymastee/gslib

max_forks_repo_name: lymastee/gslib