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 <mpi.h>
#include <sys/stat.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

#include "dnest.h"
#include "dnestvars.h"

/*
 * dnest
 * call this function to do sampling
 * ========================================================
 * arguments:
 * argc:        number of command-line options (mandatory)
 * argv:        command-line options  (mandatory)
 * fptrset:     function set pointers required  (mandatory)
 * num_params:  number of parameters (mandatory)
 * param_range: parameter ranges  (mandatory)
 * prior_type:  prior types  (optional)
 * prior_info:  prior informations  (optional)
 * sample_dir:  output directory for sampling  (mandatory)
 * optfile:     option files  (optional)
 * opts:        options struct. if optfile is empty, use this struct (optional)
 * args:        any other arguments transferred to cdnest.  (optional)
 *              useful when referring to external variables or calling external functions (optional)
 * ========================================================
 * optional arguments can be set to "NULL"
 */
double dnest(int argc, char** argv, DNestFptrSet *fptrset, int num_params,
             double *param_range, int *prior_type, double *prior_info,
             char *sample_dir, char *optfile, DNestOptions *opts, void *args)
{
  int optid;

  MPI_Comm_rank(MPI_COMM_WORLD, &dnest_thistask);
  MPI_Comm_size(MPI_COMM_WORLD, &dnest_totaltask);

  if(dnest_thistask == dnest_root)
  {
    printf("#=======================================================\n");
    printf("# Starting CDNest.\n");
    printf("# Use %d cores.\n", dnest_totaltask);
  }

  dnest_check_fptrset(fptrset);

  // cope with argv
  if(dnest_thistask == dnest_root )
  {
    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( (optid = getopt(argc, argv, "r:s:pt:clx:g:")) != -1)
    {
      switch(optid)
      {
        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;
      }
    }
  }

  MPI_Bcast(&dnest_flag_restart, 1, MPI_INT, dnest_root, MPI_COMM_WORLD);
  MPI_Bcast(&dnest_flag_postprc, 1, MPI_INT, dnest_root, MPI_COMM_WORLD);
  MPI_Bcast(&dnest_flag_sample_info, 1, MPI_INT,dnest_root, MPI_COMM_WORLD);
  MPI_Bcast(&dnest_post_temp, 1, MPI_DOUBLE, dnest_root, MPI_COMM_WORLD);
  MPI_Bcast(&dnest_flag_limits, 1, MPI_INT, dnest_root, MPI_COMM_WORLD);

  setup(argc, argv, fptrset, num_params, param_range, prior_type, prior_info, sample_dir, optfile, opts, args);

  if(dnest_flag_postprc == 1)
  {
    dnest_postprocess(dnest_post_temp, optfile, opts);
    MPI_Barrier(MPI_COMM_WORLD);
    finalise();
    return post_logz;
  }

  if(dnest_flag_sample_info == 1)
  {
    dnest_postprocess(dnest_post_temp, optfile, opts);
    finalise();
    return post_logz;
  }

  if(dnest_flag_restart==1)
    dnest_restart();

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

  dnest_postprocess(dnest_post_temp, optfile, opts);

  finalise();

  return post_logz;
}

// postprocess, calculate evidence, generate posterior sample.
void dnest_postprocess(double temperature,char *optfile, DNestOptions *opts)
{
  if(dnest_thistask == dnest_root)
  {
    options_load(optfile, opts);
    postprocess(temperature);
  }
  MPI_Bcast(&post_logz, 1, MPI_DOUBLE, dnest_root, MPI_COMM_WORLD);
}

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

  // used to gather levels' information
  if(dnest_thistask == dnest_root)
  {
    buf_size_above = malloc(dnest_totaltask * sizeof(int));
    buf_displs = malloc(dnest_totaltask * sizeof(int));
  }

  if(dnest_thistask == dnest_root)
  {
    printf("#=======================================================\n");
    printf("# Starting diffusive nested sampling.\n");
  }
  MPI_Barrier(MPI_COMM_WORLD);

  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();
    MPI_Barrier(MPI_COMM_WORLD);

    //gather levels
    MPI_Gather(levels, size_levels*sizeof(Level), MPI_BYTE,
             copies_of_levels, size_levels*sizeof(Level), MPI_BYTE, dnest_root, MPI_COMM_WORLD);

    //gather limits
    if(dnest_flag_limits == 1)
    {
      MPI_Gather(limits, size_levels*particle_offset_double*2, MPI_DOUBLE,
               copies_of_limits, size_levels*particle_offset_double*2, MPI_DOUBLE, dnest_root, MPI_COMM_WORLD );
    }

    //gather size_above
    MPI_Gather(&size_above, 1, MPI_INT, buf_size_above, 1, MPI_INT, dnest_root, MPI_COMM_WORLD);

    // task 0 responsible for updating levels
    if(dnest_thistask == dnest_root)
    {
      size_all_above_incr = 0;
      for(i = 0; i<dnest_totaltask; i++)
      {
        size_all_above_incr += buf_size_above[i];

        buf_size_above[i] *= sizeof(LikelihoodType);
      }

      // store new points following the end of all_above array.
      buf_displs[0] = size_all_above * sizeof(LikelihoodType);
      for(i=1; i<dnest_totaltask; i++)
      {
        buf_displs[i] = buf_displs[i-1] + buf_size_above[i-1];
      }

      //update size_all_above
      size_all_above += size_all_above_incr;

      if(size_all_above > options.new_level_interval*2)
      {
        printf("# Error, all above overflow.\n");
        exit(0);
      }
    }

    // gather above into all_above, stored in task 0, note that its size is different among tasks
    MPI_Gatherv(above, size_above * sizeof(LikelihoodType), MPI_BYTE,
                all_above, buf_size_above, buf_displs, MPI_BYTE, dnest_root, MPI_COMM_WORLD);

    // reset size_above for each task
    size_above = 0;

    count_mcmc_steps += options.thread_steps * dnest_totaltask;

    if(dnest_thistask == dnest_root)
    {
      //backup levels_combine
      levels_orig = malloc(size_levels_combine * sizeof(Level));
      memcpy(levels_orig, levels_combine, size_levels_combine*sizeof(Level));

      //scan over all copies of levels
      pl = copies_of_levels;
      for(i=0; i< dnest_totaltask; i++)
      {
        for(j=0; j<size_levels_combine; j++)
        {
          levels_combine[j].accepts += (pl[j].accepts - levels_orig[j].accepts);
          levels_combine[j].tries += (pl[j].tries - levels_orig[j].tries);
          levels_combine[j].visits += (pl[j].visits - levels_orig[j].visits);
          levels_combine[j].exceeds += (pl[j].exceeds - levels_orig[j].exceeds);
          //printf("%d %d\n", thistask, (pl[j].accepts - levels_orig[j].accepts) );;
        }
        pl += size_levels_combine;
      }

      free(levels_orig);

      // scan over all copies of limits
      if(dnest_flag_limits == 1)
      {
        plimits = copies_of_limits;
        for(i=0; i< dnest_totaltask; i++)
        {
          for(j=0; j < size_levels; 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 ],
                plimits[j * 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 ],
                plimits[j * particle_offset_double *2 + k*2 + 1]);
            }

          }
          plimits += size_levels * particle_offset_double * 2;
        }

        // 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();

      size_levels = size_levels_combine;
      memcpy(levels, levels_combine, size_levels * sizeof(Level));
    }

    //broadcast levels
    MPI_Bcast(&size_levels, 1, MPI_INT, dnest_root, MPI_COMM_WORLD);
    //MPI_Bcast(&count_saves, 1, MPI_INT, root, MPI_COMM_WORLD);
    MPI_Bcast(levels, size_levels * sizeof(Level), MPI_BYTE, dnest_root,  MPI_COMM_WORLD);

    if(dnest_flag_limits == 1)
      MPI_Bcast(limits, size_levels * particle_offset_double *2, MPI_DOUBLE, dnest_root, MPI_COMM_WORLD);

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

      if(dnest_thistask == dnest_root )
      {
        // save levels, limits, sync samples when running a number of steps
        if( count_saves % num_saves == 0 )
        {
          save_levels();
          if(dnest_flag_limits == 1)
            save_limits();
          fflush(fsample_info);
          fsync(fileno(fsample_info));
          fflush(fsample);
          fsync(fileno(fsample));
          printf("# Save levels, limits, and sync samples at N= %d.\n", count_saves);
        }
      }

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

  //dnest_save_restart();

  if(dnest_thistask == dnest_root)
  {
    //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);

    free(buf_size_above);
    free(buf_displs);
  }
}

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

  if(!enough_levels(levels_combine, size_levels_combine) && size_all_above >= options.new_level_interval)
  {
    // in descending order
    qsort(all_above, size_all_above, sizeof(LikelihoodType), dnest_cmp);
    int index = (int)( (1.0/compression) * size_all_above);

    Level level_tmp = {all_above[index], 0.0, 0, 0, 0, 0};
    levels_combine[size_levels_combine] = level_tmp;
    size_levels_combine++;

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

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

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

    //created_level = true;
  }

  recalculate_log_X();

/*  if(created_level)
    save_levels();

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

    //save_particle();

    if(!created_level)
      save_levels();
  }*/

}

void recalculate_log_X()
{
  int i;

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

void renormalise_visits()
{
  size_t i;

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

    if(levels_combine[i].visits >= regularisation)
    {
      levels_combine[i].exceeds = ( (double) (levels_combine[i].exceeds+1) / (double)(levels_combine[i].visits + 1) ) * regularisation;
      levels_combine[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_combine; 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);
  }
  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_combine; 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(dnest_thistask == dnest_root)
  {
    if(count_saves%1 == 0)
      printf("#[%.1f%%] Saving particle to disk. N= %d.\n", 100.0*count_saves/options.max_num_saves, count_saves);

    whichtask = gsl_rng_uniform_int(dnest_gsl_r,dnest_totaltask);
  }

  MPI_Bcast(&whichtask, 1, MPI_INT, dnest_root, MPI_COMM_WORLD);

  if(whichtask != dnest_root)
  {
    if(dnest_thistask == whichtask)
    {
      int size_message = dnest_size_of_modeltype + 2*sizeof(int) + 2*sizeof(double);
      particle_message = (void *)malloc(size_message);
      whichparticle = gsl_rng_uniform_int(dnest_gsl_r,options.num_particles);
      memcpy(particle_message, particles + whichparticle * particle_offset_size, dnest_size_of_modeltype);
      memcpy(particle_message + dnest_size_of_modeltype, &log_likelihoods[whichparticle].value, 2*sizeof(double));
      memcpy(particle_message + dnest_size_of_modeltype + 2*sizeof(double), &level_assignments[whichparticle], sizeof(int));
      memcpy(particle_message + dnest_size_of_modeltype + 2*sizeof(double) + sizeof(int), &whichparticle, sizeof(int));

      MPI_Send(particle_message, size_message, MPI_BYTE, dnest_root, 1, MPI_COMM_WORLD);

      //printf("%f %f\n", log_likelihoods[whichparticle].value, log_likelihoods[whichparticle].tiebreaker);

      free(particle_message);
    }
    if(dnest_thistask == dnest_root)
    {
      MPI_Status status;
      int size_message = dnest_size_of_modeltype + 2*sizeof(int) + 2*sizeof(double);
      int whichlevel;
      LikelihoodType logl;

      particle_message = (void *)malloc(size_message);

      MPI_Recv(particle_message, size_message, MPI_BYTE, whichtask, 1, MPI_COMM_WORLD, &status);

      memcpy(&logl, particle_message + dnest_size_of_modeltype, 2*sizeof(double) );
      memcpy(&whichlevel, particle_message + dnest_size_of_modeltype + 2*sizeof(double), sizeof(int) );
      memcpy(&whichparticle, particle_message + dnest_size_of_modeltype + 2*sizeof(double) + sizeof(int), sizeof(int) );

      //printf("%f %f\n", logl.value, logl.tiebreaker);

      print_particle(fsample, particle_message);

      fprintf(fsample_info, "%d %e %f %d\n", whichlevel,
        logl.value,
        logl.tiebreaker,
        whichtask * options.num_particles + whichparticle);

      free(particle_message);
    }
  }
  else
  {
    if(dnest_thistask == dnest_root)
    {
      whichparticle =  gsl_rng_uniform_int(dnest_gsl_r,options.num_particles);

      print_particle(fsample, particles + whichparticle * particle_offset_size);

      fprintf(fsample_info, "%d %e %f %d\n", level_assignments[whichparticle],
        log_likelihoods[whichparticle].value,
        log_likelihoods[whichparticle].tiebreaker,
        whichtask * options.num_particles + 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);

  logl_proposal.value = log_likelihoods_cal(proposal);
  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]);

  // enforce uniform exploration if levels are enough
  if(enough_levels(levels, size_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.lam;
}

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_thistask != dnest_root)
    return;

  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()
{
  if(dnest_thistask != dnest_root )
    return;

  fclose(fsample);
  fclose(fsample_info);
}

void setup(int argc, char** argv, DNestFptrSet *fptrset, int num_params,
           double *param_range, int *prior_type, double *prior_info,
           char *sample_dir, char *optfile, DNestOptions *opts, void *args)
{
  int i, j;

  // root task.
  dnest_root = 0;

  if(dnest_thistask == dnest_root)
  {
    dnest_check_directory(sample_dir);
  }
  MPI_Barrier(MPI_COMM_WORLD);

  // 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(options_file, optfile);
  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) + dnest_thistask);
#else
  gsl_rng_set(dnest_gsl_r, 9999 + dnest_thistask);
  printf("# debugging, task %d dnest random seed %d\n", dnest_thistask, 9999 + dnest_thistask);
#endif

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

  if(param_range != NULL)
  {
    dnest_param_range = malloc(num_params*2*sizeof(double));
    memcpy(dnest_param_range, param_range, num_params*2*sizeof(double));
  }
  if(prior_type != NULL)
  {
    dnest_prior_type = malloc(num_params*sizeof(int));
    memcpy(dnest_prior_type, prior_type, num_params*sizeof(int));
  }
  if(prior_info != NULL)
  {
    dnest_prior_info = malloc(num_params*2*sizeof(double));
    memcpy(dnest_prior_info, prior_info, num_params*2*sizeof(double));
  }
  if(args != NULL)
  {
    dnest_args = args;
  }

  // read options
  if(dnest_thistask == dnest_root)
    options_load(optfile, opts);
  MPI_Bcast(&options, sizeof(DNestOptions), MPI_BYTE, dnest_root, MPI_COMM_WORLD);

  //dnest_post_temp = 1.0;
  compression = exp(1.0);
  regularisation = options.new_level_interval*sqrt(options.lam);
  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
  if(dnest_thistask == dnest_root)
    all_above = (LikelihoodType *)malloc(2*options.new_level_interval * sizeof(LikelihoodType));

  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_thistask == dnest_root)
    {
      levels_combine = (Level *)malloc(options.max_num_levels * sizeof(Level));
      copies_of_levels = (Level *)malloc(dnest_totaltask * options.max_num_levels * sizeof(Level));
    }

    if(dnest_flag_limits == 1)
    {
      limits = malloc(options.max_num_levels * particle_offset_double * 2 * sizeof(double));
      if(limits == NULL)
      {
        printf("Cannot allocate memory for limits.\n"
               "This usually happens when both the numbers of parameters and levels are extremely large.\n"
               "Please do not switch on '-l' option in argv passed to dnest.\n");
        exit(EXIT_FAILURE);
      }
      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;
        }
      }

      if(dnest_thistask == dnest_root)
      {
        copies_of_limits = malloc( dnest_totaltask * options.max_num_levels * particle_offset_double * 2 * sizeof(double));
        if(copies_of_limits == NULL)
        {
          printf("Cannot allocate memory for limits.\n"
                 "This usually happens when both the numbers of parameters and levels are extremely large.\n"
                 "Please do not switch on '-l' option in argv passed to dnest.\n");
          exit(EXIT_FAILURE);
        }
      }
    }
  }
  else
  {
    levels = (Level *)malloc(LEVEL_NUM_MAX * sizeof(Level));
    if(dnest_thistask == dnest_root)
    {
      levels_combine = (Level *)malloc(LEVEL_NUM_MAX * sizeof(Level));
      copies_of_levels = (Level *)malloc(dnest_totaltask * LEVEL_NUM_MAX * sizeof(Level));
    }

    if(dnest_flag_limits == 1)
    {
      limits = malloc(LEVEL_NUM_MAX * particle_offset_double * 2 * sizeof(double));
      if(limits == NULL)
      {
        printf("Cannot allocate memory for limits.\n"
               "This usually happens when both the numbers of parameters and levels are extremely large.\n"
               "Please do not switch on '-l' option in argv passed to dnest.\n");
        exit(EXIT_FAILURE);
      }
      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;
        }
      }

      if(dnest_thistask == dnest_root)
      {
        copies_of_limits = malloc(dnest_totaltask * LEVEL_NUM_MAX * particle_offset_double * 2 * sizeof(double));
        if(copies_of_limits == NULL)
        {
          printf("Cannot allocate memory for limits.\n"
                 "This usually happens when both the numbers of parameters and levels are extremely large.\n"
                 "Please do not switch on '-l' option in argv passed to dnest.\n");
          exit(EXIT_FAILURE);
        }
      }
    }
  }

  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;
  size_all_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++;

  if(dnest_thistask == dnest_root)
  {
    size_levels_combine = 0;
    levels_combine[size_levels_combine] = level_tmp;
    size_levels_combine++;
  }

  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);
    log_likelihoods[i].value = log_likelihoods_cal_initial(particles+i*particle_offset_size);
    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);

  if(dnest_thistask == dnest_root)
  {
    free(all_above);
    free(levels_combine);
    free(copies_of_levels);
    if(dnest_flag_limits == 1)
      free(copies_of_limits);
  }
  gsl_rng_free(dnest_gsl_r);

  free(dnest_perturb_accept);

  if(dnest_param_range != NULL)
  {
    free(dnest_param_range);
  }
  if(dnest_prior_type != NULL)
  {
    free(dnest_prior_type);
  }
  if(dnest_prior_info != NULL)
  {
    free(dnest_prior_info);
  }

  if(dnest_thistask == dnest_root)
  {
    printf("# Finalizing CDNest.\n");
    printf("#=======================================================\n");
  }
}

int dnest_search_pardict(DNestPARDICT *pardict, int num_pardict, char *tag)
{
  int i;
  for(i=0; i<num_pardict; i++)
  {
    if(strcmp(pardict[i].tag, tag) == 0)
    {
      return i;
    }
  }

  fprintf(stderr, "# cdnest no match of tag %s.\n", tag);
  return num_pardict+1;
}

void options_load(char *optfile, DNestOptions *opts)
{
  if(strlen(optfile) > 0)
  {
    DNestPARDICT *pardict;
    int num_pardict;
    pardict = malloc(10 * sizeof(DNestPARDICT));
    enum TYPE {INT, DOUBLE, STRING};

    FILE *fp;
    char str[BUF_MAX_LENGTH], buf1[BUF_MAX_LENGTH], buf2[BUF_MAX_LENGTH], buf3[BUF_MAX_LENGTH];

    int i, j, nt, idx;
    nt = 0;
    strcpy(pardict[nt].tag, "NumberParticles");
    pardict[nt].addr = &options.num_particles;
    pardict[nt].isset = 0;
    pardict[nt++].id = INT;

    strcpy(pardict[nt].tag, "NewLevelIntervalFactor");
    pardict[nt].addr = &options.new_level_interval_factor;
    pardict[nt].isset = 0;
    pardict[nt++].id = DOUBLE;

    strcpy(pardict[nt].tag, "SaveIntervalFactor");
    pardict[nt].addr = &options.save_interval_factor;
    pardict[nt].isset = 0;
    pardict[nt++].id = DOUBLE;

    strcpy(pardict[nt].tag, "ThreadStepsFactor");
    pardict[nt].addr = &options.thread_steps_factor;
    pardict[nt].isset = 0;
    pardict[nt++].id = DOUBLE;

    strcpy(pardict[nt].tag, "MaxNumberLevels");
    pardict[nt].addr = &options.max_num_levels;
    pardict[nt].isset = 0;
    pardict[nt++].id = INT;

    strcpy(pardict[nt].tag, "BacktrackingLength");
    pardict[nt].addr = &options.lam;
    pardict[nt].isset = 0;
    pardict[nt++].id = DOUBLE;

    strcpy(pardict[nt].tag, "StrengthEqualPush");
    pardict[nt].addr = &options.beta;
    pardict[nt].isset = 0;
    pardict[nt++].id = DOUBLE;

    strcpy(pardict[nt].tag, "MaxNumberSaves");
    pardict[nt].addr = &options.max_num_saves;
    pardict[nt].isset = 0;
    pardict[nt++].id = INT;

    strcpy(pardict[nt].tag, "PTol");
    pardict[nt].addr = &options.max_ptol;
    pardict[nt].isset = 0;
    pardict[nt++].id = DOUBLE;

    num_pardict = nt;

    /* default values */
    options.new_level_interval_factor = 2;
    options.save_interval_factor = options.new_level_interval_factor;
    options.thread_steps_factor = 10;
    options.num_particles = 1;
    options.max_num_levels = 0;
    options.lam = 10.0;
    options.beta = 100.0;
    options.max_ptol = 0.1;
    options.max_num_saves = 10000;

    fp = fopen(options_file, "r");
    if(fp == NULL)
    {
      fprintf(stderr, "# ERROR: Cannot open options file %s.\n", options_file);
      exit(0);
    }

    while(!feof(fp))
    {
      sprintf(str,"empty");
      fgets(str, 200, fp);
      if(sscanf(str, "%s%s%s", buf1, buf2, buf3)<2)
        continue;
      if(buf1[0]=='%' || buf1[0] == '#')
        continue;
      for(i=0, j=-1; i<nt; i++)
        if(strcmp(buf1, pardict[i].tag) == 0 && pardict[i].isset == 0)
        {
          j = i;
          pardict[i].isset = 1;
          //printf("%s %s\n", buf1, buf2);
          break;
        }
      if(j >=0)
      {
        switch(pardict[j].id)
        {
          case DOUBLE:
            *((double *) pardict[j].addr) = atof(buf2);
            break;
          case STRING:
            strcpy(pardict[j].addr, buf2);
            break;
          case INT:
            *((unsigned int *)pardict[j].addr) = (unsigned int) atof(buf2);
            break;
        }
      }
      else
      {
        fprintf(stderr, "# Error in file %s: Tag '%s' is not allowed or multiple defined.\n",
                      options_file, buf1);
        exit(0);
      }
    }
    fclose(fp);

    /* check options */
    idx = dnest_search_pardict(pardict, num_pardict, "SaveIntervalFactor");
    if(pardict[idx].isset == 0)  /* if not set */
    {
      options.save_interval_factor = options.new_level_interval_factor;
    }

    free(pardict);
  }
  else
  {
    options.new_level_interval_factor = opts->new_level_interval_factor;
    options.save_interval_factor = opts->save_interval_factor;
    options.thread_steps_factor = opts->thread_steps_factor;
    options.num_particles = opts->num_particles;
    options.max_num_levels = opts->max_num_levels;
    options.lam = opts->lam;
    options.beta = opts->beta;
    options.max_ptol = opts->max_ptol;
    options.max_num_saves = opts->max_num_saves;
  }

  options.thread_steps = dnest_num_params * options.thread_steps_factor * options.num_particles;
  options.new_level_interval =  dnest_totaltask * options.thread_steps * options.new_level_interval_factor;
  options.save_interval =  dnest_totaltask * options.thread_steps * options.save_interval_factor;

  //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 < dnest_totaltask * 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);
  }

  char fname[STR_MAX_LENGTH];
  strcpy(fname, dnest_sample_dir);
  strcat(fname, "/DNEST_OPTIONS");
  FILE *fp = fopen(fname, "w");
  if(fp == NULL)
  {
    fprintf(stderr, "# ERROR: Cannot write file %s.\n", fname);
    exit(0);
  }
  fprintf(fp, "NumberParticles          %d  # Number of particles\n", options.num_particles);
  fprintf(fp, "NewLevelIntervalFactor   %.2f  # New level interval factor\n", options.new_level_interval_factor);
  fprintf(fp, "SaveIntervalFactor       %.2f  # Save interval factor\n", options.save_interval_factor);
  fprintf(fp, "ThreadStepsFactor        %.2f  # ThreadSteps factor\n", options.thread_steps_factor);
  fprintf(fp, "MaxNumberLevels          %d  # Maximum number of levels\n", options.max_num_levels);
  fprintf(fp, "BacktrackingLength       %.1f  # Backtracking scale length\n", options.lam);
  fprintf(fp, "StrengthEqualPush        %.1f  # Strength of effect to force histogram to equal push\n", options.beta);
  fprintf(fp, "MaxNumberSaves           %d    # Maximum number of saves\n", options.max_num_saves);
  fprintf(fp, "PTol                     %.1e  # Likelihood tolerance in loge\n", options.max_ptol);
  fprintf(fp, "ThreadSteps              %d   #\n", options.thread_steps);
  fprintf(fp, "SaveInterval             %d   #\n", options.save_interval);
  fprintf(fp, "NewLevelInterval         %d  #\n", options.new_level_interval);
  fprintf(fp, "SampleFile               %s  #\n", options.sample_file);
  fprintf(fp, "SampleInfoFile           %s  #\n", options.sample_info_file);
  fprintf(fp, "SamplerStateFile         %s  #\n", options.sampler_state_file);
  fprintf(fp, "PosteriorSampleFile      %s  #\n", options.posterior_sample_file);
  fprintf(fp, "PosteriorSampleInfoFile  %s  #\n", options.posterior_sample_info_file);
  fprintf(fp, "LevelsFile               %s  #\n", options.levels_file);
  if(dnest_flag_limits == 1)
    fprintf(fp, "LimitsFile               %s  #\n", options.limits_file);
  fprintf(fp, "NumberCores              %d   # Number of cores\n", dnest_totaltask);
  fprintf(fp, "NumberParameters         %d   # Number of parameters\n", dnest_num_params);
  fclose(fp);
}

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);
  }

}

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

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

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

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

inline 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);
}

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

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

inline 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;
}

inline int dnest_get_size_levels()
{
  return size_levels;
}

inline int dnest_get_which_level_update()
{
  return dnest_which_level_update;
}

inline int dnest_get_which_particle_update()
{
  return dnest_which_particle_update;
}

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

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

inline void dnest_get_posterior_sample_file(char *fname)
{
  strcpy(fname, options.posterior_sample_file);
  return;
}

inline void dnest_get_limit(int ilevel, int jparam, double *limit1, double *limit2)
{
  *limit1 = limits[ilevel*dnest_num_params*2 + jparam * 2 + 0];
  *limit2 = limits[ilevel*dnest_num_params*2 + jparam * 2 + 1];
  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 at task %d.\
      \nSet to the default function in dnest.\n", dnest_thistask);
    fptrset->from_prior = dnest_from_prior;
  }

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

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

  if(fptrset->log_likelihoods_cal == NULL)
  {
    printf("\"log_likelihoods_cal\" function is not defined at task %d.\n", dnest_thistask);
    exit(0);
  }

  if(fptrset->log_likelihoods_cal_initial == NULL)
  {
    printf("\"log_likelihoods_cal_initial\" function is not defined at task %d. \
      \nSet to the same as \"log_likelihoods_cal\" function.\n", dnest_thistask);
    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 at task %d. \
      \nSet to the same as \"log_likelihoods_cal\" function.\n", dnest_thistask);
    fptrset->log_likelihoods_cal_restart = fptrset->log_likelihoods_cal;
  }

  if(fptrset->perturb == NULL)
  {
    printf("\"perturb\" function is not defined at task %d.\
      \nSet to the default function in dnest.\n", dnest_thistask);
    fptrset->perturb = dnest_perturb;
  }

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

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

  if(fptrset->kill_action == NULL)
  {
    printf("\"kill_action\" function is not defined at task %d.\
      \nSet to the default function in dnest.\n", dnest_thistask);
    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;
}

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

void dnest_check_directory(char *sample_dir)
{
  /* check if sample_dir exists
   * if not, create it;
   * if exists, check if it is a directory;
   * if not, throw an error.*/
  struct stat st;
  int status;
  status = stat(sample_dir, &st);
  if(status != 0)
  {
    printf("================================\n"
           "Directory %s not exist! create it.\n", sample_dir);
    status = mkdir(sample_dir, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
    if(status!=0)
    {
      printf("Cannot create %s\n"
             "================================\n", sample_dir);
    }
  }
  else
  {
    if(!S_ISDIR(st.st_mode))
    {
      printf("================================\n"
             "%s is not a direcotry!\n"
             "================================\n", sample_dir);
      exit(-1);
    }
  }
  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];

  if(dnest_thistask == dnest_root)
  {
    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);
    }

    particles_all = (void *)malloc( options.num_particles *  dnest_totaltask * dnest_size_of_modeltype );

    log_likelihoods_all = (LikelihoodType *)malloc(dnest_totaltask * options.num_particles * sizeof(LikelihoodType));
    level_assignments_all = (unsigned int*)malloc(dnest_totaltask * options.num_particles * sizeof(unsigned int));
  }

  MPI_Gather(particles, options.num_particles * dnest_size_of_modeltype, MPI_BYTE,
    particles_all, options.num_particles * dnest_size_of_modeltype, MPI_BYTE, dnest_root, MPI_COMM_WORLD);

  MPI_Gather(level_assignments, options.num_particles * sizeof(unsigned int), MPI_BYTE,
    level_assignments_all, options.num_particles * sizeof(unsigned int), MPI_BYTE, dnest_root, MPI_COMM_WORLD);

  MPI_Gather(log_likelihoods, options.num_particles * sizeof(LikelihoodType), MPI_BYTE,
    log_likelihoods_all, options.num_particles * sizeof(LikelihoodType), MPI_BYTE, dnest_root, MPI_COMM_WORLD);

  if(dnest_thistask == dnest_root )
  {
    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_combine, sizeof(int), 1, fp);

    for(i=0; i<size_levels_combine; 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_combine[i], sizeof(Level), 1, fp);
    }

    for(j=0; j<dnest_totaltask; j++)
    {
      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_all[j*options.num_particles + i], sizeof(int), 1, fp);
        fwrite(&log_likelihoods_all[j*options.num_particles + i], sizeof(LikelihoodType), 1, fp);
      }
    }

    if(dnest_flag_limits == 1)
    {
      for(i=0; i<size_levels_combine; 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(j=0; j<dnest_totaltask; j++)
    {
      for(i=0; i<options.num_particles; i++)
      {
        //print_particle(fp, particles_all + (j * options.num_particles + i) * particle_offset_size);
        fwrite(particles_all + (j * options.num_particles + i) * particle_offset_size, dnest_size_of_modeltype, 1, fp);
      }
    }

    fclose(fp);
    free(particles_all);
    free(log_likelihoods_all);
    free(level_assignments_all);
  }

  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;

  if(dnest_thistask == dnest_root)
  {
    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);

    particles_all = (void *)malloc( options.num_particles *  dnest_totaltask * dnest_size_of_modeltype );
    log_likelihoods_all = (LikelihoodType *)malloc(dnest_totaltask * options.num_particles * sizeof(LikelihoodType));
    level_assignments_all = (unsigned int*)malloc(dnest_totaltask * options.num_particles * sizeof(unsigned int));

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

    /* consider that the newly input max_num_levels may be different from the saved size of levels */
    if(options.max_num_levels != 0)
    {
      if(size_levels_combine > 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_combine);
        size_levels = options.max_num_levels;
      }
      else
      {
        size_levels = size_levels_combine;
      }

    }
    else  /* not input max_num_levels, directly use the saved size of levels */
    {
      if(size_levels_combine > LEVEL_NUM_MAX)
      {
        printf("# the saved size of levels %d exceeds LEVEL_NUM_MAX %d. \n", size_levels_combine, LEVEL_NUM_MAX);
        exit(EXIT_FAILURE);
      }
      else
      {
        size_levels = size_levels_combine;
      }
    }
    // read levels
    for(i=0; i<size_levels_combine; i++)
    {
      if(i<size_levels) // not read all the levels
        fread(&levels_combine[i], sizeof(Level), 1, fp);
      else
        fseek(fp, sizeof(Level), SEEK_CUR); /* offset the file point */
    }
    memcpy(levels, levels_combine, size_levels * sizeof(Level));

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

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

    // read limits
    if(dnest_flag_limits == 1)
    {
      for(i=0; i<size_levels_combine; 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(j=0; j<dnest_totaltask; j++)
    {
      for(i=0; i<options.num_particles; i++)
      {
        particle = (particles_all + (j * options.num_particles + i) * dnest_size_of_modeltype);
        fread(particle, dnest_size_of_modeltype, 1, fp);
      }
    }

    fclose(fp);
  }

  MPI_Bcast(&count_saves, 1, MPI_INT, dnest_root, MPI_COMM_WORLD);
  MPI_Bcast(&count_mcmc_steps, 1, MPI_INT, dnest_root, MPI_COMM_WORLD);
  MPI_Bcast(&size_levels, 1, MPI_INT, dnest_root, MPI_COMM_WORLD);
  MPI_Bcast(levels, size_levels * sizeof(Level), MPI_BYTE, dnest_root,  MPI_COMM_WORLD);

  if(count_saves > options.max_num_saves)
  {
    if(dnest_thistask == dnest_root)
    {
      printf("# Number of samples already larger than the input number, exit!\n");
    }
    MPI_Barrier(MPI_COMM_WORLD);
    exit(0);
  }
  size_levels_combine = size_levels; /* reset szie_levels_combine */

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

  if(dnest_flag_limits == 1)
    MPI_Bcast(limits, size_levels * particle_offset_double * 2, MPI_DOUBLE, dnest_root, MPI_COMM_WORLD);

  MPI_Scatter(level_assignments_all, options.num_particles * sizeof(unsigned int), MPI_BYTE,
      level_assignments, options.num_particles * sizeof(unsigned int), MPI_BYTE, dnest_root, MPI_COMM_WORLD);

  MPI_Scatter(log_likelihoods_all, options.num_particles * sizeof(LikelihoodType), MPI_BYTE,
      log_likelihoods, options.num_particles * sizeof(LikelihoodType), MPI_BYTE, dnest_root, MPI_COMM_WORLD);

  MPI_Scatter(particles_all, options.num_particles * dnest_size_of_modeltype, MPI_BYTE,
    particles, options.num_particles * dnest_size_of_modeltype, MPI_BYTE, dnest_root, MPI_COMM_WORLD);

  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);
    //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]--;
    }

  }

  if(dnest_thistask == dnest_root)
  {
    free(particles_all);
    free(log_likelihoods_all);
    free(level_assignments_all);
  }
  return;
}

void dnest_from_prior(void *model)
{
  int i;
  double *pm = (double *)model;

  for(i=0; i<dnest_num_params; i++)
  {
    if(dnest_prior_type[i] == GAUSSIAN )
    {
      pm[i] = dnest_randn() * dnest_prior_info[i*2+1] + dnest_prior_info[i*2+0];
      dnest_wrap(&pm[i], dnest_param_range[i*2+0], dnest_param_range[i*2+1]);
    }
    else if(dnest_prior_type[i] == LOG)
    {
      pm[i] = log(dnest_param_range[i*2+0]) + dnest_rand()*(log(dnest_param_range[i*2+1]) - log(dnest_param_range[i*2+0]));
      pm[i] = exp(pm[i]);
    }
    else
    {
      pm[i] = dnest_param_range[i*2+0] + dnest_rand()*(dnest_param_range[i*2+1] - dnest_param_range[i*2+0]);
    }
  }
}

double dnest_perturb(void *model)
{
  double *pm = (double *)model;
  double logH = 0.0, width;
  int which;

  which = dnest_rand_int(dnest_num_params);

  width = ( dnest_param_range[which*2+1] - dnest_param_range[which*2+0] );

  if(dnest_prior_type[which] == UNIFORM)
  {
    pm[which] += dnest_randh() * width;
    dnest_wrap(&(pm[which]), dnest_param_range[which*2+0], dnest_param_range[which*2+1]);
  }
  else if(dnest_prior_type[which] == LOG)
  {
    logH -= (-log(pm[which]));
    pm[which] += dnest_randh() * width;
    dnest_wrap(&(pm[which]), dnest_param_range[which*2+0], dnest_param_range[which*2+1]);
    logH += (-log(pm[which]));
  }
  else
  {
    logH -= (-0.5*pow((pm[which] - dnest_prior_info[which*2+0])/dnest_prior_info[which*2+1], 2.0) );
    pm[which] += dnest_randh() * width;
    dnest_wrap(&pm[which], dnest_param_range[which*2+0], dnest_param_range[which*2+1]);
    logH += (-0.5*pow((pm[which] - dnest_prior_info[which*2+0])/dnest_prior_info[which*2+1], 2.0) );
  }

  return logH;
}

inline void dnest_print_particle(FILE *fp, const void *model)
{
  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;
}

inline void dnest_restart_action(int iflag)
{
  return;
}

inline void dnest_accept_action()
{
  return;
}

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

max_stars_repo_name: LiyrAstroph/CDNest


text:
#ifndef BIO_BAYESIAN_HIERARCHICAL_CLUSTERING_H_
#define BIO_BAYESIAN_HIERARCHICAL_CLUSTERING_H_

#include "bio/defs.h"

#include <boost/graph/adjacency_list.hpp>
#include <boost/io/ios_state.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <boost/numeric/ublas/vector_proxy.hpp>
#include <boost/numeric/ublas/operation.hpp>
#include <boost/numeric/ublas/triangular.hpp>
#include <boost/numeric/ublas/storage.hpp>
#include <boost/numeric/ublas/lu.hpp>
#include <boost/numeric/ublas/io.hpp>

#include <set>
#include <list>
#include <limits>
#include <iterator>
#include <algorithm>
#include <iostream>

#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_sf_exp.h>
#include <gsl/gsl_sf_log.h>
#include <gsl/gsl_math.h>

BIO_NS_START

/** @file
See http://www.gatsby.ucl.ac.uk/~heller/bhc.pdf for the motivation for these algorithms.
*/

/** Aggregates a data set with its associated score. */
template <class DataIt, class Score>
struct Cluster
{
    typedef std::set<DataIt> data_set_t;

    data_set_t data_set;
    Score score;
};

/** A Normal-Wishart prior.

Each datum should be an iterator that points to a sequence of reals.

I am not convinced this is implemented correctly...
*/
struct NormalWishartPrior
{
    typedef double real_t;
    typedef boost::numeric::ublas::vector<real_t> vector_t;
    typedef boost::numeric::ublas::matrix<real_t> matrix_t;
    matrix_t s;
    vector_t m;
    real_t r;
    real_t v;
    unsigned get_num_dims() const { return s.size1(); }

    NormalWishartPrior(unsigned num_dimensions = 0)
        : s(num_dimensions,  num_dimensions)
        , m(num_dimensions)
    {
    }

    template <class DataIt>
    real_t
    operator()(
        DataIt data_begin,
        DataIt data_end) const
    {
        using namespace boost::numeric::ublas;

        //if we have no data
        if (data_end == data_begin)
        {
            return real_t(1);
        }

        //calculate the number of data points
        unsigned n = 0; //std::count_if(data_begin, data_end, true);
        for (DataIt di = data_begin;
            data_end != di;
            ++di)
        {
            ++n;
        }

        //put the data set in a matrix
        matrix_t x(get_num_dims(), n);
        unsigned i = 0;
        for (DataIt di = data_begin;
            data_end != di;
            ++di, ++i)
        {
            for (unsigned j = 0;
                get_num_dims() != j;
                ++j)
            {
                x(j,i) = (**di)[j];
            }
        }

        //calculate the sum of the data in each dimension
        vector_t x_sum(get_num_dims());
        for (unsigned j = 0;
            get_num_dims() != j;
            ++j)
        {
            x_sum(j) = 0;
            for (unsigned i = 0;
                n != i;
                ++i)
            {
                x_sum(j) += x(j, i);
            }
        }

        //calculate s_dash - the first term is s
        matrix_t s_dash(s);
        {
            //add X.XT to s_dash
            s_dash += prod(x, trans(x));

            //add m.mT term to s_dash
            s_dash += (r * n / (n + r)) * outer_prod(m, m);

            //add x_sum.x_sumT term to s_dash
            s_dash += (1.0 / (n + r)) * outer_prod(x_sum, x_sum);

            //subtract m.x_sum term from s_dash
            s_dash -= (r / (n + r)) * (outer_prod(m, x_sum) + outer_prod(x_sum, m));
        }

        //std::cout << s << std::endl;
        //std::cout << s_dash << std::endl;
        //std::cout << prod(x, trans(x)) << std::endl;

        //calculate the determinants
        const real_t s_det = get_determinant(s);
        const real_t s_dash_det = get_determinant(s_dash);

        //v_dash
        const real_t v_dash = v + n;

        //the number of dimensions
        const real_t k = get_num_dims();

        //the number of data points
        const real_t N = n;

        //the natural logarithm of the result
        real_t ln_p_D_given_H1 = 0.0;

        ln_p_D_given_H1 -= N * k / 2 * gsl_sf_log(2 * M_PI);
        ln_p_D_given_H1 += k / 2 * gsl_sf_log(r / (N + r));
        ln_p_D_given_H1 += v / 2 * gsl_sf_log(fabs(s_det));
        ln_p_D_given_H1 -= v_dash / 2 * gsl_sf_log(fabs(s_dash_det));
        ln_p_D_given_H1 += v_dash * k / 2 * gsl_sf_log(2);
        ln_p_D_given_H1 -= v * k / 2 * gsl_sf_log(2);

        //add/subtract the pi terms
        for (unsigned d = 0; k != d; ++d)
        {
            ln_p_D_given_H1 +=
                gsl_sf_lngamma((v_dash + 1 - d) / 2)
                - gsl_sf_lngamma((v + 1 - d) / 2);

        }

        return gsl_sf_exp(ln_p_D_given_H1);
    }

    template <class M>
    static real_t
    get_determinant(const M & m)
    {
        BOOST_ASSERT(m.size1() == m.size2());

        matrix_t lu(m);
        //permutation_matrix<unsigned> p(m.size());

        lu_factorize(lu);

        real_t result = 0;
        for (unsigned j = 0;
            lu.size1() != j;
            ++j)
        {
            result += lu(j,j);
        }
        return result;
    }
};

/** A Bernoulli - Beta prior.

Each datum should be an iterator that points to a sequence of bools (or convertible types).
*/
struct BernoulliBetaPrior
{
    typedef double real_t;
    typedef std::vector<real_t> param_vec_t;
    param_vec_t alpha;
    param_vec_t beta;

    BernoulliBetaPrior(unsigned num_dimensions = 0)
        : alpha(num_dimensions)
        , beta(num_dimensions)
    {
    }

    template <class DataIt>
    real_t
    operator()(
        DataIt data_begin,
        DataIt data_end)
    {
        if (data_end == data_begin)
        {
            return real_t(1);
        }

        assert(alpha.size() == beta.size());
        assert((*data_begin)->size() == alpha.size());

        //the log of p(D|H1)
        double ln_p_D_H1 = 0.0;

        unsigned n = 0;
        for (DataIt i = data_begin; data_end != i; ++i)
        {
            ++n;
        }

        //for each dimension
        for (size_t d = 0; alpha.size() != d; ++d)
        {
            //work out m_d
            unsigned m_d = 0;
            for (DataIt i = data_begin; data_end != i; ++i)
            {
                assert((*i)->size() == alpha.size());

                if ((**i)[d])
                {
                    ++m_d;
                }
            }

            ln_p_D_H1 +=
                gsl_sf_lngamma(alpha[d] + beta[d])
                + gsl_sf_lngamma(alpha[d] + m_d)
                + gsl_sf_lngamma(beta[d] + n - m_d)
                - gsl_sf_lngamma(alpha[d])
                - gsl_sf_lngamma(beta[d])
                - gsl_sf_lngamma(alpha[d] + beta[d] + n);
        }

        return gsl_sf_exp(ln_p_D_H1);
    }
};

/** A Multinomial Dirichlet prior.

Each datum should be an iterator that points to a sequence of bools (or convertible types).
*/
struct MultinomialDirichletPrior
{
    typedef double real_t;
    typedef std::vector<real_t> param_vec_t;
    param_vec_t alpha;

    MultinomialDirichletPrior(unsigned num_dimensions = 0)
        : alpha(num_dimensions)
    {
    }

    template <class DataIt>
    real_t
    operator()(
        DataIt data_begin,
        DataIt data_end)
    {
        if (data_end == data_begin)
        {
            return real_t(1);
        }

        assert((*data_begin)->size() == alpha.size());

        //the log of p(D|H1) ie log of the result
        double ln_p_D_H1 = 0.0;

        //TODO - not implemented

        return gsl_sf_exp(ln_p_D_H1);
    }
};

/** The score of a cluster based on a prior model and the score of the two clusters
it is composed of. */
struct BayesianClusterScore
{
    typedef double real_t;

    real_t p_Dk_given_H1k;
    real_t p_Dk_given_Tk;

    BayesianClusterScore(real_t p_Dk_given_H1k = 1, real_t p_Dk_given_Tk = 1)
        : p_Dk_given_H1k(p_Dk_given_H1k)
        , p_Dk_given_Tk(p_Dk_given_Tk)
    {
    }

    bool operator<(const BayesianClusterScore & rhs) const
    {
        if (real_t(0) == p_Dk_given_Tk)
        {
            return false;
        }
        if (real_t(0) == rhs.p_Dk_given_Tk)
        {
            return false;
        }
        return operator()() < rhs();
    }

    real_t operator()() const
    {
        return p_Dk_given_H1k / p_Dk_given_Tk;
    }
};

BIO_NS_END

std::ostream &
operator<<(std::ostream & os, const BIO_NS::BayesianClusterScore & score)
{
    boost::io::ios_all_saver ias(os);
    os.precision(3);

    os << score.p_Dk_given_H1k << "/" << score.p_Dk_given_Tk << "=" << score();
    return os;
}

BIO_NS_START

/** Scores clusters according to http://www.gatsby.ucl.ac.uk/~heller/bhcnew.pdf */
template <class Model>
struct BayesianClusterScorer
{
    typedef double real_t;
    typedef Model model_t;

    model_t & model;
    BayesianClusterScorer(model_t & model)
        : model(model)
    {
    }

    typedef BayesianClusterScore score_t;
    static score_t min_score() { return score_t(-1, 1); }

    template <class Cluster>
    score_t
    operator()(
        const Cluster & cluster) const
    {
        const real_t marginal_likelihood = model(cluster.data_set.begin(), cluster.data_set.end());
        return score_t(marginal_likelihood, marginal_likelihood);
    }

    /** Returns the score of the data set resulting from merging the two data sets. */
    template <class Cluster>
    score_t
    operator()(
        const Cluster & cluster_1,
        const Cluster & cluster_2) const
    {
        //merge the two data sets
        typename Cluster::data_set_t merged_data_set;
        set_union(
            cluster_1.data_set.begin(),
            cluster_1.data_set.end(),
            cluster_2.data_set.begin(),
            cluster_2.data_set.end(),
            std::inserter(merged_data_set, merged_data_set.begin()));

        score_t result;
        const real_t pi_k = 0.5; //?????
        result.p_Dk_given_H1k =
            pi_k * model(merged_data_set.begin(), merged_data_set.end());
        result.p_Dk_given_Tk =
            result.p_Dk_given_H1k
            + (1 - pi_k) * cluster_1.score.p_Dk_given_Tk * cluster_2.score.p_Dk_given_Tk;

        return result;
    }
};

/** Defines some types for given data types and scorers. */
template <class DataIt, class Scorer>
struct ClusterTraits
{
    typedef DataIt data_it;
    typedef Scorer scorer_t;

    typedef typename data_it::value_type data_t;
    typedef std::set<data_it> data_set_t;

    typedef Cluster<data_it, typename scorer_t::score_t> cluster_t;
    typedef typename cluster_t::data_set_t::iterator cluster_data_it;

    typedef std::list<cluster_t> cluster_list_t;
    typedef typename cluster_list_t::iterator cluster_it;
};

/** Clusters the data points according to pairs' scores. */
template <class DataIt, class Scorer, class Output>
void
cluster(
    DataIt data_begin,
    DataIt data_end,
    Scorer & scorer,
    Output & output)
{
    using namespace boost;
    using namespace std;

    typedef ClusterTraits<DataIt, Scorer> traits_t;

    //initialise the list of clusters with one entry for every data point
    typename traits_t::cluster_list_t cluster_list;
    for (DataIt d = data_begin;
        data_end != d;
        ++d)
    {
        typename traits_t::cluster_t cluster;
        cluster.data_set.insert(d);
        cluster.score = scorer(cluster);
        output(cluster_list.insert(cluster_list.begin(), cluster));
    }

    //while we have at least 2 clusters to merge
    while (1 < cluster_list.size())
    {
        //look for the best pair of clusters to merge in the hierarchy
        typename std::pair<typename traits_t::cluster_it, typename traits_t::cluster_it> best_clusters =
            make_pair(cluster_list.end(), cluster_list.end());
        typename Scorer::score_t best_score = Scorer::min_score();

        //for each pair of clusters, c1 and c2
        for (typename traits_t::cluster_it c1 = cluster_list.begin();
            cluster_list.end() != c1;
            ++c1)
        {
            typename traits_t::cluster_it c2 = c1;
            ++c2;
            for ( ;
                cluster_list.end() != c2;
                ++c2)
            {
                //score the pair of clusters
                const typename Scorer::score_t score = scorer(*c1, *c2);

                //is it the best score so far?
                if (best_score < score)
                {
                    //update best
                    best_score = score;
                    best_clusters.first = c1;
                    best_clusters.second = c2;
                }
            }
        }
        assert(best_clusters.first != best_clusters.second);
        assert(cluster_list.end() != best_clusters.first);
        assert(cluster_list.end() != best_clusters.second);

        //create and add the merged cluster
        typename traits_t::cluster_t merged_cluster;
        merged_cluster.score = best_score;
        set_union(
            best_clusters.first->data_set.begin(),
            best_clusters.first->data_set.end(),
            best_clusters.second->data_set.begin(),
            best_clusters.second->data_set.end(),
            inserter(merged_cluster.data_set, merged_cluster.data_set.begin()));

        //let the output know
        output(
            best_clusters.first,
            best_clusters.second,
            cluster_list.insert(cluster_list.begin(), merged_cluster));

        //remove the old data sets and add the merged one
        cluster_list.erase(best_clusters.first);
        cluster_list.erase(best_clusters.second);
    }
}

struct LoggingClusterer
{
    std::ostream & os;

    LoggingClusterer(std::ostream & os) : os(os)
    {
    }

    template <class ClusterIt>
    void
    operator()(
        const ClusterIt & cluster)
    {
        using namespace std;

        os << "Initial cluster: ";
        for (typename ClusterIt::value_type::data_set_t::const_iterator d = cluster->data_set.begin();
            cluster->data_set.end() != d;
            ++d)
        {
            os << **d << ",";
        }
        os << " score:" << cluster->score << endl;
    }

    template <class ClusterIt>
    void
    operator()(
        ClusterIt cluster_1,
        ClusterIt cluster_2,
        ClusterIt merged_cluster)
    {
        using namespace std;

        //print what we are doing
        cout << "Clustering: ";
        for (typename ClusterIt::value_type::data_set_t::const_iterator d = cluster_1->data_set.begin();
            cluster_1->data_set.end() != d;
            ++d)
        {
            os << **d << ",";
        }
        cout << " with ";
        for (typename ClusterIt::value_type::data_set_t::const_iterator d = cluster_2->data_set.begin();
            cluster_2->data_set.end() != d;
            ++d)
        {
            os << **d << ",";
        }
        os << " score:" << merged_cluster->score << endl;
    }
};

template <class DataIt, class Scorer>
struct TreeBuildingClusterer
{
    typedef DataIt data_it;
    typedef Scorer scorer_t;
    typedef ClusterTraits<data_it, scorer_t> traits_t;

    typedef std::pair<bool, typename traits_t::data_it> vertex_property;
    typedef boost::adjacency_list<
        boost::listS,
        boost::vecS,
        boost::directedS,
        vertex_property,
        typename Scorer::score_t> tree_t;
    typedef typename boost::graph_traits<tree_t>::vertex_iterator tree_vertex_it;
    typedef typename boost::graph_traits<tree_t>::vertex_descriptor tree_vertex;
    tree_t tree;

    /** A map from clusters to tree nodes. */
    typedef std::map<typename traits_t::cluster_t *, tree_vertex> cluster_vertex_map_t;
    cluster_vertex_map_t cluster_vertex_map;

    LoggingClusterer logger;
    bool logging_on;

    TreeBuildingClusterer(bool logging_on = false) : logger(std::cout), logging_on(logging_on)
    {
    }

    void
    reset()
    {
        cluster_vertex_map.clear();
        tree.clear();
    }

    void
    operator()(
        typename traits_t::cluster_it cluster)
    {
        using namespace std;
        using namespace boost;

        if (logging_on)
        {
            logger(cluster);
        }

        cluster_vertex_map[&*cluster] = add_vertex(make_pair(true, *cluster->data_set.begin()), tree);
    }

    void
    operator()(
        typename traits_t::cluster_it cluster_1,
        typename traits_t::cluster_it cluster_2,
        typename traits_t::cluster_it merged_cluster)
    {
        using namespace std;
        using namespace boost;

        if (logging_on)
        {
            logger(cluster_1, cluster_2, merged_cluster);
        }

        //find the vertex descriptors for cluster_1 and cluster_2
        assert(cluster_vertex_map.end() != cluster_vertex_map.find(&*cluster_1));
        assert(cluster_vertex_map.end() != cluster_vertex_map.find(&*cluster_2));

        cluster_vertex_map[&*merged_cluster] = add_vertex(make_pair(false, typename traits_t::data_it()), tree);

        add_edge(cluster_vertex_map[&*merged_cluster], cluster_vertex_map[&*cluster_1], tree);
        add_edge(cluster_vertex_map[&*merged_cluster], cluster_vertex_map[&*cluster_2], tree);
    }
};

BIO_NS_END

#endif //BIO_BAYESIAN_HIERARCHICAL_CLUSTERING_H_

ext: h

max_stars_repo_name: JohnReid/biopsy

max_issues_repo_name: JohnReid/biopsy

max_forks_repo_name: JohnReid/biopsy


text:
/* specfunc/bessel_Y0.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/* Author:  G. Jungman */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include "gsl_sf_trig.h"
#include "gsl_sf_bessel.h"

#include "error.h"

#include "bessel.h"
#include "bessel_amp_phase.h"
#include "cheb_eval.c"

/*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/

/* based on SLATEC besy0, 1980 version, w. fullerton */

/* chebyshev expansions

 series for by0        on the interval  0.	    to  1.60000d+01
					with weighted error   1.20e-17
					 log weighted error  16.92
			       significant figures required  16.15
				    decimal places required  17.48
*/

static double by0_data[13] = {
  -0.011277839392865573,
  -0.128345237560420350,
  -0.104378847997942490,
   0.023662749183969695,
  -0.002090391647700486,
   0.000103975453939057,
  -0.000003369747162423,
   0.000000077293842676,
  -0.000000001324976772,
   0.000000000017648232,
  -0.000000000000188105,
   0.000000000000001641,
  -0.000000000000000011
};
static cheb_series by0_cs = {
  by0_data,
  12,
  -1, 1,
  8
};

/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/

int gsl_sf_bessel_Y0_e(const double x, gsl_sf_result * result)
{
  const double two_over_pi = 2.0/M_PI;
  const double xmax        = 1.0/GSL_DBL_EPSILON;

  /* CHECK_POINTER(result) */

  if (x <= 0.0) {
    DOMAIN_ERROR(result);
  }
  else if(x < 4.0) {
    gsl_sf_result J0;
    gsl_sf_result c;
    int stat_J0 = gsl_sf_bessel_J0_e(x, &J0);
    cheb_eval_e(&by0_cs, 0.125*x*x-1.0, &c);
    result->val = two_over_pi*(-M_LN2 + log(x))*J0.val + 0.375 + c.val;
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val) + c.err;
    return stat_J0;
  }
  else if(x < xmax) {
    /* Leading behaviour of phase is x, which is exact,
     * so the error is bounded.
     */
    const double z  = 32.0/(x*x) - 1.0;
    gsl_sf_result c1;
    gsl_sf_result c2;
    gsl_sf_result sp;
    const int stat_c1 = cheb_eval_e(&_gsl_sf_bessel_amp_phase_bm0_cs,  z, &c1);
    const int stat_c2 = cheb_eval_e(&_gsl_sf_bessel_amp_phase_bth0_cs, z, &c2);
    const int stat_sp = gsl_sf_bessel_sin_pi4_e(x, c2.val/x, &sp);
    const double sqrtx = sqrt(x);
    const double ampl  = (0.75 + c1.val) / sqrtx;
    result->val  = ampl * sp.val;
    result->err  = fabs(sp.val) * c1.err/sqrtx + fabs(ampl) * sp.err;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_ERROR_SELECT_3(stat_sp, stat_c1, stat_c2);
  }
  else {
    UNDERFLOW_ERROR(result);
  }
}

/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/

#include "eval.h"

double gsl_sf_bessel_Y0(const double x)
{
  EVAL_RESULT(gsl_sf_bessel_Y0_e(x, &result));
}

ext: c


text:
/**
 * @file tde.h
 * @brief Time delay estimation
 * @author Kenichi Kumatani
 */

#ifndef TDE_H
#define TDE_H

#include <stdio.h>
#include <assert.h>
#include <float.h>

#include <gsl/gsl_block.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_fft_real.h>
#include <gsl/gsl_fft_complex.h>
#include <common/refcount.h>
#include "common/jexception.h"

#include "stream/stream.h"
#include "feature/feature.h"
#include "modulated/modulated.h"
#include "btk.h"

// ----- definition for class `CCTDE' -----
//
/**
   @class find the time difference which provides the maximum correlation between tow signals.
   @usage
   1. constrct sample feature objects :
       samp1 =
       samp2 =
   2. construct this object and feed the sample features into it:
*/
class CCTDE : public VectorFeatureStream {
public:
  /**
     @brief
     @param
     @param
     @param
     @param
   */
  CCTDE( SampleFeaturePtr& samp1, SampleFeaturePtr& samp2, int fftLen=512, unsigned nHeldMaxCC=1, int freqLowerLimit=-1, int freqUpperLimit=-1, const String& nm= "CCTDE" );
  ~CCTDE();

  virtual const gsl_vector* next(int frame_no = -5 );
  virtual const gsl_vector* nextX( unsigned chanX = 0, int frame_no = -5);
  virtual void  reset();
  void  allsamples( int fftLen = -1 );
  void set_target_frequency_range( int freqLowerLimit, int freqUpperLimit ){
    freq_lower_limit_ = freqLowerLimit;
    freq_upper_limit_ = freqUpperLimit;
  }
  const unsigned *sample_delays() const { return sample_delays_; }
  const double *cc_values() const {return cc_values_; }

#ifdef ENABLE_LEGACY_BTK_API
  void setTargetFrequencyRange( int freqLowerLimit, int freqUpperLimit ){ set_target_frequency_range(freqLowerLimit, freqUpperLimit); }
  const unsigned *getSampleDelays(){ return sample_delays(); }
  const double *getCCValues(){ return cc_values(); }
#endif

private:
  const gsl_vector* detect_cc_peaks_( double **samples, size_t stride );

  typedef list<SampleFeaturePtr>	ChannelList_;
  typedef ChannelList_::iterator	ChannelIterator_;
  ChannelList_	channelL_;   // must be 2.
  unsigned	*sample_delays_; // sample delays
  double	*cc_values_;     // cross-correlation (CC) values
  unsigned	nHeldMaxCC_;    // how many CC values are held
  unsigned	fftLen_;
  int		samplerate_;
  gsl_vector	*window_;
  int		freq_lower_limit_;
  int		freq_upper_limit_;
  vector<unsigned > _frameCounter;
};

typedef Inherit<CCTDE, VectorFeatureStreamPtr> CCTDEPtr;

#endif


text:
#ifndef LOGSPLINEPDF
#define LOGSPLINEPDF

#include <gsl/gsl_rng.h>

#include "splinetable.h"

#ifdef __cplusplus
extern "C" {
#endif

void logsplinepdf_n_sample(double *result, int results, int burnin,
    double *coords, int dim, struct splinetable *table, int derivatives,
    double (* proposal)(void*), double (* proposal_pdf)(double, double, void*),
    void *proposal_info, const gsl_rng *rng);

void splinepdf_n_sample(double *result, int results, int burnin,
    double *coords, int dim, struct splinetable *table, int derivatives,
    double (* proposal)(void*), double (* proposal_pdf)(double, double, void*),
    void *proposal_info, const gsl_rng *rng);

#ifdef __cplusplus
}
#endif

#endif


text:
/*
 * Adapted from bits and pieces in https://code.google.com/p/eq-polar-alignment/
 * Could not find copyright statement from there, but license seems to be GPLv3...
 *
 * Copyright 2014, Kalle Vahlman, zuh@iki.fi
 * Licensed under the GNU GPL v3
 *
 */

#include <stdio.h>
#include <time.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_multiroots.h>
#include <astrometry/anwcs.h>
#include <astrometry/sip.h>
#include <astrometry/sip_qfits.h>

struct coord
{
    double ra;
    double dec;
    double x;
    double y;
};

struct state
{
    anwcs_t* h_wcs;
    anwcs_t* v_wcs;
    double pixel_scale;
    struct coord ncp;
    struct coord polaris;
    struct coord lambda;
    struct coord axis;

};

double julian_date()
{
    time_t now = time(NULL);
    struct tm *t = localtime(&now);
    int day = t->tm_mday;
    int month = t->tm_mon+1;
    int year = 1900 + t->tm_year;

    return 367*year-7*(year+(month+9)/12)/4-3*((year+(month-9)/7)/100+1)/4+275*month/9+day+1721029;
}

int calculate_ncp (double julian, struct coord *c)
{
    double t, zeta, theta;
    if (c == NULL)
        return 0;

    t = (julian-2451545)/36525;
    zeta = (2306.2181*t+0.30188*t*t+0.017998*t*t*t)/3600;
    theta = (2004.3109*t-0.42665*t*t-0.041833*t*t*t)/3600;
    c->dec = 90 - theta;
    c->ra = -zeta;

    return 1;
}

int solve_field(struct state *ps, const char *dir, const char *in, const char *out)
{
    char run[512];
    // TODO: Run the solving with code rather than system()
#define SOLVE "solve-field -B none -M none -S none -R none -U none -N none " \
              "--no-plots --temp-dir %s --dir %s --out %s " \
              "--downsample 2 --ra %.4f --dec %.4f --radius 10 %s"
    sprintf(run, SOLVE, dir, dir, out, ps->ncp.ra, ps->ncp.dec, in);
    printf("%s\n", run);
    return system(run);
}

int fvec(const gsl_vector *x, void *params, gsl_vector *f)
{
    struct state *s = (struct state *)params;
    double ra, dec, xp, yp;
    double xi = gsl_vector_get(x, 0);
    double yi = gsl_vector_get(x, 1);

    anwcs_pixelxy2radec(s->h_wcs, xi, yi, &ra, &dec);
    anwcs_radec2pixelxy(s->v_wcs, ra, dec, &xp, &yp);
    xp = xp - xi;
    yp = yp - yi;
    gsl_vector_set(f, 0, xp);
    gsl_vector_set(f, 1, yp);

    return GSL_SUCCESS;
}

int solve_axis(struct state *ps)
{
    int s;
    size_t iter=0;
    const size_t n=2;
    gsl_multiroot_fsolver *solver;
    gsl_multiroot_function f = { &fvec, n, ps};
    gsl_vector *x = gsl_vector_alloc(n);

    gsl_vector_set(x, 0, ps->polaris.x);
    gsl_vector_set(x, 1, ps->polaris.y);
    solver = gsl_multiroot_fsolver_alloc(gsl_multiroot_fsolver_hybrids, 2);
    gsl_multiroot_fsolver_set(solver, &f, x);
    do {
        iter++;
        s = gsl_multiroot_fsolver_iterate(solver);
        if (s) break;
        s = gsl_multiroot_test_residual(solver->f,1e-7);
    } while (s == GSL_CONTINUE && iter < 1000);
    ps->axis.x = gsl_vector_get(solver->x, 0);
    ps->axis.y = gsl_vector_get(solver->x, 1);

    gsl_multiroot_fsolver_free(solver);
    gsl_vector_free(x);
    return s;
}

#define MAX(a, b) a < b ? b : a

int plot(struct state *ps, char *in, char *out)
{
    char convert[2048];
    double arcmin = 60/ps->pixel_scale;
    double cropwidth = MAX(400, 120 * arcmin);
    double cropheight = MAX(400, 120 * arcmin);
    double cropx = ps->ncp.x - cropwidth/2;
    double cropy = ps->ncp.y - cropheight/2;
    int fontsize = ps->pixel_scale < 24 ? 16 : 8;

#define CONVERT "convert -fill none " \
    "-pointsize %d " \
    "-font Courier -stroke cyan " \
    "-draw 'point %.0f,%.0f' " \
    "-draw 'circle %.0f,%.0f %.0f,%.0f' " \
    "-draw 'circle %.0f,%.0f %.0f,%.0f' " \
    "-draw 'circle %.0f,%.0f %.0f,%.0f' " \
    "-draw 'circle %.0f,%.0f %.0f,%.0f' " \
    "-draw 'circle %.0f,%.0f %.0f,%.0f' " \
    "-draw \"text %.0f,%.0f '%s'\" " \
    "-draw \"text %.0f,%.0f '%s'\" " \
    "-draw \"text %.0f,%.0f '%s'\" " \
    "-draw \"text %.0f,%.0f '%s'\" " \
    "-draw \"text %.0f,%.0f '%s'\" " \
    "-stroke red " \
    "-draw 'line %.0f,%.0f %.0f,%.0f' " \
    "-draw 'line %.0f,%.0f %.0f,%.0f' " \
    "-draw 'circle %.0f,%.0f %.0f,%.0f' " \
    "-draw 'circle %.0f,%.0f %.0f,%.0f' " \
    "-stroke white " \
    "-draw 'line %.0f,%.0f %.0f,%.0f' " \
    "-draw 'line %.0f,%.0f %.0f,%.0f' " \
    "-font Symbol -fill white " \
    "-draw \"text %.0f,%.0f '%s'\" " \
    "-fill none " \
    "-stroke orange " \
    "-draw 'line %.0f,%.0f %.0f,%.0f' " \
    "-draw 'line %.0f,%.0f %.0f,%.0f' " \
    "-font Symbol -fill orange " \
    "-draw \"text %.0f,%.0f '%s'\" " \
    "-fill none " \
    "-font Courier -stroke white -fill white " \
    "-draw \"text %.0f,%.0f 'Current pixel offset from NCP: %.0f, %.0f'\" " \
    "-draw \"text %.0f,%.0f '%s %s %s %s'\" " \
    "-crop '%.0fx%.0f+%.0f+%.0f' " \
    "%s %s"

    sprintf(convert, CONVERT, fontsize,
        // Target circles
        ps->ncp.x, ps->ncp.y,
        ps->ncp.x, ps->ncp.y, ps->ncp.x+2*arcmin, ps->ncp.y,
        ps->ncp.x, ps->ncp.y, ps->ncp.x+5*arcmin, ps->ncp.y,
        ps->ncp.x, ps->ncp.y, ps->ncp.x+10*arcmin, ps->ncp.y,
        ps->ncp.x, ps->ncp.y, ps->ncp.x+20*arcmin, ps->ncp.y,
        ps->ncp.x, ps->ncp.y, ps->ncp.x+40*arcmin, ps->ncp.y,
        ps->ncp.x-2*arcmin+2, ps->ncp.y, ps->pixel_scale < 24 ? "2\\'" : "",
        ps->ncp.x-5*arcmin+2, ps->ncp.y, ps->pixel_scale < 24 ? "5\\'" : "",
        ps->ncp.x-10*arcmin+2, ps->ncp.y, ps->pixel_scale < 24 ? "10\\'" : "",
        ps->ncp.x-20*arcmin+2, ps->ncp.y, "20\\'",
        ps->ncp.x-40*arcmin+2, ps->ncp.y, "40\\'",
        // Current axis
        ps->axis.x-5*arcmin, ps->axis.y-5*arcmin, ps->axis.x+5*arcmin, ps->axis.y+5*arcmin,
        ps->axis.x-5*arcmin, ps->axis.y+5*arcmin, ps->axis.x+5*arcmin, ps->axis.y-5*arcmin,
        ps->axis.x, ps->axis.y, ps->axis.x+2*arcmin, ps->axis.y,
        ps->axis.x, ps->axis.y, ps->axis.x+5*arcmin, ps->axis.y,
        // Polaris
        ps->polaris.x, ps->polaris.y-1*arcmin,
            ps->polaris.x, ps->polaris.y+1*arcmin,
        ps->polaris.x-1*arcmin, ps->polaris.y,
            ps->polaris.x+1*arcmin, ps->polaris.y,
        ps->polaris.x+1*arcmin, ps->polaris.y-1*arcmin, "a",
        // λ UMi
        ps->lambda.x, ps->lambda.y-1*arcmin,
            ps->lambda.x, ps->lambda.y+1*arcmin,
        ps->lambda.x-1*arcmin, ps->lambda.y,
            ps->lambda.x+1*arcmin, ps->lambda.y,
        ps->lambda.x+1*arcmin, ps->lambda.y-1*arcmin, "l",
        // Offset report
        cropx + 16, cropy + 32, ps->ncp.x-ps->axis.x, ps->ncp.y-ps->axis.y,
        cropx + 16, cropy + 64,
            (ps->ncp.x == ps->axis.x && ps->ncp.y == ps->axis.y)
                ? "Well done, perfect alignment!"
                : "Adjust mount alignment to",
            (ps->ncp.x == ps->axis.x && ps->ncp.y == ps->axis.y)
                ? ""
                : ps->ncp.x < ps->axis.x ? "left" : "right",
            (ps->ncp.x == ps->axis.x && ps->ncp.y == ps->axis.y) ? "" : "and",
            (ps->ncp.x == ps->axis.x && ps->ncp.y == ps->axis.y)
                ? ""
                : ps->ncp.y < ps->axis.y ? "up" : "down",
        // Crop
        cropwidth, cropheight, cropx, cropy,
        in, out);

    printf("%s\n", convert);
    return system(convert);
}

int main (int argc, char **argv)
{
    char dir[20] = "/tmp/polar-XXXXXX";
    char h_dst[30];
    char v_dst[30];
    sip_t *sip;
    struct state ps = { NULL, NULL, 0,
                      { 0, 0, 0, 0 }, // NCP
                      { 37.9529, 89.2642, 0, 0 }, // Polaris
                      { 259.2367, 89.0378, 0, 0 }, // λ UMi
                      { 0, 0, 0, 0 } // axis
                    };

    if (argc < 3) {
        printf("E: Not enough parameters!\n");
        printf("Usage: %s <horizontal image> <vertical image> [output]\n", argv[0]);
        return -1;
    }

    if (mkdtemp(dir) == NULL) {
        printf("E: Could not create working directory!\n");
        return -1;
    }

    // Calculate current ra,dec position of North Celestial Pole
    calculate_ncp(julian_date(), &ps.ncp);

    // Solve horizontal and vertical images of Polaris and Lambda UMi
    sprintf(h_dst, "%s/h.wcs", dir);
    sprintf(v_dst, "%s/v.wcs", dir);
    solve_field(&ps, dir, argv[1], "h.wcs");
    solve_field(&ps, dir, argv[2], "v.wcs");

    // Pick pixel coordinates from horizontal image
    ps.h_wcs = anwcs_open(h_dst, 0);
    anwcs_radec2pixelxy(ps.h_wcs, ps.polaris.ra, ps.polaris.dec, &ps.polaris.x, &ps.polaris.y);
    anwcs_radec2pixelxy(ps.h_wcs, ps.lambda.ra, ps.lambda.dec, &ps.lambda.x, &ps.lambda.y);
    anwcs_radec2pixelxy(ps.h_wcs, ps.ncp.ra, ps.ncp.dec, &ps.ncp.x, &ps.ncp.y);

    // Solve rotational axis in pixel coordinates of the horizontal image
    ps.v_wcs = anwcs_open(v_dst, 0);
    solve_axis(&ps);
    anwcs_pixelxy2radec(ps.h_wcs, ps.axis.x, ps.axis.y, &ps.axis.ra, &ps.axis.dec);
    anwcs_free(ps.v_wcs);

    // Grab pixel scale from horizontal image
    sip = anwcs_get_sip(ps.h_wcs);
    ps.pixel_scale = sip_pixel_scale(sip);
    anwcs_free(ps.h_wcs);

    // Print what we found
    printf("Pixel scale: %.4f\n", ps.pixel_scale);
    printf("Polaris: %+9.4f ra, %+9.4f dec (%4.0f, %4.0f)\n",
        ps.lambda.ra, ps.lambda.dec, ps.lambda.x, ps.lambda.y);
    printf("λ UMi  : %+9.4f ra, %+9.4f dec (%4.0f, %4.0f)\n",
        ps.polaris.ra, ps.polaris.dec, ps.polaris.x, ps.polaris.y);
    printf("NCP    : %+9.4f ra, %+9.4f dec (%4.0f, %4.0f)\n",
        ps.ncp.ra, ps.ncp.dec, ps.ncp.x, ps.ncp.y);
    printf("Axis   : %+9.4f ra, %+9.4f dec (%4.0f, %4.0f)\n",
        ps.axis.ra, ps.axis.dec, ps.axis.x, ps.axis.y);
    printf("Offset : %4.0f x, %4.0f y\n",
        ps.ncp.x-ps.axis.x, ps.ncp.y-ps.axis.y);

    // Plot a chart of the polar alignment, if user asks for it
    if (argc == 4)
        plot(&ps, argv[1], argv[3]);

    return 0;
}

max_stars_repo_name: zuh/starflow

max_issues_repo_name: zuh/starflow

max_forks_repo_name: zuh/starflow