text:
/*
Copyright (c) 2015, Patrick Weltevrede
All rights reserved.

Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:

1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.

2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.

3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/

#include <stdio.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>
#include <gsl/gsl_sort.h>
#include "psrsalsa.h"
int main(int argc, char **argv)
{
  int read_log, centered_at_zero, file_column1, file_column2, colspecified, polspecified, filename, truncate, cdf, select;
  long ndata, i, j, k, nrbins, nrbinsy, *distr;
  float *cmap, *distr_float;
  double *data_x, *data_y, min_x_data, max_x_data, min_y_data, max_y_data;
  double min_x, max_x, min_y, max_y, dx, dy, x, y, s, rangex_min, rangex_max, rangey_min, rangey_max, extra_phase, select1, select2;
  int nrbins_specified, dx_specified, dy_specified, output_sigma, output_fraction, nrbins_specifiedy, twoDmode, showGraphics, rangex_set, rangey_set, plotlog, file_column2_defined, xlabelset, linewidth;
  char title[500], *filename_ptr, output_suffix[100], *oname;
  FILE *ofile;
  double min, max;
  psrsalsaApplication application;
  datafile_definition datain;
  pgplot_options_definition pgplot_options;
  pgplot_clear_options(&pgplot_options);
  initApplication(&application, "pdist", "[options] inputfile(s)");
  application.switch_verbose = 1;
  application.switch_debug = 1;
  application.switch_iformat = 1;
  application.switch_formatlist = 1;
  application.switch_device = 1;
  application.switch_cmap = 1;
  application.switch_cmaplist = 1;
  nrbins_specified = 0;
  nrbins_specifiedy = 0;
  dx_specified = 0;
  dy_specified = 0;
  output_sigma = 0;
  output_fraction = 0;
  read_log = 0;
  twoDmode = 0;
  centered_at_zero = 1;
  extra_phase = 0;
  showGraphics = 0;
  xlabelset = 0;
  sprintf(title, "Distribution");
  rangex_set = 0;
  rangey_set = 0;
  file_column1 = 1;
  file_column2 = 2;
  colspecified = 0;
  polspecified = 0;
  plotlog = 0;
  file_column2_defined = 1;
  application.cmap = PPGPLOT_HEAT;
  strcpy(output_suffix, "hist");
  filename = 0;
  truncate = 0;
  cdf = 0;
  select = 0;
  linewidth = 1;
  if(argc < 2) {
    printf("Program to generate or plot a histogram by binning data. Also a cummulative\n");
    printf("distribution can be generated. Usage:\n\n");
    printApplicationHelp(&application);
    printf("Input options:\n");
    printf("-2               Turn on 2D mode: Read in two columns of data, resulting in the\n");
    printf("                 distribution N(x,y) rathern than N(x).\n");
    printf("-col \"c1 c2\"     Specify two column numbers (counting from 1) to be used.\n");
    printf("                 Without the -2 option only one value needs to be specified.\n");
    printf("                 This option implies the input files are ascii files with each\n");
    printf("                 line having an equal number of columns. Lines starting with a #\n");
    printf("                 will be ignored.\n");
    printf("-pol \"p1 p2\"     This option implies that the input file is recognized as a\n");
    printf("                 pulsar format. Polarization p1 (and p2 in 2D mode), counting\n");
    printf("                 from zero, are used (all bins, freqs and subints). The default\n");
    printf("                 default is \"0 1\", or the integrated onpulse pulse energy for\n");
    printf("                 penergy output files in mode 1.\n");
    printf("\nOutput options:\n");
    printf("-ext             Specify suffix, default is '%s'\n", output_suffix);
    printf("-output filename Write output to this file [def=.%s extension].\n", output_suffix);
    printf("-frac            Output fraction of counts per bin rather than counts.\n");
    printf("-sigma           Generate extra column with the sqrt of the nr of counts.\n");
    printf("\nDistribution options:\n");
    printf("-cdf             The cummulative distribution is generated. Use no -dx or -n.\n");
    printf("-dx value        Specify bin width (can also use -n).\n");
    printf("-dy value        Specify bin width used for the second column if -2 is used.\n");
    printf("-log             The base-10 log of the input values is used.\n");
    printf("-n nr            The range of input values is divided in nr bins.\n");
    printf("-nx nr           Same as -n\n");
    printf("-ny nr           Similar to -nx, but now for the second column if -2 is used.\n");
    printf("-rangex \"x1 x2\"  The range of bins produced in output will cover x1..x2, which\n");
    printf("                 should include all input values (see also -trunc). It therefore\n");
    printf("                 allows more bins to be generated than strictly necessary.\n");
    printf("-rangey          As -rangex, but now for second input column.\n");
    printf("-select \"v1 v2\"  Enabling this option will no longer result in a distribution.\n");
    printf("                 Instead all samples within the range of values specified with\n");
    printf("                 v1 and v2 will be outputted as two columns: the sample number\n");
    printf("                 and the value.\n");
    printf("-trunc           Modifies behaviour of -rangex and -rangey. Values outside the\n");
    printf("                 specified range are set to the extremes of the range.\n");
    printf("-zero            By default one of the bins will be centred at zero. This option\n");
    printf("                 results in the edge of a bin to coincides with zero. In both\n");
    printf("                 cases the reported locations of the bins are their centres.\n");
    printf("                 This means that in gnuplot you want to use \"with histeps\".\n");
    printf("-zeroshift off   Put in an extra offset off to where zero would fall with\n");
    printf("                 respect to a bin. -zero is equivalent to -zeroshift 0.5. \n");
    printf("\nPlot options:\n");
    printf("-plot            Plot the distribution rather than writing output to file.\n");
    printf("-plotlog         The log10 of the counts is plotted (implies -plot and bins with\n");
    printf("                 zero counts are set to 1).\n");
    printf("-overplot        Like plot, but multiple input files are overplotted. The ranges\n");
    printf("                 are determined by first input file.\n");
    printf("-title  'title'  Set the title of the plot.\n");
    printf("-xlabel 'label'  Set label on horizontal axis.\n");
    printf("-labels           \"label_ch label_lw label_f box_ch box_lw box_f\"\n");
    printf("                  default is \"%.1f %d %d %.1f %d %d\"\n", pgplot_options.box.label_ch, pgplot_options.box.label_lw, pgplot_options.box.label_f, pgplot_options.box.box_labelsize, pgplot_options.box.box_lw, pgplot_options.box.box_f);
    printf("-lw              Set line width of the histogram.\n");
    printf("-vp              \"left right bottom top\"\n");
    printf("                 Modify the dimensions of the plot. Default is \"%.2f %.2f %.2f %.2f\"\n", pgplot_options.viewport.dxplot, pgplot_options.viewport.xsize, pgplot_options.viewport.dyplot, pgplot_options.viewport.ysize);
    printf("\n");
    printf("Please use the appropriate citation when using results of this software in your publications:\n\n");
    printf("More information about fitting distributions (in the context of pulse energies) can be found in:\n");
    printf(" - Weltevrede et al. 2006, A&A, 458, 269\n\n");
    printCitationInfo();
    terminateApplication(&application);
    return 0;
  }else {
    for(i = 1; i < argc; i++) {
      int index;
      index = i;
      if(processCommandLine(&application, argc, argv, &index)) {
 i = index;
      }else if(strcmp(argv[i], "-2") == 0) {
 twoDmode = 1;
      }else if(strcmp(argv[i], "-plot") == 0) {
 showGraphics = 1;
      }else if(strcmp(argv[i], "-overplot") == 0) {
 showGraphics = 2;
      }else if(strcmp(argv[i], "-plotlog") == 0) {
 plotlog = 1;
 if(showGraphics == 0)
   showGraphics = 1;
      }else if(strcmp(argv[i], "-xlabel") == 0) {
 xlabelset = i+1;
 i++;
      }else if(strcmp(argv[i], "-title") == 0) {
 strcpy(title, argv[i+1]);
 i++;
      }else if(strcmp(argv[i], "-output") == 0) {
 filename = i+1;
 i++;
      }else if(strcmp(argv[i], "-ext") == 0) {
 strcpy(output_suffix, argv[i+1]);
 i++;
      }else if(strcmp(argv[i], "-labels") == 0) {
 if(parse_command_string(application.verbose_state, argc, argv, i+1, 0, -1, "%f %d %d %f %d %d", &(pgplot_options.box.label_ch), &(pgplot_options.box.label_lw), &(pgplot_options.box.label_f), &(pgplot_options.box.box_labelsize), &(pgplot_options.box.box_lw), &(pgplot_options.box.box_f), NULL) == 0) {
   printerror(application.verbose_state.debug, "ERROR pdist: Cannot parse '%s' option.", argv[i]);
   return 0;
 }
 i++;
      }else if(strcmp(argv[i], "-lw") == 0) {
 if(parse_command_string(application.verbose_state, argc, argv, i+1, 0, -1, "%d", &linewidth, NULL) == 0) {
   printerror(application.verbose_state.debug, "ERROR pdist: Cannot parse '%s' option.", argv[i]);
   return 0;
 }
 i++;
      }else if(strcmp(argv[i], "-vp") == 0) {
 if(parse_command_string(application.verbose_state, argc, argv, i+1, 0, -1, "%f %f %f %f", &(pgplot_options.viewport.dxplot), &(pgplot_options.viewport.xsize), &(pgplot_options.viewport.dyplot), &(pgplot_options.viewport.ysize), NULL) == 0) {
   printerror(application.verbose_state.debug, "ERROR pdist: Cannot parse '%s' option.", argv[i]);
   return 0;
 }
 i++;
      }else if(strcmp(argv[i], "-col") == 0 || strcmp(argv[i], "-pol") == 0) {
 if(strcmp(argv[i], "-pol") == 0)
   polspecified = 1;
 else
   colspecified = 1;
 int ret;
 ret = parse_command_string(application.verbose_state, argc, argv, i+1, 0, 1, "%d %d", &file_column1, &file_column2, NULL);
 file_column2_defined = 0;
 if(ret == 2) {
   file_column2_defined = 1;
 }else if(ret == 0) {
   printerror(application.verbose_state.debug, "ERROR pdist: Cannot parse %s option, need 1 or 2 integer values.\n", argv[i]);
   return 0;
 }
 i++;
      }else if(strcmp(argv[i], "-n") == 0 || strcmp(argv[i], "-nx") == 0) {
 if(dx_specified) {
   printerror(application.verbose_state.debug, "Cannot use -n and -dx option simultaniously.\n");
   return 0;
 }
 nrbins_specified = 1;
 if(parse_command_string(application.verbose_state, argc, argv, i+1, 0, -1, "%ld", &nrbins, NULL) == 0) {
   printerror(application.verbose_state.debug, "ERROR pdist: Cannot parse '%s' option.", argv[i]);
   return 0;
 }
 i++;
      }else if(strcmp(argv[i], "-ny") == 0) {
 if(dy_specified) {
   printerror(application.verbose_state.debug, "Cannot use -ny and -dy option simultaniously.\n");
   return 0;
 }
 nrbins_specifiedy = 1;
 if(parse_command_string(application.verbose_state, argc, argv, i+1, 0, -1, "%ld", &nrbinsy, NULL) == 0) {
   printerror(application.verbose_state.debug, "ERROR pdist: Cannot parse '%s' option.", argv[i]);
   return 0;
 }
 i++;
      }else if(strcmp(argv[i], "-dx") == 0) {
 if(nrbins_specified) {
   printerror(application.verbose_state.debug, "Cannot use -n and -dx option simultaniously.\n");
   return 0;
 }
 dx_specified = 1;
 if(parse_command_string(application.verbose_state, argc, argv, i+1, 0, -1, "%lf", &dx, NULL) == 0) {
   printerror(application.verbose_state.debug, "ERROR pdist: Cannot parse '%s' option.", argv[i]);
   return 0;
 }
 i++;
      }else if(strcmp(argv[i], "-dy") == 0) {
 if(nrbins_specifiedy) {
   printerror(application.verbose_state.debug, "Cannot use -ny and -dy option simultaniously.\n");
   return 0;
 }
 dy_specified = 1;
 if(parse_command_string(application.verbose_state, argc, argv, i+1, 0, -1, "%lf", &dy, NULL) == 0) {
   printerror(application.verbose_state.debug, "ERROR pdist: Cannot parse '%s' option.", argv[i]);
   return 0;
 }
 i++;
      }else if(strcmp(argv[i], "-rangex") == 0) {
 rangex_set = 1;
 if(parse_command_string(application.verbose_state, argc, argv, i+1, 0, -1, "%lf %lf", &rangex_min, &rangex_max, NULL) == 0) {
   printerror(application.verbose_state.debug, "ERROR pdist: Cannot parse '%s' option.", argv[i]);
   return 0;
 }
 i++;
      }else if(strcmp(argv[i], "-rangey") == 0) {
 rangey_set = 1;
 if(parse_command_string(application.verbose_state, argc, argv, i+1, 0, -1, "%lf %lf", &rangey_min, &rangey_max, NULL) == 0) {
   printerror(application.verbose_state.debug, "ERROR pdist: Cannot parse '%s' option.", argv[i]);
   return 0;
 }
 i++;
      }else if(strcmp(argv[i], "-select") == 0) {
 select = 1;
 if(parse_command_string(application.verbose_state, argc, argv, i+1, 0, -1, "%lf %lf", &select1, &select2, NULL) == 0) {
   printerror(application.verbose_state.debug, "ERROR pdist: Cannot parse '%s' option.", argv[i]);
   return 0;
 }
 i++;
      }else if(strcmp(argv[i], "-trunc") == 0) {
 truncate = 1;
      }else if(strcmp(argv[i], "-sigma") == 0) {
 output_sigma = 1;
      }else if(strcmp(argv[i], "-log") == 0) {
 read_log = 1;
      }else if(strcmp(argv[i], "-zero") == 0) {
 centered_at_zero = 0;
      }else if(strcmp(argv[i], "-zeroshift") == 0) {
 if(parse_command_string(application.verbose_state, argc, argv, i+1, 0, -1, "%lf", &extra_phase, NULL) == 0) {
   printerror(application.verbose_state.debug, "ERROR pdist: Cannot parse '%s' option.", argv[i]);
   return 0;
 }
 i++;
      }else if(strcmp(argv[i], "-frac") == 0) {
 output_fraction = 1;
      }else if(strcmp(argv[i], "-cdf") == 0) {
 cdf = 1;
      }else {
 if(argv[i][0] == '-') {
   printerror(application.verbose_state.debug, "pdist: Unknown option: %s\n\nRun pdist without command line arguments to show help", argv[i]);
   terminateApplication(&application);
   return 0;
 }else {
   if(applicationAddFilename(i, application.verbose_state) == 0)
     return 0;
 }
      }
    }
  }
  if(applicationFilenameList_checkConsecutive(argv, application.verbose_state) == 0) {
    return 0;
  }
  if(numberInApplicationFilenameList(&application, argv, application.verbose_state) == 0) {
    printerror(application.verbose_state.debug, "ERROR pdist: No files specified");
    return 0;
  }
  if(select) {
    if(twoDmode) {
      printerror(application.verbose_state.debug, "ERROR pdist: The -2 option cannot be used with -select.\n");
      return 0;
    }
    if(cdf) {
      printerror(application.verbose_state.debug, "ERROR pdist: The -cdf option cannot be used with -select.\n");
      return 0;
    }
    if(showGraphics) {
      printerror(application.verbose_state.debug, "ERROR pdist: No plots can be generated when using -select.\n");
      return 0;
    }
    if(output_fraction) {
      printerror(application.verbose_state.debug, "ERROR pdist: The -frac option cannot be used with -select.\n");
      return 0;
    }
    if(output_sigma) {
      printerror(application.verbose_state.debug, "ERROR pdist: The -sigma option cannot be used with -select.\n");
      return 0;
    }
  }
  if(cdf) {
    if(nrbins_specified || dx_specified) {
      printerror(application.verbose_state.debug, "ERROR pdist: For a cdf -n and -dx should not be used.\n");
      return 0;
    }
    if(twoDmode) {
      printerror(application.verbose_state.debug, "ERROR pdist: For a cdf only one input column is expected.\n");
      return 0;
    }
    if(output_sigma) {
      printerror(application.verbose_state.debug, "ERROR pdist: For a cdf -sigma is not supported.\n");
      return 0;
    }
    if(rangex_set || rangey_set || truncate) {
      printerror(application.verbose_state.debug, "ERROR pdist: For a cdf the -rangex and -rangey options are not supported.\n");
      return 0;
    }
  }else {
    if(select == 0) {
      if(nrbins_specified == 0 && dx_specified == 0) {
 printerror(application.verbose_state.debug, "ERROR pdist: Specify resolution with the -n or -dx option.\n");
 return 0;
      }
      if(twoDmode) {
 if(nrbins_specifiedy == 0 && dy_specified == 0) {
 printerror(application.verbose_state.debug, "ERROR pdist: Specify resolution with the -ny or -dy option.\n");
 return 0;
 }
      }
    }
  }
  int initial_iformat, currentfile;
  currentfile = 1;
  initial_iformat = application.iformat;
  while((filename_ptr = getNextFilenameFromList(&application, argv, application.verbose_state)) != NULL) {
    application.iformat = initial_iformat;
    if(application.iformat <= 0 && colspecified == 0) {
      application.iformat = guessPSRData_format(filename_ptr, 1, application.verbose_state);
      if(application.iformat == -2 || application.iformat == -3)
 return 0;
    }
    if(colspecified)
      application.iformat = -1;
    if(polspecified && application.iformat <= 0) {
      printerror(application.verbose_state.debug, "ERROR pdist: Input file cannot be opened. Please check if file %s exists and otherwise specify the correct input format with the -iformat option if the format is supported, but not automatically recognized. Data is not recognized as pulsar data, but -pol was used. Use -col for ascii files.\n", filename_ptr);
      return 0;
    }
    if(application.iformat > 0) {
      verbose_definition verbose;
      cleanVerboseState(&verbose);
      copyVerboseState(application.verbose_state, &verbose);
      if(application.iformat == PSRCHIVE_ASCII_format) {
 if(verbose.debug == 0)
   verbose.verbose = 0;
      }
      i = openPSRData(&datain, filename_ptr, application.iformat, 0, 1, 1, verbose);
      if(i == 0) {
 printerror(application.verbose_state.debug, "ERROR pdist: Error opening data");
 return 0;
      }
      if(polspecified == 0) {
 file_column1 = 0;
 file_column2 = 1;
      }
      if(application.iformat == PSRCHIVE_ASCII_format && datain.gentype == GENTYPE_PENERGY) {
 if(application.verbose_state.verbose)
   printf("The file is a penergy file generated in mode 1.\n");
 if(polspecified == 0) {
   printwarning(application.verbose_state.debug, "WARNING pdist: -pol not specified, default in the integrated on-pulse energy (-pol 2)");
   file_column1 = 2;
 }
      }
      if(application.verbose_state.verbose)
 printf("Loading %ld points", datain.NrSubints*datain.NrFreqChan*datain.NrBins);
      if(twoDmode)
 printf(" times two input polarizations");
      printf("\n");
      data_x = malloc(datain.NrSubints*datain.NrFreqChan*datain.NrBins*sizeof(double));
      data_y = malloc(datain.NrSubints*datain.NrFreqChan*datain.NrBins*sizeof(double));
      if(data_x == NULL || data_y == NULL) {
 printerror(application.verbose_state.debug, "ERROR pdist: Memory allocation error.\n");
 return 0;
      }
      long subintnr, freqnr, binnr;
      long double total_x, total_y;
      ndata = 0;
      total_x = 0;
      total_y = 0;
      min_x_data = max_x_data = 0;
      for(subintnr = 0; subintnr < datain.NrSubints; subintnr++) {
 for(freqnr = 0; freqnr < datain.NrFreqChan; freqnr++) {
   for(binnr = 0; binnr < datain.NrBins; binnr++) {
     float dummy_float;
     if(readPulsePSRData(&datain, subintnr, file_column1, freqnr, binnr, 1, &dummy_float, application.verbose_state) == 0) {
       printerror(application.verbose_state.debug, "ERROR pdist: Read error, shouldn't happen.\n");
       return 0;
     }
     data_x[ndata] = dummy_float;
     if(read_log) {
       if(data_x[ndata] <= 0) {
  printerror(application.verbose_state.debug, "ERROR pdist: Cannot take logarithm of a value <= 0.\n");
  return 0;
       }
       data_x[ndata] = log10(data_x[ndata]);
     }
     if(data_x[ndata] > max_x_data || ndata == 0) {
       max_x_data = data_x[ndata];
     }
     if(data_x[ndata] < min_x_data || ndata == 0) {
       min_x_data = data_x[ndata];
     }
     total_x += data_x[ndata];
     if(twoDmode) {
       if(readPulsePSRData(&datain, subintnr, file_column2, freqnr, binnr, 1, &dummy_float, application.verbose_state) == 0) {
  printerror(application.verbose_state.debug, "ERROR pdist: Read error, shouldn't happen.\n");
  return 0;
       }
       data_y[ndata] = dummy_float;
       if(read_log) {
  if(data_y[ndata] <= 0) {
    printerror(application.verbose_state.debug, "ERROR pdist: Cannot take logarithm of a value <= 0.\n");
    return 0;
  }
  data_y[ndata] = log10(data_y[ndata]);
       }
       if(ndata == 0) {
  max_y_data = data_y[ndata];
  min_y_data = data_y[ndata];
       }
       if(data_y[ndata] > max_y_data) {
  max_y_data = data_y[ndata];
       }
       if(data_y[ndata] < min_y_data) {
  min_y_data = data_y[ndata];
       }
       total_y += data_y[ndata];
     }
     ndata++;
   }
 }
      }
      if(application.verbose_state.verbose) {
 printf("xrange = %e ... %e\n", min_x_data, max_x_data);
 printf("average = %Le\n", total_x/(long double)ndata);
 if(twoDmode) {
   printf("yrange = %e ... %e\n", min_y_data, max_y_data);
   printf("average = %Le\n", total_y/(long double)ndata);
 }
      }
      closePSRData(&datain, 0, application.verbose_state);
    }else {
      int skiplines = 0;
      if(twoDmode) {
 if(application.verbose_state.verbose)
   fprintf(stdout, "Loading x values from ascii file\n");
 if(read_ascii_column_double(filename_ptr, skiplines, '#', -1, 1, &ndata, file_column1, 1.0, read_log, &data_x, &min_x_data, &max_x_data, NULL, application.verbose_state, 1) == 0) {
   printerror(application.verbose_state.debug, "ERROR pdist: cannot load file.\n");
   if(colspecified) {
     printwarning(application.verbose_state.debug, "WARNING pdist: Using the -col option implies the input file is a simple ascii file. For penergy output (in mode 1), or an other recognized pulsar format, use the -pol option instead.\n");
   }
   return 0;
 }
 if(application.verbose_state.verbose)
   fprintf(stdout, "Loading y values from ascii file\n");
 if(file_column2_defined == 0) {
   printerror(application.verbose_state.debug, "In 2D mode, two input columns should be specified with the -col option.\n");
   return 0;
 }
 if(read_ascii_column_double(filename_ptr, skiplines, '#', -1, 1, &ndata, file_column2, 1.0, read_log, &data_y, &min_y_data, &max_y_data, NULL, application.verbose_state, 1) == 0) {
   printerror(application.verbose_state.debug, "ERROR pdist: cannot load file.\n");
   if(colspecified) {
     printwarning(application.verbose_state.debug, "WARNING pdist: Using the -col option implies the input file is a simple ascii file. For penergy output (in mode 1), or an other recognized pulsar format, use the -pol option instead.\n");
   }
   return 0;
 }
      }else {
 if(application.verbose_state.verbose)
   fprintf(stdout, "Loading values from ascii file\n");
 if(read_ascii_column_double(filename_ptr, skiplines, '#', -1, 1, &ndata, file_column1, 1.0, read_log, &data_x, &min_x_data, &max_x_data, NULL, application.verbose_state, 1) == 0) {
   printerror(application.verbose_state.debug, "ERROR pdist: cannot load file.\n");
   if(colspecified) {
     printwarning(application.verbose_state.debug, "WARNING pdist: Using the -col option implies the input file is a simple ascii file. For penergy output (in mode 1), or an other recognized pulsar format, use the -pol option instead.\n");
   }
   return 0;
 }
 data_y = malloc(ndata*sizeof(double));
 if(data_y == NULL) {
   printerror(application.verbose_state.debug, "ERROR pdist: Memory allocation error.\n");
   return 0;
 }
      }
    }
    if(cdf == 0 && select == 0) {
      switch(set_binning_histogram(min_x_data, max_x_data, rangex_set, rangex_min, rangex_max, nrbins_specified, nrbins, centered_at_zero, extra_phase, &min_x, &max_x, &dx, application.verbose_state)) {
      case 0: break;
      case 1:
      case 2: return 0;
      default: {
 printerror(application.verbose_state.debug, "ERROR pdist: Unknown return value for call to set_binning_histogram.\n");
 return 0;
      }
      }
      if(rangex_set) {
 if(min_x > min_x_data || max_x < max_x_data) {
   if(truncate == 0) {
     printerror(application.verbose_state.debug, "ERROR pdist: Values are outside specified -rangex option. You may want to use -trunc?\n");
     return 0;
   }
 }
      }
      if(twoDmode) {
 switch(set_binning_histogram(min_y_data, max_y_data, rangey_set, rangey_min, rangey_max, nrbins_specifiedy, nrbinsy, centered_at_zero, extra_phase, &min_y, &max_y, &dy, application.verbose_state)) {
 case 0: break;
 case 1:
 case 2: return 0;
 default: {
   printerror(application.verbose_state.debug, "ERROR pdist: Unknown return value for call to set_binning_histogram.\n");
   return 0;
 }
 }
 if(application.verbose_state.verbose) {
   fprintf(stdout, "Going to use binsize dx=%e and dy=%e.\n", dx, dy);
 }
      }
      if(rangey_set) {
 if(min_y > min_y_data || max_y < max_y_data) {
   if(truncate == 0) {
     printerror(application.verbose_state.debug, "ERROR pdist: Values are outside specified -rangey option. You may want to use -trunc?\n");
     return 0;
   }
 }
      }
    }
    if(cdf == 0) {
      if(truncate) {
 if(rangex_set) {
   for(i = 0; i < ndata; i++) {
     if(data_x[i] < min_x)
       data_x[i] = min_x;
     if(data_x[i] > max_x)
       data_x[i] = max_x;
   }
 }
 if(rangey_set) {
   for(i = 0; i < ndata; i++) {
     if(data_y[i] < min_y)
       data_y[i] = min_y;
     if(data_y[i] > max_y)
       data_y[i] = max_y;
   }
 }
      }
    }
    if(cdf == 0 && select == 0)
      nrbins = calculate_bin_number(max_x, dx, min_x, centered_at_zero, extra_phase) + 1;
    else
      nrbins = ndata;
    if(select == 0)
      fprintf(stdout, "Distribution will have %ld bins.\n", nrbins);
    if(twoDmode) {
      nrbinsy = calculate_bin_number(max_y, dy, min_y, centered_at_zero, extra_phase)+ 1;
      fprintf(stdout, "Distribution will have %ld y-bins.\n", nrbinsy);
    }else {
      nrbinsy = 1;
    }
    distr = (long *)malloc(nrbins*nrbinsy*sizeof(long));
    if(distr == NULL) {
      printerror(application.verbose_state.debug, "ERROR pdist: Cannot allocate memory.\n");
      return 0;
    }
    for(i = 0; i < nrbins*nrbinsy; i++)
      distr[i] = 0;
    if(showGraphics) {
      pgplot_options.box.drawtitle = 1;
      strcpy(pgplot_options.box.title, title);
      if(xlabelset) {
 strcpy(pgplot_options.box.xlabel, argv[xlabelset]);
      }else {
 if(read_log)
   strcpy(pgplot_options.box.xlabel, "log X");
 else
   strcpy(pgplot_options.box.xlabel, "X");
      }
      if(twoDmode) {
 if(read_log)
   strcpy(pgplot_options.box.ylabel, "log Y");
 else
   strcpy(pgplot_options.box.ylabel, "Y");
      }else {
 if(plotlog)
   strcpy(pgplot_options.box.ylabel, "log N");
 else
   strcpy(pgplot_options.box.ylabel, "N");
      }
      strcpy(pgplot_options.viewport.plotDevice, application.pgplotdevice);
      if(showGraphics == 2) {
 pgplot_options.viewport.dontclose = 1;
      }
      if(currentfile > 1 && showGraphics == 2) {
 pgplot_options.viewport.noclear = 1;
 pgplot_options.viewport.dontopen = 1;
 pgplot_options.box.drawbox = 0;
 pgplot_options.box.drawtitle = 0;
 pgplot_options.box.drawlabels = 0;
      }
    }else {
      if(filename != 0) {
 oname = (char *) calloc(strlen(argv[filename])+1, 1);
      }else {
 oname = (char *) calloc(strlen(filename_ptr)+strlen(output_suffix)+2, 1);
      }
      if(oname == NULL) {
 printerror(application.verbose_state.debug, "ERROR pdist: Memory allocation error.\n");
 return 0;
      }
      if(filename != 0) {
 strcpy(oname, argv[filename]);
      }else {
 sprintf(oname,"%s.%s", filename_ptr, output_suffix);
      }
      printf("Output Ascii to %s\n",oname);
      ofile = fopen(oname,"w+");
      if(ofile == NULL) {
 printerror(application.verbose_state.debug, "ERROR pdist: Cannot open '%s'", oname);
 perror("");
 return 0;
      }
    }
    if(twoDmode) {
      for(i = 0; i < ndata; i++) {
 j = calculate_bin_number(data_x[i], dx, min_x, centered_at_zero, extra_phase);
 k = calculate_bin_number(data_y[i], dy, min_y, centered_at_zero, extra_phase);
 if(j < 0 || j >= nrbins) {
   printerror(application.verbose_state.debug, "BUG!\n");
   return 0;
 }
 if(k < 0 || k >= nrbinsy) {
   printerror(application.verbose_state.debug, "BUG!\n");
   return 0;
 }
 distr[k*nrbins+j] += 1;
      }
      for(i = 0; i < nrbins; i++) {
 for(j = 0; j < nrbinsy; j++) {
   x = calculate_bin_location(i, dx, min_x, centered_at_zero, extra_phase);
   y = calculate_bin_location(j, dy, min_y, centered_at_zero, extra_phase);
   if(!showGraphics) {
       fprintf(ofile, "%e %e", x, y);
     if(!output_fraction)
       fprintf(ofile, " %ld", distr[j*nrbins+i]);
     else
       fprintf(ofile, " %e", distr[j*nrbins+i]/(double)ndata);
     if(output_sigma) {
       if(distr[j*nrbins+i] == 0)
  s = 1;
       else
  s = sqrt(distr[j*nrbins+i]);
       if(output_fraction)
  s /= (double)ndata;
  fprintf(ofile, " %e", s);
     }
   }
   if(!showGraphics) {
     fprintf(ofile, "\n");
   }
 }
      }
      if(showGraphics) {
 cmap = (float *)malloc(nrbins*nrbinsy*sizeof(float));
 if(cmap == NULL) {
   printerror(application.verbose_state.debug, "Cannot allocate memory.\n");
   return 0;
 }
 if(plotlog) {
   if(distr[0] > 0) {
     min = log10(distr[0]);
     max = log10(distr[0]);
   }else {
     min = max = 0;
   }
 }else {
   min = distr[0];
   max = distr[0];
 }
 for(i = 0; i < nrbins; i++) {
   for(j = 0; j <nrbinsy; j++) {
     if(plotlog) {
       if(distr[i*nrbinsy+j] > 0)
  cmap[i*nrbinsy+j] = log10(distr[i*nrbinsy+j]);
       else
  cmap[i*nrbinsy+j] = 0;
     }
     else
       cmap[i*nrbinsy+j] = distr[i*nrbinsy+j];
     if(cmap[i*nrbinsy+j] > max)
       max = cmap[i*nrbinsy+j];
     if(cmap[i*nrbinsy+j] < min)
       min = cmap[i*nrbinsy+j];
   }
 }
 printf("Count range: %f to %f\n", min, max);
 int showwedge = 0;
 if(pgplotMap(&pgplot_options, cmap, nrbins, nrbinsy,
       calculate_bin_location(0, dx, min_x, centered_at_zero, extra_phase),
       calculate_bin_location(0, dx, max_x, centered_at_zero, extra_phase),
       calculate_bin_location(0, dx, min_x, centered_at_zero, extra_phase)-0.5*dx,
       calculate_bin_location(0, dx, max_x, centered_at_zero, extra_phase)+0.5*dx,
       calculate_bin_location(0, dy, min_y, centered_at_zero, extra_phase),
       calculate_bin_location(0, dy, max_y, centered_at_zero, extra_phase),
       calculate_bin_location(0, dy, min_y, centered_at_zero, extra_phase)-0.5*dy,
       calculate_bin_location(0, dy, max_y, centered_at_zero, extra_phase)+0.5*dy,
       application.cmap, 0, 0, 0, NULL, 0, 0, 1.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, showwedge, 0, 0, application.verbose_state) == 0) {
   printerror(application.verbose_state.debug, "ERROR pdist: Plotting failed.\n");
   return 0;
 }
 free(cmap);
      }
    }else {
      if(select) {
 nrbins = 0;
 for(i = 0; i < ndata; i++) {
   if(data_x[i] >= select1 && data_x[i] <= select2) {
     nrbins++;
   }
 }
 j = 0;
 for(i = 0; i < ndata; i++) {
   if(data_x[i] >= select1 && data_x[i] <= select2) {
     distr[j] = i;
     data_y[j] = data_x[i];
     j++;
   }
 }
      }else {
 if(cdf == 0) {
   for(i = 0; i < ndata; i++) {
     j = calculate_bin_number(data_x[i], dx, min_x, centered_at_zero, extra_phase);
     if(j < 0 || j >= nrbins) {
       printerror(application.verbose_state.debug, "BUG! bin number %ld outside range.\n", j);
       exit(0);
     }
     distr[j] += 1;
   }
 }else {
   gsl_sort(data_x, 1, ndata);
   for(i = 0; i < ndata; i++) {
     data_y[i] = (i+1)/(double)(ndata);
   }
 }
      }
      if(!showGraphics) {
 for(i = 0; i < nrbins; i++) {
   x = i*dx;
     if(cdf == 0) {
       if(select == 0) {
  fprintf(ofile, "%e", calculate_bin_location(i, dx, min_x, centered_at_zero, extra_phase));
       }else {
  fprintf(ofile, "%ld", distr[i]);
       }
     }else {
       fprintf(ofile, "%e", data_x[i]);
     }
   if(!output_fraction) {
       if(cdf == 0) {
  if(select == 0) {
    fprintf(ofile, " %ld", distr[i]);
  }else {
    fprintf(ofile, " %e", data_y[i]);
  }
       }else {
  fprintf(ofile, " %e", data_y[i]);
       }
   }else {
       if(cdf == 0) {
  fprintf(ofile, " %e", distr[i]/(double)ndata);
       }else {
  fprintf(ofile, " %e", data_y[i]);
       }
   }
   if(output_sigma) {
     if(distr[i] == 0)
       s = 1;
     else
       s = sqrt(distr[i]);
     if(output_fraction)
       s /= (double)ndata;
     fprintf(ofile, " %e", s);
   }
   fprintf(ofile, "\n");
 }
      }else {
 float *distr_x_float;
 distr_float = malloc(nrbins*sizeof(float));
 if(distr_float == NULL) {
   printerror(application.verbose_state.debug, "ERROR pdist: Cannot allocate memory\n");
   return 0;
 }
 if(cdf == 0) {
   for(i = 0; i < nrbins; i++) {
     distr_float[i] = distr[i];
     if(output_fraction)
       distr_float[i] /= (double)ndata;
     if(plotlog) {
       if(distr[i] > 0) {
  distr_float[i] = log10(distr_float[i]);
       }else {
  if(output_fraction) {
    distr_float[i] = log10(1.0/(float)ndata);
  }else {
    distr_float[i] = 0;
  }
       }
     }
   }
 }else {
   distr_x_float = malloc(nrbins*sizeof(float));
   if(distr_x_float == NULL) {
     printerror(application.verbose_state.debug, "ERROR pdist: Cannot allocate memory\n");
     return 0;
   }
   for(i = 0; i < nrbins; i++) {
     distr_float[i] = data_y[i];
     distr_x_float[i] = data_x[i];
   }
 }
 int forceMinZero, dontsetranges, colour;
 forceMinZero = 1;
 if(plotlog)
   forceMinZero = 0;
 dontsetranges = 0;
 colour = 1;
 if(showGraphics == 2) {
   if(currentfile > 1)
     dontsetranges = 1;
   colour = currentfile;
 }
 if(cdf == 0) {
   if(pgplotGraph1(&pgplot_options, distr_float, NULL, NULL, nrbins,
     calculate_bin_location(0, dx, min_x, centered_at_zero, extra_phase),
     calculate_bin_location(nrbins-1, dx, min_x, centered_at_zero, extra_phase),
     dontsetranges,
     calculate_bin_location(0, dx, min_x, centered_at_zero, extra_phase) - dx,
     calculate_bin_location(nrbins-1, dx, min_x, centered_at_zero, extra_phase) + dx,
     0, 0, forceMinZero, 1, 0, linewidth, 0, colour, 1, NULL, -1, application.verbose_state) == 0) {
     printerror(application.verbose_state.debug, "ERROR pdist: Cannot plot graph\n");
     return 0;
   }
 }else {
   if(pgplotGraph1(&pgplot_options, distr_float, distr_x_float, NULL, nrbins,
     0,
     0,
     dontsetranges,
     0,
     0,
     0, 0, forceMinZero, 1*0, 0, linewidth, 0, colour, 1, NULL, -1, application.verbose_state) == 0) {
     printerror(application.verbose_state.debug, "ERROR pdist: Cannot plot graph\n");
     return 0;
   }
   free(distr_x_float);
 }
 free(distr_float);
      }
    }
    free(data_x);
    free(data_y);
    free(distr);
    if(showGraphics == 0) {
      fclose(ofile);
      free(oname);
    }
    if(showGraphics == 1) {
      if(currentfile != numberInApplicationFilenameList(&application, argv, application.verbose_state)) {
 i = pgplot_device_type(application.pgplotdevice, application.verbose_state);
 if(i < 3 || i > 10) {
   printf("Press a key to continue\n");
   fflush(stdout);
   pgetch();
 }
      }
    }
    currentfile += 1;
  }
  if(showGraphics == 2) {
    ppgend();
  }
  terminateApplication(&application);
  return 0;
}

max_stars_repo_name: David-McKenna/psrsalsa

max_issues_repo_name: David-McKenna/psrsalsa

max_forks_repo_name: David-McKenna/psrsalsa


text:
#include <gsl/gsl_combination.h>

#include <caml/alloc.h>
#include <caml/bigarray.h>
#include <caml/memory.h>

static void combi_of_val(gsl_combination *c, value vc)
{
    c->n = Int_val(Field(vc, 0));
    c->k = Int_val(Field(vc, 1));
    c->data = Data_bigarray_val(Field(vc, 2));
}

CAMLprim value ml_gsl_combination_init_first(value vc)
{
    gsl_combination c;
    combi_of_val(&c, vc);
    gsl_combination_init_first(&c);
    return Val_unit;
}

CAMLprim value ml_gsl_combination_init_last(value vc)
{
    gsl_combination c;
    combi_of_val(&c, vc);
    gsl_combination_init_last(&c);
    return Val_unit;
}

CAMLprim value ml_gsl_combination_valid(value vc)
{
    int r;
    gsl_combination c;
    combi_of_val(&c, vc);
    r = gsl_combination_valid(&c);
    return Val_not(Val_bool(r));
}

CAMLprim value ml_gsl_combination_next(value vc)
{
    gsl_combination c;
    combi_of_val(&c, vc);
    gsl_combination_next(&c);
    return Val_unit;
}

CAMLprim value ml_gsl_combination_prev(value vc)
{
    gsl_combination c;
    combi_of_val(&c, vc);
    gsl_combination_prev(&c);
    return Val_unit;
}

ext: c

max_stars_repo_name: nrlucaroni/pareto

max_issues_repo_name: nrlucaroni/pareto

max_forks_repo_name: nrlucaroni/pareto


text:
/*  Copyright (C) 2006 Imperial College London and others.

    Please see the AUTHORS file in the main source directory for a full list
    of copyright holders.

    Prof. C Pain
    Applied Modelling and Computation Group
    Department of Earth Science and Engineering
    Imperial College London

    amcgsoftware@imperial.ac.uk

    This library is free software; you can redistribute it and/or
    modify it under the terms of the GNU Lesser General Public
    License as published by the Free Software Foundation,
    version 2.1 of the License.

    This library 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
    Lesser General Public License for more details.

    You should have received a copy of the GNU Lesser General Public
    License along with this library; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
    USA
*/

#ifndef COMMAND_LINE_OPTIONS_H
#define COMMAND_LINE_OPTIONS_H
#include "confdefs.h"
#include "version.h"

#include "Tokenize.h"

#ifdef _AIX
#include <unistd.h>
#else
#include <getopt.h>
#endif

#ifdef HAVE_MPI
#include <mpi.h>
#endif

#ifdef HAVE_PETSC
#include <petsc.h>
#endif

#include <signal.h>

#include <cstdlib>
#include <iostream>
#include <map>
#include <sstream>
#include <string>
#include <vector>
#include "fmangle.h"

extern std::map<std::string, std::string> fl_command_line_options;

void ParseArguments(int argc, char** argv);
void PetscInit(int argc, char** argv);
void print_version(std::ostream& stream = std::cerr);
void qg_usage(char *cmd);

extern "C" {
  void set_global_debug_level_fc(int *val);
}

#endif

ext: h


text:
/* monte/gsl_monte_vegas.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Michael Booth
 *
 * 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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

/* header for the gsl "vegas" routines.  Mike Booth, May 1998 */

#ifndef __GSL_MONTE_VEGAS_H__
#define __GSL_MONTE_VEGAS_H__

#include <gsl/gsl_rng.h>
#include <gsl/gsl_monte.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

enum {GSL_VEGAS_MODE_IMPORTANCE = 1,
      GSL_VEGAS_MODE_IMPORTANCE_ONLY = 0,
      GSL_VEGAS_MODE_STRATIFIED = -1};

typedef struct {
  /* grid */
  size_t dim;
  size_t bins_max;
  unsigned int bins;
  unsigned int boxes; /* these are both counted along the axes */
  double * xi;
  double * xin;
  double * delx;
  double * weight;
  double vol;

  double * x;
  int * bin;
  int * box;

  /* distribution */
  double * d;

  /* control variables */
  double alpha;
  int mode;
  int verbose;
  unsigned int iterations;
  int stage;

  /* scratch variables preserved between calls to vegas1/2/3  */
  double jac;
  double wtd_int_sum;
  double sum_wgts;
  double chi_sum;
  double chisq;

  double result;
  double sigma;

  unsigned int it_start;
  unsigned int it_num;
  unsigned int samples;
  unsigned int calls_per_box;

  FILE * ostream;

} gsl_monte_vegas_state;

int gsl_monte_vegas_integrate(gsl_monte_function * f,
                              double xl[], double xu[],
                              size_t dim, size_t calls,
                              gsl_rng * r,
                              gsl_monte_vegas_state *state,
                              double* result, double* abserr);

gsl_monte_vegas_state* gsl_monte_vegas_alloc(size_t dim);

int gsl_monte_vegas_init(gsl_monte_vegas_state* state);

void gsl_monte_vegas_free (gsl_monte_vegas_state* state);

__END_DECLS

#endif /* __GSL_MONTE_VEGAS_H__ */

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
#include "multivar_support.h"
#include "multi_modelstruct.h"
#include "emulator_struct.h"
#include "modelstruct.h"
#include "libEmu/estimate_threaded.h"
#include <math.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>

/**
 * estimates the thetas for an allocated multvariate model and outputs the
 * info to outfp
 *
 * requries:
 * m must have been succesfully allocated and pca'd
 * outfp is an open fileptr
 */
void estimate_multi(multi_modelstruct *m, FILE* outfp)
{
	int i;
	for(i = 0; i < m->nr; i++){
		estimate_thetas_threaded(m->pca_model_array[i], m->pca_model_array[i]->options);
	}
	// dump the trained modelstruct, this maybe should go somewhere else
	dump_multi_modelstruct(outfp, m);
}

multi_emulator *alloc_multi_emulator(multi_modelstruct *m)
{
	int i;
	multi_emulator *e = (multi_emulator*)malloc(sizeof(multi_emulator));
	e->nt = m->nt;
	e->nr = m->nr;
	e->nparams = m->nparams;
	e->nmodel_points = m->nmodel_points;
	// have to read these out of the pca array
	e->nregression_fns = m->pca_model_array[0]->options->nregression_fns;
	e->nthetas = m->pca_model_array[0]->options->nthetas;
	e->model = m;

	// allocate the nr emulators we need
	e->emu_struct_array = (emulator_struct**)malloc(sizeof(emulator_struct*)*e->nr);

	for(i = 0; i < e->nr; i++)
		e->emu_struct_array[i] = alloc_emulator_struct(m->pca_model_array[i]);

	return e;
}

/**
 * free allocated memory
 */
void free_multi_emulator(multi_emulator *e)
{
	int i;
	free_multimodelstruct(e->model);
	for(i = 0; i < e->nr; i++)
		free_emulator_struct(e->emu_struct_array[i]);
	free(e->emu_struct_array);
}

/**
 * sample the multivariate emulator at the_point, returns values in the *PCA* space.
 * The returned values will be vectors of length emu->model->nr *not* emu->model->nt
 *
 * requires:
 * emu has been allocated from a multi_modelstruct, and
 * the_mean and the_variance point to allocated gsl_vectors of length (emu->model->nr)
 *
 * @argument emu: a multi_emulator structure created from an estimated multi_modelstruct with alloc_multi_emulator
 * @argument the_point: the location in parameter space (nparams) length vector
 * @argument the_mean: the emulated mean values (nr length)
 * @argument the_variance: the emulated variance values (nr length)
 */

void emulate_point_multi_pca(multi_emulator *emu, gsl_vector *the_point,
												 gsl_vector *the_mean, gsl_vector *the_variance)
{
	int i;
	// sample the nr emulators
	// the results are directly fed back  into the_mean and the_variance
	for(i = 0; i < emu->nr; i++)
		emulate_point(emu->emu_struct_array[i], the_point, gsl_vector_ptr(the_mean, i), gsl_vector_ptr(the_variance, i));
}

/**
 * emulate the values of a point in a multivariate model, returns values in the REAL Observable space.
 * This means that the values output here will be in the same space as the training data.
 *
 * requries:
 * emu has been allocated from a multi_modelstruct
 * the_mean and the_variance point to allocated gsl_vectors of length (emu->model->nt)
 *
 * @argument emu: a multi_emulator structure created from an estimated multi_modelstruct with alloc_multi_emulator
 * @argument the_point: the location in parameter space (nparams) length vector
 * @argument the_mean: the emulated mean values (t length)
 * @argument the_variance: the emulated variance values (t length)
 */
void emulate_point_multi(multi_emulator *emu, gsl_vector *the_point,
												 gsl_vector *the_mean, gsl_vector *the_variance)
{
	int i, j;
	int nr = emu->nr;
	int nt = emu->nt;
	double vec_mat_sum;
	// the mean and variance in the PCA space
	gsl_vector *mean_pca = gsl_vector_alloc(nr);
	gsl_vector *var_pca = gsl_vector_alloc(nr);

	// first we sample the nr emulators
	for(i = 0; i < emu->nr; i++)
		emulate_point(emu->emu_struct_array[i], the_point, gsl_vector_ptr(mean_pca, i), gsl_vector_ptr(var_pca, i));

	/**
	 *  now we have to project back into the REAL space from the PCA space
	 * mean_real (nt) = training_mean (nt) + pca_evec_r %*% diag(sqrt(pca_evals_R)) %*% mean_pca
	 */
	gsl_vector *mean_real = gsl_vector_alloc(nt);
	gsl_vector *temp = gsl_vector_alloc(nt);
	gsl_vector_memcpy(mean_real, emu->model->training_mean);

	vec_mat_sum = 0.0;
	for(i = 0; i < nt; i++){
		for(j = 0; j < nr; j++)
			vec_mat_sum += gsl_matrix_get(emu->model->pca_evecs_r, i, j) * sqrt(gsl_vector_get(emu->model->pca_evals_r, j))*gsl_vector_get(mean_pca, j);
		// save the sum scaled by the sqrt of the eval
		//gsl_vector_set(temp, i, vec_mat_sum*sqrt(gsl_vector_get(emu->model->pca_evals_r, i)));
		gsl_vector_set(temp, i, vec_mat_sum);
		vec_mat_sum = 0.0;
	}

	gsl_vector_add(mean_real, temp);
	gsl_vector_memcpy(the_mean, mean_real); // save the final mean

	/**
	 * project back the variance
	 * yVar_i = Sum_k=0^{nr} ( U_i_k * U_i_k * lambda_k * V_k)
	 */
	gsl_vector_set_zero(temp);
	vec_mat_sum = 0.0;
	for(i = 0; i < nt; i++){
		for(j = 0; j < nr; j++)
			vec_mat_sum += pow(gsl_matrix_get(emu->model->pca_evecs_r, i, j), 2.0) *
				gsl_vector_get(emu->model->pca_evals_r, j) * gsl_vector_get(var_pca, j);
		gsl_vector_set(the_variance, i, vec_mat_sum);
		vec_mat_sum = 0.0;
	}

	gsl_vector_free(mean_real);
	gsl_vector_free(temp);
	gsl_vector_free(mean_pca);
	gsl_vector_free(var_pca);
}

/**
 * repeats a lot of the code in emulate_point_multi, this is annoying
 * the_covar is an nt X nt matrix
 *
 * for facundo, 26.03.2013
 */

void emulate_point_multi_covar(multi_emulator *emu, gsl_vector *the_point,
															 gsl_vector *the_mean, gsl_vector *the_variance,
															 gsl_matrix *the_covar)
{
	int i, j, k;
	int nr = emu->nr;
	int nt = emu->nt;
	double covar_contrib = 0.0;

	double vec_mat_sum;
	// the mean and variance in the PCA space
	gsl_vector *mean_pca = gsl_vector_alloc(nr);
	gsl_vector *var_pca = gsl_vector_alloc(nr);

	// first we sample the nr emulators
	for(i = 0; i < emu->nr; i++)
		emulate_point(emu->emu_struct_array[i], the_point, gsl_vector_ptr(mean_pca, i), gsl_vector_ptr(var_pca, i));

	/**
	 *  now we have to project back into the REAL space from the PCA space
	 * mean_real (nt) = training_mean (nt) + pca_evec_r %*% diag(sqrt(pca_evals_R)) %*% mean_pca
	 */
	gsl_vector *mean_real = gsl_vector_alloc(nt);
	gsl_vector *temp = gsl_vector_alloc(nt);
	gsl_vector_memcpy(mean_real, emu->model->training_mean);

	vec_mat_sum = 0.0;
	for(i = 0; i < nt; i++){
		for(j = 0; j < nr; j++)
			vec_mat_sum += gsl_matrix_get(emu->model->pca_evecs_r, i, j) *
				sqrt(gsl_vector_get(emu->model->pca_evals_r, j))*
				gsl_vector_get(mean_pca, j);

		// save the sum scaled by the sqrt of the eval
		//gsl_vector_set(temp, i, vec_mat_sum*sqrt(gsl_vector_get(emu->model->pca_evals_r, i)));
		gsl_vector_set(temp, i, vec_mat_sum);
		vec_mat_sum = 0.0;
	}

	gsl_vector_add(mean_real, temp);
	gsl_vector_memcpy(the_mean, mean_real); // save the final mean

	/**
	 * project back the variance
	 * yVar_i = Sum_k=0^{nr} ( U_i_k * U_i_k * lambda_k * V_k)
	 */
	gsl_vector_set_zero(temp);
	vec_mat_sum = 0.0;
	for(i = 0; i < nt; i++){
		for(j = 0; j < nr; j++)
			vec_mat_sum += pow(gsl_matrix_get(emu->model->pca_evecs_r, i, j), 2.0) *
				gsl_vector_get(emu->model->pca_evals_r, j) * gsl_vector_get(var_pca, j);
		gsl_vector_set(the_variance, i, vec_mat_sum);
		vec_mat_sum = 0.0;
	}

	// reset the covariance
	gsl_matrix_set_zero(the_covar);

	/**
	 * loop over points in the nt x nt covariance matrix
	 */
	for(i = 0; i < nt; i++){
		for(j = 0; j < nt; j++){
			/**
			 * C_ij = sum_k^{nr} [ U_i_k * U_j_k * lambda_k * V_k ]
			 *
			 * sum contributions over k */
			covar_contrib = 0.0;
			for(k = 0; k < nr; k++){
				covar_contrib += gsl_matrix_get(emu->model->pca_evecs_r, i, k) *
					gsl_matrix_get(emu->model->pca_evecs_r, j, k) *
					gsl_vector_get(emu->model->pca_evals_r, k) *  gsl_vector_get(var_pca, k);
			}
			gsl_matrix_set(the_covar, i, j, covar_contrib);
		}
	}

	gsl_vector_free(mean_real);
	gsl_vector_free(temp);
	gsl_vector_free(mean_pca);
	gsl_vector_free(var_pca);

}

max_stars_repo_name: jackdawjackdaw/emulator

max_issues_repo_name: jackdawjackdaw/emulator

max_forks_repo_name: jackdawjackdaw/emulator


text:
//
// Created by karl_ on 2020-5-3.
//

#include <math.h>
#include <cblas.h>
#include <imgproc/color-convert.h>

gboolean imgproc_to_HSL(CoreImage *src, CoreImage **dst) {
    CoreSize *size;
    CoreColorSpace src_color_space;
    CorePixelType src_pixel_type;
    gdouble *src_data_cast, *dst_data;
    gdouble *hsl_pixel, *rgb_pixel, c_max, c_min, c;
    gsize i, j, area;

    size = core_image_get_size(src);
    area = core_size_get_area(size);
    src_color_space = core_image_get_color_space(src);
    src_pixel_type = core_image_get_pixel_type(src);

    g_return_val_if_fail(src_color_space == CORE_COLOR_SPACE_RGB, FALSE);

    if (core_pixel_is_double(src_pixel_type)) {
        src_data_cast = core_image_get_data(src);
    } else if (core_pixel_is_uint8(src_pixel_type)) {
        src_data_cast = malloc(sizeof(gdouble) * core_size_get_area(size));
        for (i = 0; i < core_size_get_area(size) * 3; ++i) {
            src_data_cast[i] = ((guint8 *) core_image_get_data(src))[i] / core_pixel_get_range(src_pixel_type);
        }
    } else {
        return FALSE;
    }

    dst_data = g_malloc(sizeof(gdouble) * core_size_get_area(size));
    for (i = 0; i < area; ++i) {
        hsl_pixel = dst_data + i * 3;
        rgb_pixel = src_data_cast + i * 3;
        c_min = 1;
        c_max = 0;
        c = c_max - c_min;
        for (j = 0; j < 3; ++j) {
            if (rgb_pixel[j] > c_max) {
                c_max = rgb_pixel[j];
            }
            if (rgb_pixel[j] < c_min) {
                c_max = rgb_pixel[j];
            }
        }

        /* Lightness */
        hsl_pixel[2] = (c_min + c_max) / 2;

        /* Saturation */
        if (!c) {
            hsl_pixel[1] = 0;
        } else {
            hsl_pixel[1] = c / (1 - fabs(hsl_pixel[2] * 2 - 1));
        }

        /* Hue */
        if (!c) {
            hsl_pixel[0] = 0;
        } else {
            if (rgb_pixel[0] == c_max) {
                hsl_pixel[0] = 1.0 / 6 * (rgb_pixel[1] - rgb_pixel[2]) / c;
            } else if (rgb_pixel[1] == c_max) {
                hsl_pixel[0] = 1.0 / 6 * (2 + (rgb_pixel[2] - rgb_pixel[0]) / c);
            } else if (rgb_pixel[2] == c_max) {
                hsl_pixel[0] = 1.0 / 6 * (4 + (rgb_pixel[0] - rgb_pixel[1]) / c);
            }
        }
    }

    if (*dst == NULL) {
        *dst = core_image_new_with_data(dst_data, CORE_COLOR_SPACE_HSL, CORE_PIXEL_D3, size, FALSE);
    } else {
        core_image_assign_data(*dst, dst_data, CORE_COLOR_SPACE_HSL, CORE_PIXEL_D3, size, FALSE);
    }

    if (core_pixel_is_uint8(src_pixel_type)) {
        free(src_data_cast);
    }
    g_object_unref(size);

    return TRUE;
}

gboolean imgproc_to_RGB(CoreImage *src, CoreImage **dst) {
    CoreSize *size;
    CoreColorSpace color_space;
    CorePixelType pixel_type;
    gdouble *src_data, *dst_data;
    gdouble *hsl_pixel, *rgb_pixel;
    gdouble c, h, h_ceil, x, m, r, g, b;
    gsize i, area;

    g_return_val_if_fail(src, FALSE);
    color_space = core_image_get_color_space(src);
    pixel_type = core_image_get_pixel_type(src);
    g_return_val_if_fail(color_space == CORE_COLOR_SPACE_HSL, FALSE);
    g_return_val_if_fail(pixel_type == CORE_PIXEL_D3, FALSE);
    size = core_image_get_size(src);
    area = core_size_get_area(size);
    src_data = core_image_get_data(src);
    dst_data = g_malloc(sizeof(gdouble) * area);
    for (i = 0; i < area; ++i) {
        hsl_pixel = src_data + i * 3;
        rgb_pixel = dst_data + i * 3;

        c = (1 - fabs(2 * hsl_pixel[2] - 1));
        h = hsl_pixel[0] * 6.0;
        x = c * (1 - fabs(h - floor(h / 2) * 2 - 1));
        h_ceil = ceil(h) + DBL_EPSILON;
        if (h_ceil >= 6) {
            r = c;
            g = 0;
            b = x;
        } else if (h_ceil >= 5) {
            r = x;
            g = 0;
            b = c;
        } else if (h_ceil >= 4) {
            r = 0;
            g = x;
            b = c;
        } else if (h_ceil >= 3) {
            r = 0;
            g = c;
            b = x;
        } else if (h_ceil >= 2) {
            r = x;
            g = c;
            b = 0;
        } else if (h_ceil >= 1) {
            r = c;
            g = x;
            b = 0;
        } else {
            r = g = b = 0;
        }
        m = hsl_pixel[2] - c / 2;
        rgb_pixel[0] = r + m;
        rgb_pixel[1] = g + m;
        rgb_pixel[2] = b + m;
    }

    if (*dst == NULL) {
        *dst = core_image_new_with_data(dst_data, CORE_COLOR_SPACE_RGB, CORE_PIXEL_D3, size, FALSE);
    } else {
        core_image_assign_data(*dst, dst_data, CORE_COLOR_SPACE_RGB, CORE_PIXEL_D3, size, FALSE);
    }

    g_object_unref(size);
    return TRUE;
}

gboolean imgproc_to_grayscale(CoreImage *src, CoreImage **dst) {
    gsize i;
    CoreSize *size;
    gdouble *dst_data_cast, *src_data_cast;
    gpointer src_data, dst_data;
    gsize channel, area;
    CorePixelType pixel_type;
    g_return_val_if_fail(src != NULL, FALSE);
    g_return_val_if_fail(dst != NULL, FALSE);

    size = core_image_get_size(src);
    area = core_size_get_area(size);
    pixel_type = core_image_get_pixel_type(src);
    channel = core_image_get_channel(src);
    src_data = core_image_get_data(src);

    if (core_pixel_is_uint8(pixel_type)) {
        src_data_cast = g_malloc(area * sizeof(gdouble) * channel);
        for (i = 0; i < area * channel; ++i) {
            src_data_cast[i] = (gdouble) ((guint8 *) src_data)[i];
        }
    } else {
        /* do not need to cast if pixel type is Dx */
        src_data_cast = src_data;
    }
    dst_data_cast = g_malloc(area * sizeof(gdouble));

    cblas_dcopy(area, src_data_cast, channel, dst_data_cast, 1);
    for (i = 1; i < channel; ++i) {
        cblas_daxpy(area, 1.0, src_data_cast + i, channel, dst_data_cast, 1);
    }
    if (channel > 1) {
        cblas_daxpy(area, 1.0 / channel - 1.0, dst_data_cast, 1, dst_data_cast, 1);
    }

    if (core_pixel_is_uint8(pixel_type)) {
        pixel_type = CORE_PIXEL_U1;
    } else if (core_pixel_is_double(pixel_type)) {
        pixel_type = CORE_PIXEL_D1;
    } else {
        g_return_val_if_fail(FALSE, FALSE);
    }

    /* cast dst_data from double back to uchar if needed */
    if (core_pixel_is_uint8(pixel_type)) {
        for (i = 0; i < area; ++i) {
            ((guint8 *) dst_data_cast)[i] = dst_data_cast[i];
        }
    }
    dst_data = dst_data_cast;
    if (core_pixel_is_uint8(pixel_type)) {
        dst_data = g_realloc(dst_data, sizeof(guint8) * area);
    }

    /* assign result */
    if (*dst == NULL) {
        *dst = core_image_new_with_data(dst_data, CORE_COLOR_SPACE_GRAY_SCALE, pixel_type, size, FALSE);
    } else {
        core_image_assign_data(*dst, dst_data, CORE_COLOR_SPACE_GRAY_SCALE, pixel_type, size, FALSE);
    }

    /* free additional allocated memory */
    if (core_pixel_is_uint8(pixel_type)) {
        g_free(src_data_cast);
    }
    g_object_unref(size);
    return TRUE;
}

gboolean imgproc_to_binary_threshold(CoreImage *src, CoreImage **dst, gdouble threshold, gboolean inverse) {
    CoreImage *gray_src = NULL;
    CoreSize *size;
    gsize i, area;
    gdouble *src_data, *dst_data;
    CoreColorSpace color_space;
    CorePixelType pixel_type;
    g_return_val_if_fail(src != NULL, FALSE);
    g_return_val_if_fail(dst != NULL, FALSE);

    if (core_image_get_channel(src) > 1) {
        imgproc_to_grayscale(src, &gray_src);
    } else {
        gray_src = g_object_ref(src);
    }
    size = core_image_get_size(gray_src);
    area = core_size_get_area(size);
    src_data = core_image_get_data(gray_src);
    dst_data = g_malloc(sizeof(gdouble) * area);
    for (i = 0; i < area; ++i) {
        dst_data[i] = ((unsigned) (src_data[i] > threshold) ^ (unsigned) inverse) * 255.0;
    }
    if (*dst == NULL) {
        core_image_new_with_data(dst_data, CORE_COLOR_SPACE_BIN, CORE_PIXEL_U1, size, FALSE);
    } else {
        core_image_assign_data(*dst, dst_data, CORE_COLOR_SPACE_BIN, CORE_PIXEL_U1, size, FALSE);
    }
    g_object_unref(size);
    g_object_unref(gray_src);
    return TRUE;
}

max_stars_repo_name: zhengbanbufengbi/EasyPhotoshop

max_issues_repo_name: zhengbanbufengbi/EasyPhotoshop

max_forks_repo_name: zhengbanbufengbi/EasyPhotoshop


text:
/*
 *  Copyright 2014 NVIDIA Corporation
 *
 *  Licensed under the Apache License, Version 2.0 (the "License");
 *  you may not use this file except in compliance with the License.
 *  You may obtain a copy of the License at
 *
 *      http://www.apache.org/licenses/LICENSE-2.0
 *
 *  Unless required by applicable law or agreed to in writing, software
 *  distributed under the License is distributed on an "AS IS" BASIS,
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *  See the License for the specific language governing permissions and
 *  limitations under the License.
 */

/* include the ATLAS headers */

extern "C"
{
#include <cblas.h>
}

/* choose precision to train and classify.  Only float and double are
 * currently suppored
 */

typedef float floatType_t;

/* macro to convert 2d coords to 1d offset */

#define INDX(row,col,ld) (((col) * (ld)) + (row))

/* macros for max/min to combine with argmin */

#define MYMAX(val,array,i,index) \
if( array[i] > val ) \
{ \
  val = array[i]; \
  index = i; \
} \

#define MYMIN(val,array,i,index) \
if( array[i] < val ) \
{ \
  val = array[i]; \
  index = i; \
} \

/* macro to clip values from min to max */

#define CLIP(val,min,max) \
if( (val) < (min) ) val = (min); \
else if( (val) > (max) ) val = (max);

/* hardcoded constants for training and test set size and feature
 * vector size
 */

#define FEATURE_VECTOR_SIZE (1899)
#define TRAINING_SET_SIZE (4000)
#define TEST_SET_SIZE (1000)

/* CUDA debugging */

#ifdef DEBUG
#define CUDA_CALL(F)  if( (F) != cudaSuccess ) \
  {printf("Error %s at %s:%d\n", cudaGetErrorString(cudaGetLastError()), \
   __FILE__,__LINE__); exit(-1);}
#define CUDA_CHECK()  if( (cudaPeekAtLastError()) != cudaSuccess ) \
  {printf("Error %s at %s:%d\n", cudaGetErrorString(cudaGetLastError()), \
   __FILE__,__LINE__-1); exit(-1);}
#else
#define CUDA_CALL(F) (F)
#define CUDA_CHECK()
#endif

/* function defs */

void readMatrixFromFile( char *, int *, const int, const int );

void calculateBI( floatType_t const *,
                  floatType_t const *,
                  floatType_t const *,
                  int ,
                  floatType_t *, floatType_t *,
                  int *, int *,
                  floatType_t const );

void svmTrain( floatType_t const *, floatType_t const *, floatType_t const,
               const int, const int,
               const floatType_t, floatType_t * );

void svmPredict( floatType_t const *, floatType_t const *,
                 int const, int const, int * );


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

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_chebyshev.h>

/*-*-*-*-*-*-*-*-*-*-*-* Allocators *-*-*-*-*-*-*-*-*-*-*-*/

gsl_cheb_series *
gsl_cheb_alloc(const size_t order)
{
  gsl_cheb_series * cs = (gsl_cheb_series *) malloc(sizeof(gsl_cheb_series));

  if(cs == 0) {
    GSL_ERROR_VAL("failed to allocate gsl_cheb_series struct", GSL_ENOMEM, 0);
  }

  cs->order    = order;
  cs->order_sp = order;

  cs->c = (double *) malloc((order+1) * sizeof(double));

  if(cs->c == 0) {
    GSL_ERROR_VAL("failed to allocate cheb coefficients", GSL_ENOMEM, 0);
  }

  cs->f = (double *) malloc((order+1) * sizeof(double));

  if(cs->f == 0) {
    GSL_ERROR_VAL("failed to allocate cheb function space", GSL_ENOMEM, 0);
  }

  return cs;
}

void gsl_cheb_free(gsl_cheb_series * cs)
{
  RETURN_IF_NULL (cs);
  free(cs->f);
  free(cs->c);
  free(cs);
}

/*-*-*-*-*-*-*-*-*-*-*-* Initializer *-*-*-*-*-*-*-*-*-*-*-*/

int gsl_cheb_init(gsl_cheb_series * cs, const gsl_function *func,
                  const double a, const double b)
{
  size_t k, j;

  if(a >= b) {
    GSL_ERROR_VAL("null function interval [a,b]", GSL_EDOM, 0);
  }
  cs->a = a;
  cs->b = b;
  /* cs->err = 0.0; */

  {
    double bma = 0.5 * (cs->b - cs->a);
    double bpa = 0.5 * (cs->b + cs->a);
    double fac = 2.0/(cs->order +1.0);

    for(k = 0; k<=cs->order; k++) {
      double y = cos(M_PI * (k+0.5)/(cs->order+1));
      cs->f[k] = GSL_FN_EVAL(func, (y*bma + bpa));
    }

    for(j = 0; j<=cs->order; j++) {
      double sum = 0.0;
      for(k = 0; k<=cs->order; k++)
        sum += cs->f[k]*cos(M_PI * j*(k+0.5)/(cs->order+1));
      cs->c[j] = fac * sum;
    }

  }
  return GSL_SUCCESS;
}

size_t
gsl_cheb_order (const gsl_cheb_series * cs)
{
  return cs->order;
}

size_t
gsl_cheb_size (const gsl_cheb_series * cs)
{
  return (cs->order + 1);
}

double *
gsl_cheb_coeffs (const gsl_cheb_series * cs)
{
  return cs->c;
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#include <stdio.h>
#include <stdarg.h>
#include <string.h>
#include <math.h>
#include <gbpLib.h>
#include <gbpRNG.h>
#include <gbpMCMC.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_fit.h>
#include <gsl/gsl_interp.h>

void read_MCMC_covariance(MCMC_info *MCMC, char *filename) {
    FILE *  fp;
    double *V;
    int     n_P, i_P, j_P;
    if(MCMC->V == NULL)
        SID_log("Initializing the covariance matrix from file {%s}...", SID_LOG_OPEN, filename);
    else
        SID_log("Updating the covariance matrix from file {%s}...", SID_LOG_OPEN, filename);
    fp = fopen(filename, "r");
    SID_fread_verify(&n_P, sizeof(int), 1, fp);
    V = (double *)SID_malloc(sizeof(double) * n_P * n_P);
    SID_fread_verify(V, sizeof(double), n_P * n_P, fp);
    fclose(fp);
    set_MCMC_covariance(MCMC, V);
    SID_free(SID_FARG V);
    SID_log("Done.", SID_LOG_CLOSE);
}

max_stars_repo_name: gbpoole/gbpCode

max_issues_repo_name: gbpoole/gbpCode

max_forks_repo_name: gbpoole/gbpCode


text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define COMPEARTH_PRIVATE_GEM3 1
#include "compearth.h"
#ifdef DEBUG
#ifdef COMPEARTH_USE_MKL
#ifdef __clang__
#pragma clang diagnostic push
#pragma clang diagnostic ignored "-Wreserved-id-macro"
#pragma clang diagnostic ignored "-Wstrict-prototypes"
#endif
#include <mkl_cblas.h>
#ifdef __clang__
#pragma clang diagnostic pop
#endif
#else
#include <cblas.h>
#endif
#endif

/*!
 * @brief Compute a rotation matrix given an axis and angle
 *
 * @param[in] n     Number of matrices to generate.
 * @param[in] v     Rotation axis.  This is an array of dimension [3 x n]
 *                  with leading dimension 3.
 * @param[in] xi    Rotation angles (degrees).  This is an array of
 *                  dimension [n].
 *
 * @param[in] U     Rotation matrices.  This is an array of dimension
 *                  [3 x 3 x n] where each leading [3 x 3] matrix is
 *                  in column major format.
 *
 * @result 0 indicates success.
 *
 * @date 2016 - Ben Baker translated Carl Tape's rotmat_gen.m to C
 *
 * @copyright MIT
 *
 */
int compearth_eulerUtil_rotmatGen(const int n,
                                  const double *__restrict__ v,
                                  const double *__restrict__ xi,
                                  double *__restrict__ U)
{
    double R1[9], R2[9], R3[9], R4[9], R5[9], R54[9], R543[9], R5432[9];
    double ele, nvph_deg, nvth_deg, rho, vph, vth, vph_deg, vth_deg;
    int i, ierr;
    const double pi180i = 180.0/M_PI;
    ierr = 0;
    for (i=0; i<n; i++)
    {
        compearth_matlab_cart2sph(1, &v[3*i+0], &v[3*i+1], &v[3*i+2],
                                  &vph, &ele, &rho);
        vth = M_PI_2 - ele;
        vth_deg = vth*pi180i;
        vph_deg = vph*pi180i;
        nvph_deg =-vph_deg;
        nvth_deg =-vth_deg;
        ierr = ierr + compearth_eulerUtil_rotmat(1, &nvph_deg, 3, R1);
        ierr = ierr + compearth_eulerUtil_rotmat(1, &nvth_deg, 2, R2);
        ierr = ierr + compearth_eulerUtil_rotmat(1, &xi[i],    3, R3);
        ierr = ierr + compearth_eulerUtil_rotmat(1, &vth_deg,  2, R4);
        ierr = ierr + compearth_eulerUtil_rotmat(1, &vph_deg,  3, R5);
        if (ierr != 0)
        {
            fprintf(stderr, "%s: Error computing rotation matrices %d\n",
                    __func__, i);
        }
        gemm3_colMajorNoTransNoTrans(R5, R4, R54);
        gemm3_colMajorNoTransNoTrans(R54, R3, R543);
        gemm3_colMajorNoTransNoTrans(R543, R2, R5432);
        gemm3_colMajorNoTransNoTrans(R5432, R1, &U[9*i]);
#ifdef DEBUG
        double U9[9];
        cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans,
                    3, 3, 3, 1.0, R5, 3, R4, 3, 0.0, R54, 3);
        cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans,
                    3, 3, 3, 1.0, R54, 3, R3, 3, 0.0, R543, 3);
        cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans,
                    3, 3, 3, 1.0, R543, 3, R2, 3, 0.0, R5432, 3);
        cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans,
                    3, 3, 3, 1.0, R5432, 3, R1, 3, 0.0, U9, 3);
        for (int k=0; k<9; k++)
        {
            if (fabs(U9[k] - U[9*i+k]) > 1.e-14)
            {
                printf("Computation failed: %f %f\n", U9[k], U[9*i+k]);
                ierr = ierr + 1;
            }
        }
#endif
    }
    return ierr;
}

max_stars_repo_name: OUCyf/mtbeach

max_issues_repo_name: carltape/mtbeach

max_forks_repo_name: carltape/mtbeach


text:
#pragma once

#include "CesiumGltf/MetadataArrayView.h"
#include "CesiumGltf/PropertyType.h"
#include "CesiumGltf/PropertyTypeTraits.h"

#include <gsl/span>

#include <cassert>
#include <cstddef>
#include <cstdint>
#include <string_view>
#include <type_traits>

namespace CesiumGltf {
/**
 * @brief Indicates the status of a property view.
 *
 * The {@link MetadataPropertyView} constructor always completes successfully. However,
 * it may not always reflect the actual content of the {@link FeatureTableProperty}, but
 * instead indicate that its {@link MetadataPropertyView::size} is 0. This enumeration
 * provides the reason.
 */
enum class MetadataPropertyViewStatus {
  /**
   * @brief This property view is valid and ready to use.
   */
  Valid,

  /**
   * @brief This property view does not exist in the FeatureTable.
   */
  InvalidPropertyNotExist,

  /**
   * @brief This property view does not have a correct type with what is
   * specified in {@link ClassProperty::type}.
   */
  InvalidTypeMismatch,

  /**
   * @brief This property view does not have a valid value buffer view index.
   */
  InvalidValueBufferViewIndex,

  /**
   * @brief This array property view does not have a valid array offset buffer
   * view index.
   */
  InvalidArrayOffsetBufferViewIndex,

  /**
   * @brief This string property view does not have a valid string offset buffer
   * view index.
   */
  InvalidStringOffsetBufferViewIndex,

  /**
   * @brief This property view has a valid value buffer view index, but buffer
   * view specifies an invalid buffer index
   */
  InvalidValueBufferIndex,

  /**
   * @brief This property view has a valid array string buffer view index, but
   * buffer view specifies an invalid buffer index
   */
  InvalidArrayOffsetBufferIndex,

  /**
   * @brief This property view has a valid string offset buffer view index, but
   * buffer view specifies an invalid buffer index
   */
  InvalidStringOffsetBufferIndex,

  /**
   * @brief This property view has buffer view's offset not aligned by 8 bytes
   */
  InvalidBufferViewNotAligned8Bytes,

  /**
   * @brief This property view has an out-of-bound buffer view
   */
  InvalidBufferViewOutOfBound,

  /**
   * @brief This property view has an invalid buffer view's length which is not
   * a multiple of the size of its type or offset type
   */
  InvalidBufferViewSizeNotDivisibleByTypeSize,

  /**
   * @brief This property view has an invalid buffer view's length which cannot
   * fit all the instances of the feature table
   */
  InvalidBufferViewSizeNotFitInstanceCount,

  /**
   * @brief This array property view has both component count and offset buffer
   * view
   */
  InvalidArrayComponentCountAndOffsetBufferCoexist,

  /**
   * @brief This array property view doesn't have either component count or
   * offset buffer view
   */
  InvalidArrayComponentCountOrOffsetBufferNotExist,

  /**
   * @brief This property view have an unknown offset type
   */
  InvalidOffsetType,

  /**
   * @brief This property view has offset values not sorted ascendingly
   */
  InvalidOffsetValuesNotSortedAscending,

  /**
   * @brief This property view has an offset point to an out of bound value
   */
  InvalidOffsetValuePointsToOutOfBoundBuffer
};

/**
 * @brief A view on the data of the FeatureTableProperty
 *
 * It provides utility to retrieve the actual data stored in the
 * {@link FeatureTableProperty::bufferView} like an array of elements.
 * Data of each instance can be accessed through the {@link get(int64_t instance)} method
 *
 * @param ElementType must be uin8_t, int8_t, uint16_t, int16_t,
 * uint32_t, int32_t, uint64_t, int64_t, float, double, bool, std::string_view,
 * and MetadataArrayView<T> with T must be one of the types mentioned above
 */
template <typename ElementType> class MetadataPropertyView {
public:
  /**
   * @brief Constructs a new instance viewing a non-existent property.
   */
  MetadataPropertyView()
      : _status{MetadataPropertyViewStatus::InvalidPropertyNotExist},
        _valueBuffer{},
        _arrayOffsetBuffer{},
        _stringOffsetBuffer{},
        _offsetType{},
        _offsetSize{},
        _componentCount{},
        _instanceCount{},
        _normalized{} {}

  /**
   * @brief Construct a new instance pointing to the data specified by
   * FeatureTableProperty
   * @param valueBuffer The raw buffer specified by {@link FeatureTableProperty::bufferView}
   * @param arrayOffsetBuffer The raw buffer specified by {@link FeatureTableProperty::arrayOffsetBufferView}
   * @param stringOffsetBuffer The raw buffer specified by {@link FeatureTableProperty::stringOffsetBufferView}
   * @param offsetType The offset type of the arrayOffsetBuffer and stringOffsetBuffer that is specified by {@link FeatureTableProperty::offsetType}
   * @param componentCount The number of elements for fixed array value which is specified by {@link FeatureTableProperty::componentCount}
   * @param instanceCount The number of instances specified by {@link FeatureTable::count}
   * @param normalized Whether this property has a normalized integer type.
   */
  MetadataPropertyView(
      MetadataPropertyViewStatus status,
      gsl::span<const std::byte> valueBuffer,
      gsl::span<const std::byte> arrayOffsetBuffer,
      gsl::span<const std::byte> stringOffsetBuffer,
      PropertyType offsetType,
      int64_t componentCount,
      int64_t instanceCount,
      bool normalized) noexcept
      : _status{status},
        _valueBuffer{valueBuffer},
        _arrayOffsetBuffer{arrayOffsetBuffer},
        _stringOffsetBuffer{stringOffsetBuffer},
        _offsetType{offsetType},
        _offsetSize{getOffsetSize(offsetType)},
        _componentCount{componentCount},
        _instanceCount{instanceCount},
        _normalized{normalized} {}

  /**
   * @brief Gets the status of this property view.
   *
   * Indicates whether the view accurately reflects the property's data, or
   * whether an error occurred.
   */
  MetadataPropertyViewStatus status() const noexcept { return _status; }

  /**
   * @brief Get the value of an instance of the FeatureTable.
   * @param instance The instance index
   * @return The value of the instance
   */
  ElementType get(int64_t instance) const noexcept {
    assert(
        _status == MetadataPropertyViewStatus::Valid &&
        "Check the status() first to make sure view is valid");
    assert(
        size() > 0 &&
        "Check the size() of the view to make sure it's not empty");
    assert(instance >= 0 && "instance index must be positive");

    if constexpr (IsMetadataNumeric<ElementType>::value) {
      return getNumeric(instance);
    }

    if constexpr (IsMetadataBoolean<ElementType>::value) {
      return getBoolean(instance);
    }

    if constexpr (IsMetadataString<ElementType>::value) {
      return getString(instance);
    }

    if constexpr (IsMetadataNumericArray<ElementType>::value) {
      return getNumericArray<typename MetadataArrayType<ElementType>::type>(
          instance);
    }

    if constexpr (IsMetadataBooleanArray<ElementType>::value) {
      return getBooleanArray(instance);
    }

    if constexpr (IsMetadataStringArray<ElementType>::value) {
      return getStringArray(instance);
    }
  }

  /**
   * @brief Get the number of instances in the FeatureTable
   * @return The number of instances in the FeatureTable
   */
  int64_t size() const noexcept { return _instanceCount; }

  /**
   * @brief Get the component count of this property. Only applicable when the
   * property is an array type.
   *
   * @return The component count of this property.
   */
  int64_t getComponentCount() const noexcept { return _componentCount; }

  /**
   * @brief Whether this property has a normalized integer type.
   *
   * @return Whether this property has a normalized integer type.
   */
  bool isNormalized() const noexcept { return _normalized; }

private:
  ElementType getNumeric(int64_t instance) const noexcept {
    return reinterpret_cast<const ElementType*>(_valueBuffer.data())[instance];
  }

  bool getBoolean(int64_t instance) const noexcept {
    const int64_t byteIndex = instance / 8;
    const int64_t bitIndex = instance % 8;
    const int bitValue =
        static_cast<int>(_valueBuffer[byteIndex] >> bitIndex) & 1;
    return bitValue == 1;
  }

  std::string_view getString(int64_t instance) const noexcept {
    const size_t currentOffset =
        getOffsetFromOffsetBuffer(instance, _stringOffsetBuffer, _offsetType);
    const size_t nextOffset = getOffsetFromOffsetBuffer(
        instance + 1,
        _stringOffsetBuffer,
        _offsetType);
    return std::string_view(
        reinterpret_cast<const char*>(_valueBuffer.data() + currentOffset),
        (nextOffset - currentOffset));
  }

  template <typename T>
  MetadataArrayView<T> getNumericArray(int64_t instance) const noexcept {
    if (_componentCount > 0) {
      const gsl::span<const std::byte> vals(
          _valueBuffer.data() + instance * _componentCount * sizeof(T),
          _componentCount * sizeof(T));
      return MetadataArrayView<T>{vals};
    }

    const size_t currentOffset =
        getOffsetFromOffsetBuffer(instance, _arrayOffsetBuffer, _offsetType);
    const size_t nextOffset = getOffsetFromOffsetBuffer(
        instance + 1,
        _arrayOffsetBuffer,
        _offsetType);
    const gsl::span<const std::byte> vals(
        _valueBuffer.data() + currentOffset,
        (nextOffset - currentOffset));
    return MetadataArrayView<T>{vals};
  }

  MetadataArrayView<std::string_view>
  getStringArray(int64_t instance) const noexcept {
    if (_componentCount > 0) {
      const gsl::span<const std::byte> offsetVals(
          _stringOffsetBuffer.data() + instance * _componentCount * _offsetSize,
          (_componentCount + 1) * _offsetSize);
      return MetadataArrayView<std::string_view>(
          _valueBuffer,
          offsetVals,
          _offsetType,
          _componentCount);
    }

    const size_t currentOffset =
        getOffsetFromOffsetBuffer(instance, _arrayOffsetBuffer, _offsetType);
    const size_t nextOffset = getOffsetFromOffsetBuffer(
        instance + 1,
        _arrayOffsetBuffer,
        _offsetType);
    const gsl::span<const std::byte> offsetVals(
        _stringOffsetBuffer.data() + currentOffset,
        (nextOffset - currentOffset + _offsetSize));
    return MetadataArrayView<std::string_view>(
        _valueBuffer,
        offsetVals,
        _offsetType,
        (nextOffset - currentOffset) / _offsetSize);
  }

  MetadataArrayView<bool> getBooleanArray(int64_t instance) const noexcept {
    if (_componentCount > 0) {
      const size_t offsetBits = _componentCount * instance;
      const size_t nextOffsetBits = _componentCount * (instance + 1);
      const gsl::span<const std::byte> buffer(
          _valueBuffer.data() + offsetBits / 8,
          (nextOffsetBits / 8 - offsetBits / 8 + 1));
      return MetadataArrayView<bool>(buffer, offsetBits % 8, _componentCount);
    }

    const size_t currentOffset =
        getOffsetFromOffsetBuffer(instance, _arrayOffsetBuffer, _offsetType);
    const size_t nextOffset = getOffsetFromOffsetBuffer(
        instance + 1,
        _arrayOffsetBuffer,
        _offsetType);

    const size_t totalBits = nextOffset - currentOffset;
    const gsl::span<const std::byte> buffer(
        _valueBuffer.data() + currentOffset / 8,
        (nextOffset / 8 - currentOffset / 8 + 1));
    return MetadataArrayView<bool>(buffer, currentOffset % 8, totalBits);
  }

  static int64_t getOffsetSize(PropertyType offsetType) noexcept {
    switch (offsetType) {
    case CesiumGltf::PropertyType::Uint8:
      return sizeof(uint8_t);
    case CesiumGltf::PropertyType::Uint16:
      return sizeof(uint16_t);
    case CesiumGltf::PropertyType::Uint32:
      return sizeof(uint32_t);
    case CesiumGltf::PropertyType::Uint64:
      return sizeof(uint64_t);
    default:
      return 0;
    }
  }

  static size_t getOffsetFromOffsetBuffer(
      size_t instance,
      const gsl::span<const std::byte>& offsetBuffer,
      PropertyType offsetType) noexcept {
    switch (offsetType) {
    case PropertyType::Uint8: {
      assert(instance < offsetBuffer.size() / sizeof(uint8_t));
      const uint8_t offset = *reinterpret_cast<const uint8_t*>(
          offsetBuffer.data() + instance * sizeof(uint8_t));
      return static_cast<size_t>(offset);
    }
    case PropertyType::Uint16: {
      assert(instance < offsetBuffer.size() / sizeof(uint16_t));
      const uint16_t offset = *reinterpret_cast<const uint16_t*>(
          offsetBuffer.data() + instance * sizeof(uint16_t));
      return static_cast<size_t>(offset);
    }
    case PropertyType::Uint32: {
      assert(instance < offsetBuffer.size() / sizeof(uint32_t));
      const uint32_t offset = *reinterpret_cast<const uint32_t*>(
          offsetBuffer.data() + instance * sizeof(uint32_t));
      return static_cast<size_t>(offset);
    }
    case PropertyType::Uint64: {
      assert(instance < offsetBuffer.size() / sizeof(uint64_t));
      const uint64_t offset = *reinterpret_cast<const uint64_t*>(
          offsetBuffer.data() + instance * sizeof(uint64_t));
      return static_cast<size_t>(offset);
    }
    default:
      assert(false && "Offset type has unknown type");
      return 0;
    }
  }

  MetadataPropertyViewStatus _status;
  gsl::span<const std::byte> _valueBuffer;
  gsl::span<const std::byte> _arrayOffsetBuffer;
  gsl::span<const std::byte> _stringOffsetBuffer;
  PropertyType _offsetType;
  int64_t _offsetSize;
  int64_t _componentCount;
  int64_t _instanceCount;
  bool _normalized;
};
} // namespace CesiumGltf

max_stars_repo_name: yieryi/cesium-native

max_issues_repo_name: yieryi/cesium-native

max_forks_repo_name: yieryi/cesium-native


text:
/* Sound display/edit/etc
 *
 * originally intended as a re-implementation of my much-missed dpysnd -- the Foonly/SAIL/E/Mus10/Grnlib sound editor from ca 1983.
 */

#include "snd.h"

snd_state *ss = NULL;

static bool ignore_mus_error(int type, char *msg)
{
  XEN result = XEN_FALSE;

  if (XEN_HOOKED(ss->mus_error_hook))
    result = run_or_hook(ss->mus_error_hook,
			 XEN_LIST_2(C_TO_XEN_INT(type),
				    C_TO_XEN_STRING(msg)),
			 S_mus_error_hook);
  return(XEN_NOT_FALSE_P(result));
}

#if HAVE_SETJMP_H
  void top_level_catch(int ignore);
#endif

void mus_error_to_snd(int type, char *msg)
{
  if (!ss)
    {
      fprintf(stderr, "%s", msg);
      return;
    }

  if (!(ignore_mus_error(type, msg)))
    {
#if HAVE_EXTENSION_LANGUAGE
      if (msg == NULL)
	XEN_ERROR(XEN_ERROR_TYPE("mus-error"),
		  XEN_LIST_1(C_TO_XEN_STRING((char *)mus_error_type_to_string(type))));
      else XEN_ERROR(XEN_ERROR_TYPE("mus-error"),
		     XEN_LIST_1(C_TO_XEN_STRING(msg)));
#endif
      snd_error("%s: %s", mus_error_type_to_string(type), msg);
#if HAVE_SETJMP_H
      ss->jump_ok = true;
      top_level_catch(1); /* sigh -- try to keep going */
#endif
    }
}

static void mus_print_to_snd(char *msg)
{
  if (!ss)
    {
      fprintf(stderr, "%s", msg);
      return;
    }
  if (!(ignore_mus_error(MUS_NO_ERROR, msg)))
    if (msg)
      {
	int i, len;
	listener_append(";");
	len = strlen(msg);

	for (i = 1; i < len - 1; i++)
	  if ((msg[i] == '\n') && (msg[i + 1] == ' '))
	    msg[i + 1] = ';';

	if (msg[0] == '\n')
	  listener_append((char *)(msg + 1));
	else listener_append(msg);

	if (msg[strlen(msg) - 1] != '\n')
	  listener_append("\n");
      }
}

static void initialize_load_path(void)
{
  /* look for SND_PATH env var, add dirs to %load-path or load_path */
  char *path;
  path = getenv("SND_PATH");
  if (path)
    {
      /* colon-separated list of directory names, pushed on load-path in reverse order (hopefully = search order) */
      int i, len, dirs = 1, curdir = 0, start = 0;
      char **dirnames;

      len = strlen(path);
      for (i = 0; i < len; i++)
	if (path[i] == ':')
	  dirs++;

      dirnames = (char **)calloc(dirs, sizeof(char *));
      for (i = 0; i < len; i++)
	{
	  if ((path[i] == ':') ||
	      (i == len - 1))
	    {
	      if (i > start)
		{
		  int j, lim;
		  char *tmp;

		  if (i == (len - 1))
		    lim = i + 1;
		  else lim = i;

		  tmp = (char *)calloc(lim - start + 1, sizeof(char));
		  for (j = start; j < lim; j++)
		    tmp[j - start] = path[j];

		  dirnames[curdir++] = mus_expand_filename(tmp);
		  start = i + 1;
		  free(tmp);
		}
	    }
	}

      for (i = curdir - 1; i >= 0; i--)
	{
	  XEN_ADD_TO_LOAD_PATH(dirnames[i]);
	  free(dirnames[i]);
	}
      free(dirnames);
    }
}

void snd_set_global_defaults(bool need_cleanup)
{
  if (need_cleanup)
    {
      if (ss->HTML_Program) {free(ss->HTML_Program); ss->HTML_Program = NULL;}
      if (ss->HTML_Dir) {free(ss->HTML_Dir); ss->HTML_Dir = NULL;}
      if (ss->Temp_Dir) {free(ss->Temp_Dir); ss->Temp_Dir = NULL;}
      if (ss->Save_Dir) {free(ss->Save_Dir); ss->Save_Dir = NULL;}
      if (ss->Ladspa_Dir) {free(ss->Ladspa_Dir); ss->Ladspa_Dir = NULL;}
      if (ss->Save_State_File) {free(ss->Save_State_File); ss->Save_State_File = NULL;}
      if (ss->Eps_File) {free(ss->Eps_File); ss->Eps_File = NULL;}
      if (ss->Listener_Prompt) {free(ss->Listener_Prompt); ss->Listener_Prompt = NULL;}
      if (ss->Open_File_Dialog_Directory) {free(ss->Open_File_Dialog_Directory); ss->Open_File_Dialog_Directory = NULL;}

      /* not sure about the next two... */
      if ((cursor_style(ss) == CURSOR_PROC) && (XEN_PROCEDURE_P(ss->cursor_proc)))
	snd_unprotect_at(ss->cursor_proc_loc);
      if ((zoom_focus_style(ss) == ZOOM_FOCUS_PROC) && (XEN_PROCEDURE_P(ss->zoom_focus_proc)))
	snd_unprotect_at(ss->zoom_focus_proc_loc);
    }

  ss->Transform_Size =              DEFAULT_TRANSFORM_SIZE;
  ss->Fft_Window =                  DEFAULT_FFT_WINDOW;
  ss->Fft_Window_Alpha =            DEFAULT_FFT_WINDOW_ALPHA;
  ss->Fft_Window_Beta =             DEFAULT_FFT_WINDOW_BETA;
  ss->Transform_Graph_Type =        DEFAULT_TRANSFORM_GRAPH_TYPE;
  ss->Sinc_Width =                  DEFAULT_SINC_WIDTH;
  ss->Zero_Pad =                    DEFAULT_ZERO_PAD;
  ss->Wavelet_Type =                DEFAULT_WAVELET_TYPE;
  ss->Transform_Type =              DEFAULT_TRANSFORM_TYPE;
  ss->Transform_Normalization =     DEFAULT_TRANSFORM_NORMALIZATION;
  ss->Show_Transform_Peaks =        DEFAULT_SHOW_TRANSFORM_PEAKS;
  ss->Show_Sonogram_Cursor =        DEFAULT_SHOW_SONOGRAM_CURSOR;
  ss->Fft_Log_Magnitude =           DEFAULT_FFT_LOG_MAGNITUDE;
  ss->Fft_Log_Frequency =           DEFAULT_FFT_LOG_FREQUENCY;
  ss->Fft_With_Phases =             DEFAULT_FFT_WITH_PHASES;
  ss->Max_Transform_Peaks =         DEFAULT_MAX_TRANSFORM_PEAKS;
  ss->Log_Freq_Start =              DEFAULT_LOG_FREQ_START;
  ss->Min_dB =                      DEFAULT_MIN_DB;
  ss->lin_dB =                      pow(10.0, DEFAULT_MIN_DB * 0.05);
  ss->Show_Selection_Transform =    DEFAULT_SHOW_SELECTION_TRANSFORM;
  ss->Default_Output_Chans =        DEFAULT_OUTPUT_CHANS;
  ss->Default_Output_Srate =        DEFAULT_OUTPUT_SRATE;
  ss->Default_Output_Header_Type =  DEFAULT_OUTPUT_HEADER_TYPE;
  ss->Default_Output_Data_Format =  DEFAULT_OUTPUT_DATA_FORMAT;
  ss->Audio_Input_Device =          DEFAULT_AUDIO_INPUT_DEVICE;
  ss->Audio_Output_Device =         DEFAULT_AUDIO_OUTPUT_DEVICE;
  ss->Dac_Size =                    DEFAULT_DAC_SIZE;
  ss->Dac_Combines_Channels =       DEFAULT_DAC_COMBINES_CHANNELS;
  ss->Auto_Resize =                 DEFAULT_AUTO_RESIZE;
  ss->Auto_Update =                 DEFAULT_AUTO_UPDATE;
  ss->Auto_Update_Interval =        DEFAULT_AUTO_UPDATE_INTERVAL;
  ss->Ask_Before_Overwrite =        DEFAULT_ASK_BEFORE_OVERWRITE;
  ss->With_Toolbar =                DEFAULT_WITH_TOOLBAR;
  ss->With_Tooltips =               DEFAULT_WITH_TOOLTIPS;
  ss->Remember_Sound_State =        DEFAULT_REMEMBER_SOUND_STATE;
  ss->Ask_About_Unsaved_Edits =     DEFAULT_ASK_ABOUT_UNSAVED_EDITS;
  ss->Save_As_Dialog_Src =          DEFAULT_SAVE_AS_DIALOG_SRC;
  ss->Save_As_Dialog_Auto_Comment = DEFAULT_SAVE_AS_DIALOG_AUTO_COMMENT;
  ss->Show_Full_Duration =          DEFAULT_SHOW_FULL_DURATION;
  ss->Show_Full_Range =             DEFAULT_SHOW_FULL_RANGE;
  ss->Initial_Beg =                 DEFAULT_INITIAL_BEG;
  ss->Initial_Dur =                 DEFAULT_INITIAL_DUR;
  ss->With_Background_Processes =   DEFAULT_WITH_BACKGROUND_PROCESSES;
  ss->With_File_Monitor =           DEFAULT_WITH_FILE_MONITOR;
  ss->Selection_Creates_Region =    DEFAULT_SELECTION_CREATES_REGION;
  ss->Channel_Style =               DEFAULT_CHANNEL_STYLE;
  ss->Sound_Style =                 DEFAULT_SOUND_STYLE;
  ss->Graphs_Horizontal =           DEFAULT_GRAPHS_HORIZONTAL;
  ss->Graph_Style =                 DEFAULT_GRAPH_STYLE;
  ss->Region_Graph_Style =          DEFAULT_GRAPH_STYLE;
  ss->Time_Graph_Type =             DEFAULT_TIME_GRAPH_TYPE;
  ss->X_Axis_Style =                DEFAULT_X_AXIS_STYLE;
  ss->Beats_Per_Minute =            DEFAULT_BEATS_PER_MINUTE;
  ss->Beats_Per_Measure =           DEFAULT_BEATS_PER_MEASURE;
  ss->With_Relative_Panes =         DEFAULT_WITH_RELATIVE_PANES;
  ss->With_GL =                     DEFAULT_WITH_GL;
  ss->Dot_Size =                    DEFAULT_DOT_SIZE;
  ss->Grid_Density =                DEFAULT_GRID_DENSITY;
  ss->Zoom_Focus_Style =            DEFAULT_ZOOM_FOCUS_STYLE;
  ss->zoom_focus_proc =             XEN_UNDEFINED;
  ss->zoom_focus_proc_loc =         NOT_A_GC_LOC;
  ss->Max_Regions =                 DEFAULT_MAX_REGIONS;
  ss->Show_Y_Zero =                 DEFAULT_SHOW_Y_ZERO;
  ss->Show_Grid =                   DEFAULT_SHOW_GRID;
  ss->Show_Axes =                   DEFAULT_SHOW_AXES;
  ss->Show_Indices =                DEFAULT_SHOW_INDICES;
  ss->Show_Backtrace =              DEFAULT_SHOW_BACKTRACE;
  ss->With_Inset_Graph =            DEFAULT_WITH_INSET_GRAPH;
  ss->With_Interrupts =             DEFAULT_WITH_INTERRUPTS;
  ss->With_Menu_Icons =             DEFAULT_WITH_MENU_ICONS;
  ss->With_Smpte_Label =            DEFAULT_WITH_SMPTE_LABEL;
  ss->With_Pointer_Focus =          DEFAULT_WITH_POINTER_FOCUS;
  ss->Play_Arrow_Size =             DEFAULT_PLAY_ARROW_SIZE;
  ss->Sync_Style =                  DEFAULT_SYNC_STYLE;
  ss->Listener_Prompt =             mus_strdup(DEFAULT_LISTENER_PROMPT);
  ss->listener_prompt_length =      mus_strlen(ss->Listener_Prompt);
  ss->Minibuffer_History_Length =   DEFAULT_MINIBUFFER_HISTORY_LENGTH;
  ss->Clipping =                    DEFAULT_CLIPPING;
  ss->Optimization =                DEFAULT_OPTIMIZATION;
  ss->Print_Length =                DEFAULT_PRINT_LENGTH;
  ss->View_Files_Sort =             DEFAULT_VIEW_FILES_SORT;
  ss->Just_Sounds =                 DEFAULT_JUST_SOUNDS;
  ss->Open_File_Dialog_Directory =  NULL;
  ss->HTML_Dir =                    mus_strdup(DEFAULT_HTML_DIR);
  ss->HTML_Program =                mus_strdup(DEFAULT_HTML_PROGRAM);
  ss->Cursor_Size =                 DEFAULT_CURSOR_SIZE;
  ss->Cursor_Style =                DEFAULT_CURSOR_STYLE;
  ss->Tracking_Cursor_Style =       DEFAULT_TRACKING_CURSOR_STYLE;
  ss->With_Tracking_Cursor =        DEFAULT_WITH_TRACKING_CURSOR;
  ss->cursor_proc =                 XEN_UNDEFINED;
  ss->cursor_proc_loc =             NOT_A_GC_LOC;
  ss->Verbose_Cursor =              DEFAULT_VERBOSE_CURSOR;
  ss->Cursor_Update_Interval =      DEFAULT_CURSOR_UPDATE_INTERVAL;
  ss->Cursor_Location_Offset =      DEFAULT_CURSOR_LOCATION_OFFSET;
  ss->Show_Mix_Waveforms =          DEFAULT_SHOW_MIX_WAVEFORMS;
  ss->Mix_Waveform_Height =         DEFAULT_MIX_WAVEFORM_HEIGHT;
  ss->Mix_Tag_Width =               DEFAULT_MIX_TAG_WIDTH;
  ss->Mix_Tag_Height =              DEFAULT_MIX_TAG_HEIGHT;
  ss->With_Mix_Tags =               DEFAULT_WITH_MIX_TAGS;
  ss->Mark_Tag_Width =              DEFAULT_MARK_TAG_WIDTH;
  ss->Mark_Tag_Height =             DEFAULT_MARK_TAG_HEIGHT;
  ss->Show_Marks =                  DEFAULT_SHOW_MARKS;
  ss->Color_Map =                   DEFAULT_COLOR_MAP;
  ss->Color_Map_Size =              DEFAULT_COLOR_MAP_SIZE;
  ss->Color_Cutoff =                DEFAULT_COLOR_CUTOFF;
  ss->Color_Scale =                 DEFAULT_COLOR_SCALE;
  ss->Color_Inverted =              DEFAULT_COLOR_INVERTED;
  ss->Color_Map =                   DEFAULT_COLOR_MAP;
  ss->Wavo_Hop =                    DEFAULT_WAVO_HOP;
  ss->Wavo_Trace =                  DEFAULT_WAVO_TRACE;
  ss->Spectro_Hop =                 DEFAULT_SPECTRO_HOP;
  ss->Spectro_X_Scale =             DEFAULT_SPECTRO_X_SCALE;
  ss->Spectro_Y_Scale =             DEFAULT_SPECTRO_Y_SCALE;
  ss->Spectro_Z_Scale =             DEFAULT_SPECTRO_Z_SCALE;
  ss->Spectro_Z_Angle =             DEFAULT_SPECTRO_Z_ANGLE;
  ss->Spectro_X_Angle =             DEFAULT_SPECTRO_X_ANGLE;
  ss->Spectro_Y_Angle =             DEFAULT_SPECTRO_Y_ANGLE;
  ss->Spectrum_End =                DEFAULT_SPECTRUM_END;
  ss->Spectrum_Start =              DEFAULT_SPECTRUM_START;
  ss->Enved_Base =                  DEFAULT_ENVED_BASE;
  ss->Enved_Power =                 DEFAULT_ENVED_POWER;
  ss->Enved_Wave_p =                DEFAULT_ENVED_WAVE_P;
  ss->Enved_Style =                 DEFAULT_ENVED_STYLE;
  ss->Enved_Target =                DEFAULT_ENVED_TARGET;
  ss->Enved_Filter_Order =          DEFAULT_ENVED_FILTER_ORDER;
  ss->Eps_Bottom_Margin =           DEFAULT_EPS_BOTTOM_MARGIN;
  ss->Eps_Left_Margin =             DEFAULT_EPS_LEFT_MARGIN;
  ss->Eps_Size =                    DEFAULT_EPS_SIZE;
  ss->Expand_Control_Min =          DEFAULT_EXPAND_CONTROL_MIN;
  ss->Expand_Control_Max =          DEFAULT_EXPAND_CONTROL_MAX;
  ss->Amp_Control_Min =             DEFAULT_AMP_CONTROL_MIN;
  ss->Amp_Control_Max =             DEFAULT_AMP_CONTROL_MAX;
  ss->Speed_Control_Min =           DEFAULT_SPEED_CONTROL_MIN;
  ss->Speed_Control_Max =           DEFAULT_SPEED_CONTROL_MAX;
  ss->Contrast_Control_Min =        DEFAULT_CONTRAST_CONTROL_MIN;
  ss->Contrast_Control_Max =        DEFAULT_CONTRAST_CONTROL_MAX;
  ss->Contrast_Control_Amp =        DEFAULT_CONTRAST_CONTROL_AMP;
  ss->Expand_Control_Length =       DEFAULT_EXPAND_CONTROL_LENGTH;
  ss->Expand_Control_Ramp =         DEFAULT_EXPAND_CONTROL_RAMP;
  ss->Expand_Control_Hop =          DEFAULT_EXPAND_CONTROL_HOP;
  ss->Expand_Control_Jitter =       DEFAULT_EXPAND_CONTROL_JITTER;
  ss->Reverb_Control_Feedback =     DEFAULT_REVERB_CONTROL_FEEDBACK;
  ss->Reverb_Control_Lowpass =      DEFAULT_REVERB_CONTROL_LOWPASS;
  ss->Reverb_Control_Scale_Min =    DEFAULT_REVERB_CONTROL_SCALE_MIN;
  ss->Reverb_Control_Scale_Max =    DEFAULT_REVERB_CONTROL_SCALE_MAX;
  ss->Reverb_Control_Decay =        DEFAULT_REVERB_CONTROL_DECAY;
  ss->Speed_Control_Tones =         DEFAULT_SPEED_CONTROL_TONES;
  ss->Speed_Control_Style =         DEFAULT_SPEED_CONTROL_STYLE;
  ss->Reverb_Control_Length_Min =   DEFAULT_REVERB_CONTROL_LENGTH_MIN;
  ss->Reverb_Control_Length_Max =   DEFAULT_REVERB_CONTROL_LENGTH_MAX;
  ss->Filter_Control_Order =        DEFAULT_FILTER_CONTROL_ORDER;
  ss->Filter_Control_In_Db =        DEFAULT_FILTER_CONTROL_IN_DB;
  ss->Filter_Control_In_Hz =        DEFAULT_FILTER_CONTROL_IN_HZ;
  ss->Show_Controls =               DEFAULT_SHOW_CONTROLS;

  if (MUS_DEFAULT_TEMP_DIR != (char *)NULL)
    ss->Temp_Dir = mus_strdup(MUS_DEFAULT_TEMP_DIR);
  else ss->Temp_Dir = NULL;

  if (MUS_DEFAULT_SAVE_DIR != (char *)NULL)
    ss->Save_Dir = mus_strdup(MUS_DEFAULT_SAVE_DIR);
  else ss->Save_Dir = NULL;

  if (DEFAULT_LADSPA_DIR != (char *)NULL)
    ss->Ladspa_Dir = mus_strdup(DEFAULT_LADSPA_DIR);
  else ss->Ladspa_Dir = NULL;

  if (DEFAULT_SAVE_STATE_FILE != (char *)NULL)
    ss->Save_State_File = mus_strdup(DEFAULT_SAVE_STATE_FILE);
  else ss->Save_State_File = NULL;

  if (DEFAULT_PEAK_ENV_DIR != (char *)NULL)
    ss->Peak_Env_Dir = mus_strdup(DEFAULT_PEAK_ENV_DIR);
  else ss->Peak_Env_Dir = NULL;

  if (DEFAULT_EPS_FILE != (char *)NULL)
    ss->Eps_File = mus_strdup(DEFAULT_EPS_FILE);
  else ss->Eps_File = NULL;
}

#if HAVE_SETJMP_H && HAVE_SCHEME
static void jump_to_top_level(void)
{
  top_level_catch(1);
}
#endif

#if HAVE_GSL
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_errno.h>

/* default gsl error handler apparently aborts main program! */

static void snd_gsl_error(const char *reason, const char *file, int line, int gsl_errno)
{
  XEN_ERROR(XEN_ERROR_TYPE("gsl-error"),
	    XEN_LIST_6(C_TO_XEN_STRING("GSL: ~A, ~A in ~A line ~A, gsl err: ~A"),
		       C_TO_XEN_STRING(gsl_strerror(gsl_errno)),
		       C_TO_XEN_STRING(reason),
		       C_TO_XEN_STRING(file),
		       C_TO_XEN_INT(line),
		       C_TO_XEN_INT(gsl_errno)));
}
#endif

#if SND_AS_WIDGET
  snd_state *snd_main(int argc, char **argv)
#else
  int main(int argc, char **argv)
#endif
{
  int i;

#if HAVE_GSL
  /* if HAVE_GSL and the environment variable GSL_IEEE_MODE exists, use it */
  /* GSL_IEEE_MODE=double-precision,mask-underflow,mask-denormalized */
  if (getenv("GSL_IEEE_MODE") != NULL)
    gsl_ieee_env_setup();
  gsl_set_error_handler(snd_gsl_error);
#endif

  ss = (snd_state *)calloc(1, sizeof(snd_state)); /* not calloc! */
  ss->fam_ok = false;
  ss->startup_errors = NULL;

#if HAVE_GTK_3
  g_type_init();
#endif

  mus_sound_initialize(); /* has to precede version check (mus_audio_moniker needs to be setup in Alsa/Oss) */
  xen_initialize();

#if HAVE_SCHEME && HAVE_SETJMP_H
  s7_set_error_exiter(s7, jump_to_top_level);
#endif

  for (i = 1; i < argc; i++)
    {
      if (strcmp(argv[i], "--version") == 0)
	{
	  fprintf(stdout, "%s", version_info());
	  snd_exit(0);
	}
      else
	{
	  if (strcmp(argv[i], "--help") == 0)
	    {
	      fprintf(stdout, "%s", "Snd is a sound editor; see http://ccrma.stanford.edu/software/snd/.\n");
	      fprintf(stdout, "%s", version_info());
	      snd_exit(0);
	    }
	}
    }

  initialize_format_lists();
  snd_set_global_defaults(false);

#if MUS_DEBUGGING
  ss->Trap_Segfault = false;
#else
  ss->Trap_Segfault = DEFAULT_TRAP_SEGFAULT;
#endif
  ss->jump_ok = false;
  allocate_regions(max_regions(ss));
  ss->init_window_x = DEFAULT_INIT_WINDOW_X;
  ss->init_window_y = DEFAULT_INIT_WINDOW_Y;
  ss->init_window_width = DEFAULT_INIT_WINDOW_WIDTH;
  ss->init_window_height = DEFAULT_INIT_WINDOW_HEIGHT;
  ss->click_time = 100;
  init_sound_file_extensions();

  ss->max_sounds = 4;                 /* expands to accommodate any number of files */
  ss->sound_sync_max = 0;
  ss->stopped_explicitly = false;     /* C-g sets this flag so that we can interrupt various loops */
  ss->checking_explicitly = false;
  ss->selection_play_stop = false;
  ss->reloading_updated_file = 0;
  ss->selected_sound = NO_SELECTION;
  ss->sounds = (snd_info **)calloc(ss->max_sounds, sizeof(snd_info *));
  ss->print_choice = PRINT_SND;
  ss->graph_hook_active = false;
  ss->lisp_graph_hook_active = false;
  ss->exiting = false;
  ss->deferred_regions = 0;
  ss->fam_connection = NULL;
  ss->snd_error_data = NULL;
  ss->snd_error_handler = NULL;
  ss->snd_warning_data = NULL;
  ss->snd_warning_handler = NULL;
  ss->xen_error_data = NULL;
  ss->xen_error_handler = NULL;
  ss->update_sound_channel_style = NOT_A_CHANNEL_STYLE;

#if HAVE_GL && WITH_GL2PS
  ss->gl_printing = false;
#endif
  g_xen_initialize();
  ss->search_proc = XEN_UNDEFINED;
  ss->search_expr = NULL;
  ss->search_tree = NULL;
  mus_error_set_handler(mus_error_to_snd);
  mus_print_set_handler(mus_print_to_snd);

  initialize_load_path(); /* merge SND_PATH entries into the load-path */

#ifdef SND_AS_WIDGET
  return(ss);
#else
  snd_doit(argc, argv);
  return(0);
#endif
}

void g_init_base(void)
{
  #define H_mus_error_hook S_mus_error_hook " (error-type error-message):  called upon mus_error. \
If it returns " PROC_TRUE ", Snd ignores the error (it assumes you've handled it via the hook)."

  ss->mus_error_hook = XEN_DEFINE_HOOK(S_mus_error_hook, 2, H_mus_error_hook);       /* arg = error-type error-message */
}