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 <string.h>
#include <time.h>
#include <unistd.h>
#include <math.h>
#include <float.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

#include "dnestvars.h"

double dnest(int argc, char** argv, DNestFptrSet *fptrset, int num_params,
             char *sample_dir, int max_num_saves, double ptol, const void *arg)
{
  int opt;

  dnest_arg = arg;

  dnest_check_fptrset(fptrset);

  // cope with argv

  dnest_post_temp = 1.0;
  dnest_flag_restart = 0;
  dnest_flag_postprc = 0;
  dnest_flag_sample_info = 0;
  dnest_flag_limits = 0;

  strcpy(file_save_restart, "restart_dnest.txt");
  strcpy(dnest_sample_postfix, "\0");
  strcpy(dnest_sample_tag, "\0");

  opterr = 0;
  optind = 0;
  while( (opt = getopt(argc, argv, "r:s:pt:clx:g:")) != -1)
  {
    switch(opt)
    {
      case 'r':
        dnest_flag_restart = 1;
        strcpy(file_restart, optarg);
        printf("# Dnest restarts.\n");
        break;
      case 's':
        strcpy(file_save_restart, optarg);
        //printf("# Dnest sets restart file %s.\n", file_save_restart);
        break;
      case 'p':
        dnest_flag_postprc = 1;
        dnest_post_temp = 1.0;
        printf("# Dnest does postprocess.\n");
        break;
      case 't':
        dnest_post_temp = atof(optarg);
        printf("# Dnest sets a temperature %f.\n", dnest_post_temp);
        if(dnest_post_temp == 0.0)
        {
          printf("# Dnest incorrect option -t %s.\n", optarg);
          exit(0);
        }
        if(dnest_post_temp < 1.0)
        {
          printf("# Dnest temperature should >= 1.0\n");
          exit(0);
        }
        break;
      case 'c':
        dnest_flag_sample_info = 1;
        printf("# Dnest recalculates sample information.\n");
        break;
      case 'l':
        dnest_flag_limits = 1;
        printf("# Dnest level-dependent sampling.\n");
        break;
      case 'x':
        strcpy(dnest_sample_postfix, optarg);
        printf("# Dnest sets sample postfix %s.\n", dnest_sample_postfix);
        break;
      case 'g':
        strcpy(dnest_sample_tag, optarg);
        printf("# Dnest sets sample tag %s.\n", dnest_sample_tag);
        break;
      case '?':
        printf("# Dnest incorrect option -%c %s.\n", optopt, optarg);
        exit(0);
        break;
      default:
        break;
    }
  }

  setup(argc, argv, fptrset, num_params, sample_dir, max_num_saves, ptol);

  if(dnest_flag_postprc == 1)
  {
    dnest_postprocess(dnest_post_temp, max_num_saves, ptol);
    finalise();
    return post_logz;
  }

  if(dnest_flag_sample_info == 1)
  {
    dnest_postprocess(dnest_post_temp, max_num_saves, ptol);
    finalise();
    return post_logz;
  }

  if(dnest_flag_restart==1)
    dnest_restart();

  initialize_output_file();
  dnest_run();
  close_output_file();

  dnest_postprocess(dnest_post_temp, max_num_saves, ptol);

  finalise();

  return post_logz;
}

// postprocess, calculate evidence, generate posterior sample.
void dnest_postprocess(double temperature, int max_num_saves, double ptol)
{
  options_load(max_num_saves, ptol);
  postprocess(temperature);
}

void dnest_run()
{
  int i, j, k, size_all_above_incr;
  Level *pl, *levels_orig;
  int *buf_size_above, *buf_displs;
  double *plimits;

  printf("# Start diffusive nested sampling.\n");

  while(true)
  {
    //check for termination
    if(options.max_num_saves !=0 &&
        count_saves != 0 && (count_saves%options.max_num_saves == 0))
      break;

    dnest_mcmc_run();

    count_mcmc_steps += options.thread_steps;

    if(dnest_flag_limits == 1)
    {
      // limits of smaller levels should be larger than those of higher levels
      for(j=size_levels-2; j >= 0; j--)
        for(k=0; k<particle_offset_double; k++)
        {
          limits[ j * particle_offset_double *2 + k*2 ] = fmin( limits[ j * particle_offset_double *2 + k*2 ],
                  limits[ (j+1) * particle_offset_double *2 + k*2 ] );
          limits[ j * particle_offset_double *2 + k*2 + 1] = fmax( limits[ j * particle_offset_double *2 + k*2 +1 ],
                  limits[ (j+1) * particle_offset_double *2 + k*2 + 1 ] );
        }
    }

    do_bookkeeping();

    if(count_mcmc_steps >= (count_saves + 1)*options.save_interval)
    {
      save_particle();

      // save levels, limits, sync samples when running a number of steps
      if( count_saves % num_saves == 0 )
      {
        if(size_levels <= options.max_num_levels)
        {
          save_levels();

          printf("# Save levels at N= %d.\n", count_saves);
        }
        if(dnest_flag_limits == 1)
          save_limits();
        fflush(fsample_info);
        fsync(fileno(fsample_info));
        fflush(fsample);
        fsync(fileno(fsample));
        printf("# Save limits, and sync samples at N= %d.\n", count_saves);
      }

      //if( count_saves % num_saves_restart == 0 )
      //{
      //  dnest_save_restart();
      //}
    }
  }

  //dnest_save_restart();

  //save levels
  save_levels();
  if(dnest_flag_limits == 1)
    save_limits();

  /* output state of sampler */
  FILE *fp;
  fp = fopen(options.sampler_state_file, "w");
  fprintf(fp, "%d %d\n", size_levels, count_saves);
  fclose(fp);
}

void do_bookkeeping()
{
  int i;
  //bool created_level = false;

  if(!enough_levels(levels, size_levels) && size_above >= options.new_level_interval)
  {
    // in descending order
    qsort(above, size_above, sizeof(LikelihoodType), dnest_cmp);
    int index = (int)( (1.0/compression) * size_above);

    Level level_tmp = {above[index], 0.0, 0, 0, 0, 0};
    levels[size_levels] = level_tmp;
    size_levels++;

    printf("# Creating level %d with log likelihood = %e.\n",
               size_levels-1, levels[size_levels-1].log_likelihood.value);

    // clear out the last index records
    for(i=index; i<size_above; i++)
    {
      above[i].value = 0.0;
      above[i].tiebreaker = 0.0;
    }
    size_above = index;

    if(enough_levels(levels, size_levels))
    {
      renormalise_visits();
      options.max_num_levels = size_levels;
      printf("# Done creating levles.\n");
    }
    else
    {
      kill_lagging_particles();
    }
  }
  recalculate_log_X();
}

void recalculate_log_X()
{
  int i;

  levels[0].log_X = 0.0;
  for(i=1; i<size_levels; i++)
  {
    levels[i].log_X = levels[i-1].log_X
    + log( (double)( (levels[i-1].exceeds + 1.0/compression * regularisation)
                    /(levels[i-1].visits + regularisation)  ) );
  }
}

void renormalise_visits()
{
  size_t i;

  for(i=0; i<size_levels; i++)
  {
    if(levels[i].tries >= regularisation)
    {
      levels[i].accepts = ((double)(levels[i].accepts+1) / (double)(levels[i].tries+1)) * regularisation;
      levels[i].tries = regularisation;
    }

    if(levels[i].visits >= regularisation)
    {
      levels[i].exceeds = ( (double) (levels[i].exceeds+1) / (double)(levels[i].visits + 1) ) * regularisation;
      levels[i].visits = regularisation;
    }
  }
}

void kill_lagging_particles()
{
  static unsigned int deletions = 0;

  bool *good;
  good = (bool *)malloc(options.num_particles * sizeof(bool));

  double max_log_push = -DBL_MAX;

  double kill_probability = 0.0;
  unsigned int num_bad = 0;
  size_t i;

  for(i=0; i<options.num_particles; i++)good[i] = true;

  for(i=0; i<options.num_particles; i++)
  {
    if( log_push(level_assignments[i]) > max_log_push)
      max_log_push = log_push(level_assignments[i]);

    kill_probability = pow(1.0 - 1.0/(1.0 + exp(-log_push(level_assignments[i]) - 4.0)), 3);
    if(gsl_rng_uniform(dnest_gsl_r) <= kill_probability)
    {
      good[i] = false;
      ++num_bad;
    }
  }

  if(num_bad < options.num_particles)
  {
    for(i=0; i< options.num_particles; i++)
    {
      if(!good[i])
      {
        int i_copy;
        do
        {
          i_copy = gsl_rng_uniform_int(dnest_gsl_r, options.num_particles);
        }while(!good[i_copy] || gsl_rng_uniform(dnest_gsl_r) >= exp(log_push(level_assignments[i_copy]) - max_log_push));

        memcpy(particles+i*particle_offset_size, particles + i_copy*particle_offset_size, dnest_size_of_modeltype);
        log_likelihoods[i] = log_likelihoods[i_copy];
        level_assignments[i] = level_assignments[i_copy];

        kill_action(i, i_copy);

        deletions++;

        printf("# Replacing lagging particle.\n");
        printf("# This has happened %d times.\n", deletions);
      }
    }
  }
  else
    printf("# Warning: all particles lagging!.\n");

  free(good);
}

/* save levels */
void save_levels()
{
  if(!save_to_disk)
    return;

  int i;
  FILE *fp;

  fp = fopen(options.levels_file, "w");
  fprintf(fp, "# log_X, log_likelihood, tiebreaker, accepts, tries, exceeds, visits\n");
  for(i=0; i<size_levels; i++)
  {
    fprintf(fp, "%14.12g %14.12g %f %llu %llu %llu %llu\n", levels[i].log_X, levels[i].log_likelihood.value,
      levels[i].log_likelihood.tiebreaker, levels[i].accepts,
      levels[i].tries, levels[i].exceeds, levels[i].visits);
  }
  fclose(fp);

  /* update state of sampler */
  fp = fopen(options.sampler_state_file, "w");
  fprintf(fp, "%d %d\n", size_levels, count_saves);
  fclose(fp);
}

void save_limits()
{
  int i, j;
  FILE *fp;

  fp = fopen(options.limits_file, "w");
  for(i=0; i<size_levels; i++)
  {
    fprintf(fp, "%d  ", i);
    for(j=0; j<particle_offset_double; j++)
      fprintf(fp, "%f  %f  ", limits[i*2*particle_offset_double+j*2], limits[i*2*particle_offset_double+j*2+1]);

    fprintf(fp, "\n");
  }
  fclose(fp);
}

/* save particle */
void save_particle()
{
  count_saves++;

  if(!save_to_disk)
    return;

  int whichparticle, whichtask;
  void *particle_message;

  if(count_saves%10 == 0)
    printf("#[%.1f%%] Saving sample N= %d.\n", 100.0*count_saves/options.max_num_saves, count_saves);

  whichparticle =  gsl_rng_uniform_int(dnest_gsl_r,options.num_particles);

  print_particle(fsample, particles + whichparticle * particle_offset_size, dnest_arg);

  fprintf(fsample_info, "%d %e %f %d\n", level_assignments[whichparticle],
        log_likelihoods[whichparticle].value,
        log_likelihoods[whichparticle].tiebreaker,
        whichparticle);
}

void dnest_mcmc_run()
{
  unsigned int which;
  unsigned int i;

  for(i = 0; i<options.thread_steps; i++)
  {

    /* randomly select out one particle to update */
    which = gsl_rng_uniform_int(dnest_gsl_r, options.num_particles);

    dnest_which_particle_update = which;

    //if(count_mcmc_steps >= 10000)printf("FFFF\n");
    //printf("%d\n", which);
    //printf("%f %f %f\n", particles[which].param[0], particles[which].param[1], particles[which].param[2]);
    //printf("level:%d\n", level_assignments[which]);
    //printf("%e\n", log_likelihoods[which].value);

    if(gsl_rng_uniform(dnest_gsl_r) <= 0.5)
    {
      update_particle(which);
      update_level_assignment(which);
    }
    else
    {
      update_level_assignment(which);
      update_particle(which);
    }

    if( !enough_levels(levels, size_levels)  && levels[size_levels-1].log_likelihood.value < log_likelihoods[which].value)
    {
      above[size_above] = log_likelihoods[which];
      size_above++;
    }
  }
}

void update_particle(unsigned int which)
{
  void *particle = particles+ which*particle_offset_size;
  LikelihoodType *logl = &(log_likelihoods[which]);

  Level *level = &(levels[level_assignments[which]]);

  void *proposal = (void *)malloc(dnest_size_of_modeltype);
  LikelihoodType logl_proposal;
  double log_H;

  memcpy(proposal, particle, dnest_size_of_modeltype);
  dnest_which_level_update = level_assignments[which];

  log_H = perturb(proposal, dnest_arg);

  logl_proposal.value = log_likelihoods_cal(proposal, dnest_arg);
  logl_proposal.tiebreaker =  (*logl).tiebreaker + gsl_rng_uniform(dnest_gsl_r);
  dnest_wrap(&logl_proposal.tiebreaker, 0.0, 1.0);

  if(log_H > 0.0)
    log_H = 0.0;

  dnest_perturb_accept[which] = 0;
  if( gsl_rng_uniform(dnest_gsl_r) <= exp(log_H) && level->log_likelihood.value < logl_proposal.value)
  {
    memcpy(particle, proposal, dnest_size_of_modeltype);
    memcpy(logl, &logl_proposal, sizeof(LikelihoodType));
    level->accepts++;

    dnest_perturb_accept[which] = 1;
    accept_action();
    account_unaccepts[which] = 0; /* reset the number of unaccepted perturb */
  }
  else
  {
    account_unaccepts[which] += 1; /* number of unaccepted perturb */
  }
  level->tries++;

  unsigned int current_level = level_assignments[which];
  for(; current_level < size_levels-1; ++current_level)
  {
    levels[current_level].visits++;
    if(levels[current_level+1].log_likelihood.value < log_likelihoods[which].value)
      levels[current_level].exceeds++;
    else
      break; // exit the loop if it does not satify higher levels
  }
  free(proposal);
}

void update_level_assignment(unsigned int which)
{
  int i;

  int proposal = level_assignments[which]
                 + (int)( pow(10.0, 2*gsl_rng_uniform(dnest_gsl_r))*gsl_ran_ugaussian(dnest_gsl_r));

  if(proposal == level_assignments[which])
    proposal =  ((gsl_rng_uniform(dnest_gsl_r) < 0.5)?(proposal-1):(proposal+1));

  proposal=mod_int(proposal, size_levels);

  double log_A = -levels[proposal].log_X + levels[level_assignments[which]].log_X;

  log_A += log_push(proposal) - log_push(level_assignments[which]);

  if(size_levels == options.max_num_levels)
    log_A += options.beta*log( (double)(levels[level_assignments[which]].tries +1)/ (levels[proposal].tries +1) );

  if(log_A > 0.0)
    log_A = 0.0;

  if( gsl_rng_uniform(dnest_gsl_r) <= exp(log_A) && levels[proposal].log_likelihood.value < log_likelihoods[which].value)
  {
    level_assignments[which] = proposal;

// update the limits of the level
    if(dnest_flag_limits == 1)
    {
      double *particle = (double *) (particles+ which*particle_offset_size);
      for(i=0; i<particle_offset_double; i++)
      {
        limits[proposal * 2 * particle_offset_double +  i*2] =
            fmin(limits[proposal * 2* particle_offset_double +  i*2], particle[i]);
        limits[proposal * 2 * particle_offset_double +  i*2+1] =
            fmax(limits[proposal * 2 * particle_offset_double +  i*2+1], particle[i]);
      }
    }

  }

}

double log_push(unsigned int which_level)
{
  if(which_level > size_levels)
  {
    printf("level overflow %d %d.\n", which_level, size_levels);
    exit(0);
  }
  if(enough_levels(levels, size_levels))
    return 0.0;

  int i = which_level - (size_levels - 1);
  return ((double)i)/options.lambda;
}

bool enough_levels(Level *l, int size_l)
{
  int i;

  if(options.max_num_levels == 0)
  {
    if(size_l >= LEVEL_NUM_MAX)
      return true;

    if(size_l < 10)
      return false;

    int num_levels_to_check = 20;
    if(size_l > 80)
      num_levels_to_check = (int)(sqrt(20) * sqrt(0.25*size_l));

    int k = size_l - 1, kc = 0;
    double tot = 0.0;
    double max = -DBL_MAX;
    double diff;

    for(i= 0; i<num_levels_to_check; i++)
    {
      diff = l[k].log_likelihood.value - l[k-1].log_likelihood.value;
      tot += diff;
      if(diff > max)
        max = diff;

      k--;
      kc++;
      if( k < 1 )
        break;
    }
    if(tot/kc < options.max_ptol && max < options.max_ptol*1.1)
      return true;
    else
      return false;
  }
  return (size_l >= options.max_num_levels);
}

void initialize_output_file()
{
  if(dnest_flag_restart !=1)
    fsample = fopen(options.sample_file, "w");
  else
    fsample = fopen(options.sample_file, "a");

  if(fsample==NULL)
  {
    fprintf(stderr, "# Cannot open file sample.txt.\n");
    exit(0);
  }
  if(dnest_flag_restart != 1)
    fprintf(fsample, "# \n");

  if(dnest_flag_restart != 1)
    fsample_info = fopen(options.sample_info_file, "w");
  else
    fsample_info = fopen(options.sample_info_file, "a");

  if(fsample_info==NULL)
  {
    fprintf(stderr, "# Cannot open file %s.\n", options.sample_info_file);
    exit(0);
  }
  if(dnest_flag_restart != 1)
    fprintf(fsample_info, "# level assignment, log likelihood, tiebreaker, ID.\n");
}

void close_output_file()
{
  fclose(fsample);
  fclose(fsample_info);
}

void setup(int argc, char** argv, DNestFptrSet *fptrset, int num_params, char *sample_dir, int max_num_saves, double ptol)
{
  int i, j;

  // root task.
  dnest_root = 0;

  // setup function pointers
  from_prior = fptrset->from_prior;
  log_likelihoods_cal = fptrset->log_likelihoods_cal;
  log_likelihoods_cal_initial = fptrset->log_likelihoods_cal_initial;
  log_likelihoods_cal_restart = fptrset->log_likelihoods_cal_restart;
  perturb = fptrset->perturb;
  print_particle = fptrset->print_particle;
  read_particle = fptrset->read_particle;
  restart_action = fptrset->restart_action;
  accept_action = fptrset->accept_action;
  kill_action = fptrset->kill_action;
  strcpy(dnest_sample_dir, sample_dir);

  // random number generator
  dnest_gsl_T = (gsl_rng_type *) gsl_rng_default;
  dnest_gsl_r = gsl_rng_alloc (dnest_gsl_T);
#ifndef Debug
  gsl_rng_set(dnest_gsl_r, time(NULL));
#else
  gsl_rng_set(dnest_gsl_r, 9999);
  printf("# debugging, dnest random seed %d\n", 9999);
#endif

  dnest_num_params = num_params;
  dnest_size_of_modeltype = dnest_num_params * sizeof(double);

  // read options
  options_load(max_num_saves, ptol);

  //dnest_post_temp = 1.0;
  compression = exp(1.0);
  regularisation = options.new_level_interval*sqrt(options.lambda);
  save_to_disk = true;

  // particles
  particle_offset_size = dnest_size_of_modeltype/sizeof(void);
  particle_offset_double = dnest_size_of_modeltype/sizeof(double);
  particles = (void *)malloc(options.num_particles*dnest_size_of_modeltype);

  // initialise sampler
  above = (LikelihoodType *)malloc(2*options.new_level_interval * sizeof(LikelihoodType));

  log_likelihoods = (LikelihoodType *)malloc(2*options.num_particles * sizeof(LikelihoodType));
  level_assignments = (unsigned int*)malloc(options.num_particles * sizeof(unsigned int));

  account_unaccepts = (unsigned int *)malloc(options.num_particles * sizeof(unsigned int));
  for(i=0; i<options.num_particles; i++)
  {
    account_unaccepts[i] = 0;
  }

  if(options.max_num_levels != 0)
  {
    levels = (Level *)malloc(options.max_num_levels * sizeof(Level));
    if(dnest_flag_limits == 1)
    {
      limits = malloc(options.max_num_levels * particle_offset_double * 2 * sizeof(double));
      for(i=0; i<options.max_num_levels; i++)
      {
        for(j=0; j<particle_offset_double; j++)
        {
          limits[i*2*particle_offset_double+ j*2] = DBL_MAX;
          limits[i*2*particle_offset_double + j*2 + 1] = -DBL_MAX;
        }
      }
    }
  }
  else
  {
    levels = (Level *)malloc(LEVEL_NUM_MAX * sizeof(Level));

    if(dnest_flag_limits == 1)
    {
      limits = malloc(LEVEL_NUM_MAX * particle_offset_double * 2 * sizeof(double));
      for(i=0; i<LEVEL_NUM_MAX; i++)
      {
        for(j=0; j<particle_offset_double; j++)
        {
          limits[i*2*particle_offset_double + j*2] = DBL_MAX;
          limits[i*2*particle_offset_double + j*2 + 1] = -DBL_MAX;
        }
      }
    }
  }

  dnest_perturb_accept = malloc(options.num_particles * sizeof(int));
  for(i=0; i<options.num_particles; i++)
  {
    dnest_perturb_accept[i] = 0;
  }

  count_mcmc_steps = 0;
  count_saves = 0;
  num_saves = (int)fmax(0.02*options.max_num_saves, 1.0);
  num_saves_restart = (int)fmax(0.2 * options.max_num_saves, 1.0);

// first level
  size_levels = 0;
  size_above = 0;
  LikelihoodType like_tmp = {-DBL_MAX, gsl_rng_uniform(dnest_gsl_r)};
  Level level_tmp = {like_tmp, 0.0, 0, 0, 0, 0};
  levels[size_levels] = level_tmp;
  size_levels++;

  for(i=0; i<options.num_particles; i++)
  {
    dnest_which_particle_update = i;
    dnest_which_level_update = 0;
    from_prior(particles+i*particle_offset_size, dnest_arg);
    log_likelihoods[i].value = log_likelihoods_cal_initial(particles+i*particle_offset_size, dnest_arg);
    log_likelihoods[i].tiebreaker = dnest_rand();
    level_assignments[i] = 0;
  }

  /*ModelType proposal;
  printf("%f %f %f \n", particles[0].param[0], particles[0].param[1], particles[0].param[2] );
  proposal = particles[0];
  printf("%f %f %f \n", proposal.param[0], proposal.param[1], proposal.param[2] );*/
}

void finalise()
{
  free(particles);
  free(above);
  free(log_likelihoods);
  free(level_assignments);
  free(levels);

  free(account_unaccepts);

  if(dnest_flag_limits == 1)
    free(limits);

  gsl_rng_free(dnest_gsl_r);

  free(dnest_perturb_accept);

  printf("# Finalizing dnest.\n");
}

void options_load(int max_num_saves, double ptol)
{
  //sscanf(buf, "%d", &options.num_particles);
  options.num_particles = 2;

  //fgets(buf, BUF_MAX_LENGTH, fp);
  //sscanf(buf, "%d", &options.new_level_interval);
  options.new_level_interval = options.num_particles * dnest_num_params*10;

  //fgets(buf, BUF_MAX_LENGTH, fp);
  //sscanf(buf, "%d", &options.save_interval);
  options.save_interval = options.new_level_interval;

  //fgets(buf, BUF_MAX_LENGTH, fp);
  //sscanf(buf, "%d", &options.thread_steps);
  options.thread_steps = options.new_level_interval;

  //fgets(buf, BUF_MAX_LENGTH, fp);
  //sscanf(buf, "%d", &options.max_num_levels);
  options.max_num_levels = 0;

  //fgets(buf, BUF_MAX_LENGTH, fp);
  //sscanf(buf, "%lf", &options.lambda);
  options.lambda = 10.0;

  //fgets(buf, BUF_MAX_LENGTH, fp);
  //sscanf(buf, "%lf", &options.beta);
  options.beta = 100.0;

  //fgets(buf, BUF_MAX_LENGTH, fp);
  //sscanf(buf, "%d", &options.max_num_saves);
  options.max_num_saves = max_num_saves;

  options.max_ptol = ptol;

  //fgets(buf, BUF_MAX_LENGTH, fp);
  //sscanf(buf, "%s", options.sample_file);
  strcpy(options.sample_file, dnest_sample_dir);
  strcat(options.sample_file,"/sample");
  strcat(options.sample_file, dnest_sample_tag);
  strcat(options.sample_file, ".txt");
  strcat(options.sample_file, dnest_sample_postfix);

  //fgets(buf, BUF_MAX_LENGTH, fp);
  //sscanf(buf, "%s", options.sample_info_file);
  strcpy(options.sample_info_file, dnest_sample_dir);
  strcat(options.sample_info_file,"/sample_info");
  strcat(options.sample_info_file, dnest_sample_tag);
  strcat(options.sample_info_file, ".txt");
  strcat(options.sample_info_file, dnest_sample_postfix);

  //fgets(buf, BUF_MAX_LENGTH, fp);
  //sscanf(buf, "%s", options.levels_file);
  strcpy(options.levels_file, dnest_sample_dir);
  strcat(options.levels_file,"/levels");
  strcat(options.levels_file, dnest_sample_tag);
  strcat(options.levels_file, ".txt");
  strcat(options.levels_file, dnest_sample_postfix);

  //fgets(buf, BUF_MAX_LENGTH, fp);
  //sscanf(buf, "%s", options.sampler_state_file);
  strcpy(options.sampler_state_file, dnest_sample_dir);
  strcat(options.sampler_state_file,"/sampler_state");
  strcat(options.sampler_state_file, dnest_sample_tag);
  strcat(options.sampler_state_file, ".txt");
  strcat(options.sampler_state_file, dnest_sample_postfix);

  //fgets(buf, BUF_MAX_LENGTH, fp);
  //sscanf(buf, "%s", options.posterior_sample_file);
  strcpy(options.posterior_sample_file, dnest_sample_dir);
  strcat(options.posterior_sample_file,"/posterior_sample");
  strcat(options.posterior_sample_file, dnest_sample_tag);
  strcat(options.posterior_sample_file, ".txt");
  strcat(options.posterior_sample_file, dnest_sample_postfix);

  //fgets(buf, BUF_MAX_LENGTH, fp);
  //sscanf(buf, "%s", options.posterior_sample_info_file);
  strcpy(options.posterior_sample_info_file, dnest_sample_dir);
  strcat(options.posterior_sample_info_file,"/posterior_sample_info");
  strcat(options.posterior_sample_info_file, dnest_sample_tag);
  strcat(options.posterior_sample_info_file, ".txt");
  strcat(options.posterior_sample_info_file, dnest_sample_postfix);

  //fgets(buf, BUF_MAX_LENGTH, fp);
  //sscanf(buf, "%s", options.limits_file);
  strcpy(options.limits_file, dnest_sample_dir);
  strcat(options.limits_file,"/limits");
  strcat(options.limits_file, dnest_sample_tag);
  strcat(options.limits_file, ".txt");
  strcat(options.limits_file, dnest_sample_postfix);

  // check options.

  if(options.new_level_interval < options.thread_steps)
  {
    printf("# incorrect options:\n");
    printf("# new level interval should be equal to or larger than");
    printf("  totaltask * thread step.\n");
    exit(0);
  }

/*  strcpy(options.sample_file, "sample.txt");
  strcpy(options.sample_info_file, "sample_info.txt");
  strcpy(options.levels_file, "levels.txt");
  strcpy(options.sampler_state_file, "sampler_state.txt");*/
}

double mod(double y, double x)
{
  if(x > 0.0)
  {
    return (y/x - floor(y/x))*x;
  }
  else if(x == 0.0)
  {
    return 0.0;
  }
  else
  {
    printf("Warning in mod(double, double) %e\n", x);
    exit(0);
  }

}

void dnest_wrap(double *x, double min, double max)
{
  *x = mod(*x - min, max - min) + min;
}

void wrap_limit(double *x, double min, double max)
{

  *x = fmax(fmin(*x, max), min);
}

int mod_int(int y, int x)
{
  if(y >= 0)
    return y - (y/x)*x;
  else
    return (x-1) - mod_int(-y-1, x);
}

double dnest_randh()
{
  return pow(10.0, 1.5 - 3.0*fabs(gsl_ran_tdist(dnest_gsl_r, 2))) * gsl_ran_ugaussian(dnest_gsl_r);
}

double dnest_rand()
{
  return gsl_rng_uniform(dnest_gsl_r);
}

int dnest_rand_int(int size)
{
  return gsl_rng_uniform_int(dnest_gsl_r, size);
}

double dnest_randn()
{
  return gsl_ran_ugaussian(dnest_gsl_r);
}

int dnest_cmp(const void *pa, const void *pb)
{
  LikelihoodType *a = (LikelihoodType *)pa;
  LikelihoodType *b = (LikelihoodType *)pb;

  // in decesending order
  if(a->value > b->value)
    return false;
  if( a->value == b->value && a->tiebreaker > b->tiebreaker)
    return false;

  return true;
}

int dnest_get_size_levels()
{
  return size_levels;
}

int dnest_get_which_level_update()
{
  return dnest_which_level_update;
}

int dnest_get_which_particle_update()
{
  return dnest_which_particle_update;
}

unsigned int dnest_get_which_num_saves()
{
  return num_saves;
}
unsigned int dnest_get_count_saves()
{
  return count_saves;
}

unsigned long long int dnest_get_count_mcmc_steps()
{
  return count_mcmc_steps;
}

void dnest_get_posterior_sample_file(char *fname)
{
  strcpy(fname, options.posterior_sample_file);
  return;
}
/*
 * version check
 *
 *  1: greater
 *  0: equal
 * -1: lower
 */
int dnest_check_version(char *version_str)
{
  int major, minor, patch;

  sscanf(version_str, "%d.%d.%d", &major, &minor, &patch);

  if(major > DNEST_MAJOR_VERSION)
    return 1;
  if(major < DNEST_MAJOR_VERSION)
    return -1;

  if(minor > DNEST_MINOR_VERSION)
    return 1;
  if(minor < DNEST_MINOR_VERSION)
    return -1;

  if(patch > DNEST_PATCH_VERSION)
    return 1;
  if(patch > DNEST_PATCH_VERSION)
    return -1;

  return 0;
}

void dnest_check_fptrset(DNestFptrSet *fptrset)
{
  if(fptrset->from_prior == NULL)
  {
    printf("\"from_prior\" function is not defined.\n");
    exit(0);
  }

  if(fptrset->print_particle == NULL)
  {
    //printf("\"print_particle\" function is not defined. \
    //  \nSet to be default function in dnest.\n");
    fptrset->print_particle = dnest_print_particle;
  }

  if(fptrset->read_particle == NULL)
  {
    //printf("\"read_particle\" function is not defined. \
    //  \nSet to be default function in dnest.\n");
    fptrset->read_particle = dnest_read_particle;
  }

  if(fptrset->log_likelihoods_cal == NULL)
  {
    printf("\"log_likelihoods_cal\" function is not defined.\n");
    exit(0);
  }

  if(fptrset->log_likelihoods_cal_initial == NULL)
  {
    //printf("\"log_likelihoods_cal_initial\" function is not defined. \
    //  \nSet to the same as \"log_likelihoods_cal\" function.\n");
    fptrset->log_likelihoods_cal_initial = fptrset->log_likelihoods_cal;
  }

  if(fptrset->log_likelihoods_cal_restart == NULL)
  {
    //printf("\"log_likelihoods_cal_restart\" function is not defined. \
    //  \nSet to the same as \"log_likelihoods_cal\" function.\n");
    fptrset->log_likelihoods_cal_restart = fptrset->log_likelihoods_cal;
  }

  if(fptrset->perturb == NULL)
  {
    printf("\"perturb\" function is not defined.\n");
    exit(0);
  }

  if(fptrset->restart_action == NULL)
  {
    //printf("\"restart_action\" function is not defined.\
    //  \nSet to the default function in dnest.\n");
    fptrset->restart_action = dnest_restart_action;
  }

  if(fptrset->accept_action == NULL)
  {
    //printf("\"accept_action\" function is not defined.\
    //  \nSet to the default function in dnest.\n");
    fptrset->accept_action = dnest_accept_action;
  }

  if(fptrset->kill_action == NULL)
  {
    //printf("\"kill_action\" function is not defined.\
    //  \nSet to the default function in dnest.\n");
    fptrset->kill_action = dnest_kill_action;
  }

  return;
}

DNestFptrSet * dnest_malloc_fptrset()
{
  DNestFptrSet * fptrset;
  fptrset = (DNestFptrSet *)malloc(sizeof(DNestFptrSet));

  fptrset->from_prior = NULL;
  fptrset->log_likelihoods_cal = NULL;
  fptrset->log_likelihoods_cal_initial = NULL;
  fptrset->log_likelihoods_cal_restart = NULL;
  fptrset->perturb = NULL;
  fptrset->print_particle = NULL;
  fptrset->read_particle = NULL;
  fptrset->restart_action = NULL;
  fptrset->accept_action = NULL;
  fptrset->kill_action = NULL;
  return fptrset;
}

void dnest_free_fptrset(DNestFptrSet * fptrset)
{
  free(fptrset);
  return;
}

/*!
 *  Save sampler state for later restart.
 */
void dnest_save_restart()
{
  FILE *fp;
  int i, j;
  void *particles_all;
  LikelihoodType *log_likelihoods_all;
  unsigned int *level_assignments_all;
  char str[200];

  sprintf(str, "%s_%d", file_save_restart, count_saves);
  fp = fopen(str, "wb");
  if(fp == NULL)
  {
    fprintf(stderr, "# Error: Cannot open file %s. \n", file_save_restart);
    exit(0);
  }

  printf("# Save restart data to file %s.\n", str);

    //fprintf(fp, "%d %d\n", count_saves, count_mcmc_steps);
    //fprintf(fp, "%d\n", size_levels_combine);

  fwrite(&count_saves, sizeof(int), 1, fp);
  fwrite(&count_mcmc_steps, sizeof(int), 1, fp);
  fwrite(&size_levels, sizeof(int), 1, fp);

  for(i=0; i<size_levels; i++)
  {
    //fprintf(fp, "%14.12g %14.12g %f %llu %llu %llu %llu\n", levels_combine[i].log_X, levels_combine[i].log_likelihood.value,
    //  levels_combine[i].log_likelihood.tiebreaker, levels_combine[i].accepts,
    //  levels_combine[i].tries, levels[i].exceeds, levels_combine[i].visits);

    fwrite(&levels[i], sizeof(Level), 1, fp);
  }

  for(i=0; i<options.num_particles; i++)
  {
    //fprintf(fp, "%d %e %f\n", level_assignments_all[j*options.num_particles + i],
    //  log_likelihoods_all[j*options.num_particles + i].value,
    //  log_likelihoods_all[j*options.num_particles + i].tiebreaker);
    fwrite(&level_assignments[j*options.num_particles + i], sizeof(int), 1, fp);
    fwrite(&log_likelihoods[j*options.num_particles + i], sizeof(LikelihoodType), 1, fp);
  }

  if(dnest_flag_limits == 1)
  {
    for(i=0; i<size_levels; i++)
    {
      //fprintf(fp, "%d  ", i);
      for(j=0; j<particle_offset_double; j++)
      {
        //fprintf(fp, "%f  %f  ", limits[i*2*particle_offset_double+j*2], limits[i*2*particle_offset_double+j*2+1]);
        fwrite(&limits[i*2*particle_offset_double+j*2], sizeof(double), 2, fp);
      }

        //fprintf(fp, "\n");
    }
  }

  for(i=0; i<options.num_particles; i++)
  {
    //print_particle(fp, particles_all + (j * options.num_particles + i) * particle_offset_size);
    fwrite(particles + (j * options.num_particles + i) * particle_offset_size, dnest_size_of_modeltype, 1, fp);
  }

  fclose(fp);

  restart_action(0);
}

void dnest_restart()
{
  FILE *fp;
  int i, j;
  void *particles_all;
  unsigned int *level_assignments_all;
  LikelihoodType *log_likelihoods_all;
  void *particle;

  fp = fopen(file_restart, "rb");
  if(fp == NULL)
  {
    fprintf(stderr, "# Error: Cannot open file %s. \n", file_restart);
    exit(0);
  }

  printf("# Reading %s\n", file_restart);

  fread(&count_saves, sizeof(int), 1, fp);
  fread(&count_mcmc_steps, sizeof(int), 1, fp);
  fread(&size_levels, sizeof(int), 1, fp);

  /* consider that the newly input max_num_levels may be different from the save one */
  if(options.max_num_levels != 0)
  {
    if(size_levels > options.max_num_levels)
    {
      printf("# input max_num_levels %d smaller than the one in restart data %d.\n", options.max_num_levels, size_levels);
      size_levels = options.max_num_levels;
    }
  }
    // read levels
  for(i=0; i<size_levels; i++)
  {
    if(i<size_levels) // not read all the levels
      fread(&levels[i], sizeof(Level), 1, fp);
    else
      fseek(fp, sizeof(Level), SEEK_CUR); /* offset the file point */
  }
    memcpy(levels, levels, size_levels * sizeof(Level));

    // read level assignment
  for(i=0; i<options.num_particles; i++)
  {
    fread(&level_assignments[j*options.num_particles + i], sizeof(int), 1, fp);
    fread(&log_likelihoods[j*options.num_particles + i], sizeof(LikelihoodType), 1, fp);

    /* reset the level assignment that exceeds the present maximum level */
    if(level_assignments[j*options.num_particles + i] > size_levels -1)
    {
      level_assignments[j*options.num_particles + i] = size_levels - 1;
    }
  }

  // read limits
  if(dnest_flag_limits == 1)
  {
    for(i=0; i<size_levels; i++)
    {
      if(i < size_levels)
      {
        for(j=0; j<particle_offset_double; j++)
        {
          fread(&limits[i*2*particle_offset_double+j*2], sizeof(double), 2, fp);
        }
      }
      else
      {
        fseek(fp, sizeof(double) * 2 * particle_offset_double, SEEK_CUR); /* offset the file point */
      }
    }
  }

  // read particles
  for(i=0; i<options.num_particles; i++)
  {
    particle = (particles + (j * options.num_particles + i) * dnest_size_of_modeltype);
    fread(particle, dnest_size_of_modeltype, 1, fp);
  }

  fclose(fp);

  if(count_saves > options.max_num_saves)
  {
    printf("# Number of samples already larger than the input number, exit!\n");
    exit(0);
  }

  num_saves = (int)fmax(0.02*(options.max_num_saves-count_saves), 1.0); /* reset num_saves */
  num_saves_restart = (int)fmax(0.2 * (options.max_num_saves-count_saves), 1.0); /* reset num_saves_restart */

  restart_action(1);

  for(i=0; i<options.num_particles; i++)
  {
    dnest_which_particle_update = i;
    dnest_which_level_update = level_assignments[i];
    //printf("%d %d %f\n", thistask, i, log_likelihoods[i].value);
    log_likelihoods[i].value = log_likelihoods_cal_restart(particles+i*particle_offset_size, dnest_arg);
    //printf("%d %d %f\n", thistask, i, log_likelihoods[i].value);
    //due to randomness, the original level assignment may be incorrect. re-asign the level
    while(log_likelihoods[i].value < levels[level_assignments[i]].log_likelihood.value)
    {
      printf("# level assignment decrease %d %f %f %d.\n", i, log_likelihoods[i].value,
        levels[level_assignments[i]].log_likelihood.value, level_assignments[i]);
      level_assignments[i]--;
    }

  }
  return;
}

void dnest_print_particle(FILE *fp, const void *model, const void *arg)
{
  int i;
  double *pm = (double *)model;

  for(i=0; i<dnest_num_params; i++)
  {
    fprintf(fp, "%e ", pm[i] );
  }
  fprintf(fp, "\n");
  return;
}

void dnest_read_particle(FILE *fp, void *model)
{
  int j;
  double *psample = (double *)model;

  for(j=0; j < dnest_num_params; j++)
  {
    if(fscanf(fp, "%lf", psample+j) < 1)
    {
      printf("%f\n", *psample);
      fprintf(stderr, "#Error: Cannot read file %s.\n", options.sample_file);
      exit(0);
    }
  }
  return;
}

void dnest_restart_action(int iflag)
{
  return;
}

void dnest_accept_action()
{
  return;
}

void dnest_kill_action(int i, int i_copy)
{
  return;
}

max_stars_repo_name: LiyrAstroph/PyCALI

max_issues_repo_name: LiyrAstroph/pyCALI

max_forks_repo_name: LiyrAstroph/pyCALI


text:
/** @file */

#ifndef __CCL_F2D_H_INCLUDED__
#define __CCL_F2D_H_INCLUDED__

#include <gsl/gsl_spline.h>
#include <gsl/gsl_interp2d.h>
#include <gsl/gsl_spline2d.h>

CCL_BEGIN_DECLS

//f2d extrapolation types for early times
typedef enum ccl_f2d_extrap_growth_t
{
  ccl_f2d_cclgrowth = 401, //Use CCL's linear growth
  ccl_f2d_customgrowth = 402, //Use a custom growth function
  ccl_f2d_constantgrowth = 403, //Use a constant growth factor
  ccl_f2d_no_extrapol = 404, //Do not extrapolate, just throw an exception
} ccl_f2d_extrap_growth_t;

//f2d interpolation types
typedef enum ccl_f2d_interp_t
{
  ccl_f2d_3 = 303, //Bicubic interpolation
} ccl_f2d_interp_t;

/**
 * Struct containing a 2D power spectrum
 */
typedef struct {
  double lkmin,lkmax; /**< Edges in log(k)*/
  double amin,amax; /**< Edges in a*/
  int is_factorizable; /**< Is this factorizable into k- and a-dependent functions? */
  int is_k_constant; /**< no k-dependence, just return 1*/
  int is_a_constant; /**< no a-dependence, just return 1*/
  int extrap_order_lok; /**< Order of extrapolating polynomial in log(k) for low k (0, 1 or 2)*/
  int extrap_order_hik; /**< Order of extrapolating polynomial in log(k) for high k (0, 1 or 2)*/
  ccl_f2d_extrap_growth_t extrap_linear_growth;  /**< Extrapolation type at high redshifts*/
  int is_log; /**< Do I hold the values of log(f(k,a))?*/
  double (*growth)(double); /**< Custom extrapolating growth function*/
  double growth_factor_0; /**< Constant extrapolating growth factor*/
  int growth_exponent; /**< Power to which growth should be exponentiated*/
  gsl_spline *fk; /**< Spline holding the values of the k-dependent factor*/
  gsl_spline *fa; /**< Spline holding the values of the a-dependent factor*/
  gsl_spline2d *fka; /**< Spline holding the values of f(k,a)*/
} ccl_f2d_t;

/**
 * Create a ccl_f2d_t structure.
 * @param na number of elements in a_arr.
 * @param a_arr array of scale factor values at which the function is defined. The array should be ordered.
 * @param nk number of elements of lk_arr.
 * @param lk_arr array of logarithmic wavenumbers at which the function is defined (i.e. this array contains ln(k), NOT k). The array should be ordered.
 * @param fka_arr array of size na * nk containing the 2D function. The 2D ordering is such that fka_arr[ia*nk+ik] = f(k=exp(lk_arr[ik]),a=a_arr[ia]).
 * @param fk_arr array of size nk containing the k-dependent part of the function. Only relevant if is_factorizable is true.
 * @param fa_arr array of size na containing the a-dependent part of the function. Only relevant if is_factorizable is true.
 * @param is_factorizable if not 0, fk_arr and fa_arr will be used as 1-D arrays to construct a factorizable 2D function.
 * @param extrap_order_lok Order of the polynomial that extrapolates on wavenumbers smaller than the minimum of lk_arr. Allowed values: 0 (constant), 1 (linear extrapolation) and 2 (quadratic extrapolation). Extrapolation happens in ln(k).
 * @param extrap_order_hik Order of the polynomial that extrapolates on wavenumbers larger than the maximum of lk_arr. Allowed values: 0 (constant), 1 (linear extrapolation) and 2 (quadratic extrapolation). Extrapolation happens in ln(k).
 * @param extrap_linear_growth: ccl_f2d_extrap_growth_t value defining how the function with scale factors below the interpolation range. Allowed values: ccl_f2d_cclgrowth (scale with the CCL linear growth factor), ccl_f2d_customgrowth (scale with a custom function of redshift passed through `growth`), ccl_f2d_constantgrowth (scale by multiplying the function at the earliest available scale factor by a constant number, defined by `growth_factor_0`), ccl_f2d_no_extrapol (throw an error if the function is ever evaluated outside the interpolation range in a). Note that, above the interpolation range (i.e. for low redshifts), the function will be assumed constant.
 * @param is_fka_log: if not zero, `fka_arr` contains ln(f(k,a)) instead of f(k,a). If the function is factorizable, then `fk_arr` holds ln(K(k)) and `fa_arr` holds ln(A(a)), where f(k,a)=K(k)*A(a).
 * @param growth: custom growth function. Irrelevant if extrap_linear_growth!=ccl_f2d_customgrowth.
 * @param growth_factor_0: custom growth function. Irrelevant if extrap_linear_growth!=ccl_f2d_constantgrowth.
 * @param growth_exponent: power to which the extrapolating growth factor should be exponentiated when extrapolating (e.g. usually 2 for linear power spectra).
 * @param interp_type: 2D interpolation method. Currently only ccl_f2d_3 is implemented (bicubic interpolation).
 * @param status Status flag. 0 if there are no errors, nonzero otherwise.
 */
ccl_f2d_t *ccl_f2d_t_new(int na,double *a_arr,
			 int nk,double *lk_arr,
			 double *fka_arr,
			 double *fk_arr,
			 double *fa_arr,
			 int is_factorizable,
			 int extrap_order_lok,
			 int extrap_order_hik,
			 ccl_f2d_extrap_growth_t extrap_linear_growth,
			 int is_fka_log,
			 double (*growth)(double),
			 double growth_factor_0,
			 int growth_exponent,
			 ccl_f2d_interp_t interp_type,
			 int *status);

/**
 * Evaluate 2D function of k and a defined by ccl_f2d_t structure.
 * @param fka ccl_f2d_t structure defining f(k,a).
 * @param lk Natural logarithm of the wavenumber.
 * @param a Scale factor.
 * @param cosmo ccl_cosmology structure, only needed if evaluating f(k,a) at small scale factors outside the interpolation range, and if fka was initialized with extrap_linear_growth = ccl_f2d_cclgrowth.
 * @param status Status flag. 0 if there are no errors, nonzero otherwise.
 */
double ccl_f2d_t_eval(ccl_f2d_t *fka,double lk,double a,void *cosmo,
		      int *status);

/**
 * F2D structure destructor.
 * Frees up all memory associated with a f2d structure.
 * @param fka Structure to be freed.
 */
void ccl_f2d_t_free(ccl_f2d_t *fka);

CCL_END_DECLS

#endif

ext: h

max_stars_repo_name: benediktdiemer/CCL

max_issues_repo_name: benediktdiemer/CCL

max_forks_repo_name: benediktdiemer/CCL


text:
#pragma once

#pragma warning(disable:4201) // nonstandard extension used : nameless struct/union
#pragma warning(disable:4238) // nonstandard extension used : class rvalue used as lvalue
#pragma warning(disable:4324) // structure was padded due to __declspec(align())

#ifndef WIN32_LEAN_AND_MEAN
    #define WIN32_LEAN_AND_MEAN
#endif
#ifndef NOMINMAX
    #define NOMINMAX
#endif
#include <windows.h>

#include <d3d12.h>

#pragma comment(lib, "d3d12.lib")
#pragma comment(lib, "dxgi.lib")

#define MY_IID_PPV_ARGS IID_PPV_ARGS
#define D3D12_GPU_VIRTUAL_ADDRESS_NULL      ((D3D12_GPU_VIRTUAL_ADDRESS)0)
#define D3D12_GPU_VIRTUAL_ADDRESS_UNKNOWN   ((D3D12_GPU_VIRTUAL_ADDRESS)-1)

#pragma region
#include <d3dx12.h>
#include <d3d12shader.h>
#include <d3dcompiler.h>
#pragma endregion

#pragma region

#pragma comment (lib, "d3dcompiler.lib")

#pragma endregion

#include <cstdint>
#include <cstdio>
#include <cstdarg>
#include <vector>
#include <memory>
#include <string>
#include <exception>
#include <numeric>
#include <wrl.h>
#include <ppltasks.h>
#include <gsl\gsl>

#include "Utility.h"
#include "VectorMath.h"
#include "EngineTuning.h"
#include "EngineProfiling.h"

template< typename T >
using NotNull = gsl::not_null< T >;

#define SHADER_ARGS(shader) (shader), \
								 sizeof(shader)

template<typename T, std::size_t N>
struct Array :public std::array<T, N>
{
    static inline std::function<T(T, T)> multiply = [](T priview, T current) {
        return priview * current;
    };
public:
    T product()
    {
        return std::accumulate(begin(), end(), 1, multiply);
    }
};
using U8 = uint8_t;
using U32 = uint32_t;
using U32x1 = U32;
using U32x2 = Array<U32, 2>;
using U32x3 = Array<U32, 3>;
using U32x4 = Array<U32, 4>;
using F32 = float;
using F32x2 = Array<F32, 2>;
using F32x3 = Array<F32, 3>;
using F32x4 = Array<F32, 4>;

#pragma region voxel
#define VOXEL_RESOLUTION 128
#define CLIP_REGION_COUNT 6
#define FACE_COUNT 6
#pragma endregion

template <typename T, size_t N > Array<T, N> DivideByMultiple(Array <T, N> value, Array <T, N> alignment)
{
    Array <T, N> result;
    for (size_t i = 0; i < N; i++)
    {
        result[i] = (T)((value[i] + alignment[i] - 1) / alignment[i]);
    }
    return result;
}


text:
#include <string.h>
#include <alloca.h>
#include <math.h>
#include <mpi.h>

#include <kdcount/kdtree.h>

#include <fastpm/libfastpm.h>
#include <gsl/gsl_rng.h>

#include <fastpm/prof.h>
#include <fastpm/logging.h>
#include <fastpm/store.h>

#include <fastpm/fof.h>

#include "pmpfft.h"
#include "pm2lpt.h"
#include "pmghosts.h"
#include "vpm.h"

#include <fastpm/io.h>
#include <fastpm/string.h>
#include <bigfile-mpi.h>

//#define FASTPM_FOF_DEBUG

struct FastPMFOFFinderPrivate {
    int ThisTask;
    int NTask;
    double * boxsize;
    MPI_Comm comm;
};

/* creating a kdtree struct
 * for store with np particles starting from start.
 * */
struct KDTreeNodeBuffer {
    void * mem;
    void * base;
    char * ptr;
    char * end;
    struct KDTreeNodeBuffer * prev;
};

static void *
_kdtree_buffered_malloc(void * userdata, size_t size)
{
    struct KDTreeNodeBuffer ** pbuffer = (struct KDTreeNodeBuffer**) userdata;

    struct KDTreeNodeBuffer * buffer = *pbuffer;

    if(buffer->base == NULL || buffer->ptr + size >= buffer->end) {
        struct KDTreeNodeBuffer * newbuffer = malloc(sizeof(newbuffer[0]));
        newbuffer->mem = buffer->mem;
        size_t newsize = 1024 * 1024 * 4; /* 4 MB for each block */
        if(newsize < size) {
            newsize = size;
        }
        newbuffer->base = fastpm_memory_alloc(buffer->mem, "KDTreeBase", newsize, FASTPM_MEMORY_STACK);
        newbuffer->ptr = newbuffer->base;
        newbuffer->end = newbuffer->base + newsize;
        newbuffer->prev = buffer;

        *pbuffer = newbuffer;
        buffer = newbuffer;
    }

    void * r = buffer->ptr;
    buffer->ptr += size;
    return r;
}

static void
_kdtree_buffered_free(void * userdata, size_t size, void * ptr)
{
    /* do nothing; */
}

static
KDNode *
_create_kdtree (KDTree * tree, int thresh,
    FastPMStore ** stores, int nstore,
    double boxsize[])
{
    /* if boxsize is NULL the tree will be non-periodic. */
    /* the allocator; started empty. */
    struct KDTreeNodeBuffer ** pbuffer = malloc(sizeof(void*));
    struct KDTreeNodeBuffer * headbuffer = malloc(sizeof(headbuffer[0]));
    *pbuffer = headbuffer;

    headbuffer->mem = stores[0]->mem;
    headbuffer->base = NULL;
    headbuffer->prev = NULL;

    int s;
    ptrdiff_t i;

    tree->userdata = pbuffer;

    tree->input.dims[0] = 0;

    for(s = 0; s < nstore; s ++) {
        tree->input.dims[0] += stores[s]->np;
    }

    tree->input.dims[1] = 3;

    if(tree->input.dims[0] < stores[0]->np_upper) {
        /* if the first store is big enough, use it for the tree */
        tree->input.buffer = (void*) &stores[0]->x[0][0];
    } else {
        /* otherwise, allocate a big buffer and make a copy */
        tree->input.buffer = _kdtree_buffered_malloc(pbuffer,
                    tree->input.dims[0] * sizeof(stores[0]->x[0]));
        memcpy(tree->input.buffer, &stores[0]->x[0][0], stores[0]->np * sizeof(stores[0]->x[0]));
    }

    i = stores[0]->np;

    /* copy the other positions to the base pointer. */
    for(s = 1; s < nstore; s ++) {
        memcpy(((char*) tree->input.buffer) + i * sizeof(stores[0]->x[0]),
                &stores[s]->x[0][0],
                stores[s]->np * sizeof(stores[0]->x[0]));

        i = i + stores[s]->np;
    }

    tree->input.strides[0] = sizeof(stores[0]->x[0]);
    tree->input.strides[1] = sizeof(stores[0]->x[0][0]);
    tree->input.elsize = sizeof(stores[0]->x[0][0]);
    tree->input.cast = NULL;

    tree->ind = _kdtree_buffered_malloc(pbuffer, tree->input.dims[0] * sizeof(tree->ind[0]));
    for(i = 0; i < tree->input.dims[0]; i ++) {
        tree->ind[i] = i;
    }
    tree->ind_size = tree->input.dims[0];

    tree->malloc = _kdtree_buffered_malloc;
    tree->free = _kdtree_buffered_free;

    tree->thresh = thresh;

    tree->boxsize = boxsize;

    KDNode * root = kd_build(tree);
    fastpm_info("Creating KDTree with %td nodes for %td particles\n", tree->size, tree->ind_size);
    return root;
}

void
_free_kdtree (KDTree * tree, KDNode * root)
{
    kd_free(root);
    struct KDTreeNodeBuffer * buffer, *q, **pbuffer = tree->userdata;

    for(buffer = *pbuffer; buffer; buffer = q) {
        if(buffer->base)
            fastpm_memory_free(buffer->mem, buffer->base);
        q = buffer->prev;
        free(buffer);
    }
    free(pbuffer);
}

void
fastpm_fof_init(FastPMFOFFinder * finder, FastPMStore * store, PM * pm)
{
    finder->priv = malloc(sizeof(FastPMFOFFinderPrivate));
    finder->p = store;
    finder->pm = pm;

    finder->event_handlers = NULL;

    if (finder->periodic)
        finder->priv->boxsize = pm_boxsize(pm);
    else
        finder->priv->boxsize = NULL;

    finder->priv->comm = pm_comm(pm);
    MPI_Comm comm = finder->priv->comm;
    MPI_Comm_rank(comm, &finder->priv->ThisTask);
    MPI_Comm_size(comm, &finder->priv->NTask);
}

static void
_fof_local_find(FastPMFOFFinder * finder,
            FastPMStore * p,
            PMGhostData * pgd,
            ptrdiff_t * head, double linkinglength)
{
    /* local find of p and the the ghosts */
    KDTree tree;

    FastPMStore * stores[2] = {p, pgd->p};

    KDNode * root = _create_kdtree(&tree, finder->kdtree_thresh, stores, 2, finder->priv->boxsize);

    kd_fof(root, linkinglength, head);

    _free_kdtree(&tree, root);
}

static int
_merge(uint64_t * src, ptrdiff_t isrc, uint64_t * dest, ptrdiff_t idest, ptrdiff_t * head)
{
    int merge = 0;
    ptrdiff_t j = head[idest];
    if(src[isrc] < dest[j]) {
        merge = 1;
    }

    if(merge) {
        dest[j] = src[isrc];
    }
    return merge;
}

struct reduce_fof_data {
    ptrdiff_t * head;
    size_t nmerged;
};

static void
FastPMReduceFOF(FastPMStore * src, ptrdiff_t isrc, FastPMStore * dest, ptrdiff_t idest, int ci, void * userdata)
{

    struct reduce_fof_data * data = userdata;

    data->nmerged += _merge(src->minid, isrc, dest->minid, idest, data->head);
}

static void
_fof_global_merge(
    FastPMFOFFinder * finder,
    FastPMStore * p,
    PMGhostData * pgd,
    uint64_t * minid,
    ptrdiff_t * head
)
{
    ptrdiff_t i;

    MPI_Comm comm = finder->priv->comm;

    size_t npmax = p->np;

    MPI_Allreduce(MPI_IN_PLACE, &npmax, 1, MPI_LONG, MPI_MAX, comm);

    /* initialize minid, used as a global tag of groups as we merge */
    for(i = 0; i < p->np; i ++) {
        /* assign unique ID to each particle; could use a better scheme with true offsets */
        minid[i] = i + finder->priv->ThisTask * npmax;
        #ifdef FASTPM_FOF_DEBUG
        /* for debugging, overwrite the previous unique ID with the true ID of particles */
        minid[i] = p->id[i];
        #endif
    }

    /* send minid */
    p->minid = minid; /* only send up to p->np */
    pm_ghosts_send(pgd, COLUMN_MINID);
    p->minid = NULL;

    /* copy over to minid for storage; FIXME: allow overriding  */
    for(i = 0; i < pgd->p->np; i ++) {
        minid[i + p->np] = pgd->p->minid[i];
    }

    /* reduce the minid of the head items according to the local connection. */

    while(1) {
        size_t nmerge = 0;
        for(i = 0; i < p->np + pgd->p->np; i ++) {
            nmerge += _merge(minid, i, minid, i, head);
        }
        if(nmerge == 0) break;
    }

#ifdef FASTPM_FOF_DEBUG
    {
        FILE * fp = fopen(fastpm_strdup_printf("dump-pos-%d.f8", finder->priv->ThisTask), "w");
        fwrite(p->x, p->np, sizeof(double) * 3, fp);
        fwrite(pgd->p->x, pgd->p->np, sizeof(double) * 3, fp);
        fclose(fp);
    }
    {
        FILE * fp = fopen(fastpm_strdup_printf("dump-id-%d.f8", finder->priv->ThisTask), "w");
        fwrite(p->id, p->np, sizeof(int64_t), fp);
        fwrite(pgd->p->id, pgd->p->np, sizeof(int64_t) * 3, fp);
        fclose(fp);
    }

#endif
    int iter = 0;

    while(1) {

        /* prepare the communication buffer, every ghost has
         * the minid and task of the current head. such that
         * they will connect to the other ranks correctly */

        for(i = 0; i < pgd->p->np; i ++) {
            pgd->p->minid[i] = minid[head[i + p->np]];
        }

        /* at this point all items on ghosts have local minid and task, ready to reduce */

        struct reduce_fof_data data = {
            .head = head,
            .nmerged = 0,
        };

        /* merge ghosts into the p, reducing the MINID on p */

        p->minid = minid; /* only update up to p->np */
        pm_ghosts_reduce(pgd, COLUMN_MINID, FastPMReduceFOF, &data);
        p->minid = NULL;

        size_t nmerged = data.nmerged;

        /* at this point all items on p->fof with head[i] = i have present minid and task */

        MPI_Allreduce(MPI_IN_PLACE, &nmerged, 1, MPI_LONG, MPI_SUM, comm);

        MPI_Barrier(comm);

        fastpm_info("FOF reduction iteration %d : merged %td crosslinks\n", iter, nmerged);

        if(nmerged == 0) break;

        for(i = 0; i < p->np + pgd->p->np; i ++) {
            if(minid[i] < minid[head[i]]) {
                fastpm_raise(-1, "p->fof invariance is broken i = %td np = %td\n", i, p->np);
            }
        }

        iter++;
    }

    /* previous loop only updated head[i]; now make sure every particles has the correct minid. */
    for(i = 0; i < p->np + pgd->p->np; i ++) {
        minid[i] = minid[head[i]];
    }

    #ifdef FASTPM_FOF_DEBUG
    {
    for(i = 0; i < p->np + pgd->p->np ; i ++) {
        uint64_t id = i <p->np?p->id[i]:pgd->p->id[i - p->np];
        if(minid[i] == 88) {
            fastpm_ilog(INFO, "%d MINID == %ld ID = %ld headminid == %ld i = %td / %td head=%td", finder->priv->ThisTask,
                minid[i], id, minid[head[i]], i, p->np, head[i]);
        }
        if(id == 88 || id == 96 || id == 152 || id == 160) {
            fastpm_ilog(INFO, "%d MINID == %ld ID = %ld headminid == %ld i = %td / %td head=%td", finder->priv->ThisTask,
                minid[i], id, minid[head[i]], i, p->np, head[i]);
        }
    }
    }
    #endif
}

/* set head[i] to hid*/
static size_t
_assign_halo_attr(FastPMFOFFinder * finder, PMGhostData * pgd, ptrdiff_t * head, size_t np, size_t np_ghosts, int nmin)
{
    ptrdiff_t * offset = fastpm_memory_alloc(finder->p->mem, "FOFOffset", sizeof(offset[0]) * (np + np_ghosts), FASTPM_MEMORY_STACK);
    uint8_t * has_remote = fastpm_memory_alloc(finder->p->mem, "FOFHasRemote", sizeof(has_remote[0]) * (np + np_ghosts), FASTPM_MEMORY_STACK);
    uint8_t * has_local = fastpm_memory_alloc(finder->p->mem, "FOFHasLocal", sizeof(has_local[0]) * (np + np_ghosts), FASTPM_MEMORY_STACK);

    ptrdiff_t i;
    for(i = 0; i < np + np_ghosts; i ++) {
        offset[i] = 0;
        has_remote[i] = 0;
        has_local[i] = 0;
    }

    /* set offset to number of particles in the halo */
    for(i = 0; i < np + np_ghosts; i ++) {
        offset[head[i]] ++;
    }

    for(i = 0; i < np; i ++) {
        has_local[head[i]] = 1;
    }

    /* if the group is connected to a remote component */
    /* if the particle has a ghost, then */
    pm_ghosts_has_ghosts(pgd, has_remote);

    /* if connected to a particle that has ghost */
    for(i = 0; i < np; i ++) {
        if(has_remote[i]) has_remote[head[i]] = 1;
    }

    /* if connected to a ghost */
    for(i = np; i < np + np_ghosts; i ++) {
        has_remote[head[i]] = 1;
    }

    size_t it = 0;

    /* assign attr index for groups at least contain 1 local particle */
    for(i = 0; i < np + np_ghosts; i ++) {
        if(has_local[i] && (has_remote[i] || offset[i] >= nmin)) {
            offset[i] = it;
            it ++;
        } else {
            offset[i] = -1;
        }
    }

    size_t nhalos = it;

    /* update head [i] to offset[head[i]], which stores the index in the halo store for this particle. */
    for(i = 0; i < np + np_ghosts; i ++) {
        head[i] = offset[head[i]];
        /* this will not happen if nmin == 1 */
        if(head[i] != -1 && head[i] > nhalos) {
            fastpm_raise(-1, "head[i] (%td) > nhalos (%td) This shall not happen.\n", head[i], nhalos);
        }
    }

    fastpm_memory_free(finder->p->mem, has_local);
    fastpm_memory_free(finder->p->mem, has_remote);
    fastpm_memory_free(finder->p->mem, offset);

    return it;
}
static double
periodic_add(double x, double wx, double y, double wy, double L)
{
    if(wx > 0) {
        while(y - x > L / 2) y -= L;
        while(y - x < -L / 2) y += L;

        return wx * x + wy * y;
    } else {
        return wy * y;
    }
}

/*
 * apply mask to a halo storage.
 * relable head, and halo->id
 * */
static void
fastpm_fof_subsample(FastPMFOFFinder * finder, FastPMStore * halos, FastPMParticleMaskType * mask, ptrdiff_t * head)
{
    /* mapping goes from the old head[i] value to the new head[i] value */
    ptrdiff_t * mapping = fastpm_memory_alloc(finder->p->mem, "FOFMapping", sizeof(mapping[0]) * halos->np, FASTPM_MEMORY_STACK);

    ptrdiff_t i;
    for(i = 0; i < halos->np; i ++) {
        mapping[i] = -1;
        halos->id[i] = i;
    }

    /* remove non-contributing halo segments */
    fastpm_store_subsample(halos, mask, halos);

    for(i = 0; i < halos->np; i ++) {
        mapping[halos->id[i]] = i;
        halos->id[i] = i;
    }

    /* adjust head[i] */

    for(i = 0; i < finder->p->np; i ++) {
        if(head[i] >= 0) {
            head[i] = mapping[head[i]];
        }
    }

    fastpm_memory_free(finder->p->mem, mapping);
}

static void
fastpm_fof_apply_length_cut(FastPMFOFFinder * finder, FastPMStore * halos, ptrdiff_t * head)
{
    MPI_Comm comm = finder->priv->comm;

    FastPMParticleMaskType * mask = fastpm_memory_alloc(finder->p->mem, "LengthMask", sizeof(mask[0]) * halos->np, FASTPM_MEMORY_STACK);
    ptrdiff_t i;
    for(i = 0; i < halos->np; i ++) {
        /* remove halos that are shorter than nmin */
        if(halos->length[i] < finder->nmin) {
            mask[i] = 0;
        } else {
            mask[i] = 1;
        }
    }
    fastpm_fof_subsample(finder, halos, mask, head);
    fastpm_memory_free(finder->p->mem, mask);

    fastpm_info("After length cut we have %td halos (%td including ghost halos).\n",
            fastpm_store_get_mask_sum(halos, comm),
            fastpm_store_get_np_total(halos, comm));
}

/* first every undecided halo; then the real halo with particles */
static int
FastPMLocalSortByMinID(const int i1,
                    const int i2,
                    FastPMStore * p)
{
    int v1 = (p->minid[i1] < p->minid[i2]);
    int v2 = (p->minid[i1] > p->minid[i2]);

    return v2 - v1;
}

static int
FastPMTargetMinID(FastPMStore * store, ptrdiff_t i, void * userdata)
{
    FastPMFOFFinder * finder = userdata;

    const uint32_t GOLDEN32 = 2654435761ul;
    /* const uint64_t GOLDEN64 = 11400714819323198549; */
    /* may over flow, but should be okay here as the periodicity is ggt NTask */
    int key = (store->minid[i] * GOLDEN32) % (unsigned) finder->priv->NTask;
    return key;
}

static int
FastPMTargetTask(FastPMStore * store, ptrdiff_t i, void * userdata)
{
    return store->task[i];
}

/* for debugging, move particles to a spatially unrelated rank */
static int
FastPMTargetFOF(FastPMStore * store, ptrdiff_t i, void * userdata)
{
#ifdef FASTPM_FOF_DEBUG
    PM * pm = userdata;
    int NTask;
    MPI_Comm_size(pm_comm(pm), &NTask);
    const uint32_t GOLDEN32 = 2654435761ul;
    /* const uint64_t GOLDEN64 = 11400714819323198549; */
    /* may over flow, but should be okay here as the periodicity is ggt NTask */
    int key = (store->id[i] * GOLDEN32) % (unsigned) NTask;
    return key;
#else
    return FastPMTargetPM(store, i, userdata);
#endif
}
/*
 * This function group halos by halos->minid[i].
 *
 * All halos with the same minid will have the same attribute values afterwards, except
 * a few book keeping items named in this function (see comments inside)
 *
 * */
static void
fastpm_fof_reduce_halo_attrs(FastPMFOFFinder * finder, FastPMStore * halos,
        void (* add_func) (FastPMFOFFinder * finder, FastPMStore * halo1, ptrdiff_t i1, FastPMStore * halo2, ptrdiff_t i2 ),
        void (* reduce_func)(FastPMFOFFinder * finder, FastPMStore * halo1, ptrdiff_t i1)
)
{

    ptrdiff_t i;

    ptrdiff_t first = -1;
    uint64_t lastminid = 0;

    /* ind is the array to use to replicate items */
    int * ind = fastpm_memory_alloc(finder->p->mem, "HaloPermutation", sizeof(ind[0]) * halos->np, FASTPM_MEMORY_STACK);

    /* the following items will have mask[i] == 1, but we will mark some to 0 if
     * they are not the principle (first) halos segment with this minid */
    for(i = 0; i < halos->np + 1; i++) {
        /* we use i == halos->np to terminate the last segment */
        if(first == -1 || i == halos->np || lastminid != halos->minid[i]) {
            if (first >= 0) {
                /* a segment ended */
                ptrdiff_t j;
                for(j = first; j < i; j ++) {
                    ind[j] = first;
                }
            }
            if (i < halos->np) {
                /* a segment started */
                halos->mask[i] = 1;
                lastminid = halos->minid[i];
            }
            /* starting a segment of the same halos */
            first = i;
        } else {
            /* inside segment, simply add the ith halo to the first of the segment */
            halos->mask[i] = 0;
            add_func(finder, halos, first, halos, i);
        }
    }

    /* use permute to replicate the first halo attr to the rest:
     *
     * we do not want to replicate
     * - fof, as fof.task is the original mpi rank of the halo
     * - id, as it is the original location of the halo on the original mpi rank.
     * - mask, whether it is primary or not
     *
     * we need to replicate because otherwise when we return the head array on the
     * original ranks will be violated.
     *   */

    FastPMStore save[1];

    /* FIXME: add a method for this! */
    memcpy(save->columns, halos->columns, sizeof(save->columns));

    halos->minid = NULL;
    halos->task = NULL;
    halos->id = NULL;
    halos->mask = NULL;

    fastpm_store_permute(halos, ind);

    memcpy(halos->columns, save->columns, sizeof(save->columns));

    for(i = 0; i < halos->np; i++) {
        reduce_func(finder, halos, i);
    }

    fastpm_memory_free(finder->p->mem, ind);
}

static void
fastpm_fof_compute_halo_attrs(FastPMFOFFinder * finder, FastPMStore * halos,
            ptrdiff_t * head,
            void (*convert_func)(FastPMFOFFinder * finder, FastPMStore * p, ptrdiff_t i, FastPMStore * halos),
            void (*add_func)(FastPMFOFFinder * finder, FastPMStore * halos, ptrdiff_t i1, FastPMStore * halos2, ptrdiff_t i2),
            void (*reduce_func)(FastPMFOFFinder * finder, FastPMStore * halos, ptrdiff_t i)
)
{
    MPI_Comm comm = finder->priv->comm;

    FastPMStore h1[1];
    fastpm_store_init(h1, "FOF", 1, halos->attributes, FASTPM_MEMORY_HEAP);
    ptrdiff_t i;

    for(i = 0; i < finder->p->np; i++) {
        ptrdiff_t hid = head[i];
        if(hid < 0) continue;

        if(hid >= halos->np) {
            fastpm_raise(-1, "halo of a particle out of bounds (%td > %td)\n", hid, halos->np);
        }

        /* initialize h1 with existing halo attributes for this particle */
        fastpm_store_take(halos, hid, h1, 0);

        convert_func(finder, finder->p, i, h1);

        add_func(finder, halos, hid, h1, 0);
    }

    fastpm_store_destroy(h1);
    /* decompose halos by minid (gather); if all halo segments of the same minid are on the same rank, we can combine
     * these into a single entry, then replicate for each particle to look up;
     * halo segments that have no local particles are never exchanged to another rank. */
    if(0 != fastpm_store_decompose(halos,
            (fastpm_store_target_func) FastPMTargetMinID, finder, comm)) {

        fastpm_raise(-1, "out of space sending halos by MinID.\n");
    }

    /* now head[i] is no longer the halo attribute of particle i. */

    /* to combine, first local sort by minid; those without local particles are moved to the beginning
     * so we can easily skip them. */
    fastpm_store_sort(halos, FastPMLocalSortByMinID);

    /* reduce and update properties */
    fastpm_fof_reduce_halo_attrs(finder, halos, add_func, reduce_func);

    /* decompose halos by task (return) */
    if (0 != fastpm_store_decompose(halos,
            (fastpm_store_target_func) FastPMTargetTask, finder, comm)) {
        fastpm_raise(-1, "out of space for gathering halos this shall never happen.\n");
    }
    /* local sort by id (restore the order) */
    fastpm_store_sort(halos, FastPMLocalSortByID);

    /* now head[i] is again the halo attribute of particle i. */
}

/*
 * compute the attrs of the local halo segments based on local particles.
 * head : the hid of each particle;
 * fofsave : the minid of each particle (unique label of each halo)
 *
 * if a halo has no local particles,  halos->mask[i] is set to 0.
 * if a halo has any local particles, and halos->mask[i] is set to 1.
 * */
static void
fastpm_fof_remove_empty_halos(FastPMFOFFinder * finder, FastPMStore * halos, uint64_t * minid, ptrdiff_t * head)
{

    /* */
    ptrdiff_t i;

    /* set minid and task of the halo; all of the halo particles of the same minid needs to be reduced */
    for(i = 0; i < finder->p->np; i++) {
        /* set the minid of the halo; need to take care of the ghosts too.
         * even though we do not add them to the attributes */
        ptrdiff_t hid = head[i];

        if(hid < 0) continue;

        if(halos->mask[hid] == 0) {
            /* halo will be reduced by minid */
            halos->minid[hid] = minid[i];
            halos->mask[hid] = 1;
        } else {
            if(halos->minid[hid] != minid[i]) {
                fastpm_raise(-1, "Consistency check failed after FOF global merge.\n");
            }
        }
    }

    /* halo will be returned to this task;
     * we save ThisTask here.*/
    for(i = 0; i < halos->np; i ++) {
        halos->task[i] = finder->priv->ThisTask;
    }

    fastpm_fof_subsample(finder, halos, halos->mask, head);
}

static void
_add_basic_halo_attrs(FastPMFOFFinder * finder, FastPMStore * halos, ptrdiff_t hid, FastPMStore * h1, ptrdiff_t i)
{
    double * boxsize = finder->priv->boxsize;

    if(halos->aemit)
        halos->aemit[hid] += h1->aemit[i];

    int d;

    for(d = 0; d < 3; d++) {
        if(halos->v)
            halos->v[hid][d] += h1->v[i][d];
        if(halos->dx1)
            halos->dx1[hid][d] += h1->dx1[i][d];
        if(halos->dx2)
            halos->dx2[hid][d] += h1->dx2[i][d];
    }

    if(halos->x) {
        for(d = 0; d < 3; d++) {
            if(boxsize) {
                halos->x[hid][d] = periodic_add(
                    halos->x[hid][d] / halos->length[hid], halos->length[hid],
                    h1->x[i][d] / h1->length[i], h1->length[i], boxsize[d]);
            } else {
                halos->x[hid][d] += h1->x[i][d];
            }
        }
    }

    if(halos->q) {
        for(d = 0; d < 3; d ++) {
            if (boxsize) {
                halos->q[hid][d] = periodic_add(
                    halos->q[hid][d] / halos->length[hid], halos->length[hid],
                    h1->q[i][d] / h1->length[i], h1->length[i], boxsize[d]);
            } else {
                halos->q[hid][d] += h1->q[i][d];
            }
        }
    }
    /* do this after the loop because x depends on the old length. */
    halos->length[hid] += h1->length[i];
}

/* convert a particle to the first halo in the halo store */
static void
_convert_basic_halo_attrs(FastPMFOFFinder * finder, FastPMStore * p, ptrdiff_t i, FastPMStore * halos)
{
    int hid = 0;

    int d;

    double q[3];

    if(halos->q && fastpm_store_has_q(p)) {
        fastpm_store_get_q_from_id(p, p->id[i], q);
    }
    halos->length[hid] = 1;

    for(d = 0; d < 3; d++) {
        if(halos->x)
            halos->x[hid][d] = p->x[i][d];
        if(halos->v)
            halos->v[hid][d] = p->v[i][d];
        if(halos->dx1)
            halos->dx1[hid][d] = p->dx1[i][d];
        if(halos->dx2)
            halos->dx2[hid][d] = p->dx2[i][d];
        if(halos->q && fastpm_store_has_q(p)) {
            halos->q[hid][d] = q[d];
        }
    }

    if(halos->aemit)
        halos->aemit[hid] = p->aemit[i];
}
static void
_reduce_basic_halo_attrs(FastPMFOFFinder * finder, FastPMStore * halos, ptrdiff_t i)
{
    int d;
    double n = halos->length[i];

    for(d = 0; d < 3; d++) {
        if(halos->x)
            halos->x[i][d] /= n;
        if(halos->v)
            halos->v[i][d] /= n;
        if(halos->dx1)
            halos->dx1[i][d] /= n;
        if(halos->dx2)
            halos->dx2[i][d] /= n;
        if(halos->q)
            halos->q[i][d] /= n;
    }
    if(halos->aemit)
        halos->aemit[i] /= n;
}

static void
_add_extended_halo_attrs(FastPMFOFFinder * finder, FastPMStore * h1, ptrdiff_t i1, FastPMStore * h2, ptrdiff_t i2)
{
    int d;
    if(h1->rvdisp) {
        for(d = 0; d < 9; d ++) {
            h1->rvdisp[i1][d] += h2->rvdisp[i2][d];
        }
    }
    if(h1->vdisp) {
        for(d = 0; d < 6; d ++) {
            h1->vdisp[i1][d] += h2->vdisp[i2][d];
        }
    }
    if(h1->rdisp) {
        for(d = 0; d < 6; d ++) {
            h1->rdisp[i1][d] += h2->rdisp[i2][d];
        }
    }
}

static void
_convert_extended_halo_attrs(FastPMFOFFinder * finder, FastPMStore * p, ptrdiff_t i, FastPMStore * halos)
{
    int hid = 0;

    int d;
    double rrel[3];

    for(d = 0; d < 3; d ++) {
        rrel[d] = p->x[i][d] - halos->x[hid][d];

        if(finder->priv->boxsize) {
            double L = finder->priv->boxsize[d];
            while(rrel[d] > L / 2) rrel[d] -= L;
            while(rrel[d] < -L / 2) rrel[d] += L;
        }
    }

    if(halos->vdisp) {
        /* FIXME: add hubble expansion term based on aemit and hubble function? needs to modify Finder object */

        double vrel[3];
        for(d = 0; d < 3; d ++) {
            vrel[d] = p->v[i][d] - halos->v[hid][d];
        }
        for(d = 0; d < 3; d ++) {
            halos->vdisp[hid][d] = vrel[d] * vrel[d];
            halos->vdisp[hid][d + 3] = vrel[d] * vrel[(d + 1) % 3];
        }
    }
    if(halos->rvdisp) {
        /* FIXME: add hubble expansion term based on aemit and hubble function? needs to modify Finder object */

        double vrel[3];
        for(d = 0; d < 3; d ++) {
            vrel[d] = p->v[i][d] - halos->v[hid][d];
        }
        for(d = 0; d < 3; d ++) {
            halos->rvdisp[hid][d] = rrel[d] * vrel[d];
            halos->rvdisp[hid][d + 3] = rrel[d] * vrel[(d + 1) % 3];
            halos->rvdisp[hid][d + 6] = rrel[d] * vrel[(d + 2) % 3];
        }
    }
    if(halos->rdisp) {
        for(d = 0; d < 3; d ++) {
            halos->rdisp[hid][d] = rrel[d] * rrel[d];
            halos->rdisp[hid][d + 3] = rrel[d] * rrel[(d + 1) % 3];
        }
    }
}
static void
_reduce_extended_halo_attrs(FastPMFOFFinder * finder, FastPMStore * halos, ptrdiff_t hid)
{
    double n = halos->length[hid];
    int d;
    if(halos->rvdisp) {
        for(d = 0; d < 9; d ++) {
            halos->rvdisp[hid][d] /= n;
        }
    }
    if(halos->vdisp) {
        for(d = 0; d < 6; d ++) {
            halos->vdisp[hid][d] /= n;
        }
    }
    if(halos->rdisp) {
        for(d = 0; d < 3; d ++) {
            halos->rdisp[hid][d] /= n;
        }
    }
}

/* This function creates the storage object for halo segments that are local on this
 * rank. We store mnay attributes. We only allow a flucutation of 2 around avg_halos.
 * this should be OK, since we will only redistribute by the MinID, which are supposed
 * to be very uniform.
 * */
static void
fastpm_fof_create_local_halos(FastPMFOFFinder * finder, FastPMStore * halos, size_t nhalos)
{

    MPI_Comm comm = finder->priv->comm;

    FastPMColumnTags attributes = finder->p->attributes;
    attributes |= COLUMN_MASK;
    attributes |= COLUMN_LENGTH | COLUMN_MINID | COLUMN_TASK;
    attributes |= COLUMN_RDISP | COLUMN_VDISP | COLUMN_RVDISP;
    attributes |= COLUMN_ACC; /* ACC used as the first particle position offset */
    attributes &= ~COLUMN_POTENTIAL;
    attributes &= ~COLUMN_DENSITY;
    attributes &= ~COLUMN_TIDAL;

    /* store initial position only for periodic case. non-periodic suggests light cone and
     * we cannot infer q from ID sensibly. (crashes there) */
    if(finder->priv->boxsize) {
        attributes |= COLUMN_Q;
    } else {
        attributes &= ~COLUMN_Q;
    }

    double avg_halos;
    double max_halos;
    /* + 1 to ensure avg_halos > 0 */
    MPIU_stats(comm, nhalos + 1, "->", &avg_halos, &max_halos);

    fastpm_info("Allocating %d halos per rank for final catalog.\n", (size_t) max_halos * 2);

    /* give it enough space for rebalancing. */
    fastpm_store_init(halos, NULL, (size_t) (max_halos * 2),
            attributes,
            FASTPM_MEMORY_HEAP);

    halos->np = nhalos;
    halos->meta = finder->p->meta;

    ptrdiff_t i;
    for(i = 0; i < halos->np; i++) {
        halos->id[i] = i;

        halos->mask[i] = 0; /* unselect the halos ; will turn this only if any particle is used */
        /* everthing should have been set to zero already by fastpm_store_init */
    }
}

void
fastpm_fof_execute(FastPMFOFFinder * finder, FastPMStore * halos)
{
    /* initial decompose -- reduce number of ghosts */
    FastPMStore * p = finder->p;
    PM * pm = finder->pm;
    MPI_Comm comm = finder->priv->comm;

    /* only do wrapping for periodic data */
    if(finder->priv->boxsize)
        fastpm_store_wrap(p, finder->priv->boxsize);

    double npmax, npmin, npstd, npmean;

    MPIU_stats(comm, p->np, "<->s", &npmin, &npmean, &npmax, &npstd);

    fastpm_info("load balance before fof decompose : min = %g max = %g mean = %g std = %g\n",
        npmin, npmax, npmean, npstd
        );

#if 1
    /* still route particles to the pm pencils as if they are periodic. */
    /* should still work (albeit use crazy memory) if we skip this. */
    if(0 != fastpm_store_decompose(p,
                (fastpm_store_target_func) FastPMTargetFOF, pm, comm)
    ) {
        fastpm_raise(-1, "out of storage space decomposing for FOF\n");
    }
#endif
    MPIU_stats(comm, p->np, "<->s", &npmin, &npmean, &npmax, &npstd);

    fastpm_info("load balance after fof decompose : min = %g max = %g mean = %g std = %g\n",
        npmin, npmax, npmean, npstd
        );

    /* create ghosts mesh size is usually > ll so we are OK here. */
    double below[3], above[3];

    int d;
    for(d = 0; d < 3; d ++) {
        /* bigger padding reduces number of iterations */
        below[d] = -finder->linkinglength * 1;
        above[d] = finder->linkinglength * 1;
    }

    PMGhostData * pgd = pm_ghosts_create_full(pm, p,
            COLUMN_POS | COLUMN_ID | COLUMN_MINID,
            below, above
        );

    pm_ghosts_send(pgd, COLUMN_POS);
    pm_ghosts_send(pgd, COLUMN_ID);

    size_t np_and_ghosts = p->np + pgd->p->np;

    #ifdef FASTPM_FOF_DEBUG
    fastpm_ilog(INFO, "Rank %d has %td particles including ghost\n", finder->priv->ThisTask, np_and_ghosts);
    #endif

    ptrdiff_t * head = fastpm_memory_alloc(p->mem, "FOFHead",
                    sizeof(head[0]) * np_and_ghosts, FASTPM_MEMORY_STACK);

    FastPMStore savebuff[1];
    fastpm_store_init(savebuff, p->name, np_and_ghosts, COLUMN_MINID, FASTPM_MEMORY_STACK);

    _fof_local_find(finder, p, pgd, head, finder->linkinglength);

    _fof_global_merge (finder, p, pgd, savebuff->minid, head);

    /* assign halo attr entries. This will keep only candidates that can possibly reach to nmin */
    size_t nsegments = _assign_halo_attr(finder, pgd, head, p->np, pgd->p->np, finder->nmin);

    fastpm_info("Found %td halos segments >= %d particles; or cross linked. \n", nsegments, finder->nmin);
    pm_ghosts_free(pgd);

    /* create local halos and modify head to index the local halos */
    fastpm_fof_create_local_halos(finder, halos, nsegments);
    /* remove halos without any local particles */
    fastpm_fof_remove_empty_halos(finder, halos, savebuff->minid, head);

    fastpm_store_destroy(savebuff);

    /* reduce the primary halo attrs */
    fastpm_fof_compute_halo_attrs(finder, halos, head, _convert_basic_halo_attrs, _add_basic_halo_attrs, _reduce_basic_halo_attrs);

    #ifdef FASTPM_FOF_DEBUG
    {
        int i;
        for(i  = 0; i < halos->np; i ++) {
            fastpm_ilog(INFO, "Task = %d, Halo[%d] = %d mask=%d MINID=%ld\n", finder->priv->ThisTask, i, halos->length[i], halos->mask[i], halos->minid[i]);
        }
    }
    #endif

    /* apply length cut */
    fastpm_fof_apply_length_cut(finder, halos, head);

    /* reduce the primary halo attrs */
    fastpm_fof_compute_halo_attrs(finder, halos, head, _convert_extended_halo_attrs, _add_extended_halo_attrs, _reduce_extended_halo_attrs);

    /* the event is called with full halos, only those where mask==1 are primary
     * the others are ghosts with the correct properties but shall not show up in the
     * catalog.
     * */
    FastPMHaloEvent event[1];
    event->halos = halos;
    event->p = finder->p;
    event->ihalo = head;

    fastpm_emit_event(finder->event_handlers, FASTPM_EVENT_HALO,
                    FASTPM_EVENT_STAGE_AFTER, (FastPMEvent*) event, finder);

    fastpm_memory_free(finder->p->mem, head);

    fastpm_store_subsample(halos, halos->mask, halos);

    fastpm_info("After event: %td halos.\n", fastpm_store_get_np_total(halos, comm));
}

void
fastpm_fof_destroy(FastPMFOFFinder * finder)
{
    fastpm_destroy_event_handlers(&finder->event_handlers);
    free(finder->priv);
}

ext: c

max_stars_repo_name: sbird/FastPMRunner

max_issues_repo_name: sbird/FastPMRunner

max_forks_repo_name: sbird/FastPMRunner


text:
#pragma once

#include <gsl/gsl>
#include <string>
#include <memory>
#include <vector>
#include <string_view>

#include "aas/aas_math.h"

namespace aas {

	struct MeshVertex {
		static constexpr size_t sMaxBoneAttachments = 6;
		DX::XMFLOAT4 pos;
		DX::XMFLOAT4 normal;
		DX::XMFLOAT2 uv;
		uint16_t padding;
		uint16_t attachmentCount;
		uint16_t attachmentBone[sMaxBoneAttachments];
		float attachmentWeight[sMaxBoneAttachments];

		bool IsAttachedTo(int bone_idx) const {
			for (int i = 0; i < attachmentCount; i++) {
				if (attachmentBone[i] == bone_idx) {
					return true;
				}
			}
			return false;
		}
	};

	struct MeshFace {
		static constexpr size_t sVertexCount = 3;
		int16_t material_idx;
		int16_t vertices[sVertexCount];
	};

	struct MeshBone {
		int16_t flags;
		int16_t parent_id;
		char name[48];
		Matrix3x4 full_world_inverse;
	};

	class Mesh {
	public:
		Mesh(std::string filename, std::vector<uint8_t> data);

		const std::vector<std::string_view> &GetMaterials() const {
			return materials_;
		}

		std::string_view GetMaterial(int materialIdx) const {
			return materials_[materialIdx];
		}

		gsl::span<MeshBone> GetBones() const {
			return bones_;
		}

		const MeshBone &GetBone(int boneIdx) const {
			Expects(boneIdx >= 0 && boneIdx < bones_.size());
			return bones_[boneIdx];
		}

		gsl::span<MeshVertex> GetVertices() const {
			return vertices_;
		}

		const MeshVertex &GetVertex(int vertexIdx) const {
			return vertices_[vertexIdx];
		}

		const size_t GetVertexCount() const {
			return vertices_.size();
		}

		MeshVertex &GetVertex(int vertexIdx) {
			return vertices_[vertexIdx];
		}

		gsl::span<MeshFace> GetFaces() const {
			return faces_;
		}

		const MeshFace &GetFace(int faceIdx) const {
			return faces_[faceIdx];
		}

		void RenormalizeClothVertices(int clothBoneId);

		const std::string &GetFilename() const {
			return filename_;
		}

	private:
		std::string filename_;
		std::vector<uint8_t> data_;
		std::vector<std::string_view> materials_;

		// These are views into the data buffer
		gsl::span<MeshBone> bones_;
		gsl::span<MeshVertex> vertices_;
		gsl::span<MeshFace> faces_;
	};

	std::unique_ptr<Mesh> LoadMeshFile(std::string_view filename);

}

max_stars_repo_name: edoipi/TemplePlus

max_issues_repo_name: edoipi/TemplePlus

max_forks_repo_name: edoipi/TemplePlus


text:
#if !defined(fvSupport_h)
#define fvSupport_h

#include <petsc.h>

#define MAX_FVM_RHS_FUNCTION_FIELDS 4

typedef PetscErrorCode (*FVMRHSFluxFunction)(PetscInt dim, const PetscFVFaceGeom *fg, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar fieldL[], const PetscScalar fieldR[],
                                             const PetscScalar gradL[], const PetscScalar gradR[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar auxL[], const PetscScalar auxR[],
                                             const PetscScalar gradAuxL[], const PetscScalar gradAuxR[], PetscScalar flux[], void *ctx);

typedef PetscErrorCode (*FVMRHSPointFunction)(PetscInt dim, const PetscFVCellGeom *cg, const PetscInt uOff[], const PetscScalar u[], const PetscInt aOff[], const PetscScalar a[], PetscScalar f[], void *ctx);

typedef PetscErrorCode (*FVAuxFieldUpdateFunction)(PetscReal time, PetscInt dim, const PetscFVCellGeom *cellGeom, const PetscScalar *conservedValues, PetscScalar *auxField, void *ctx);

/**
 * struct to describe how to compute RHS finite volume flux source terms
 */
struct _FVMRHSFluxFunctionDescription {
    FVMRHSFluxFunction function;
    void *context;

    PetscInt field;
    PetscInt inputFields[MAX_FVM_RHS_FUNCTION_FIELDS];
    PetscInt numberInputFields;

    PetscInt auxFields[MAX_FVM_RHS_FUNCTION_FIELDS];
    PetscInt numberAuxFields;
};

typedef struct _FVMRHSFluxFunctionDescription FVMRHSFluxFunctionDescription;

/**
 * struct to describe how to compute RHS finite volume point source terms
 */
struct _FVMRHSPointFunctionDescription {
    FVMRHSPointFunction function;
    void *context;

    PetscInt fields[MAX_FVM_RHS_FUNCTION_FIELDS];
    PetscInt numberFields;

    PetscInt inputFields[MAX_FVM_RHS_FUNCTION_FIELDS];
    PetscInt numberInputFields;

    PetscInt auxFields[MAX_FVM_RHS_FUNCTION_FIELDS];
    PetscInt numberAuxFields;
};

typedef struct _FVMRHSPointFunctionDescription FVMRHSPointFunctionDescription;

/**
  DMPlexTSComputeRHSFunctionFVM - Form the local forcing F from the local input X using flux and pointfunctions specified by the user

  Input Parameters:
+ dm - The mesh
. t - The time
. locX  - Local solution
- user - The user context

  Output Parameter:
. F  - Global output vector

  Level: developer

.seealso: DMPlexComputeJacobianActionFEM()
**/
PETSC_EXTERN PetscErrorCode ABLATE_DMPlexComputeRHSFunctionFVM(FVMRHSFluxFunctionDescription *fluxFunctionDescriptions, PetscInt numberFluxFunctionDescription,
                                                               FVMRHSPointFunctionDescription *pointFunctionDescriptions, PetscInt numberPointFunctionDescription, DM dm, PetscReal time, Vec locX, Vec F);

/**
 * Populate the boundary with gradient information
 * @param dm
 * @param auxDM
 * @param time
 * @param locXVec
 * @param locAuxField
 * @param numberUpdateFunctions
 * @param updateFunctions
 * @param data
 * @return
 */

/**
 * Takes all local vector
 * @return
 */
PETSC_EXTERN PetscErrorCode ABLATE_DMPlexComputeFluxResidual_Internal(FVMRHSFluxFunctionDescription functionDescription[], PetscInt numberFunctionDescription, DM, IS, PetscReal, Vec, Vec, PetscReal, Vec);

/**
   Form the local forcing F from the local input X using pointwise functions specified by the user

  Input Parameters:
+ dm - The mesh
. t - The time
. locX  - Local solution

  Output Parameter:
. F  - local output vector

**/
PETSC_EXTERN PetscErrorCode ABLATE_DMPlexComputePointResidual_Internal(FVMRHSPointFunctionDescription *functionDescription, PetscInt numberFunctionDescription, DM, IS, PetscReal, Vec, Vec, PetscReal, Vec);

/**
 * reproduces the petsc call with grad fixes for multiple fields
 * @param dm
 * @param fvm
 * @param locX
 * @param grad
 * @return
 */
PETSC_EXTERN PetscErrorCode DMPlexReconstructGradientsFVM_MulfiField(DM dm, PetscFV fvm, Vec locX, Vec grad);

/**
 * reproduces the petsc call with grad fixes for multiple fields
 * @param dm
 * @param fvm
 * @param locX
 * @param grad
 * @return
 */
PETSC_EXTERN PetscErrorCode DMPlexGetDataFVM_MulfiField(DM dm, PetscFV fv, Vec *cellgeom, Vec *facegeom, DM *gradDM);

/**
 * Function to update all cells.  This should be merged into other update calls
 * @param dm
 * @param auxDM
 * @param time
 * @param locXVec
 * @param locAuxField
 * @param numberUpdateFunctions
 * @param updateFunctions
 * @param ctx
 * @return
 */
PETSC_EXTERN PetscErrorCode FVFlowUpdateAuxFieldsFV(DM dm, DM auxDM, PetscReal time, Vec locXVec, Vec locAuxField, PetscInt numberUpdateFunctions, FVAuxFieldUpdateFunction *updateFunctions, void **ctx);

#endif

max_stars_repo_name: mschulwitz/ablate

max_issues_repo_name: mschulwitz/ablate

max_forks_repo_name: mschulwitz/ablate


text:
#pragma once
#include "GradUtil.h"

#ifndef _NOGSL
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>
#else
#include "CustomSolver.h"
#endif

#include "BooleanDAG.h"

class ActualEvaluators;
class SmoothEvaluators;

class Predicate {
public:
    virtual string print() = 0;
    virtual double evaluate(ActualEvaluators* eval) = 0;
    virtual double evaluate(gsl_vector* grad, ActualEvaluators* eval) = 0;
    virtual double evaluate(SmoothEvaluators* eval) = 0;
    virtual double evaluate(gsl_vector* grad, SmoothEvaluators* eval) = 0;
    virtual double grad_norm(ActualEvaluators* eval) = 0;
    virtual const gsl_vector* grad(ActualEvaluators* eval) = 0;
    virtual bool isPure() = 0;
    virtual void makeImpure() = 0;
    virtual bool isBasic() = 0;
    virtual bool isDiff() = 0;
};

class BasicPredicate: public Predicate { // TODO: override equal
    bool pure;
public:
    int nid;
    BooleanDAG* dag;
    Predicate* oldPred;

    BasicPredicate(BooleanDAG* _dag, int _nid): dag(_dag), nid(_nid), oldPred(NULL) { pure = true;}
    BasicPredicate(BooleanDAG* _dag, int _nid, Predicate* _oldPred): dag(_dag), nid(_nid), oldPred(_oldPred) { pure = true;}
    virtual void makeImpure() {
        pure = false;
    }
    virtual bool isPure() {
        return pure;
    }
    virtual bool isBasic() { return true; }
    virtual bool isDiff() { return false; }
    virtual string print() {
        stringstream s;
        //s << (*dag)[nid]->lprint();
        s << nid;
        if (oldPred != NULL) {
            s << "(" << oldPred->print() << ")";
        }
        return s.str();
    }
    virtual double evaluate(ActualEvaluators* eval);
    virtual double evaluate(gsl_vector* grad, ActualEvaluators* eval);
    virtual double evaluate(SmoothEvaluators* eval);
    virtual double evaluate(gsl_vector* grad, SmoothEvaluators* eval);
    virtual double grad_norm(ActualEvaluators* eval);
    virtual const gsl_vector* grad(ActualEvaluators* eval);

};

class DiffPredicate: public Predicate {
    bool pure;
public:
    Predicate* p1;
    Predicate* p2;
    gsl_vector* tmp_p1; // TODO: Not sure where is the best place to create these tmp vectors for gradient calculation
    gsl_vector* tmp_p2;
    DiffPredicate(Predicate* _p1, Predicate* _p2, int ncontrols): p1(_p1), p2(_p2) {
        tmp_p1 = gsl_vector_alloc(ncontrols);
        tmp_p2 = gsl_vector_alloc(ncontrols);
        pure = true;
    }
    ~DiffPredicate() {
        gsl_vector_free(tmp_p1);
        gsl_vector_free(tmp_p2);
    }
    virtual void makeImpure() {
        pure = false;
    }
    virtual bool isPure() {
        return pure;
    }
    virtual bool isBasic() { return false; }
    virtual bool isDiff() { return true; }
    virtual string print() {
        stringstream s;
        s << p1->print() << " - " << p2->print() << endl;
        return s.str();
    }
    virtual double evaluate(ActualEvaluators* eval);
    virtual double evaluate(gsl_vector* grad, ActualEvaluators* eval);
    virtual double evaluate(SmoothEvaluators* eval);
    virtual double evaluate(gsl_vector* grad, SmoothEvaluators* eval);
    virtual double grad_norm(ActualEvaluators* eval);
    virtual const gsl_vector* grad(ActualEvaluators* eval);

};

max_stars_repo_name: natebragg/sketch-backend

max_issues_repo_name: natebragg/sketch-backend

max_forks_repo_name: natebragg/sketch-backend


text:
#pragma once

#include <gsl\gsl>
#include <winrt\base.h>
#include <d3d11.h>
#include "DrawableGameComponent.h"
#include "PointSpriteMaterial.h"

namespace Rendering
{
	class PointSpriteDemo final : public Library::DrawableGameComponent
	{
	public:
		PointSpriteDemo(Library::Game& game, const std::shared_ptr<Library::Camera>& camera);
		PointSpriteDemo(const PointSpriteDemo&) = delete;
		PointSpriteDemo(PointSpriteDemo&&) = default;
		PointSpriteDemo& operator=(const PointSpriteDemo&) = default;
		PointSpriteDemo& operator=(PointSpriteDemo&&) = default;
		~PointSpriteDemo();

		virtual void Initialize() override;
		virtual void Draw(const Library::GameTime& gameTime) override;

	private:
		void InitializeRandomPoints();

		PointSpriteMaterial mMaterial;
		winrt::com_ptr<ID3D11Buffer> mVertexBuffer;
		std::size_t mVertexCount{ 0 };
		bool mUpdateMaterial{ true };
	};
}


text:
/**
  * Copyright 2019 José Manuel Abuín Mosquera <josemanuel.abuin@usc.es>
  *
  * This file is part of Matrix Market Suite.
  *
  * Matrix Market Suite 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.
  *
  * Matrix Market Suite 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 Matrix Market Suite. If not, see <http://www.gnu.org/licenses/>.
  */

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

#include <cblas.h>

#include "DMxV-MPI.h"

void usageDMxVMPI(){

	fprintf(stderr, "\n");
	fprintf(stderr, "Usage: MM-Suite DMxV [options] <input-matrix> <input-vector>\n");
	fprintf(stderr, "\nInput/output options:\n\n");
	fprintf(stderr, "       -o STR        Output file name. Default: stdout\n");
	fprintf(stderr, "       -r            Input format is row per line. Default: False\n");
	fprintf(stderr, "\nParameters options:\n\n");
	fprintf(stderr, "       -a DOUBLE     Alpha. Default: 1.0\n");
	fprintf(stderr, "       -b DOUBLE     Beta. Default: 0.0\n");
	fprintf(stderr, "\nPerformance options:\n\n");
	fprintf(stderr, "       -t INT        Number of threads to use in OpenBLAS. Default: 1\n");
	fprintf(stderr, "\n");

}

int DMxVMPI(int argc, char *argv[], int numProcs, int myid) {

	int 			ret_code = 1;
	int			option;

	unsigned long 		*II;
	unsigned long 		*J;
	double 			*values;

	unsigned long 		M;
	unsigned long 		local_M;
	unsigned long 		N;
	unsigned long 		long nz;

	double 			*vectorValues;
	unsigned long 		M_Vector;
	unsigned long 		N_Vector;
	unsigned long long 	nz_vector;

	char			*outputFileName = NULL;

	char			*inputMatrixFile = NULL;
	char			*inputVectorFile = NULL;
	char			*outputVectorFile = NULL;

	int			inputFormatRow = 0;

	double			alpha = 1.0;
	double			beta = 0.0;

	int 		numThreads = 1;

	int			i, j;

	while ((option = getopt(argc, argv,"ro:a:b:t:")) >= 0) {
		switch (option) {
			case 'o' :
				//free(outputFileName);

				outputFileName = (char *) malloc(sizeof(char)*strlen(optarg)+1);
				strcpy(outputFileName,optarg);

				break;

			case 'r':
				inputFormatRow = 1;
				break;

			case 'b':
				beta = atof(optarg);
				break;

			case 'a':
				alpha = atof(optarg);
				break;

			case 't':
				numThreads = atoi(optarg);
				break;

			default: break;
		}

	}

	if ((optind + 3 != argc) && (optind + 2 != argc)) {
		if (myid == 0) {
			//fprintf(stderr,"[%s] Argc: %d, optind: %d\n",__func__, argc, optind);
			usageDMxVMPI();
		}
		return 0;
	}

	openblas_set_num_threads(numThreads);

	if(optind + 3 == argc) { //We have an output vector

		outputVectorFile = (char *)malloc(sizeof(char)*strlen(argv[optind+2])+1);
		strcpy(outputVectorFile,argv[optind+2]);
	}

	if(outputFileName == NULL) {
		outputFileName = (char *) malloc(sizeof(char)*6);
		sprintf(outputFileName,"stdout");
	}

	inputMatrixFile = (char *)malloc(sizeof(char)*strlen(argv[optind])+1);
	inputVectorFile = (char *)malloc(sizeof(char)*strlen(argv[optind+1])+1);

	strcpy(inputMatrixFile,argv[optind]);
	strcpy(inputVectorFile,argv[optind+1]);

	//Read matrix
	if(inputFormatRow) {
		if(!readDenseCoordinateMatrixMPIRowLine(inputMatrixFile,&II,&J,&values,&M,&local_M,&N,&nz,myid, numProcs)){
			fprintf(stderr, "[%s] Can not read Matrix\n",__func__);
			return 0;
		}
	}
	else{
		if(!readDenseCoordinateMatrixMPI(inputMatrixFile,&II,&J,&values,&M,&local_M,&N,&nz,myid, numProcs)){
			fprintf(stderr, "[%s] Can not read Matrix\n",__func__);
			return 0;
		}
	}

	//Read input vector
	if(!readDenseVector(inputVectorFile, &vectorValues,&M_Vector,&N_Vector,&nz_vector)){
		fprintf(stderr, "[%s] Can not read Vector\n",__func__);
		return 0;
	}

	/*
	void cblas_dgemv(const enum CBLAS_ORDER order,
                 const enum CBLAS_TRANSPOSE TransA, const int M, const int N,
                 const double alpha, const double *A, const int lda,
                 const double *X, const int incX, const double beta,
                 double *Y, const int incY);
                 */

        double *partial_result=(double *) calloc(local_M,sizeof(double));
        double* y = (double*)calloc(N,sizeof(double));

        if(y == NULL){
		fprintf(stderr,"[%s] Error reserving memory for final result vector in processor %d\n",__func__,myid);
		return 0;
	}

        //Read output vector if any
	if(outputVectorFile != NULL) {
		if(!readDenseVector(outputVectorFile, &y,&M_Vector,&N_Vector,&nz_vector)){
			fprintf(stderr, "[%s] Can not read Vector %s\n",__func__, outputVectorFile);
			return 0;
		}

		for( i = (local_M * myid), j = 0; i< (local_M * myid + local_M) && j< local_M; i++, j++) {
        		partial_result[j] = y [i];
        	}
	}

        double t_real = realtime();

        //y := alpha * A * x + beta * y

	//cblas_dgemv(CblasRowMajor,CblasNoTrans,local_M,N,1.0,values,N,vectorValues,1,0.0,partial_result,1);
	cblas_dgemv(CblasRowMajor,CblasNoTrans,local_M,N,alpha,values,N,vectorValues,1,beta,partial_result,1);

	MPI_Allgather (partial_result,local_M,MPI_DOUBLE,y,local_M,MPI_DOUBLE,MPI_COMM_WORLD);
	MPI_Barrier(MPI_COMM_WORLD);

	fprintf(stderr, "\n[%s] Time spent in DMxV: %.6f sec\n", __func__, realtime() - t_real);

	if (myid == 0){

		writeDenseVector(outputFileName, y,M_Vector,N_Vector,nz_vector);

	}
	return ret_code;
}


text:
//
//  global.h
//  Agile_RRT
//
//  Created by Timothy M. Caldwell.
//  Copyright (c) 2015 Timothy M. Caldwell. Refer to license.txt
//

#ifndef __Agile_RRT__global__
#define __Agile_RRT__global__

#include <fstream>
//#include <string>
//#include <sstream>

#include <iostream>
#include <iomanip>

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

//boost libraries
#include <boost/numeric/odeint.hpp>

//gsl libraries
#include <gsl/gsl_spline.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_errno.h>

//Eigen libraries
#include <eigen3/Eigen/Dense>
#include <eigen3/Eigen/Core>
#include <eigen3/Eigen/Geometry>

//namespaces
using namespace std;
namespace odeint = boost::numeric::odeint;

extern const int NI;                              // Constant; Number of inputs to the dynamic system
extern const int NS;                              // Constant; Number of states to the dynamic system

extern const double kLARGENUM;                    // Constant; A sufficiently large number

// Constants; For integration using odeint
extern const int kMAX_NUM_INTEGRATION_STEPS;
extern const double kINTEGRATION_DT_START;

//integrate state types
extern const int kFREEDYNAMICS_TYPE;
extern const int kFEEDFORWARD_TYPE;
extern const int kTRACKPOINT_TYPE;
extern const int kTRACKTRAJECTORY_TYPE;
extern const int kLINEARSTEERINGPROJECTION_TYPE;
extern const int kLINEARSTEERINGPROJECTIONWCOST_TYPE;

//integrate linear state and ricatti types
extern const int kLTI_TYPE;
extern const int kLTV_TYPE;

//STMSlot2 types
extern const int kLTI_OPENLOOP_TYPE;
extern const int kLTI_CLOSEDLOOP_TYPE;
extern const int kLTV_OPENLOOP_TYPE;
extern const int kLTV_CLOSEDLOOP_TYPE;

//control type for linear state equation
extern const int kFREE_CONTROLTYPE;
extern const int kFEEDFORWARD_CONTROLTYPE;
extern const int kMINENERGY_OPENLOOP_CONTROLTYPE;
extern const int kMINENERGY_CLOSEDLOOP_CONTROLTYPE;

//Gramian types
extern const int kWW_GRAMIANTYPE;
extern const int kWWK_GRAMIANTYPE;
extern const int kSSK_GRAMIANTYPE;

typedef vector<double> ode_state_type;            // odeint integrator state_type

///////////////////////////////
//Currently implemented as global. Needs to be updated for NS and NI
///////////////////////////////
// Eigen Matrix types of size NS or NI by NS or NI
typedef Eigen::Matrix< double , 8 , 8 > NSxNS_type;
typedef Eigen::Matrix< double , 8 , 1 > NSxNI_type;
typedef Eigen::Matrix< double , 1 , 1 > NIxNI_type;
typedef Eigen::Matrix< double , 1 , 8 > NIxNS_type;
typedef Eigen::Matrix< double , 8 , 1 > NSx1_type;
typedef Eigen::Matrix< double , 1 , 1 > NIx1_type;

struct constraints_struct {
  Eigen::Matrix<double, 8, 2> xbounds;              // list of upper and lower bounds for each state
  Eigen::Matrix<double, 1, 2> ubounds;              // list of upper and lower bounds for each input

  // list of circle obstacles of radius obs_radii[.] at point (obs_Xs[.], obs_Ys[.])
  vector< double > obs_radii;
  vector< double > obs_Xs;
  vector< double > obs_Ys;
};
// list of upper and lower bounds for each state and input. Should be the same as those set in constraints_struct
struct sampling_bounds_struct {
  Eigen::Matrix<double, 8, 2> xbounds;
  Eigen::Matrix<double, 1, 2> ubounds;
};

// stepper type used by odeint for integration
//typedef odeint::controlled_runge_kutta< odeint::runge_kutta_dopri5<ode_state_type> > stepper_type; // Alternative stepper choice
typedef odeint::controlled_runge_kutta< odeint::runge_kutta_cash_karp54<ode_state_type> > stepper_type;

// Printing parameters
extern const int kPRINT_PRECISION;                // Printing precision number
extern const double kPRINT_CHOP_BOUND;            // Lower bound for chopping (i.e. if abs(x)<kPRINT_CHOP_BOUND, then 0 is returned)

// prints a vec to stream in a form easily readable by Mathematica.
extern void PrintVec(const vector<double> & vec, ostream * stream = &(std::cout));
extern void PrintVec(const vector<double> & vec, const string & name, ostream * stream);

#endif // __Agile_RRT__global__