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 <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_version.h>
#include <gsl/gsl_multifit_nlin.h>
#include <gsl/gsl_blas.h>
#include "psrsalsa.h"
typedef struct {
  size_t n;
  double *x;
  double *y;
  double *sigma;
  fitfunc_collection_type *fitfunction;
  size_t iter;
}levmar_internal_fitter_data_def;
static levmar_internal_fitter_data_def levmar_internal_fitter_data;
typedef struct {
  gsl_vector_view func_params;
  const gsl_multifit_fdfsolver_type *solver_type;
  gsl_multifit_function_fdf user_itteration_funcs;
  gsl_multifit_fdfsolver *solver;
  gsl_matrix *covar;
}levmar_internal_gsl_info;
typedef struct {
  double *covar;
  double *alpha;
  double *beta;
  double *alpha_tmp;
  double *beta_tmp;
  double *func_params_laststep;
  double *derivatives;
  double alambda;
  double chisq;
  long nrfitparams;
  levmar_internal_fitter_data_def *data;
}levmar_internal_psrsalsa_info;
void print_gsl_version_used(FILE *stream)
{
  int major, minor;
  major = GSL_VERSION_NUMBER/100.0;
  minor = GSL_VERSION_NUMBER-100*major;
  fprintf(stream, "%s (header) %d.%d (specified during compilation)", GSL_VERSION, major, minor);
}
double evaluate_fitfunc_collection(fitfunc_collection_type *function, double x, verbose_definition verbose)
{
  int i;
  double y;
  y = 0;
  for(i = 0; i < function->nrfuncs; i++) {
    if(function->func[i].type == FUNC_POLYNOMAL) {
      y += function->func[i].value[0]*pow(x, function->func[i].param[0]);
    }else {
      fflush(stdout);
      printerror(verbose.debug, "ERROR evaluate_fitfunction_collection: Unknown funtional type in specified function.");
      exit(0);
    }
  }
  return y;
}
double evaluate_fitfunc_collection_deriv_param(fitfunc_collection_type *function, int functionnr, int paramnr, double x, verbose_definition verbose)
{
  if(function->func[functionnr].type == FUNC_POLYNOMAL) {
    if(paramnr == 0) {
      return pow(x, function->func[functionnr].param[0]);
    }
  }else {
    fflush(stdout);
    printerror(verbose.debug, "ERROR evaluate_fitfunc_collection_deriv_param: Unknown funtional type in specified function.");
    exit(0);
  }
  fflush(stdout);
  printerror(verbose.debug, "ERROR evaluate_fitfunc_collection_deriv_param: Unknown paramer number specified for given funtional type.");
  exit(0);
}
int print_fitfunctions(fitfunc_collection_type *function, int novalue, int showerror, int index, verbose_definition verbose)
{
  int i, j, n;
  if(showerror)
    showerror = 1;
  for(n = 0; n <= showerror; n++) {
    if(n == 0)
      printf("y = ");
    j = 1;
    for(i = 0; i < function->nrfuncs; i++) {
      if(index == i || index < 0) {
 if(function->func[i].type == FUNC_POLYNOMAL) {
   if(novalue || showerror) {
     if(showerror && n == 1) {
       printf("    a%d = %13e +- %13e   (%lf +- %lf)\n", j, function->func[i].value[0], function->func[i].error[0], function->func[i].value[0], function->func[i].error[0]);
     }else {
       printf("a%d*x**%lf", j, function->func[i].param[0]);
     }
   }else {
     printf("%lf*x**%lf", function->func[i].value[0], function->func[i].param[0]);
   }
   j += 1;
 }else {
   fflush(stdout);
   printerror(verbose.debug, "ERROR print_fitfunctions: Unrecognized function type.");
   return 0;
 }
 if(i == function->nrfuncs - 1)
   printf("\n");
 else if(n == 0)
   printf(" + ");
      }
    }
  }
  return 1;
}
void countnrparameters_fitfunction(fitfunc_collection_type *function, int *nrfitparameters)
{
  int n, fnr, nrvalues;
  *nrfitparameters = 0;
  for(fnr = 0; fnr < function->nrfuncs; fnr++) {
    if(function->func[fnr].type == FUNC_POLYNOMAL) {
      nrvalues = 1;
    }else {
      fflush(stdout);
      printerror(0, "ERROR countnrparameters_fitfunction: Unknown funtional type.");
      exit(0);
    }
    for(n = 0; n < nrvalues; n++) {
      if(function->func[fnr].fit_flag[n]) {
 *nrfitparameters += 1;
      }
    }
  }
}
int set_fitted_parameters_fitfunc_collection(fitfunc_collection_type *function, double *fitparameters, verbose_definition verbose)
{
  int fnr, nrvalues, valnr, fitparamnr;
  fitparamnr = 0;
  for(fnr = 0; fnr < function->nrfuncs; fnr++) {
    if(function->func[fnr].type == FUNC_POLYNOMAL) {
      nrvalues = 1;
    }else {
      fflush(stdout);
      printerror(verbose.debug, "ERROR set_fitted_parameters_fitfunc_collection: Unknown funtional type.");
      return 0;
    }
    for(valnr = 0; valnr < nrvalues; valnr++) {
      if(function->func[fnr].fit_flag[valnr]) {
 function->func[fnr].value[valnr] = fitparameters[fitparamnr++];
      }
    }
  }
  return 1;
}
int levmar_itteration_calc_function_internal_gsl(const gsl_vector *fitparams, void *data, gsl_vector *f)
{
  int n;
  size_t npts = ((levmar_internal_fitter_data_def *)data)->n;
  double *xdata = ((levmar_internal_fitter_data_def *)data)->x;
  double *ydata = ((levmar_internal_fitter_data_def *)data)->y;
  double *sigma = ((levmar_internal_fitter_data_def *)data)->sigma;
  fitfunc_collection_type *fitfunction_internal = ((levmar_internal_fitter_data_def *)data)->fitfunction;
  double value;
  verbose_definition noverbose;
  cleanVerboseState(&noverbose);
  if(fitparams->stride != 1) {
    printerror(0, "ERROR levmar_itteration_calc_function_internal_gsl: It is assumed the gsl vectors have a stride of 1.");
    exit(0);
  }
  fitfunc_collection_type newfunction;
  memcpy(&newfunction, fitfunction_internal, sizeof(fitfunc_collection_type));
  if(set_fitted_parameters_fitfunc_collection(&newfunction, fitparams->data, noverbose) == 0) {
    printerror(0, "ERROR levmar_itteration_calc_function_internal_gsl: Copying new fit parameters failed.");
    exit(0);
  }
  for(n = 0; n < npts; n++) {
    value = evaluate_fitfunc_collection(&newfunction, xdata[n], noverbose);
    gsl_vector_set(f, n, (value - ydata[n])/sigma[n]);
  }
  return GSL_SUCCESS;
}
int levmar_itteration_calc_deriv_internal_gsl(const gsl_vector *fitparams, void *data, gsl_matrix *J)
{
  int fnr, j, n, nrvalues, paramnr;
  size_t npts = ((levmar_internal_fitter_data_def *)data)->n;
  double *xdata = ((levmar_internal_fitter_data_def *)data)->x;
  double *sigma = ((levmar_internal_fitter_data_def *)data)->sigma;
  fitfunc_collection_type *fitfunction_internal = ((levmar_internal_fitter_data_def *)data)->fitfunction;
  double deriv;
  verbose_definition noverbose;
  cleanVerboseState(&noverbose);
  if(fitparams->stride != 1) {
    printerror(0, "ERROR levmar_itteration_calc_function_internal_gsl: It is assumed the gsl vectors have a stride of 1.");
    exit(0);
  }
  fitfunc_collection_type newfunction;
  memcpy(&newfunction, fitfunction_internal, sizeof(fitfunc_collection_type));
  if(set_fitted_parameters_fitfunc_collection(&newfunction, fitparams->data, noverbose) == 0) {
    printerror(0, "ERROR levmar_itteration_calc_function_internal_gsl: Copying new fit parameters failed.");
    exit(0);
  }
  j = 0;
  for(fnr = 0; fnr < newfunction.nrfuncs; fnr++) {
    if(newfunction.func[fnr].type == FUNC_POLYNOMAL) {
      nrvalues = 1;
    }else {
      fflush(stdout);
      printerror(noverbose.debug, "ERROR levmar_itteration_calc_deriv_internal_gsl: Unknown funtional type.");
      return 0;
    }
    for(paramnr = 0; paramnr < nrvalues; paramnr++) {
      if(newfunction.func[fnr].fit_flag[paramnr]) {
 for(n = 0; n < npts; n++) {
   deriv = evaluate_fitfunc_collection_deriv_param(&newfunction, fnr, paramnr, xdata[n], noverbose);
   gsl_matrix_set(J, n, j, deriv/sigma[n]);
 }
 j++;
      }
    }
  }
  return GSL_SUCCESS;
}
int levmar_itteration_calc_func_and_deriv_internal_gsl(const gsl_vector *x, void *data, gsl_vector *f, gsl_matrix *J)
{
  levmar_itteration_calc_function_internal_gsl(x, data, f);
  levmar_itteration_calc_deriv_internal_gsl(x, data, J);
  return GSL_SUCCESS;
}
void levmar_itteration_calc_alpha_beta_chi2_internal_psrsalsa(double *params, void *info)
{
  long curDataPoint, nrDataPoints, curparam, curparam2, nrfitparams, fnr, nrvalues, paramnr;
  double *alpha, *beta, *derivatives, *x, *y, *sigma, chi2, ypred, tmp1, tmp2, tmp3, tmp4;
  verbose_definition noverbose;
  cleanVerboseState(&noverbose);
  nrDataPoints = ((levmar_internal_psrsalsa_info *)info)->data->n;
  x = ((levmar_internal_psrsalsa_info *)info)->data->x;
  y = ((levmar_internal_psrsalsa_info *)info)->data->y;
  sigma = ((levmar_internal_psrsalsa_info *)info)->data->sigma;
  alpha = ((levmar_internal_psrsalsa_info *)info)->alpha;
  beta = ((levmar_internal_psrsalsa_info *)info)->beta;
  derivatives = ((levmar_internal_psrsalsa_info *)info)->derivatives;
  nrfitparams = ((levmar_internal_psrsalsa_info *)info)->nrfitparams;
  fitfunc_collection_type newfunction;
  memcpy(&newfunction, levmar_internal_fitter_data.fitfunction, sizeof(fitfunc_collection_type));
  if(set_fitted_parameters_fitfunc_collection(&newfunction, params, noverbose) == 0) {
    printerror(0, "ERROR levmar_itteration_calc_func_and_deriv_internal_nr: Copying new fit parameters failed.");
    exit(0);
  }
  for(curparam = 0; curparam < nrfitparams; curparam++) {
    beta[curparam] = 0;
    for(curparam2 = curparam; curparam2 < nrfitparams; curparam2++) {
      alpha[curparam*nrfitparams+curparam2] = 0;
    }
  }
  chi2 = 0;
  for(curDataPoint = 0; curDataPoint < nrDataPoints; curDataPoint++) {
    ypred = evaluate_fitfunc_collection(&newfunction, x[curDataPoint], noverbose);
    curparam = 0;
    for(fnr = 0; fnr < newfunction.nrfuncs; fnr++) {
      if(newfunction.func[fnr].type == FUNC_POLYNOMAL) {
 nrvalues = 1;
      }else {
 fflush(stdout);
 printerror(noverbose.debug, "ERROR levmar_itteration_calc_alpha_beta_chi2_internal_psrsalsa: Unknown funtional type.");
 exit(0);
      }
      for(paramnr = 0; paramnr < nrvalues; paramnr++) {
 if(newfunction.func[fnr].fit_flag[paramnr]) {
   derivatives[curparam++] = evaluate_fitfunc_collection_deriv_param(&newfunction, fnr, paramnr, x[curDataPoint], noverbose);
 }
      }
    }
    tmp1 = 1.0/(sigma[curDataPoint]*sigma[curDataPoint]);
    tmp2 = (y[curDataPoint] - ypred);
    tmp3 = tmp1*tmp2;
    chi2 += tmp2*tmp3;
    for(curparam = 0; curparam < nrfitparams; curparam++) {
      tmp4 = derivatives[curparam]*tmp3;
      beta[curparam] += tmp4;
      for(curparam2 = curparam; curparam2 < nrfitparams; curparam2++) {
 alpha[curparam*nrfitparams+curparam2] += derivatives[curparam]*derivatives[curparam2]*tmp1;
      }
    }
  }
  for(curparam = 0; curparam < nrfitparams; curparam++) {
    for(curparam2 = curparam+1; curparam2 < nrfitparams; curparam2++) {
      alpha[curparam2*nrfitparams+curparam] = alpha[curparam*nrfitparams+curparam2];
    }
  }
  ((levmar_internal_psrsalsa_info *)info)->chisq = chi2;
}
int initialize_levmar_check_start_and_current_are_the_same(fitfunc_collection_type *function, verbose_definition verbose)
{
  int fnr, nrvalues, paramnr;
  for(fnr = 0; fnr < function->nrfuncs; fnr++) {
    if(function->func[fnr].type == FUNC_POLYNOMAL) {
      nrvalues = 1;
    }else {
      fflush(stdout);
      printerror(verbose.debug, "ERROR levmar_itteration_calc_func_and_deriv_internal_nr: Unknown funtional type.");
      return 0;
    }
    for(paramnr = 0; paramnr < nrvalues; paramnr++) {
      if(function->func[fnr].fit_flag[paramnr]) {
 if(function->func[fnr].value[paramnr] != function->func[fnr].start[paramnr]) {
   fflush(stdout);
   printerror(verbose.debug, "ERROR initialize_levmar_check_start_and_current_are_the_same: Start and current value appear to be different (%lf != %lf)", function->func[fnr].value[paramnr], function->func[fnr].start[paramnr]);
   return 0;
 }
      }
    }
  }
  return 1;
}
int levmar_initialize(fitfunc_collection_type *function, double **func_params, double **func_params_err, int *nrfitparameters, double *datax, double *datay, double *datasigma, long nrdatapoints, levmar_internal_gsl_info *gsl_params, levmar_internal_psrsalsa_info *psrsalsa_params,
#ifdef NRAVAIL
levmar_internal_nr_info *nr_params,
#endif
int algorithm, int mode, int showcovariance, verbose_definition verbose)
{
  int fnr, nrvalues, paramnr, n, n2;
  if(mode == 0) {
    countnrparameters_fitfunction(function, nrfitparameters);
    if(verbose.verbose) printf("  Initializing fitter with %d fit parameters\n", *nrfitparameters);
  }else if(mode == 1) {
    levmar_internal_fitter_data.n = nrdatapoints;
    levmar_internal_fitter_data.x = datax;
    levmar_internal_fitter_data.y = datay;
    levmar_internal_fitter_data.sigma = datasigma;
    levmar_internal_fitter_data.fitfunction = function;
    *func_params = malloc(*nrfitparameters*sizeof(double));
    *func_params_err = malloc(*nrfitparameters*sizeof(double));
    if(*func_params == NULL || *func_params_err == NULL) {
      fflush(stdout);
      printerror(verbose.debug, "ERROR levmar_initialize: Cannot allocate memory");
      return 0;
    }
    paramnr = 0;
    for(fnr = 0; fnr < function->nrfuncs; fnr++) {
      if(function->func[fnr].type == FUNC_POLYNOMAL) {
 nrvalues = 1;
      }else {
 fflush(stdout);
 printerror(verbose.debug, "ERROR levmar_initialize: Unknown funtional type.");
 return 0;
      }
      for(n = 0; n < nrvalues; n++) {
 if(function->func[fnr].fit_flag[n]) {
   (*func_params)[paramnr++] = function->func[fnr].value[n];
 }
 function->func[fnr].error[n] = sqrt(-1);
      }
    }
    if(algorithm == 1) {
      gsl_params->func_params = gsl_vector_view_array(*func_params, *nrfitparameters);
      gsl_params->covar = gsl_matrix_alloc(*nrfitparameters, *nrfitparameters);
      gsl_params->user_itteration_funcs.f = &levmar_itteration_calc_function_internal_gsl;
      gsl_params->user_itteration_funcs.df = &levmar_itteration_calc_deriv_internal_gsl;
      gsl_params->user_itteration_funcs.fdf = &levmar_itteration_calc_func_and_deriv_internal_gsl;
      gsl_params->user_itteration_funcs.n = nrdatapoints;
      gsl_params->user_itteration_funcs.p = *nrfitparameters;
      gsl_params->user_itteration_funcs.params = &levmar_internal_fitter_data;
      gsl_params->solver_type = gsl_multifit_fdfsolver_lmsder;
      gsl_params->solver = gsl_multifit_fdfsolver_alloc (gsl_params->solver_type, nrdatapoints, *nrfitparameters);
      gsl_multifit_fdfsolver_set(gsl_params->solver, &(gsl_params->user_itteration_funcs), &(gsl_params->func_params.vector));
    }else if(algorithm == 2) {
      psrsalsa_params->func_params_laststep = malloc(sizeof(double)*(*nrfitparameters));
      psrsalsa_params->beta = malloc(sizeof(double)*(*nrfitparameters));
      psrsalsa_params->derivatives = malloc(sizeof(double)*(*nrfitparameters));
      psrsalsa_params->covar = malloc(sizeof(double)*(*nrfitparameters)*(*nrfitparameters));
      psrsalsa_params->alpha = malloc(sizeof(double)*(*nrfitparameters)*(*nrfitparameters));
      psrsalsa_params->alpha_tmp = malloc(sizeof(double)*(*nrfitparameters)*(*nrfitparameters));
      psrsalsa_params->beta_tmp = malloc(sizeof(double)*(*nrfitparameters));
      if(psrsalsa_params->func_params_laststep == NULL || psrsalsa_params->beta == NULL || psrsalsa_params->derivatives == NULL || psrsalsa_params->alpha == NULL || psrsalsa_params->covar == NULL || psrsalsa_params->beta_tmp == NULL || psrsalsa_params->alpha_tmp == NULL) {
 printerror(verbose.debug, "ERROR levmar_initialize: Memory allocation error");
 return 0;
      }
      psrsalsa_params->nrfitparams = *nrfitparameters;
      psrsalsa_params->data = &levmar_internal_fitter_data;
      psrsalsa_params->alambda = 0.001;
      levmar_itteration_calc_alpha_beta_chi2_internal_psrsalsa(*func_params, psrsalsa_params);
    }else if(algorithm == 3) {
#ifdef NRAVAIL
      nr_params->ia = malloc(sizeof(int)*(*nrfitparameters));
      nr_params->func_params_laststep = malloc(sizeof(double)*(*nrfitparameters));
      if(nr_params->ia == NULL || nr_params->func_params_laststep == NULL) {
 printerror(verbose.debug, "ERROR levmar_initialize: Memory allocation error");
 return 0;
      }
      for(n = 0; n < *nrfitparameters; n++)
 nr_params->ia[n] = 1;
      double **matrix_nr_d(long nrl, long nrh, long ncl, long nch);
      nr_params->covar = matrix_nr_d(1, *nrfitparameters, 1, *nrfitparameters);
      nr_params->alpha = matrix_nr_d(1, *nrfitparameters, 1, *nrfitparameters);
      void mrqmin_nr_d(double x[], double y[], double sig[], int ndata, double a[], int ia[],
     int ma, double **covar, double **alpha, double *chisq,
     void (*funcs)(double, double [], double *, double [], int), double *alamda);
      nr_params->alambda = -1;
      mrqmin_nr_d(levmar_internal_fitter_data.x-1, levmar_internal_fitter_data.y-1, levmar_internal_fitter_data.sigma-1, nrdatapoints, (*func_params)-1, nr_params->ia-1, *nrfitparameters, nr_params->covar, nr_params->alpha, &nr_params->chisq, levmar_itteration_calc_func_and_deriv_internal_nr, &(nr_params->alambda));
#else
      printerror(verbose.debug, "ERROR levmar_initialize: Code is not compiled with NR support");
      return 0;
#endif
    }
  }else if(mode == 2) {
    paramnr = 0;
    for(fnr = 0; fnr < function->nrfuncs; fnr++) {
      if(function->func[fnr].type == FUNC_POLYNOMAL) {
 nrvalues = 1;
      }else {
 fflush(stdout);
 printerror(0, "ERROR levmar_initialize: Unknown funtional type.");
 return 0;
      }
      for(n = 0; n < nrvalues; n++) {
 if(function->func[fnr].fit_flag[n]) {
   function->func[fnr].value[n] = (*func_params)[paramnr];
   function->func[fnr].error[n] = (*func_params_err)[paramnr];
   paramnr++;
 }
      }
    }
    if(algorithm == 1) {
      if(showcovariance) {
 printf("Covariance matrix: \n");
 for(n = 0; n < *nrfitparameters; n++) {
   for(n2 = 0; n2 < *nrfitparameters; n2++) {
     printf("%lf ", gsl_matrix_get(gsl_params->covar,n,n2));
   }
   printf("\n");
 }
      }
      gsl_multifit_fdfsolver_free(gsl_params->solver);
      gsl_matrix_free(gsl_params->covar);
    }else if(algorithm == 2) {
      if(showcovariance) {
 printf("Covariance matrix: \n");
 for(n = 0; n < *nrfitparameters; n++) {
   for(n2 = 0; n2 < *nrfitparameters; n2++) {
     printf("%lf ", psrsalsa_params->covar[n*(*nrfitparameters)+n2]);
   }
   printf("\n");
 }
      }
      free(psrsalsa_params->beta);
      free(psrsalsa_params->derivatives);
      free(psrsalsa_params->alpha);
      free(psrsalsa_params->covar);
      free(psrsalsa_params->func_params_laststep);
      free(psrsalsa_params->alpha_tmp);
      free(psrsalsa_params->beta_tmp);
    }else if(algorithm == 3) {
#ifdef NRAVAIL
      if(showcovariance) {
 printf("Covariance matrix: \n");
 for(n = 0; n < *nrfitparameters; n++) {
   for(n2 = 0; n2 < *nrfitparameters; n2++) {
     printf("%lf ", nr_params->covar[n+1][n2+1]);
   }
   printf("\n");
 }
      }
      void free_matrix_nr_d(double **m, long nrl, long nrh, long ncl, long nch);
      free_matrix_nr_d(nr_params->covar, 1, *nrfitparameters, 1, *nrfitparameters);
      free_matrix_nr_d(nr_params->alpha, 1, *nrfitparameters, 1, *nrfitparameters);
      free(nr_params->ia);
      free(nr_params->func_params_laststep);
#else
      printerror(verbose.debug, "ERROR levmar_initialize: Code is not compiled with NR support");
      return 0;
#endif
    }
    free(*func_params);
    free(*func_params_err);
  }
  return 1;
}
int fit_levmar_internal(int algorithm, int nrdatapoints, int nrfitparameters, double *func_params, double *func_params_err, levmar_internal_gsl_info *gsl_params, levmar_internal_psrsalsa_info *psrsalsa_params,
#ifdef NRAVAIL
   levmar_internal_nr_info *nr_params,
#endif
   double epsabs, double epsrel, int maxiter, verbose_definition verbose)
{
  int status, i, converged;
  long iter;
  double chisq_last_nr;
#ifdef NRAVAIL
  void mrqmin_nr_d(double x[], double y[], double sig[], int ndata, double a[], int ia[],
     int ma, double **covar, double **alpha, double *chisq,
     void (*funcs)(double, double [], double *, double [], int), double *alamda);
#else
  if(algorithm == 3) {
    printerror(verbose.debug, "ERROR fit_levmar_internal: Code is not compiled with NR support");
    return 0;
  }
#endif
  if(algorithm == 2) {
    memcpy(psrsalsa_params->func_params_laststep, func_params, sizeof(double)*nrfitparameters);
    chisq_last_nr = psrsalsa_params->chisq;
#ifdef NRAVAIL
  }else if(algorithm == 3) {
    memcpy(nr_params->func_params_laststep, func_params, sizeof(double)*nrfitparameters);
    chisq_last_nr = nr_params->chisq;
#endif
  }
  iter = 0;
  do {
    iter++;
    levmar_internal_fitter_data.iter += 1;
    if(algorithm == 1) {
      status = gsl_multifit_fdfsolver_iterate(gsl_params->solver);
      if(status == GSL_ETOLF || status == GSL_ETOLX || status == GSL_ETOLG) {
 int status2;
 status2 = gsl_multifit_test_delta(gsl_params->solver->dx, gsl_params->solver->x, epsabs, epsrel);
 if(status2 == GSL_SUCCESS) {
   status = 0;
 }else {
   if(verbose.debug == 0) {
     printwarning(verbose.debug, "Machine precision reached, but required precision tolerance conditions are not met.");
   }else {
     printwarning(verbose.debug, "Machine precision reached, but required precision tolerance conditions are not met:");
     for(i = 0; i < nrfitparameters; i++) {
       printwarning(verbose.debug, "parameter %d: value=%e stepsize=%e", i+1, gsl_vector_get(gsl_params->solver->x, i), gsl_vector_get(gsl_params->solver->dx, i));
     }
   }
 }
      }
      if(verbose.debug) {
 printf("fit_levmar_internal: itteration %ld - status = %s\n", iter, gsl_strerror(status));
      }
      if(status) {
 break;
      }else {
 status = gsl_multifit_test_delta(gsl_params->solver->dx, gsl_params->solver->x, epsabs, epsrel);
      }
    }else if(algorithm == 2) {
      for(i = 0; i < nrfitparameters*nrfitparameters; i++)
 psrsalsa_params->alpha_tmp[i] = psrsalsa_params->alpha[i];
      for(i = 0; i < nrfitparameters; i++)
 psrsalsa_params->beta_tmp[i] = psrsalsa_params->beta[i];
      for (i = 0; i < nrfitparameters; i++) {
 psrsalsa_params->alpha[i*nrfitparameters+i] *= 1.0+psrsalsa_params->alambda;
      }
      int linalg_solve_matrix_eq_gauss_jordan(double *matrixa, double *matrixb, int n, int m, verbose_definition verbose);
      if(linalg_solve_matrix_eq_gauss_jordan(psrsalsa_params->alpha, psrsalsa_params->beta, nrfitparameters, 1, verbose) != 0) {
 printerror(verbose.debug, "ERROR fit_levmar_internal: Solving matrix equation failed.");
 return 9;
      }
      for (i = 0; i < nrfitparameters; i++)
 func_params[i] += psrsalsa_params->beta[i];
      levmar_itteration_calc_alpha_beta_chi2_internal_psrsalsa(func_params, psrsalsa_params);
      status = GSL_SUCCESS;
      if(chisq_last_nr < psrsalsa_params->chisq) {
 status = GSL_ETOLF;
 memcpy(func_params, psrsalsa_params->func_params_laststep, sizeof(double)*nrfitparameters);
 for(i = 0; i < nrfitparameters*nrfitparameters; i++)
   psrsalsa_params->alpha[i] = psrsalsa_params->alpha_tmp[i];
 for(i = 0; i < nrfitparameters; i++)
   psrsalsa_params->beta[i] = psrsalsa_params->beta_tmp[i];
 psrsalsa_params->alambda *= 10.1;
      }
      converged = 1;
      for(i = 0; i < nrfitparameters; i++) {
 if(fabs(func_params[i] - psrsalsa_params->func_params_laststep[i]) >= 1.2e-16+epsabs+epsrel*fabs(psrsalsa_params->func_params_laststep[i])) {
   converged = 0;
   break;
 }
      }
      if(converged && status == GSL_ETOLF)
 status = GSL_SUCCESS;
      if(status == GSL_ETOLF) {
       printwarning(verbose.debug, "Machine precision reached, but required precision tolerance conditions are not met.");
      }else if(converged == 0) {
 status = GSL_CONTINUE;
      }
      if(chisq_last_nr > psrsalsa_params->chisq) {
 memcpy(psrsalsa_params->func_params_laststep, func_params, sizeof(double)*nrfitparameters);
 chisq_last_nr = psrsalsa_params->chisq;
 psrsalsa_params->alambda *= 0.1;
      }
      if(verbose.debug) {
 printf("fit_levmar_internal: itteration %ld - status = %s\n", iter, gsl_strerror(status));
      }
#ifdef NRAVAIL
    }else if(algorithm == 3) {
      mrqmin_nr_d(levmar_internal_fitter_data.x-1, levmar_internal_fitter_data.y-1, levmar_internal_fitter_data.sigma-1, nrdatapoints, func_params-1, nr_params->ia-1, nrfitparameters, nr_params->covar, nr_params->alpha, &(nr_params->chisq), levmar_itteration_calc_func_and_deriv_internal_nr, &(nr_params->alambda));
      status = GSL_SUCCESS;
      for(i = 0; i < nrfitparameters; i++) {
 if(fabs(func_params[i] - nr_params->func_params_laststep[i]) >= epsabs+epsrel*fabs(nr_params->func_params_laststep[i])) {
   status = GSL_CONTINUE;
   break;
 }
      }
      if(verbose.debug) {
 printf("fit_levmar_internal: itteration %ld - status = %s\n", iter, gsl_strerror(status));
      }
      if(status == GSL_CONTINUE) {
 if(chisq_last_nr <= nr_params->chisq) {
   status = GSL_ETOLF;
   printwarning(verbose.debug, "Machine precision reached, but required precision tolerance conditions are not met.");
 }else {
   memcpy(nr_params->func_params_laststep, func_params, sizeof(double)*nrfitparameters);
   chisq_last_nr = nr_params->chisq;
 }
      }
#endif
    }
  }while(status == GSL_CONTINUE && iter < maxiter);
  if(algorithm == 1) {
#if GSL_VERSION_NUMBER >= 200
    if(verbose.debug) {
      printf("fit_levmar_internal: Entering gsl >= 2.0 specific part of the code\n");
    }
    gsl_matrix *J = gsl_matrix_alloc(nrdatapoints, nrfitparameters);
    gsl_multifit_fdfsolver_jac(gsl_params->solver, J);
    gsl_multifit_covar(J, 0.0, gsl_params->covar);
#else
    if(verbose.debug) {
      printf("fit_levmar_internal: Entering gsl < 2.0 specific part of the code\n");
    }
    gsl_multifit_covar(gsl_params->solver->J, 0.0, gsl_params->covar);
#endif
    levmar_internal_fitter_data.fitfunction->chi2 = gsl_blas_dnrm2(gsl_params->solver->f);
    levmar_internal_fitter_data.fitfunction->chi2_red = (levmar_internal_fitter_data.fitfunction->chi2)/sqrt(nrdatapoints - nrfitparameters);
    levmar_internal_fitter_data.fitfunction->chi2_red *= levmar_internal_fitter_data.fitfunction->chi2_red;
    levmar_internal_fitter_data.fitfunction->chi2 *= levmar_internal_fitter_data.fitfunction->chi2;
  }else if(algorithm == 2) {
    int linalg_solve_matrix_eq_gauss_jordan(double *matrixa, double *matrixb, int n, int m, verbose_definition verbose);
    if(linalg_solve_matrix_eq_gauss_jordan(psrsalsa_params->alpha, psrsalsa_params->beta, nrfitparameters, 1, verbose) != 0) {
      printerror(verbose.debug, "ERROR fit_levmar_internal: Solving matrix equation failed.");
      return 9;
    }
    for(i = 0; i < nrfitparameters*nrfitparameters; i++)
      psrsalsa_params->covar[i] = psrsalsa_params->alpha[i];
    levmar_internal_fitter_data.fitfunction->chi2 = psrsalsa_params->chisq;
    levmar_internal_fitter_data.fitfunction->chi2_red = psrsalsa_params->chisq/(double)(nrdatapoints - nrfitparameters);
#ifdef NRAVAIL
  }else if(algorithm == 3) {
    nr_params->alambda = 0;
    mrqmin_nr_d(levmar_internal_fitter_data.x-1, levmar_internal_fitter_data.y-1, levmar_internal_fitter_data.sigma-1, nrdatapoints, func_params-1, nr_params->ia-1, nrfitparameters, nr_params->covar, nr_params->alpha, &nr_params->chisq, levmar_itteration_calc_func_and_deriv_internal_nr, &nr_params->alambda);
    levmar_internal_fitter_data.fitfunction->chi2 = nr_params->chisq;
    levmar_internal_fitter_data.fitfunction->chi2_red = nr_params->chisq/(double)(nrdatapoints - nrfitparameters);
#endif
  }
  if(verbose.debug) {
    printf("chisq/dof = %g after %ld iterations\n", levmar_internal_fitter_data.fitfunction->chi2_red, iter);
  }
  if(algorithm == 1) {
    for(i = 0; i < nrfitparameters; i++) {
      func_params[i] = gsl_vector_get(gsl_params->solver->x, i);
      func_params_err[i] = sqrt(gsl_matrix_get(gsl_params->covar,i,i));
    }
  }else if(algorithm == 2) {
    for(i = 0; i < nrfitparameters; i++) {
      func_params_err[i] = sqrt(psrsalsa_params->covar[i*nrfitparameters+i]);
    }
#ifdef NRAVAIL
  }else if(algorithm == 3) {
    for(i = 0; i < nrfitparameters; i++) {
      func_params_err[i] = sqrt(nr_params->covar[i+1][i+1]);
    }
#endif
  }
  if(levmar_internal_fitter_data.iter == maxiter)
    return 2;
  if(algorithm == 1) {
    if(status == 0)
      return 0;
    if(status == GSL_ETOLF || status == GSL_ETOLX || status == GSL_ETOLG) {
      return 3;
    }
    if(status == GSL_ENOPROG) {
      fflush(stdout);
      printerror(verbose.debug, "fit_levmar_internal: Cannot determine suitable trial step. Continue itterating might help");
      return 4;
    }
  }else if(algorithm == 2 || algorithm == 3) {
    return 0;
  }
  fflush(stdout);
  printerror(verbose.debug, "fit_levmar_internal: UNKNOWN RETURN CODE");
  return 10;
}
int fit_levmar(int algorithm, fitfunc_collection_type *function, double *data_x, double *data_y, double *data_sigma, long ndata, int oneatatime, int force_chi2_1, double epsabs, double epsrel, int maxiter, int *status, int showresults, int showcovariance, verbose_definition verbose)
{
  int loopnr, improved, i, j, k, fnr, valuenr, nrfitparameters_one, nrfitparameters, free_data_sigma, nrvalues;
  double old_chi2, sig_scale, *func_params, *func_params_err;
  verbose_definition noverbose;
  levmar_internal_gsl_info gsl_params;
  levmar_internal_psrsalsa_info psrsalsa_params;
#ifndef NRAVAIL
  if(algorithm == 3) {
    printerror(verbose.debug, "ERROR fit_levmar: Code is not compiled with NR support");
    return 0;
  }
#else
  levmar_internal_nr_info nr_params;
#endif
  cleanVerboseState(&noverbose);
  noverbose.nocounters = 1;
  if(verbose.verbose) {
    printf("Levenberg-Marquardt algorithm:\n");
    if(verbose.debug) {
      printf("  Fit function: ");
      print_fitfunctions(function, 1, 0, -1, verbose);
    }
  }
  if(initialize_levmar_check_start_and_current_are_the_same(function, verbose) == 0) {
    fflush(stdout);
    printerror(verbose.debug, "ERROR fit_levmar: Start and current value appear to be different");
    return 1;
  }
  free_data_sigma = 0;
  if(data_sigma == NULL) {
    data_sigma = malloc(ndata*sizeof(double));
    for(i = 0; i < ndata; i++)
      data_sigma[i] = 1;
    free_data_sigma = 1;
    force_chi2_1 = 1;
  }
  if(levmar_initialize(function, &func_params, &func_params_err, &nrfitparameters, data_x, data_y, data_sigma, ndata, &gsl_params, &psrsalsa_params,
#ifdef NRAVAIL
         &nr_params,
#endif
         algorithm, 0, showcovariance, verbose) == 0) {
    return 1;
  }
  levmar_internal_fitter_data.iter = 0;
  if(oneatatime == 0) {
    if(levmar_initialize(function, &func_params, &func_params_err, &nrfitparameters, data_x, data_y, data_sigma, ndata, &gsl_params, &psrsalsa_params,
#ifdef NRAVAIL
    &nr_params,
#endif
    algorithm, 1, showcovariance, verbose) == 0) {
      return 1;
    }
    *status = fit_levmar_internal(algorithm, ndata, nrfitparameters, func_params, func_params_err, &gsl_params, &psrsalsa_params,
#ifdef NRAVAIL
      &nr_params,
#endif
      epsabs, epsrel, maxiter, verbose);
  }else {
    fitfunc_collection_type *original_function;
    original_function = malloc(sizeof(fitfunc_collection_type));
    if(original_function == NULL) {
      printerror(verbose.debug, "ERROR fit_levmar: Memory allocation error.");
      return 1;
    }
    memcpy(original_function, function, sizeof(fitfunc_collection_type));
    old_chi2 = -1;
    loopnr = 0;
    do {
      improved = 0;
      i = 0;
      for(fnr = 0; fnr < function->nrfuncs; fnr++) {
 if(function->func[fnr].type == FUNC_POLYNOMAL) {
   nrvalues = 1;
 }else {
   fflush(stdout);
   printerror(verbose.debug, "ERROR fit_levmar: Unknown funtional type.");
   return 1;
 }
 for(valuenr = 0; valuenr < nrvalues; valuenr++) {
   if(original_function->func[fnr].fit_flag[valuenr]) {
     for(j = 0; j < function->nrfuncs; j++) {
       for(k = 0; k < MaxNrFitParameters; k++)
  function->func[j].fit_flag[k] = 0;
     }
     function->func[fnr].fit_flag[valuenr] = 1;
     if(verbose.debug) {
       printf("  Loop %d: chi2 = %f", loopnr, old_chi2);
       printf("\n");
     }
     if(levmar_initialize(function, &func_params, &func_params_err, &nrfitparameters_one, data_x, data_y, data_sigma, ndata, &gsl_params, &psrsalsa_params,
#ifdef NRAVAIL
     &nr_params,
#endif
     algorithm, 0, showcovariance, noverbose) == 0) {
       free(original_function);
       return 1;
     }
     if(levmar_initialize(function, &func_params, &func_params_err, &nrfitparameters_one, data_x, data_y, data_sigma, ndata, &gsl_params, &psrsalsa_params,
#ifdef NRAVAIL
     &nr_params,
#endif
     algorithm, 1, showcovariance, verbose) == 0) {
       free(original_function);
       return 1;
     }
     *status = fit_levmar_internal(algorithm, ndata, nrfitparameters_one, func_params, func_params_err, &gsl_params, &psrsalsa_params,
#ifdef NRAVAIL
       &nr_params,
#endif
       epsabs, epsrel, maxiter, verbose);
     if(levmar_initialize(function, &func_params, &func_params_err, &nrfitparameters_one, data_x, data_y, data_sigma, ndata, &gsl_params, &psrsalsa_params,
#ifdef NRAVAIL
     &nr_params,
#endif
     algorithm, 2, showcovariance, verbose) == 0) {
       free(original_function);
       return 1;
     }
     if(function->chi2 < old_chi2 || old_chi2 < 0) {
       old_chi2 = function->chi2;
       improved = 1;
     }
     loopnr++;
     i++;
   }
 }
      }
    }while(improved);
    for(j = 0; j < function->nrfuncs; j++) {
      for(k = 0; k < MaxNrFitParameters; k++)
 function->func[j].fit_flag[k] = original_function->func[j].fit_flag[k];
    }
    if(levmar_initialize(function, &func_params, &func_params_err, &nrfitparameters, data_x, data_y, data_sigma, ndata, &gsl_params, &psrsalsa_params,
#ifdef NRAVAIL
    &nr_params,
#endif
    algorithm, 0, showcovariance, verbose) == 0) {
      free(original_function);
      return 1;
    }
    if(levmar_initialize(function, &func_params, &func_params_err, &nrfitparameters, data_x, data_y, data_sigma, ndata, &gsl_params, &psrsalsa_params,
#ifdef NRAVAIL
    &nr_params,
#endif
    algorithm, 1, showcovariance, verbose) == 0) {
      free(original_function);
      return 1;
    }
    *status = fit_levmar_internal(algorithm, ndata, nrfitparameters, func_params, func_params_err, &gsl_params, &psrsalsa_params,
#ifdef NRAVAIL
      &nr_params,
#endif
      epsabs, epsrel, maxiter, verbose);
    free(original_function);
  }
  if(verbose.verbose) {
    printf("  After %ld iterations found chisq/dof = %g\n", levmar_internal_fitter_data.iter, levmar_internal_fitter_data.fitfunction->chi2_red);
  }
  if(levmar_initialize(function, &func_params, &func_params_err, &nrfitparameters, data_x, data_y, data_sigma, ndata, &gsl_params, &psrsalsa_params,
#ifdef NRAVAIL
         &nr_params,
#endif
         algorithm, 2, showcovariance, verbose) == 0) {
    return 1;
  }
  if(force_chi2_1) {
    sig_scale = sqrt(function->chi2_red);
    for(i = 0; i < ndata; i++) {
      data_sigma[i] *= sig_scale;
    }
    if(levmar_initialize(function, &func_params, &func_params_err, &nrfitparameters, data_x, data_y, data_sigma, ndata, &gsl_params, &psrsalsa_params,
#ifdef NRAVAIL
    &nr_params,
#endif
    algorithm, 0, showcovariance, verbose) == 0) {
      return 1;
    }
    if(levmar_initialize(function, &func_params, &func_params_err, &nrfitparameters, data_x, data_y, data_sigma, ndata, &gsl_params, &psrsalsa_params,
#ifdef NRAVAIL
    &nr_params,
#endif
    algorithm, 1, showcovariance, verbose) == 0) {
      return 1;
    }
    *status = fit_levmar_internal(algorithm, ndata, nrfitparameters, func_params, func_params_err, &gsl_params, &psrsalsa_params,
#ifdef NRAVAIL
      &nr_params,
#endif
      epsabs, epsrel, maxiter, verbose);
    if(levmar_initialize(function, &func_params, &func_params_err, &nrfitparameters, data_x, data_y, data_sigma, ndata, &gsl_params, &psrsalsa_params,
#ifdef NRAVAIL
    &nr_params,
#endif
    algorithm, 2, showcovariance, verbose) == 0) {
      return 1;
    }
  }
  if(verbose.verbose || showresults) {
    printf("  Fit function: ");
    print_fitfunctions(function, 0, 1, -1, verbose);
    printf("  ");
    print_fitfunctions(function, 0, 0, -1, verbose);
    printf("\n  chi2 = %f\n  reduced chi2 = %f (%ld points - %d params)\n", function->chi2, function->chi2_red, ndata, nrfitparameters);
  }
  if(free_data_sigma) {
    free(data_sigma);
    data_sigma = NULL;
  }
  return *status;
}

max_stars_repo_name: David-McKenna/psrsalsa

max_issues_repo_name: David-McKenna/psrsalsa

max_forks_repo_name: David-McKenna/psrsalsa


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

#ifndef GENETIC_ALGORITHM_LIBRARY_ROULETTE_WHEEL_SELECTION_H
#define GENETIC_ALGORITHM_LIBRARY_ROULETTE_WHEEL_SELECTION_H

#include <algorithm>
#include <cstddef>
#include <limits>
#include <random>
#include <utility>
#include <type_traits>
#include <vector>
#include <gsl/gsl_assert>
#include "identity.h"

template<class UniformRandomBitGenerator, class Function = identity>
class roulette_wheel_selection {
public:
	explicit roulette_wheel_selection(UniformRandomBitGenerator& g, const Function& f = Function()) noexcept(noexcept(Function(f)));
	template<class Specimen>
	void operator()(std::vector<Specimen>& specimens, std::size_t n);
private:
	UniformRandomBitGenerator& rand;
	Function probability_function;
};

template<class UniformRandomBitGenerator, class Function>
inline roulette_wheel_selection<UniformRandomBitGenerator, Function>::roulette_wheel_selection(UniformRandomBitGenerator& g, const Function& f) noexcept(noexcept(Function(f)))
	: rand(g), probability_function(f) {}

template<class UniformRandomBitGenerator, class Function>
template<class Specimen>
inline void roulette_wheel_selection<UniformRandomBitGenerator, Function>::operator()(std::vector<Specimen>& specimens, std::size_t n) {
	Expects(specimens.size() >= n);
	using sample_type = std::common_type_t<double, decltype(probability_function(std::declval<typename Specimen::rating_type>()))>;
	using iterator_type = typename std::vector<Specimen>::iterator;
	using pair_type = std::pair<sample_type, iterator_type>;
	std::vector<pair_type> samples;
	samples.reserve(specimens.size());
	for (iterator_type it = specimens.begin(); it != specimens.end(); ++it) {
		const sample_type probability = probability_function(it->rating());
		sample_type sample = std::numeric_limits<sample_type>::max();
		if (probability > 0.0) {
			const std::exponential_distribution<sample_type> distribution(probability);
			sample = distribution(rand);
		}
		samples.emplace_back(sample, it);
	}
	std::nth_element(samples.begin(), samples.begin() + n, samples.end());
	samples.resize(n);
	iterator_type it = specimens.begin();
	for (const auto& sample : samples) {
		std::iter_swap(it++, sample.second);
	}
	specimens.resize(n);
}

#endif

ext: h

max_stars_repo_name: SirEmentaler/Eugenics-Wars

max_issues_repo_name: SirEmentaler/Eugenics-Wars

max_forks_repo_name: SirEmentaler/Eugenics-Wars


text:
#include "quantum_gates.h"
#include "quac_p.h"
#include <stdlib.h>
#include <stdio.h>
#include <petsc.h>
#include <stdarg.h>
#include "qasm_parser.h"
#include "error_correction.h"
#include "dm_utilities.h"
#include <ctype.h>

void quil_read(char filename[],PetscInt *num_qubits,circuit *circ){
  FILE *fp;
  int ch = 0,lines=0,finished_allocate=0,num_pragma=0,this_qubit,max_qubit=-1;
  char *line = NULL;
  size_t len = 0;
  ssize_t read;
  PetscReal time=1.0;

  fp = fopen(filename,"r");

  if (fp == NULL){
    if (nid==0){
      printf("ERROR! File not found in projectq_qasm_read!\n");
    }
  }

  //Count number of lines
  while(!feof(fp)){
    ch = fgetc(fp);
    if(ch == '\n'){
      lines++;
    }
  }
  //Rewind file
  rewind(fp);

  *num_qubits = 0;
  num_pragma  = 0;
  this_qubit  = -1;
  while ((read = getline(&line, &len, fp)) != -1){
    if (strstr(line,"PRAGMA")){
      //We will skip pragma lines
      num_pragma = num_pragma + 1;
    } else {
      this_qubit = line[strlen(line) - 2] - '0';
      if (this_qubit>max_qubit){
        max_qubit = this_qubit;
      }
    }
  }
  *num_qubits = max_qubit+1;

  lines = lines - num_pragma;

  create_circuit(circ,lines);
  rewind(fp); //Rewind the file (again)

  while ((read = getline(&line, &len, fp)) != -1){
    if (strstr(line,"PRAGMA")){
      //We will skip pragma lines
    } else if (strlen(line)>1){
      _quil_add_gate(line,circ,time);
      time = time + 1.0;
    }
  }

  fclose(fp);
  if (line) free(line);
  return;
}

void _quil_get_angle_pi(char angle_pi[32],PetscReal *angle){
  char numerator[32],denominator[32];
  int found_denom=1,i_n,i_d,i,denom,numer,factor;
  //Gate was printed with 'pi' or 'pi/2'
  angle_pi[strlen(angle_pi)-1] = 0;
  factor = 1;
  //Search for digits
  i_d = 0;   i_n = 0;
  denom = 1; numer = 1;
  for (i=0;angle_pi[i] != '\0'; i++){
    if (angle_pi[i]=='-'){
      //Found negative sign
      factor = -1;
    } else if (angle_pi[i]=="/") {
      // Found division; all digits after this
      // belong to the denominator
      found_denom = 1;
    } if (isdigit(angle_pi[i])){
      if (found_denom==0){
        //Numerator
        numerator[i_n] = angle_pi[i];
        i_n = i_n + 1;
      } else {
        //Denominator
        denominator[i_n] = angle_pi[i];
        i_d = i_d + 1;
      }
    }
  }
  if (i_d>0){
    denom = atoi(denominator);
  }
  if (i_n>0){
    numer = atoi(numerator);
  }
  *angle = factor * numer * PETSC_PI / denom;
  return;
}

void _quil_add_gate(char *line,circuit *circ,PetscReal time){
  char *token=NULL,*ptr=NULL;
  const char s[2] = " ";
  char angle_pi[32];
  int qubit1=-1,qubit2=-1;
  PetscReal angle;
  gate_type my_gate_type;
  // Split string on " " to separate the gate and the qubits
  while (token=strsep(&line," ")) {
    size_t i, j;
    for (i=0, j=0; token[j]=token[i]; j+=!isspace(token[i++]));
    //FIXME: Not exhaustive
    if (isdigit(token[0])){
      // check qubit numbers
      if (qubit1<0){
        qubit1 = atoi(token);
      } else {
        qubit2 = atoi(token);
      }
    } else {
      // gate types
      if (strcmp(token,"CNOT")==0){
        my_gate_type = CNOT;
      } else if (strcmp(token,"I")==0){
        my_gate_type = EYE;
      } else if (strcmp(token,"CZ")==0){
        my_gate_type = CZ;;
      } else if (strcmp(token,"H")==0) {
        my_gate_type = HADAMARD;
      } else if (strcmp(token,"Z")==0) {
        my_gate_type = SIGMAZ;
      } else if (strcmp(token,"X")==0) {
        my_gate_type = SIGMAX;
      } else if (strcmp(token,"Y")==0) {
        my_gate_type = SIGMAY;
      } else if (strstr(token,"RX")) {
        my_gate_type = RX;
        if (strstr(token,"pi")){
          sscanf(token,"RX(%s)",&angle_pi);
          _quil_get_angle_pi(angle_pi,&angle);
        } else {
          //Not sure if quil prints like this
          sscanf(token,"RX(%lf)",&angle);
        }
      } else if (strstr(token,"RY")) {
        my_gate_type = RY;
        if (strstr(token,"pi")){
          sscanf(token,"RY(%s)",&angle_pi);
          _quil_get_angle_pi(angle_pi,&angle);
        } else {
          //Not sure if quil prints like this
          sscanf(token,"RY(%lf)",&angle);
        }
      } else if (strstr(token,"RZ")) {
        if (strstr(token,"pi")){
          sscanf(token,"RZ(%s)",&angle_pi);
          _quil_get_angle_pi(angle_pi,&angle);
        } else {
          //Not sure if quil prints like this
          sscanf(token,"RZ(%lf)",&angle);
        }
        my_gate_type = RZ;
      }
    }
  }

  if (qubit2>0){
    //Multiqubit gate
    add_gate_to_circuit(circ,time,my_gate_type,qubit1,qubit2);
  } else {
    //Single qubit gate
    if (my_gate_type==6||my_gate_type==7||my_gate_type==8){
      //Rotation gate
      add_gate_to_circuit(circ,time,my_gate_type,qubit1,angle);
    } else {
      add_gate_to_circuit(circ,time,my_gate_type,qubit1);
    }
  }
  return;
}

void projectq_qasm_read(char filename[],PetscInt *num_qubits,circuit *circ){
  FILE *fp;
  int ch = 0,lines=0,finished_allocate=0;
  char *line = NULL;
  size_t len = 0;
  ssize_t read;
  PetscReal time=1.0;

  fp = fopen(filename,"r");

  if (fp == NULL){
    if (nid==0){
      printf("ERROR! File not found in projectq_qasm_read!\n");
    }
  }

  //Count number of lines
  while(!feof(fp)){
    ch = fgetc(fp);
    if(ch == '\n')
      {
        lines++;
      }
  }
  //Rewind file
  rewind(fp);

  *num_qubits = 0;
  while ((read = getline(&line, &len, fp)) != -1){
    if (!finished_allocate){
      if (strstr(line,"Allocate")){
        //Get number of qubits by reading number of lines with 'Allocate'
        *num_qubits = *num_qubits+1;
      } else {
        /*Subtract off the allocate lines to get the total number
         * of gates.
         * Factor of 2 because the end has DEALLOCATE statements
         * WARNING! Comments will be included, too - not too big of a problem,
         * just will overallocate a bit */
        lines = lines - 2*(*num_qubits);
        //Allocate the circuit
        create_circuit(circ,lines);
        finished_allocate = 1;
        //Add the first gate to list
        _projectq_qasm_add_gate(line,circ,time);
        time = time + 1.0;
        //        gate_list[0]
      }
    } else if (strstr(line,"Deallocate")){
      //Breakout since we've finished going through the circuit
      break;
    } else {
      _projectq_qasm_add_gate(line,circ,time);
      time = time + 1.0;
    }

  }

  fclose(fp);
  if (line) free(line);
  return;
}

void _projectq_qasm_add_gate(char *line,circuit *circ,PetscReal time){
  char *token=NULL;
  int qubit1,qubit2;
  PetscReal angle;
  gate_type my_gate_type;
  // Split string on | to separate gate type and qubits
  while (token=strsep(&line,"|")) {
    //Strip whitespace
    size_t i, j;
    for (i=0, j=0; token[j]=token[i]; j+=!isspace(token[i++]));
    //Do direct strcmp for some, strstr for others
    //FIXME: Not exhaustive
    if (strstr(token,"Qureg")){
      // qubit numbers
      if (my_gate_type<0){
        //Multiqubit gate
        sscanf(token,"(Qureg[%d],Qureg[%d])",&qubit1,&qubit2);
        add_gate_to_circuit(circ,time,my_gate_type,qubit1,qubit2);
      } else {
        //Single qubit gate
        sscanf(token,"Qureg[%d]",&qubit1);
        if (my_gate_type==6||my_gate_type==7||my_gate_type==8){
          //Rotation gate
          add_gate_to_circuit(circ,time,my_gate_type,qubit1,angle);
        } else {
          add_gate_to_circuit(circ,time,my_gate_type,qubit1);
        }
      }
    } else {
      // gate types
      if (strcmp(token,"CX")==0){
        my_gate_type = CNOT;
      } else if (strcmp(token,"H")==0) {
        my_gate_type = HADAMARD;
      } else if (strcmp(token,"Z")==0) {
        my_gate_type = SIGMAZ;
      } else if (strcmp(token,"X")==0) {
        my_gate_type = SIGMAX;
      } else if (strcmp(token,"Y")==0) {
        my_gate_type = SIGMAY;
      } else if (strstr(token,"Rx")) {
        my_gate_type = RX;
        sscanf(token,"Rx(%lf)",&angle);
      } else if (strstr(token,"Ry")) {
        my_gate_type = RY;
        sscanf(token,"Ry(%lf)",&angle);
      } else if (strstr(token,"Rz")) {
        sscanf(token,"Rz(%lf)",&angle);
        my_gate_type = RZ;
      }
    }
  }
}

void projectq_vqe_get_expectation(char filename[],Vec rho,PetscScalar *trace_val){
  FILE *fp;
  char *token=NULL,*token2=NULL;
  char *line = NULL,gate_char;
  size_t len = 0;
  ssize_t read;
  int token_number,num_ops,qubit_number;
  operator ops[100];
  PetscReal scalar_multiply;
  PetscScalar temp_trace_val;
  fp = fopen(filename,"r");
  *trace_val = 0.0;
  if (fp == NULL){
    if (nid==0){
      printf("ERROR! File not found in projectq_vqe_get_expectation!\n");
    }
  }

  while ((read = getline(&line, &len, fp)) != -1){
    token_number = 0;
    while (token=strsep(&line,"[")) {
      if(token_number==0){
        //Strip whitespace
        size_t i, j;
        for (i=0, j=0; token[j]=token[i]; j+=!isspace(token[i++]));
        scalar_multiply = atof(token);
        token_number = 1;
      } else {
        token2=strsep(&token,"]");
        num_ops = 0;
        while (token=strsep(&token2," ")) {
          if (strcmp(token,"")==0){
            *trace_val = *trace_val + scalar_multiply;
          } else {
            //Assume qubit number in file is global system number
            //FIXME: Put logical->physical qubit mapping here
            sscanf(token,"%c%d",&gate_char,&qubit_number);
            if (gate_char=='X'){
              ops[num_ops] = subsystem_list[qubit_number]->sig_x;
            } else if (gate_char=='Y'){
              ops[num_ops] = subsystem_list[qubit_number]->sig_y;
            } else if (gate_char=='Z'){
              ops[num_ops] = subsystem_list[qubit_number]->sig_z;
            }
            num_ops = num_ops + 1;
          }
        }
        if (num_ops!=0){
          //This is a hack where I pass many ops, even though many of them
          //may not exist or be valid; they won't be accessed, at least.
          //Consider passing the array instead, and iterating inside?
          get_expectation_value(rho,&temp_trace_val,num_ops,
                                ops[0],ops[1],ops[2],ops[3],
                                ops[4],ops[5],ops[6],ops[7],
                                ops[8],ops[9],ops[10],ops[11],
                                ops[12],ops[13],ops[14],ops[15],
                                ops[16],ops[17],ops[18],ops[19]);
          temp_trace_val = temp_trace_val * scalar_multiply;
          *trace_val = *trace_val + temp_trace_val;
        }
      }
    }
  }
  return;
}

void projectq_vqe_get_expectation_squared(char filename[],Vec rho,PetscScalar *trace_val){
  FILE *fp;
  char *token=NULL,*token2=NULL;
  char *line = NULL,gate_char;
  size_t len = 0;
  ssize_t read;
  int token_number,num_ops,qubit_number;
  operator ops[100];
  PetscReal scalar_multiply;
  PetscScalar temp_trace_val;
  fp = fopen(filename,"r");
  *trace_val = 0.0;
  if (fp == NULL){
    if (nid==0){
      printf("ERROR! File not found in projectq_vqe_get_expectation!\n");
    }
  }

  while ((read = getline(&line, &len, fp)) != -1){
    token_number = 0;
    while (token=strsep(&line,"[")) {
      if(token_number==0){
        //Strip whitespace
        size_t i, j;
        for (i=0, j=0; token[j]=token[i]; j+=!isspace(token[i++]));
        scalar_multiply = atof(token);
        token_number = 1;
      } else {
        token2=strsep(&token,"]");
        num_ops = 0;
        while (token=strsep(&token2," ")) {
          if (strcmp(token,"")==0){
            *trace_val = *trace_val + scalar_multiply;
          } else {
            //Assume qubit number in file is global system number
            //FIXME: Put logical->physical qubit mapping here
            sscanf(token,"%c%d",&gate_char,&qubit_number);
            if (gate_char=='X'){
              ops[num_ops] = subsystem_list[qubit_number]->sig_x;
            } else if (gate_char=='Y'){
              ops[num_ops] = subsystem_list[qubit_number]->sig_y;
            } else if (gate_char=='Z'){
              ops[num_ops] = subsystem_list[qubit_number]->sig_z;
            }
            num_ops = num_ops + 1;
          }
        }
        if (num_ops!=0){
          //This is a hack where I pass many ops, even though many of them
          //may not exist or be valid; they won't be accessed, at least.
          //Consider passing the array instead, and iterating inside?
          get_expectation_value(rho,&temp_trace_val,num_ops,
                                ops[0],ops[1],ops[2],ops[3],
                                ops[4],ops[5],ops[6],ops[7],
                                ops[8],ops[9],ops[10],ops[11],
                                ops[12],ops[13],ops[14],ops[15],
                                ops[16],ops[17],ops[18],ops[19]);
          temp_trace_val = temp_trace_val * scalar_multiply;
          *trace_val = *trace_val + temp_trace_val;
        }
      }
    }
  }
  return;
}

void projectq_vqe_get_expectation_encoded(char filename[],Vec rho,PetscScalar *trace_val,
                                          PetscInt num_encoders,...){
  FILE *fp;
  char *token=NULL,*token2=NULL;
  char *line = NULL,gate_char;
  size_t len = 0;
  ssize_t read;
  int token_number,num_ops,qubit_number;
  va_list ap;
  operator ops[100];
  PetscInt i;
  PetscReal scalar_multiply;
  PetscScalar temp_trace_val;
  encoded_qubit encoders[50]; //50 is more systems than we will be able to do

  va_start(ap,num_encoders);
  for (i=0;i<num_encoders;i++){
    encoders[i] = va_arg(ap,encoded_qubit);
  }

  fp = fopen(filename,"r");
  *trace_val = 0.0;
  if (fp == NULL){
    if (nid==0){
      printf("ERROR! File not found in projectq_vqe_get_expectation!\n");
    }
  }

  while ((read = getline(&line, &len, fp)) != -1){
    token_number = 0;
    while (token=strsep(&line,"[")) {
      if(token_number==0){
        //Strip whitespace
        size_t i, j;
        for (i=0, j=0; token[j]=token[i]; j+=!isspace(token[i++]));
        scalar_multiply = atof(token);
        token_number = 1;
      } else {
        token2=strsep(&token,"]");
        num_ops = 0;
        while (token=strsep(&token2," ")) {
          if (strcmp(token,"")==0){
            *trace_val = *trace_val + scalar_multiply;
          } else {
            //Assume qubit number in file is global system number
            //FIXME: Put logical->physical qubit mapping here
            sscanf(token,"%c%d",&gate_char,&qubit_number);
            qubit_number = encoders[qubit_number].qubits[0];
            if (gate_char=='X'){
              ops[num_ops] = subsystem_list[qubit_number]->sig_x;
            } else if (gate_char=='Y'){
              ops[num_ops] = subsystem_list[qubit_number]->sig_y;
            } else if (gate_char=='Z'){
              ops[num_ops] = subsystem_list[qubit_number]->sig_z;
            }
            num_ops = num_ops + 1;
          }
        }
        if (num_ops!=0){
          if (num_ops>4){
            if (nid==0) {
              printf("ERROR! vqe_get_expectation only supports 4 ops for now\n");
              exit(0);
            }
          }
          get_expectation_value(rho,&temp_trace_val,num_ops,ops[0],ops[1],ops[2],ops[3]);
          temp_trace_val = temp_trace_val * scalar_multiply;
          *trace_val = *trace_val + temp_trace_val;
        }
      }
    }
  }
  return;
}

ext: c

max_stars_repo_name: bcwu/QuaC

max_issues_repo_name: bcwu/QuaC

max_forks_repo_name: bcwu/QuaC


text:
#include <stdio.h>
#include <stdlib.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_filter.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_vector.h>

int
main(void)
{
  const size_t N = 1000;                               /* length of time series */
  const size_t K = 25;                                 /* window size */
  const double t = 4.0;                                /* number of scale factors for outlier detection */
  gsl_vector *x = gsl_vector_alloc(N);                 /* input vector */
  gsl_vector *y = gsl_vector_alloc(N);                 /* output (filtered) vector */
  gsl_vector *xmedian = gsl_vector_alloc(N);           /* window medians */
  gsl_vector *xsigma = gsl_vector_alloc(N);            /* window scale estimates */
  gsl_vector_int *ioutlier = gsl_vector_int_alloc(N);  /* outlier detected? */
  gsl_filter_impulse_workspace * w = gsl_filter_impulse_alloc(K);
  gsl_rng *r = gsl_rng_alloc(gsl_rng_default);
  size_t noutlier;
  size_t i;

  /* generate input signal */
  for (i = 0; i < N; ++i)
    {
      double xi = 10.0 * sin(2.0 * M_PI * i / (double) N);
      double ei = gsl_ran_gaussian(r, 2.0);
      double u = gsl_rng_uniform(r);
      double outlier = (u < 0.01) ? 15.0*GSL_SIGN(ei) : 0.0;

      gsl_vector_set(x, i, xi + ei + outlier);
    }

  /* apply impulse detection filter */
  gsl_filter_impulse(GSL_FILTER_END_TRUNCATE, GSL_FILTER_SCALE_QN, t, x, y,
                     xmedian, xsigma, &noutlier, ioutlier, w);

  /* print results */
  for (i = 0; i < N; ++i)
    {
      double xi = gsl_vector_get(x, i);
      double yi = gsl_vector_get(y, i);
      double xmedi = gsl_vector_get(xmedian, i);
      double xsigmai = gsl_vector_get(xsigma, i);
      int outlier = gsl_vector_int_get(ioutlier, i);

      printf("%zu %f %f %f %f %d\n",
             i,
             xi,
             yi,
             xmedi + t * xsigmai,
             xmedi - t * xsigmai,
             outlier);
    }

  gsl_vector_free(x);
  gsl_vector_free(y);
  gsl_vector_free(xmedian);
  gsl_vector_free(xsigma);
  gsl_vector_int_free(ioutlier);
  gsl_filter_impulse_free(w);
  gsl_rng_free(r);

  return 0;
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
#include <stdio.h>
#include <stddef.h>
#include <stdbool.h>
#include <stdlib.h>
#include <limits.h>
#include <string.h>

#include <cpuid.h>

#include <perf_counter.h>

#if defined(USE_MKL)
	#include <mkl_cblas.h>
#elif defined(USE_OPENBLAS)
	#include <cblas.h>
#elif defined(USE_BLIS)
	#include <blis/blis.h>
#else
	#error No BLAS library enabled
#endif

extern unsigned long long median(unsigned long long array[], size_t length);
extern void read_memory(const void* memory, size_t length);

enum gemm_type {
	gemm_type_unknown,
	gemm_type_cgemm,
	gemm_type_sgemm,
};

unsigned long long benchmark_gemm(
	const void* memory, size_t cache_size,
	enum gemm_type type,
	size_t m, size_t n, size_t k,
	const float a[], const float b[], float c[],
	size_t max_iterations)
{
	unsigned long long computation_time[max_iterations];
	size_t computation_samples = 0;
	for (size_t iteration = 0; iteration < max_iterations; iteration++) {
		read_memory(memory, cache_size);

		unsigned long long start_time, end_time;
		if (!read_timer(&start_time))
			continue;

		switch (type) {
			case gemm_type_sgemm:
			{
#if defined(USE_MKL) || defined(USE_OPENBLAS)
				cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasTrans,
					m, n, k,
					1.0f, a, k, b, k, 0.0f, c, n);
#elif defined(USE_BLIS)
				float alpha = 1.0f;
				float beta = 0.0f;
				bli_sgemm(BLIS_NO_TRANSPOSE, BLIS_TRANSPOSE,
					m, n, k,
					&alpha,
					(float*) a, k, 1,
					(float*) b, k, 1,
					&beta,
					c, n, 1);
#endif
				break;
			}
			case gemm_type_cgemm:
			{
#if defined(USE_MKL) || defined(USE_OPENBLAS)
				float alpha[2] = { 1.0f, 0.0f };
				float beta[2] = { 0.0f, 0.0f };
				cblas_cgemm(CblasRowMajor, CblasNoTrans, CblasTrans,
					m, n, k,
					alpha, a, k, b, k, beta, c, n);
#elif defined(USE_BLIS)
				scomplex alpha = { 1.0f, 0.0f };
				scomplex beta = { 0.0f, 0.0f };
				bli_cgemm(BLIS_NO_TRANSPOSE, BLIS_TRANSPOSE,
					m, n, k,
					&alpha,
					(scomplex*) a, k, 1,
					(scomplex*) b, k, 1,
					&beta,
					(scomplex*) c, n, 1);
#endif
				break;
			}
		}

		if (!read_timer(&end_time))
			continue;

		computation_time[computation_samples++] = end_time - start_time;
	}

	return median(computation_time, computation_samples);
}

struct options {
	enum gemm_type type;
	size_t m;
	size_t n;
	size_t k;
	size_t iterations;
};

static void print_options_help(const char* program_name) {
	printf(
"%s parameters...\n"
"Required parameters:\n"
"  -g   --gemm               The type of GEMM operation (sgemm, cgemm)\n"
"  -m                        The M dimension\n"
"  -n                        The N dimension\n"
"  -k                        The K dimension\n"
"Optional parameters:\n"
"  -i   --iterations         # iterations (default: 151)\n",
		program_name);
}

static struct options parse_options(int argc, char** argv) {
	struct options options = {
		.type = gemm_type_unknown,
		.k = 0,
		.m = 0,
		.n = 0,
		.iterations = 151
	};
	for (int argi = 1; argi < argc; argi += 1) {
		if ((strcmp(argv[argi], "-g") == 0) || (strcmp(argv[argi], "--gemm") == 0)) {
			if (argi + 1 == argc) {
				fprintf(stderr, "Error: expected k value\n");
				exit(EXIT_FAILURE);
			}
			if (strcmp(argv[argi + 1], "sgemm") == 0) {
				options.type = gemm_type_sgemm;
			} else if (strcmp(argv[argi + 1], "cgemm") == 0) {
				options.type = gemm_type_cgemm;
			} else {
				fprintf(stderr, "Error: invalid value %s for the gemm type: expected \"sgemm\", \"cgemm\"\n", argv[argi + 1]);
				exit(EXIT_FAILURE);
			}
			argi += 1;
		} else if (strcmp(argv[argi], "-k") == 0) {
			if (argi + 1 == argc) {
				fprintf(stderr, "Error: expected K value\n");
				exit(EXIT_FAILURE);
			}
			if (sscanf(argv[argi + 1], "%zu", &options.k) != 1) {
				fprintf(stderr, "Error: can not parse %s as an unsigned integer\n", argv[argi + 1]);
				exit(EXIT_FAILURE);
			}
			if (options.k == 0) {
				fprintf(stderr, "Error: invalid value %s for the K parameter: positive value expected\n", argv[argi + 1]);
				exit(EXIT_FAILURE);
			}
			argi += 1;
		} else if (strcmp(argv[argi], "-m") == 0) {
			if (argi + 1 == argc) {
				fprintf(stderr, "Error: expected M value\n");
				exit(EXIT_FAILURE);
			}
			if (sscanf(argv[argi + 1], "%zu", &options.m) != 1) {
				fprintf(stderr, "Error: can not parse %s as an unsigned integer\n", argv[argi + 1]);
				exit(EXIT_FAILURE);
			}
			if (options.m == 0) {
				fprintf(stderr, "Error: invalid value %s for the M parameter: positive value expected\n", argv[argi + 1]);
				exit(EXIT_FAILURE);
			}
			argi += 1;
		} else if (strcmp(argv[argi], "-n") == 0) {
			if (argi + 1 == argc) {
				fprintf(stderr, "Error: expected N value\n");
				exit(EXIT_FAILURE);
			}
			if (sscanf(argv[argi + 1], "%zu", &options.n) != 1) {
				fprintf(stderr, "Error: can not parse %s as an unsigned integer\n", argv[argi + 1]);
				exit(EXIT_FAILURE);
			}
			if (options.n == 0) {
				fprintf(stderr, "Error: invalid value %s for the N parameter: positive value expected\n", argv[argi + 1]);
				exit(EXIT_FAILURE);
			}
			argi += 1;
		} else if ((strcmp(argv[argi], "--iterations") == 0) || (strcmp(argv[argi], "-i") == 0)) {
			if (argi + 1 == argc) {
				fprintf(stderr, "Error: expected iterations value\n");
				exit(EXIT_FAILURE);
			}
			if (sscanf(argv[argi + 1], "%zu", &options.iterations) != 1) {
				fprintf(stderr, "Error: can not parse %s as an unsigned integer\n", argv[argi + 1]);
				exit(EXIT_FAILURE);
			}
			if (options.iterations == 0) {
				fprintf(stderr, "Error: invalid value %s for the number of iterations: positive value expected\n", argv[argi + 1]);
				exit(EXIT_FAILURE);
			}
			argi += 1;
		} else if ((strcmp(argv[argi], "--help") == 0) || (strcmp(argv[argi], "-h") == 0)) {
			print_options_help(argv[0]);
			exit(EXIT_SUCCESS);
		} else {
			fprintf(stderr, "Error: unknown argument '%s'\n", argv[argi]);
			print_options_help(argv[0]);
			exit(EXIT_FAILURE);
		}
	}
	if (options.type == gemm_type_unknown) {
		fprintf(stderr, "Error: gemm kernel type is not specified\n");
		print_options_help(argv[0]);
		exit(EXIT_FAILURE);
	}
	if (options.m == 0) {
		fprintf(stderr, "Error: M is not specified\n");
		print_options_help(argv[0]);
		exit(EXIT_FAILURE);
	}
	if (options.n == 0) {
		fprintf(stderr, "Error: N is not specified\n");
		print_options_help(argv[0]);
		exit(EXIT_FAILURE);
	}
	if (options.k == 0) {
		fprintf(stderr, "Error: K is not specified\n");
		print_options_help(argv[0]);
		exit(EXIT_FAILURE);
	}
	return options;
}

int main(int argc, char** argv) {
	const struct options options = parse_options(argc, argv);
	size_t components = 0;
	switch (options.type) {
		case gemm_type_sgemm:
			components = 1;
			break;
		case gemm_type_cgemm:
			components = 2;
			break;
		default:
			__builtin_unreachable();
	}

	const size_t cache_size = 128 * 1024 * 1024;
	void* memory = valloc(cache_size);

#ifdef USE_BLIS
	err_t blis_status = bli_init();
	if (blis_status != BLIS_SUCCESS) {
		fprintf(stderr, "BLIS initialization failed: error code %d\n", blis_status);
		exit(EXIT_FAILURE);
	}
#endif

	float* c = (float*) valloc(options.m * options.n * components * sizeof(float));
	float* a = (float*) valloc(options.m * options.k * components * sizeof(float));
	float* b = (float*) valloc(options.k * options.n * components * sizeof(float));

	memset(c, 0, options.m * options.n * components * sizeof(float));
	memset(a, 0, options.m * options.k * components * sizeof(float));
	memset(b, 0, options.k * options.n * components * sizeof(float));

	{
		printf("Iterations: %zu\n", options.iterations);

		const unsigned long long gemm_nanoseconds = benchmark_gemm(
			memory, cache_size,
			options.type,
			options.m, options.n, options.k,
			a, b, c,
			options.iterations);
		const double gemm_gflops = 2.0 * components * components * options.m * options.n * options.k / ((double) gemm_nanoseconds);

		printf("Time: %5.3lf ms\n", gemm_nanoseconds * 1.0e-6);
		printf("Performance: %5.3lf GFLOPs/s\n", gemm_gflops);
	}
}


text:
#include <stdio.h>
#include <stdlib.h>
#include "cs.h"
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include "adkCSparse.h"
#include "adkGSL.h"

//gslMatrixToCSsparse -- returns new cs_sparse object. compressed
struct cs_di_sparse *gslMatrixToCSsparse(gsl_matrix *gm, int numNonZero){
	int i, j;
	struct cs_di_sparse *triplet, *compmat;
	double tmp;

	//now allocate new cs_sparse obj
	triplet = cs_spalloc(gm->size1, gm->size2, numNonZero, 1, 1);
	for(i=0;i<gm->size1;i++){
		for(j=0;j<gm->size2;j++){
			tmp = gsl_matrix_get(gm,i,j);
			if(tmp != 0.0) cs_entry(triplet,i,j,tmp);
		}
	}
	compmat = cs_compress(triplet);
	cs_spfree(triplet);
	return(compmat);
}

//gslMatrixToCSsparse -- returns new cs_sparse object. compressed
struct cs_di_sparse *arraysToCSsparse(int rank, int numNonZero, int *dim1, int *dim2,double *x){
	int i;
	struct cs_di_sparse *triplet, *compmat;

	//now allocate new cs_sparse obj
	triplet = cs_spalloc(rank, rank, numNonZero, 1, 1);
	for(i=0;i<numNonZero;i++) cs_entry(triplet,dim1[i],dim2[i],x[i]);
	compmat = cs_compress(triplet);
	cs_spfree(triplet);
	return(compmat);
}

//cs_max -- returns the maximum value of a cs_sparse object
double cs_max(struct cs_di_sparse *mat){
	int j,p;
	double max=0.0;
	double val;

	if (mat->nz < 0){
		for (j = 0 ; j < mat->n ; j++){
			for(p = mat->p[j] ; p < mat->p[j+1] ; p++){
				val = mat->x[p];
				if(val >max){
					max=val;
				}
			}
		}
	}
	else{ //below is untested
		for(j=0;j<mat->nzmax;j++){
			if(mat->x[j] > max){
				max = mat->x[j];
			}
		}
	}
	return(max);
}

//cs_scale -- scales each value by some constant
void cs_scale(struct cs_di_sparse *mat, double c){
	int i;
	for(i=0;i<mat->nzmax;i++){
		mat->x[i] *= c;
	}
}

//cs_scale -- scales each value by some constant
void cs_print_adk(struct cs_di_sparse *mat){
	int i,j,col,row,p;
	double val;
	gsl_matrix *gm;
	gm = gsl_matrix_alloc(mat->m,mat->n);
	gsl_matrix_set_zero(gm);

	if (mat->nz < 0){
		for (j = 0 ; j < mat->n ; j++){
			col=j;
			for(p = mat->p[j] ; p < mat->p[j+1] ; p++){
				row = mat->i [p];
				val = mat->x[p];
				gsl_matrix_set(gm,row,col,val);
			}
		}
	}
	else{

		for(i=0;i<mat->nz;i++){
			gsl_matrix_set(gm,mat->p[i],mat->i[i],mat->x[i]);
		}
	}
	gsl_matrix_prettyPrint(gm);
	gsl_matrix_free(gm);

}

max_stars_repo_name: dortegadelv/IMaDNA

max_issues_repo_name: dortegadelv/IMaDNA


text:
//
// Created by Tobias on 10/2/2020.
//

#ifndef SD3D_ORBITCAMERACONTROLLER_H
#define SD3D_ORBITCAMERACONTROLLER_H

#include <GLFW/glfw3.h>
#include <gsl-lite/gsl-lite.hpp>

#pragma warning(push, 0)

#include <glm/glm.hpp>
#include <glm/gtc/matrix_transform.hpp>

#pragma warning(pop)

#include "../graphics/Camera.h"

class OrbitCameraController {
public:
	struct OrbitCamSettings {
		glm::vec3 target{};

		float dist{};
		float minDist{};
		float maxDist{};
		float scrollSpeed{};

		float speed{};
		float deceleration{};

		OrbitCamSettings(const glm::vec3 &target, float dist, float speed,
						 float deceleration = 0.05f) :
			OrbitCamSettings{target, dist,  dist,        dist,
							 0.0f,   speed, deceleration} {}

		OrbitCamSettings(const glm::vec3 &target, float dist, float minDist,
						 float maxDist, float scrollSpeed, float speed,
						 float deceleration = 0.05f) :
			target{target},
			dist{dist},
			minDist{minDist},
			maxDist{maxDist},
			scrollSpeed{scrollSpeed},
			speed{speed},
			deceleration{deceleration} {}
	};

private:
	// split mouse vector in 2 directions
	// x -> theta
	// y -> phi

	glm::vec3 m_target{};

	float m_dist{};
	float m_minDist{};
	float m_maxDist{};
	float m_scrollSpeed{};

	float m_yaw{-90.0f};
	float m_pitch{};

	glm::vec2 m_moveVec{};

	float m_deceleration{0.05f};
	float m_speed{};

	Camera m_cam;

	// for settings window
private:
	float m_fovTemp;

	float m_fovSaved;
	OrbitCamSettings m_settingsSaved;

public:
	OrbitCameraController(Camera &&cam, const OrbitCamSettings &settings) :
		m_target{settings.target},
		m_dist{settings.dist},
		m_minDist{settings.minDist},
		m_maxDist{settings.maxDist},
		m_scrollSpeed{settings.scrollSpeed},
		m_deceleration{settings.deceleration},
		m_speed{settings.speed},
		m_cam{cam},
		m_fovTemp{gsl::narrow_cast<float>(m_cam.fov())},
		m_fovSaved{m_fovTemp},
		m_settingsSaved{settings} {}

	void update(float deltaTime);

	[[nodiscard]] const Camera &ccam() const { return m_cam; }
	[[nodiscard]] Camera &cam() { return m_cam; };

	void rotate(double dX, double dY);

	void zoom(double dScroll);

	[[nodiscard]] float min_dist() const { return m_minDist; }

	[[nodiscard]] float max_dist() const { return m_maxDist; }

	void set_min_dist(float val) { m_minDist = val; }

	void set_max_dist(float val) { m_maxDist = val; }

	void settings_gui(bool &show, glm::vec3 &clearCol);

	void reset_settings() {
		m_target = m_settingsSaved.target;
		m_dist = m_settingsSaved.dist;
		m_minDist = m_settingsSaved.minDist;
		m_maxDist = m_settingsSaved.maxDist;
		m_scrollSpeed = m_settingsSaved.scrollSpeed;
		m_deceleration = m_settingsSaved.deceleration;
		m_speed = m_settingsSaved.speed;

		m_fovTemp = m_fovSaved;
		m_cam.set_fov(m_fovTemp);
	}
};

#endif // SD3D_ORBITCAMERACONTROLLER_H


text:
//  © Shahram Talei @ 2020 The University of Alabama
//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.
//You should have received a copy of the GNU General Public License
//along with this program.  If not, see <http://www.gnu.org/licenses/>.
// This code is distributed as is and there is no warranty or technical support
// Reading stored data in a sage file
// How to run:
//first run:
//$make
//Then you can run"
// $./ExportForTracking > AllGals838.ascii
// If you have an error you may need to give permission
// $ chmod +x ExportForTracking
//
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <unistd.h>
//#include <gsl/gsl_rng.h>
//#include <ctype.h>

////////////////////////////////////////////////////////////
//// Definitions and variables
char SageDir[500];
struct Path_Names
{
    char paths[100];
}*SageFilesPath;

int SageFilesCount;
int NumGalaxies;

// This structune holds all of the information from the Sage files
struct SageGalaxies
{
  int   Type;
  int   FileNr;
  long long   GalaxyIndex;
  int   HaloIndex;
  int   FOFHaloIndex; //SubhaloIndex;
  int   TreeIndex;

  int   SnapNum;
  int   CentralGal; //GalIndex;
  float CentralMvir;

   //properties of subhalo at the last time this galaxy was a central galaaxy
  float Pos[3];
  float Vel[3];
  float Spin[3];
  int   Len;
  float Mvir;
  float Rvir;
  float Vvir;
  float Vmax;
  float VelDisp;

   //baryonic reservoirs
  float ColdGas;
  float StellarMass;
  float BulgeMass;
  float HotGas;
  float EjectedMass;
  float BlackHoleMass;
  float ICS;
  //metals
 float MetalsColdGas;
 float MetalsStellarMass;
 float MetalsBulgeMass;
 float MetalsHotGas;
 float MetalsEjectedMass;
 float MetalsICS;

  //misc
 float Sfr;
 float SfrBulge;
 float SfrICS;
 float DiskScaleRadius;
 float Cooling;
 float Heating;
 float LastMajorMerger;
 float OutflowRate;

 float infallMvir;  //infall properties
 float infallVvir;
 float infallVmax;
 float r_heat;
 float Stars;

}*SageOutput;
//typedef struct SageGalaxies

////////////////////////////////////////////////////////////
//// Functions

int CountSageFiles(int snap)
{
         //int i=0;
         int n=0;
         char line[10000];
         char file1[10000];
         FILE *fd;
         sprintf(file1, "%s/sagefile_%03d.txt",SageDir, snap);
         fd = fopen(file1, "r");
         while(fgets(line, sizeof(line), fd)!=NULL)
         {
            n++;// 1;
         }
         fclose(fd);
         return n;
}
void ReadSageFNames(int snap,struct Path_Names *SageFile)
{
        //printf("Reading sage file-name(s) for snap:%d\n",snap);
    int i=0;
    char line[10000];
    char file1[10000];
    char ABSpath[10000];//2
    char model[500];
    FILE *fd;
    sprintf(file1, "%s/sagefile_%03d.txt", SageDir, snap);
    fd = fopen(file1, "r");
    while(EOF != fscanf(fd, "%s%*[^\n]", line))
    {
       strcpy(ABSpath,line);//2
       strcpy(model,strrchr(ABSpath,'/'));
       //strcpy(SageFilesPath[i].paths, line);//1
       if (model[0] == '/')
         memmove(model, model+1, strlen(model));
       sprintf(SageFile[i].paths,"%s/%s",SageDir,model);
       i++;
    }
   fclose(fd);
   return;
 }

 int ReadSageHeader(int FileCount,struct Path_Names *SageFile)
 {
      int Ntrees;
          int NtotGal;
          int totmal;
          int i;
          FILE *fd;
          char file1[1000];
          totmal = 0;
          for(i=0; i<FileCount; i++)
             {
                  sprintf(file1, "%s", SageFile[i].paths);
                  fd = fopen(file1, "rb");
                  if(NULL == fd)
                  {
                     printf("Cannot open sage file");
                     return(-1);
                  }
                  fread(&Ntrees, sizeof(int), 1, fd);
                  fread(&NtotGal, sizeof(int), 1, fd);
                  totmal = totmal + NtotGal;
                  fclose(fd);
             }
 return totmal;
 }

void LoadSageFiles(int snap)
{
        SageFilesCount = CountSageFiles(snap);
        SageFilesPath = (struct Path_Names*)malloc(SageFilesCount * sizeof(struct Path_Names));
        ReadSageFNames(snap,SageFilesPath);
        return;
    }

    void ReadSageModel(int FileCount,struct Path_Names *SageFile,struct SageGalaxies *Output)
    {
    int Ntrees;
               int NtotGal;
               int *galpertree;
               char file1[1000];
               int offset;
               int i;
               FILE *fd;
               offset = 0;
               for(i=0; i<FileCount; i++)
               {
                  sprintf(file1, "%s", SageFile[i].paths);
                  fd = fopen(file1, "rb");
                        if(NULL == fd)
                        {
                           printf("Cannot open sage for reading");
                           return;
                        }
                        fread(&Ntrees, sizeof(int), 1, fd);
                        fread(&NtotGal, sizeof(int), 1, fd);
                        galpertree = (int*)malloc(Ntrees*sizeof(int));
                        fread(&galpertree[0], sizeof(int), Ntrees, fd);
                        fread(&Output[offset], sizeof(struct SageGalaxies), NtotGal, fd);
                        fclose(fd);
                        offset = offset + NtotGal;
                }
                        free(galpertree);
                        return;
    }

void ReadSage(int snap)
{

//printf("Extracting information from the sage file.\n");
LoadSageFiles(snap);
fflush(stdout);
NumGalaxies = ReadSageHeader(SageFilesCount,SageFilesPath);

if(( SageOutput = (struct SageGalaxies*)malloc(NumGalaxies * sizeof(struct SageGalaxies))) == NULL)
    {
        printf("Failed to allocate for target sage...");
        return;
    }
ReadSageModel(SageFilesCount, SageFilesPath, SageOutput);
return;
}

void PrintGalaxyInfo(struct SageGalaxies *Output,int index)
{
  int i;
  i=index;
printf("~~~~~~~~~~~~~~~~~~~~~~∮∮∮∮∮∮∮∮∮∮∮∮~~~~~~~~~~~~~~~~~~~~~~\n");
printf("Galaxy Information- galaxy:%d, Type:%d, FileNr:%d\n",i,Output[i].Type,Output[i].FileNr);
printf("GalIndex:%lld,HaloIndex:%d,TreeIndex:%d,CentralGal:%d,CentralMvir:%f\n",Output[i].GalaxyIndex,Output[i].HaloIndex,
Output[i].TreeIndex,Output[i].CentralGal,Output[i].CentralMvir);
printf("Pos(x,y,z):(%f,%f,%f)\nVel(vx,vy,vz):(%f,%f,%f)\nSpin:(%f,%f,%f),Len:%d\n",Output[i].Pos[0],Output[i].Pos[1],
Output[i].Pos[2],Output[i].Vel[0],Output[i].Vel[1],Output[i].Vel[2],Output[i].Spin[0],Output[i].Spin[1],Output[i].Spin[2],Output[i].Len);

printf("Mvir:%f,Rvir:%f,Vvir:%f,Vmax:%f,VelDisp:%f\n",Output[i].Mvir,Output[i].Rvir,Output[i].Vvir,Output[i].Vmax,Output[i].VelDisp);

printf("ColdGas:%f,StellarMass:%f,BulgeMass:%f,HotGas:%f\nEjectedMass:%f,BHMass:%f,ICS:%f\n",Output[i].ColdGas,Output[i].StellarMass,
Output[i].BulgeMass,Output[i].HotGas,Output[i].EjectedMass,Output[i].BlackHoleMass,Output[i].ICS);

printf("\nMetals:\nColdGas:%f,StellarMass:%f,BulgeMass:%f,HotGas:%f\nEjectedMass:%f,ICS:%f\n\n",Output[i].MetalsColdGas,
Output[i].MetalsStellarMass,Output[i].MetalsBulgeMass,Output[i].MetalsHotGas,Output[i].MetalsEjectedMass,Output[i].MetalsICS);
printf("Sfr:%f,SfrBulge:%f,SfrICS:%f\nRd:%f,Cooling:%f,Heating:%f\nLastMajorMerger:%f,OutflowRate:%f\n",Output[i].Sfr,
Output[i].SfrBulge,Output[i].SfrICS,Output[i].DiskScaleRadius,Output[i].Cooling,Output[i].Heating,Output[i].LastMajorMerger,Output[i].OutflowRate);

printf("infallMvir:%f,infallVvir:%f,infallVmax:%f,r_heat:%f\n",Output[i].infallMvir,Output[i].infallVvir,Output[i].infallVmax,Output[i].r_heat);
//printf("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n");
//printf("∮∮∮∮∮∮ ∮ ∮ ∮ ∮ ∮ ∮ ∮ ∮ ∮ ∮ ∮ ∮ ∮ ∮ ∮ ∮ ∮ ∮ ∮ ∮ ∮ ∮ ∮ ∮ ∮ ∮ ∮ ∮ ∮∮∮∮∮∮\n");
printf("~~~~~~~~~~~~~~~~~~~~~~∮∮∮∮∮∮∮∮∮∮∮∮~~~~~~~~~~~~~~~~~~~~~~\n");

return;

}
void ExportGalaxy(int snap,struct SageGalaxies *Output, int id)
{
  //snap=161;
  printf("%g,%g,%g,%g,%g,%d,%g,%d,%d,%d,%g\n",Output[id].Pos[0], Output[id].Pos[1],Output[id].Pos[2],Output[id].Mvir*(1.0e10),Output[id].Rvir,snap,Output[id].StellarMass*1.0e10,Output[id].GalaxyIndex,Output[id].HaloIndex,Output[id].CentralGal, Output[id].MetalsStellarMass*(1.0e10));//,Output[id].Stars, Output[id].DiskScaleRadius); //should remove stars for current tagging

}

int main(int argc, char * argv[])
{
//sprintf(SageDir,"/home/shahram/Desktop/Research/3_Tagging/ReadSage/sage_out");
//838: 49-264
//717: 36-264 -> new run error on 149
//16406:40-264
//19880: 41-264
//32251: 46-264
//High5M:
//  /media/shahram/SD/Sample100Mpc/32251/High5M/r2/sage_out
//  /media/shahram/SD/Sample100Mpc/m12f/sage_out

sprintf(SageDir,"/media/shahram/SD/Sample100Mpc/m12i/sage_out"); //49-264
int i,s, LastSnap,FirstSnap;
LastSnap= atoi(argv[2]);
FirstSnap=atoi(argv[1]);
for(s=FirstSnap;s<=LastSnap;s++)
{
  ReadSage(s);
  //printf("snap:%d\n",s);
for(i=0;i<NumGalaxies;i++)
  ExportGalaxy(s,SageOutput,i);
}
  return 0;
}

max_stars_repo_name: stalei/TagMix

max_issues_repo_name: stalei/TagMix

max_forks_repo_name: stalei/TagMix


text:
#pragma once

#include <math.h>
#include <setjmp.h>
#include <stdarg.h>
#include <stddef.h>
#include <stdint.h>

#include <cmocka.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>

/* A very small number */
#define EPSILON 1.0e-16

#ifdef DOXYGEN
/**
 * @brief Assert that the two given doubles are equal to within tolerance.
 *
 * The function prints an error message to standard error and terminates the
 * test by calling fail() if the doubles are not equal to within tolerance.
 *
 * @param[in]  a    The first double to compare.
 *
 * @param[in]  b    The double to compare against the first one.
 *
 * @param[in]  eps  The double tolerance value.
 */
void assert_double_equal(double a, double b, double eps);
#else
#define assert_double_equal(a, b, eps) \
    _assert_double_equal(a, b, eps, __FILE__, __LINE__)
#endif

#ifdef DOXYGEN
/**
 * @brief Assert that the two given doubles are not equal to within tolerance.
 *
 * The function prints an error message to standard error and terminates the
 * test by calling fail() if the doubles are equal to within tolerance.
 *
 * @param[in]  a    The first double to compare.
 *
 * @param[in]  b    The double to compare against the first one.
 *
 * @param[in]  eps  The double tolerance value.
 *
 * @see assert_double_equal()
 */
void assert_double_not_equal(double a, double b, double eps);
#else
#define assert_double_not_equal(a, b, eps) \
    _assert_double_not_equal(a, b, eps, __FILE__, __LINE__)
#endif

#ifdef DOXYGEN
/**
 * @brief Assert that the two given GSL vectors are equal.
 *
 * The function prints an error message to standard error and terminates the
 * test by calling fail() if the vectors are not equal to within tolerance.
 *
 * @param[in]  a    The first GSL vector to compare.
 *
 * @param[in]  b    The GSL vector to compare against the first one.
 *
 * @param[in]  eps  The double tolerance value.
 *
 * @see assert_gsl_vector_equal()
 */
void assert_gsl_vector_equal(gsl_vector *a, gsl_vector *b, double eps);
#else
#define assert_gsl_vector_equal(a, b, eps) \
    _assert_gsl_vector_equal(a, b, eps, __FILE__, __LINE__)
#endif

#ifdef DOXYGEN
/**
 * @brief Assert that the two given GSL matrices are equal.
 *
 * The function prints an error message to standard error and terminates the
 * test by calling fail() if the matrices are not equal to within tolerance.
 *
 * @param[in]  a    The first GSL matrix to compare.
 *
 * @param[in]  b    The GSL matrix to compare against the first one.
 *
 * @param[in]  eps  The double tolerance value.
 *
 * @see assert_gsl_matrix_equal()
 */
void assert_gsl_matrix_equal(gsl_matrix *a, gsl_matrix *b, double eps);
#else
#define assert_gsl_matrix_equal(a, b, eps) \
    _assert_gsl_matrix_equal(a, b, eps, __FILE__, __LINE__)
#endif

void _assert_double_equal(const double a, const double b, const double eps,
                          const char* const file, const int line);

void _assert_double_not_equal(const double a, const double b, const double eps,
                              const char* const file, const int line);

void _assert_gsl_vector_equal(const gsl_vector *a, const gsl_vector *b, double eps,
			      const char* const file, const int line);

void _assert_gsl_matrix_equal(const gsl_matrix *a, const gsl_matrix *b, double eps,
			      const char* const file, const int line);

max_stars_repo_name: mygulamali/neural-net

max_issues_repo_name: mygulamali/neural-net

max_forks_repo_name: mygulamali/neural-net