text:
/*
 */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "aXe_grism.h"
#include "aXe_utils.h"
#include "inout_aper.h"
#include "aper_conf.h"
#include <gsl/gsl_multifit_nlin.h>
#include "trfit_utils.h"

#define AXE_IMAGE_PATH "AXE_IMAGE_PATH"
#define AXE_OUTPUT_PATH "AXE_OUTPUT_PATH"
#define AXE_CONFIG_PATH "AXE_CONFIG_PATH"

int
main(int argc, char *argv[])
{

  char           *opt;

  char            grism_image[MAXCHAR];
  char            grism_image_path[MAXCHAR];

  char            conf_file[MAXCHAR];
  char            conf_file_path[MAXCHAR];

  char            aper_file[MAXCHAR];
  char            aper_file_path[MAXCHAR];

  int             i, j;
  int             num;

  object        **oblist;
  observation    *obs;

  aperture_conf  *conf;

  double exptime;

  gsl_vector *fit_result;

  if ((argc < 1) || (opt = get_online_option("help", argc, argv))) {
    fprintf(stdout,
	    "aXe_TFIT Version %s:\n"
	    "\n", RELEASE);
    exit(1);
  }
  fprintf(stdout, "aXe_TFIT: Starting...\n");

  // Get the data file name
  strcpy(grism_image, argv[1]);

  // Get the configuration file name
  strcpy(conf_file, argv[2]);
  build_path(AXE_CONFIG_PATH, conf_file, conf_file_path);

  // Read the configuration file
  conf = get_aperture_descriptor(conf_file_path);

  // Determine where the various extensions are in the FITS file
  build_path(AXE_IMAGE_PATH, grism_image, grism_image_path);
  get_extension_numbers(grism_image_path, conf, conf->optkey1, conf->optval1);

  // Get or set up the name of the output Aperture File
  if ((opt = get_online_option ("in_AF", argc, argv)))
    {
      /* get it */
      strcpy (aper_file, opt);
      strcpy (aper_file_path, opt);
    }
  else {
    // Build aperture file name
    replace_file_extension (grism_image, aper_file, ".fits",
			    ".OAF", conf->science_numext);
    build_path (AXE_OUTPUT_PATH, aper_file, aper_file_path);
  }

  fprintf(stdout,
	  "aXe_TFIT: Input configuration file name:   %s\n",
	  conf_file_path);
  fprintf(stdout,
	  "aXe_TFIT: Input data file name:            %s\n",
	  grism_image_path);
  fprintf(stdout,
	  "aXe_TFIT: SCI extension number:            %d\n",
	  conf->science_numext);
  fprintf(stdout,
	  "aXe_TFIT: ERR extension number:            %d\n",
	  conf->errors_numext);
  fprintf(stdout,
	  "aXe_TFIT: DQ extension number:             %d\n",
	  conf->dq_numext);
  fprintf(stdout,
	  "aXe_TFIT: DQ mask:                         %d\n",
	  conf->dqmask);
  fprintf(stdout,
	  "aXe_TFIT: Input aperture file name:        %s\n",
	  aper_file_path);
  fprintf(stdout, "\n\n");

  /* Loading the observation data */
  fprintf(stdout, "aXe_TFIT: ");

  //
  // try to get the descriptor 'exptime' from the 'sci'-extension
  //
  exptime = (double)get_float_from_keyword(grism_image_path,
					   conf->science_numext,
					   conf->exptimekey);
  if (isnan(exptime))
    exptime = (double)get_float_from_keyword(grism_image_path, 1,
					     conf->exptimekey);
  if (isnan(exptime))
    exptime = 1.0;

  // Load an image from AXE_DATA_PATH
  obs = load_image_t(grism_image_path, conf->science_numext,
		     conf->errors_numext, conf->dq_numext,
		     conf->dqmask, exptime, conf->rdnoise);

  // Loading the object list
  fprintf(stdout, "aXe_TFIT: Loading object list...");
  oblist = file_to_object_list_seq(aper_file_path, obs);
  fprintf(stdout, "%d objects loaded.\n", object_list_size(oblist));

  // initialize the counter
  i = 0;

  // check whether there are apertures
  if (oblist != NULL) {

    // go over all apertures
    while (oblist[i] != NULL) {

      // print the current aperture ID
      fprintf(stdout, "aXe_TFIT: Fitting object ID:%d", oblist[i]->ID);

      // go over all beams
      for (j = 0; j < oblist[i]->nbeams; j++) {

	// skip beam if ignore flag for thisbeam is set
	if (oblist[i]->beams[j].ignore !=0 || oblist[i]->beams[j].ID != 0)
	  {
	    /*fprintf(stdout,", %c ( Ignored )",
		    BEAM(oblist[i]->beams[j].ID));*/
	    continue;
	  }

	// print the current beam ID
	/*fprintf(stdout, ", %c (Fitted)", BEAM(oblist[i]->beams[j].ID));*/

	// fit the trace of the beam, modify the
	// beam definition with the new values
	// oblist[i]->beams[j].refpoint.y = fit_beamtrace(conf, obs, j, oblist[i]->beams[j]);
	fit_result = fit_beamtrace(conf, obs, j, oblist[i]->beams[j]);
	oblist[i]->beams[j].refpoint.x = gsl_vector_get(fit_result,0);
	oblist[i]->beams[j].refpoint.y = gsl_vector_get(fit_result,1);
	oblist[i]->beams[j].width = gsl_vector_get(fit_result,4);

	gsl_vector_free(fit_result);
      }

      // give feedback on the progress,
      // enhance the counter
      fprintf(stdout, " Done.\n");
      i++;
    }
  }

  // write the modified beam definitions
  // to an OAF file
  fprintf (stdout, "aXe_TFIT: Writing new aperture file...");fflush(stdout);
  num = object_list_to_file (oblist, aper_file_path, 1);
  fprintf (stdout, "%d beams written.\n", num);

  // release the memory
  free_observation(obs);
  if (oblist != NULL){
	free_oblist(oblist);}

  // print a last status message and exit
  fprintf(stdout, "aXe_TFIT: Done...\n");
  exit(0);
}

max_stars_repo_name: sosey/pyaxe

max_issues_repo_name: sosey/pyaxe

max_forks_repo_name: sosey/pyaxe


text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "compearth.h"
#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_lapacke.h>
#ifdef __clang__
#pragma clang diagnostic pop
#endif
#else
#include <lapacke.h>
#endif

/*!
 * @brief Convert from moment tnesor to scalar seismic moment.
 *
 * @param[in] nm     Number of moment tensors.
 * @param[in] im0    =1 for Silver and Jordan (1982) formula. \n
 *                   =2 for GCMT formula. \n
 *                   =3 for old `Caltech Way.'
 * @param[in] M      [nm x 6] array of moment tensors packed
 *                   {M11, M22, M33, M12, M13, M23} with leading
 *                   dimension 6.  The Mij element should be of units
 *                   of Newton-meters (although it should not affect
 *                   the function here).
 *
 * @param[out] M0    Scalar moments with units of Newton-meters.  This
 *                   is an array of dimension [nm].
 *
 * @result 0 indicates success.
 *
 * @author Carl Tape and translated to C by Ben Baker
 *
 * @copyright MIT
 *
 */
int compearth_CMT2m0(const int nm, const int im0,
                     const double *__restrict__ M,
                     double *__restrict__ M0)
{
    double Mt9[9], lams[3], work[100], M11, M22, M33, M12, M13, M23;
    const double sqrt2i = 1.0/sqrt(2.0);
    const int lwork = 100;
    int i, ierr, info;
    ierr = 0;
    if (im0 == 2 || im0 == 3)
    {
        for (i=0; i<nm; i++)
        {
            // Fill upper-triangle; column 1
            Mt9[0] = M[6*i];
            Mt9[1] = 0.0;
            Mt9[2] = 0.0;
            // column 2
            Mt9[3] = M[6*i+3];
            Mt9[4] = M[6*i+1];
            Mt9[5] = 0.0;
            // column 3
            Mt9[6] = M[6*i+4];
            Mt9[7] = M[6*i+5];
            Mt9[8] = M[6*i+2];
            // get eigenvalues in ascending order
            info = LAPACKE_dsyev_work(LAPACK_COL_MAJOR, 'N', 'U', 3,
                                      Mt9, 3, lams, work, lwork);
            if (info != 0)
            {
                fprintf(stderr, "%s: Error computing eigenvalues\n", __func__);
                return -1;
            }
            // GCMT double couple formulation
            if (im0 == 2)
            {
                M0[i] = 0.5*(fabs(lams[0]) + fabs(lams[2]));
            }
            else
            {
                M0[i] = fabs(lams[2]);
                if (fabs(lams[0]) > fabs(lams[2])){M0[i] = fabs(lams[0]);}
            }
        }
    }
    // Silver and Jordan
    else if (im0 == 1)
    {
        for (i=0; i<nm; i++)
        {
            M11 = M[6*i];
            M22 = M[6*i+1];
            M33 = M[6*i+2];
            M12 = M[6*i+3];
            M13 = M[6*i+4];
            M23 = M[6*i+5];
            M0[i] = sqrt2i*sqrt((      M11*M11 + M22*M22 + M33*M33
                                + 2.0*(M12*M12 + M13*M13 + M23*M23)));
        }
    }
    else
    {
        fprintf(stderr, "%s: Invalid magnitude type %d\n", __func__, im0);
        for (i=0; i<nm; i++){M0[i] = 0.0;}
        return -1;
    }
    return ierr;
}

ext: c

max_stars_repo_name: OUCyf/mtbeach

max_issues_repo_name: carltape/mtbeach

max_forks_repo_name: carltape/mtbeach


text:
/*****************
 * @file leastsq.h
 * @date 4/1/08
 * @author Douglas Applegate
 *
 * @brief Provides a wrapper around the GSL libraries for fitting functions
 *****************/

//CVSID $Id: leastsq.h,v 1.1 2008-04-04 01:43:12 mallen Exp $

////////////////////////////////////////////////////

#ifndef LEASTSQ_H
#define LEASTSQ_H

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

/**
 * Function prototype for fitting function and its derivative
 * @param double x position at which to evaluate function
 * @param const double* parameter values to use in evaluation
 */
typedef double (*FitFunction_t)(double, const double*);

/**
 * Provides non-linear fitting functionality
 */
class Leastsq {

 public:

  /**
   * Initializes object with basic fit information
   * @param fitFunc function to fit to data
   * @param dfitFunc array of fit function derivatives, one for each param
   * @param nparams number of parameters that fitFunc and dfitFunc accept
   */
  Leastsq(FitFunction_t fitFunc,
	  FitFunction_t* dfitFunc,
	  int nparams);

  /**
   * Default constructor
   */
  virtual ~Leastsq();

  /**
   * Fits function to the given data, starting at initial guess
   * @param x x values of data
   * @param y y values of data
   * @param ndata number of data points
   * @param guess initial parameter values for fit
   *
   * @returns status status code for fit. 1 == fit ran fine
   */
  int doFit(const double* x, const double* y, int ndata, double* guess);

  /////GETTERS

  double numParams() const {return _nparams;};

  const double** covarMatrix() const {
    return const_cast<const double**>(_covar);
  };

  double chisq() const {return _chisq;};

  const double* parameters() const {
    return const_cast<const double*>(_parameters);
  };

 protected:

  /////INSTANCE METHODS

  int gslFitFunction(const gsl_vector* x,
		     void* params,
		     gsl_vector* result);

  int gslDFitFunction(const gsl_vector* x,
		      void* params,
		      gsl_matrix* J);

  /////Instance Variables

  FitFunction_t _fitFunc; ///< Function to fit to data
  FitFunction_t* _dfitFunc; ///< Analytical derivative of fitFunc
  int _nparams; ///< Number of parameters in fitFunc and dfitFunc
  double** _covar; ///< Covariance Matrix
  double _chisq; ///< Chisq or residual from fit
  double* _parameters; ///< Best fit parameters
  const double* _x; ///< x values of data
  const double* _y; ///< y values of data
  int _ndata; ///< number of data
};

#endif

ext: h

max_stars_repo_name: deapplegate/wtgpipeline

max_issues_repo_name: deapplegate/wtgpipeline

max_forks_repo_name: deapplegate/wtgpipeline


text:
/*! @copyright (c) 2017 King Abdullah University of Science and
 *                      Technology (KAUST). All rights reserved.
 *
 * STARS-H is a software package, provided by King Abdullah
 *             University of Science and Technology (KAUST)
 *
 * @file testing/mpi_cauchy.c
 * @version 1.3.0
 * @author Aleksandr Mikhalev
 * @date 2017-11-07
 * */

#ifdef MKL
    #include <mkl.h>
#else
    #include <cblas.h>
    #include <lapacke.h>
#endif
#include <stdio.h>
#include <stdlib.h>
#include <mpi.h>
#include <starsh.h>
#include <starsh-cauchy.h>

int main(int argc, char **argv)
{
    MPI_Init(&argc, &argv);
    int mpi_size, mpi_rank;
    MPI_Comm_size(MPI_COMM_WORLD, &mpi_size);
    MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank);
    if(argc < 5)
    {
        if(mpi_rank == 0)
        {
            printf("%d arguments provided, but 4 are needed\n", argc-1);
            printf("mpi_cauchy N block_size maxrank tol\n");
        }
        MPI_Finalize();
        return 1;
    }
    int N = atoi(argv[1]), block_size = atoi(argv[2]);
    int maxrank = atoi(argv[3]);
    double tol = atof(argv[4]);
    int onfly = 0;
    char dtype = 'd', symm = 'N';
    int ndim = 2;
    STARSH_int shape[2] = {N, N};
    int info;
    srand(0);
    // Init STARS-H
    info = starsh_init();
    if(info != 0)
    {
        MPI_Finalize();
        return 1;
    }
    // Generate data for spatial statistics problem
    STARSH_cauchy *data;
    STARSH_kernel *kernel;
    info = starsh_application((void **)&data, &kernel, N, dtype, STARSH_CAUCHY,
            STARSH_CAUCHY_KERNEL1, 0);
    if(info != 0)
    {
        MPI_Finalize();
        return 1;
    }
    // Init problem with given data and kernel and print short info
    STARSH_problem *P;
    info = starsh_problem_new(&P, ndim, shape, symm, dtype, data, data,
            kernel, "Cauchy example");
    if(info != 0)
    {
        MPI_Finalize();
        return 1;
    }
    if(mpi_rank == 0)
        starsh_problem_info(P);
    // Init plain clusterization and print info
    STARSH_cluster *C;
    info = starsh_cluster_new_plain(&C, data, N, block_size);
    if(info != 0)
    {
        MPI_Finalize();
        return 1;
    }
    if(mpi_rank == 0)
        starsh_cluster_info(C);
    // Init tlr division into admissible blocks and print short info
    STARSH_blrf *F;
    STARSH_blrm *M;
    info = starsh_blrf_new_tlr_mpi(&F, P, symm, C, C);
    if(info != 0)
    {
        MPI_Finalize();
        return 1;
    }
    if(mpi_rank == 0)
        starsh_blrf_info(F);
    // Approximate each admissible block
    MPI_Barrier(MPI_COMM_WORLD);
    double time1 = MPI_Wtime();
    info = starsh_blrm_approximate(&M, F, maxrank, tol, onfly);
    if(info != 0)
    {
        MPI_Finalize();
        return 1;
    }
    MPI_Barrier(MPI_COMM_WORLD);
    time1 = MPI_Wtime()-time1;
    if(mpi_rank == 0)
    {
        starsh_blrf_info(F);
        starsh_blrm_info(M);
    }
    if(mpi_rank == 0)
        printf("TIME TO APPROXIMATE: %e secs\n", time1);
    // Measure approximation error
    MPI_Barrier(MPI_COMM_WORLD);
    time1 = MPI_Wtime();
    double rel_err = starsh_blrm__dfe_mpi(M);
    MPI_Barrier(MPI_COMM_WORLD);
    time1 = MPI_Wtime()-time1;
    if(mpi_rank == 0)
    {
        printf("TIME TO MEASURE ERROR: %e secs\nRELATIVE ERROR: %e\n",
                time1, rel_err);
        if(rel_err/tol > 10.)
        {
            printf("Resulting relative error is too big\n");
            MPI_Finalize();
            return 1;
        }
    }
    if(rel_err/tol > 10.)
    {
        MPI_Finalize();
        return 1;
    }
    // Measure time for 10 BLRM matvecs and for 10 BLRM TLR matvecs
    double *x, *y, *y_tlr;
    int nrhs = 1;
    x = malloc(N*nrhs*sizeof(*x));
    y = malloc(N*nrhs*sizeof(*y));
    y_tlr = malloc(N*nrhs*sizeof(*y_tlr));
    if(mpi_rank == 0)
    {
        int iseed[4] = {0, 0, 0, 1};
        LAPACKE_dlarnv_work(3, iseed, N*nrhs, x);
        cblas_dscal(N*nrhs, 0.0, y, 1);
        cblas_dscal(N*nrhs, 0.0, y_tlr, 1);
    }
    MPI_Barrier(MPI_COMM_WORLD);
    time1 = MPI_Wtime();
    for(int i = 0; i < 10; i++)
        starsh_blrm__dmml_mpi(M, nrhs, 1.0, x, N, 0.0, y, N);
    MPI_Barrier(MPI_COMM_WORLD);
    time1 = MPI_Wtime()-time1;
    if(mpi_rank == 0)
    {
        printf("TIME FOR 10 BLRM MATVECS: %e secs\n", time1);
    }
    MPI_Barrier(MPI_COMM_WORLD);
    time1 = MPI_Wtime();
    for(int i = 0; i < 10; i++)
        starsh_blrm__dmml_mpi_tlr(M, nrhs, 1.0, x, N, 0.0, y_tlr, N);
    MPI_Barrier(MPI_COMM_WORLD);
    time1 = MPI_Wtime()-time1;
    if(mpi_rank == 0)
    {
        cblas_daxpy(N, -1.0, y, 1, y_tlr, 1);
        printf("TIME FOR 10 TLR MATVECS: %e secs\n", time1);
        printf("MATVEC DIFF: %e\n", cblas_dnrm2(N, y_tlr, 1)
                /cblas_dnrm2(N, y, 1));
    }
    MPI_Finalize();
    return 0;
}

max_stars_repo_name: wawando/stars-h

max_issues_repo_name: wawando/stars-h

max_forks_repo_name: wawando/stars-h


text:
#ifndef OPENMC_TALLIES_FILTER_DELAYEDGROUP_H
#define OPENMC_TALLIES_FILTER_DELAYEDGROUP_H

#include <vector>

#include <gsl/gsl>

#include "openmc/tallies/filter.h"

namespace openmc {

//==============================================================================
//! Bins outgoing fission neutrons in their delayed groups.
//!
//! The get_all_bins functionality is not actually used.  The bins are manually
//! iterated over in the scoring subroutines.
//==============================================================================

class DelayedGroupFilter : public Filter
{
public:
  //----------------------------------------------------------------------------
  // Constructors, destructors

  ~DelayedGroupFilter() = default;

  //----------------------------------------------------------------------------
  // Methods

  std::string type() const override {return "delayedgroup";}

  void from_xml(pugi::xml_node node) override;

  void get_all_bins(const Particle& p, TallyEstimator estimator, FilterMatch& match)
  const override;

  void to_statepoint(hid_t filter_group) const override;

  std::string text_label(int bin) const override;

  //----------------------------------------------------------------------------
  // Accessors

  const std::vector<int>& groups() const { return groups_; }

  void set_groups(gsl::span<int> groups);

private:
  //----------------------------------------------------------------------------
  // Data members

  std::vector<int> groups_;
};

} // namespace openmc
#endif // OPENMC_TALLIES_FILTER_DELAYEDGROUP_H

max_stars_repo_name: Hit-Weixg/openmc

max_issues_repo_name: mehmeturkmen/openmc

max_forks_repo_name: mehmeturkmen/openmc


text:
/*
	Copyright (C) 2020 David Cornu
	for the Convolutional Interactive Artificial
	Neural Networks by/for Astrophysicists (CIANNA) Code
	(https://github.com/Deyht/CIANNA)

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

#ifndef DEFS_H
#define DEFS_H

#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <tgmath.h>
#include <string.h>
#include <sys/time.h>
#include <unistd.h>

#ifdef comp_CUDA
#ifdef CUDA
#include <cuda_runtime.h>
#include <cublas_v2.h>
#if defined(GEN_VOLTA) || defined(GEN_AMPERE)
#include <cuda_fp16.h>
#endif

#if defined(GEN_AMPERE)
#include <cuda_bf16.h>
#endif

#include <curand.h>
#include <curand_kernel.h>
#endif
#endif

#ifdef BLAS
#include <cblas.h>
#endif

#ifdef OPEN_MP
#include <omp.h>
#endif

static const double two_pi = 2.0*3.14159265358979323846;

#endif //DEFS_H

max_stars_repo_name: Deyht/CIANNA

max_issues_repo_name: Deyht/CIANNA

max_forks_repo_name: Deyht/CIANNA


text:
/*****************************************************************************
*
* Rokko: Integrated Interface for libraries of eigenvalue decomposition
*
* Copyright (C) 2012-2014 by Synge Todo <wistaria@phys.s.u-tokyo.ac.jp>
*
* Distributed under the Boost Software License, Version 1.0. (See accompanying
* file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
*
*****************************************************************************/

#include <stdlib.h>
#include <stdio.h>
#include <lapacke.h>

#define N 5

int main() {
  int i, j, info;
  double w[N];
  double a[N][N] = {
    5, 4, 3, 2, 1,
    4, 4, 3, 2, 1,
    3, 3, 3, 2, 1,
    2, 2, 2, 2, 1,
    1, 1, 1, 1, 1
  };
  info = LAPACKE_dsyev(LAPACK_ROW_MAJOR, 'V', 'U', N, a[0], N, w);
  printf("Eigenvalues of Frank matrix: ");
  for (i = 0; i < N; ++i) {
    printf(" %6.2f", w[i]);
  }
  printf("\n");
  printf("Eigenvectors of Frank matrix (column by column):\n");
  for (i = 0; i < N; ++i) {
    for (j = 0; j < N; ++j) {
      printf(" %6.2f", a[i][j]);
    }
    printf("\n");
  }
  return 0;
}

max_stars_repo_name: wistaria/rokko

max_issues_repo_name: wistaria/rokko

max_forks_repo_name: wistaria/rokko


text:
//------------------------------------------------------------------------------
// Util.h
// Various utility functions and basic types used throughout the library.
//
// File is under the MIT license; see LICENSE for details.
//------------------------------------------------------------------------------
#pragma once

#include <climits>     // for type size macros
#include <cstddef>     // for std::byte
#include <cstdint>     // for sized integer types
#include <new>         // for placement new
#include <optional>    // for std::optional
#include <string_view> // for std::string_view
#include <utility>     // for many random utility functions
#include <variant>     // for std::variant

using std::byte;
using std::int16_t;
using std::int32_t;
using std::int64_t;
using std::int8_t;
using std::intptr_t;
using std::nullptr_t;
using std::optional;
using std::ptrdiff_t;
using std::size_t;
using std::string_view;
using std::uint16_t;
using std::uint32_t;
using std::uint64_t;
using std::uint8_t;
using std::uintptr_t;

using namespace std::literals;

#if !defined(ASSERT_ENABLED)
#    if !defined(NDEBUG)
#        define ASSERT_ENABLED 1
#    endif
#endif

#if ASSERT_ENABLED
#    if defined(__GNUC__) || defined(__clang__)
#        define ASSERT_FUNCTION __PRETTY_FUNCTION__
#    elif defined(_MSC_VER)
#        define ASSERT_FUNCTION __FUNCSIG__
#    elif defined(__SUNPRO_CC)
#        define ASSERT_FUNCTION __func__
#    else
#        define ASSERT_FUNCTION __FUNCTION__
#    endif
#    define ASSERT(cond)                                                                 \
        do {                                                                             \
            if (!(cond))                                                                 \
                slang::assert::assertFailed(#cond, __FILE__, __LINE__, ASSERT_FUNCTION); \
        } while (false)

#else
#    define ASSERT(cond)        \
        do {                    \
            (void)sizeof(cond); \
        } while (false)
#endif

#define THROW_UNREACHABLE                                                                  \
    throw std::logic_error(std::string(__FILE__) + ":" + std::to_string(__LINE__) + ": " + \
                           "Default case should be unreachable!")

#define __ENUM_ELEMENT(x) x,
#define __ENUM_STRING(x) #x,
#define ENUM(name, elements)                                          \
    enum class name { elements(__ENUM_ELEMENT) };                     \
    inline string_view toString(name e) {                             \
        static const char* strings[] = { elements(__ENUM_STRING) };   \
        return strings[static_cast<std::underlying_type_t<name>>(e)]; \
    }                                                                 \
    inline std::ostream& operator<<(std::ostream& os, name e) { return os << toString(e); }

#define ENUM_MEMBER(name, elements)                                   \
    enum name { elements(__ENUM_ELEMENT) };                           \
    friend string_view toString(name e) {                             \
        static const char* strings[] = { elements(__ENUM_STRING) };   \
        return strings[static_cast<std::underlying_type_t<name>>(e)]; \
    }

#include <gsl/gsl>

using gsl::finally;
using gsl::make_span;
using gsl::not_null;
using gsl::span;

// Compiler-specific macros for warnings and suppressions
#ifdef __clang__
#    define NO_SANITIZE(warningName) __attribute__((no_sanitize(warningName)))
#else
#    define NO_SANITIZE(warningName)
#endif

#include <bitmask.hpp>
using bitmask_lib::bitmask;

#include <nlohmann/json_fwd.hpp>
using json = nlohmann::json;

#define HAS_METHOD_TRAIT(name)                                                               \
    template<typename, typename T>                                                           \
    struct has_##name {                                                                      \
        static_assert(always_false<T>::value,                                                \
                      "Second template parameter needs to be of function type.");            \
    };                                                                                \
    template<typename C, typename Ret, typename... Args>                                     \
    struct has_##name<C, Ret(Args...)> {                                                     \
    private:                                                                                \
        template<typename T>                                                                 \
        static constexpr auto check(T*) ->                                                   \
            typename std::is_same<decltype(std::declval<T>().name(std::declval<Args>()...)), \
                                  Ret>::type;                                                \
        template<typename>                                                                   \
        static constexpr std::false_type check(...);                                         \
        typedef decltype(check<C>(0)) type;                                                  \
                                                                                \
    public:                                                                                \
        static constexpr bool value = type::value;                                           \
    };                                                                                \
    template<typename C, typename Ret, typename... Args>                                     \
    static constexpr bool has_##name##_v = has_##name<C, Ret(Args...)>::value

namespace slang {

template<typename T>
struct always_false : std::false_type {};

/// Converts a span of characters into a string_view.
inline string_view to_string_view(span<char> text) {
    return string_view(text.data(), (size_t)text.size());
}

inline void hash_combine(size_t&) {
}

/// Hash combining function, based on the function from Boost.
template<typename T, typename... Rest>
inline void hash_combine(size_t& seed, const T& v, Rest... rest) {
    std::hash<T> hasher;
    seed ^= hasher(v) + 0x9e3779b9 + (seed << 6) + (seed >> 2);
    hash_combine(seed, rest...);
}

namespace assert {

class AssertionException : public std::logic_error {
public:
    AssertionException(const std::string& message) : std::logic_error(message) {}
};

[[noreturn]] void assertFailed(const char* expr, const char* file, int line, const char* func);

} // namespace assert

} // namespace slang

namespace detail {

template<typename Tuple, size_t Index = std::tuple_size<Tuple>::value - 1>
struct HashValueImpl {
    static void apply(size_t& seed, const Tuple& tuple) {
        HashValueImpl<Tuple, Index - 1>::apply(seed, tuple);
        slang::hash_combine(seed, std::get<Index>(tuple));
    }
};

template<typename Tuple>
struct HashValueImpl<Tuple, 0> {
    static void apply(size_t& seed, const Tuple& tuple) {
        slang::hash_combine(seed, std::get<0>(tuple));
    }
};

} // namespace detail

namespace std {

template<typename... TT>
struct hash<std::tuple<TT...>> {
    size_t operator()(const std::tuple<TT...>& tt) const {
        size_t seed = 0;
        ::detail::HashValueImpl<std::tuple<TT...>>::apply(seed, tt);
        return seed;
    }
};

} // namespace std

max_stars_repo_name: jrrk/slang

max_issues_repo_name: jrrk/slang

max_forks_repo_name: jrrk/slang


text:
#ifndef RECTREE_H
#define RECTREE_H
#include "tree.h"
#include "recedge.h"
#include "weakarg.h"
#include "mpiutils.h"
#include "wargxml.h"

#include <gsl/gsl_randist.h>
#define MAXNODES 100000

using namespace std;

namespace weakarg
{

class Param;

/**
    @brief Recombinant tree using Tree class
*/

class RecTree:public Tree
{
    friend class RecTreeAux;
protected:
// Variables
    std::vector<std::vector<int> > tabSons;///<Contains the branching order for the current local tree
    std::vector<std::vector<double> > tabSonsDist;///<Contains the distances for the current local tree
    std::vector<int> tabNode;///<Used only in makeLocalTree (made global for speed)
    std::vector<int> tabRec;///<Used only in makeLocalTree (made global for speed)
    std::vector<int> tabFather;///<Used only in makeLocalTree (made global for speed)
    std::vector<double> age;///<Used only in makeLocalTree (made global for speed)
    std::vector<int> tabEdge;///<Used only in makeLocalTree (made global for speed)
    std::vector<bool> sameLTasPrev;///<Indicates for each site whether the local tree is the same as that of the previous site
    unsigned int L;///< Number of sites in the sequences
    std::vector<RecEdge *> edge;///< The recombinant edges
// Functions
    int randomActiveNode(std::vector<int> nodelist);///<Returns a random node id from a vector of indicators for alive (-1 for not alive at current time)
    int getEdgeCoal(std::vector<int> nl,int k,double* time);///<Returns the node the recombinant edge is from and sets the time

    int moveClonalFixFrom(std::vector<int> ornodes=std::vector<int>(0));///< Fixes all from times to be in to valid nodes and reinstates **the vector ornodes. returns the number of recedges moved to younger nodes - number moved to older nodes (or -1 if this is not calculated)

    void updateEdgeTimes(int e,double reldist);///moves an edge when its node changes by proportion reldist
    void scaleEdges(int which, double dist);///<Scales the recedges "to" on edge which by an amount dist
    void fixFromTimes(int affedge,double dist);///<Fixes from times when we've moved edge affedge an amount dist
    void fixToTimes(int which, double dist);///<Fixes "to" times when we've changed an node time implicitly

    void swapEdge(int a, int b);///<Swaps two edges in the edge list (ordering only)
    void swapEdgeTo(int a, int b);///<Swaps the edges going to *two specified --nodes--*

    void orderEdges(int which);///<Checks that the edge which is in the correct place in the list (from 0.. which)

    inline int orderNodes(int which, double dist)
    {
	return(orderNodes(which,dist,-1));
    }///<places a node in a new place in the list (updated to account for recedges);
    int orderNodes(int which, double dist, int oldwhich);///<places a node in a new place in the list (updated to account for recedges); a known position is a positive oldwhich...
    void swapNode(int a, int b);///<Swaps two nodes in the list (updated to maintain recedges going to them)
    int lastCommonAncestor(int s1,int s2);///< USES TABNODES AND TABFATHER!  Returns the last common ancestor index of the tabnodes ASSUMING the current local tree.
    void makeLocalTreeKeepingFathers(unsigned int site);///< Makes a version of the local tree that has the fathers but is unoptimised
    double pairwiseDist(int s1,int s2);///< Returns the pairwise distance between two sequences accounting for recombination.  SAME WARNINGS AS lastCommonAncestor!
public:

// Constructors and destructors:
    RecTree(unsigned int numsites,WargXml *infile,bool addedges=true,bool forceages=true,bool isfinal=true,std::streampos itstart=-1);///<Creates a tree from an output file

    RecTree(unsigned int numsites,std::string newick,bool isFilename=true,bool forceages=true);///<Creates a tree from a Newick file
    RecTree(int n,double rho,double delta,std::vector<int> blocks);///<Creates a Coalescent tree with recombinant edges by simulation
    void dropEdges(double rho,double delta,vector<int> blocks);
    RecTree(Data * data,double rho,double delta,vector<int> blocks);///<Creates UPGMA tree
    RecTree(RecTree *intree,string newick,bool isFilename=false,bool forceages=true);/// Copies a rectree
    RecTree(const RecTree& rt);///< Copy constructor
    void assign(const RecTree& rt);
    ~RecTree();
// Part of the simulation initialisation:
    void setBlock(unsigned int* gstart,unsigned int* gend,double delta,std::vector<int> *blocks);///<Sets the start and end positions of a recombinant block based on geometric distribution and independent blocks
    int getEdgeCoal(double* time);///<Returns the node the recombinant edge is from and sets the time
// Part of the likelihood/prior calculation:
    double priorEdge(int e,Param * param) const;///<Returns the LOG prior of a given edge number in rectree
    double priorEdge(RecEdge* edge0,Param * param) const;///FMA_CHANGES: Computes prior of edge not in rectree
    double priorEdge(double tFrom,double tTo) const;///<Returns the UNLOGGED prior of a given edge
    double prior(Param * param) const;///<Returns the value of the prior function

    void makeLocalTree(unsigned int site,double thetaPerSite);///<Evaluates the local tree at the given site
// Tests:
    void testEdges() const;///<Tests recedges arrive in between the ages of the nodes they go to
    void testTree();///<Checks ages increase in the tree for all sites.
    inline bool sameLocalTreeAsPrev(unsigned int site) const
    {
        return sameLTasPrev[site];
    }///<Returns true if the local tree of site is the same as that of site-1
    void calcSameLTasPrev(unsigned int site);///<Estimate if the LT is the same as the previous one
// Information functions:
    inline Node * getRoot(){return(root);}
    inline double getEdgeTimeAbsFrom(int e) const
    {
        return edge[e]->getTimeFrom()+nodes[edge[e]->getEdgeFrom()]->getAge();
    }
    ;///<Returns the absolute age of the origin of a recombinant edge
    inline double getEdgeTimeAbsTo  (int e) const
    {
        return edge[e]->getTimeTo  ()+nodes[edge[e]->getEdgeTo  ()]->getAge();
    }
    ;///<Returns the absolute age of the destination of a recombinant edge
    inline int getL() const
    {
        return L;
    }///<Returns the total size of the alignment

    inline RecEdge*getRecEdge(int w) const
    {
        return edge[w];
    }///<Returns the w-th recombinant edge
    inline int getSon (int i, unsigned int site,bool isLeft,double *dist) const
    {
        site=0;
        *dist=tabSonsDist[i][isLeft];
        return tabSons[i][isLeft];
    };
    ///<Returns the left daughter node for a given local site and node (NULL if no left daughter)
    inline int affecting(unsigned int site) const
    {
        int a=0;
        for (unsigned int i=0;i<edge.size();i++)
            if (edge[i]->affectsSite(site))
                a++;
        return a;
    }
    ;///<Returns the number of edges affecting a site
    inline int numRecEdge() const
    {
        return edge.size();
    }///<Returns the number of recombinant edges
    int numRecEdgeOnBranch(int b) const
    {
        int r=0;
        for (unsigned int i=0;i<edge.size();i++)
            if (edge[i]->getEdgeTo()==b)
                r++;
        return r;
    }///<Returns the number of recedges on a given branch
    std::vector<int> alive(double t) const;///<Returns the branches alive at a given time
    bool isThere(int start,int end,int e,double time1,double time2) const;///<Returns yes if there is an edge with departure or arrival on branch e and between time1 and time2, and which affects material that intersects [start;end]
    inline RecEdge* getEdge(int i) const
    {
        return(edge[i]);
    }///<Returns the edge with index i
    double getEdgeTreeTime(int i) const;///<Returns the time an edge spans in the tree
    vector<int> getAffEdges(vector<int> * samplespace);///<Calculates the recedges relevent to a set of clonal edges
    void updateList(int which, int newloc,vector<int> * list);///<Updates a list of edges when which was moved to newloc
    int sampleEdge(vector<int> * samplespace=NULL);///<Samples a random edge that (optionally) interacts with a specific set of clonal edges
// Modification functions:
    void addEdgesFromFile(WargXml *infile,int siteoffset=0);///< adds edges to the rectree from the specified warg XML file (which is set to an iteration). Adds a "siteoffset" to add edges from different runs on different genes
    int addRecEdge(double tfrom,double tto,unsigned int gstart,unsigned int gend,int edgefrom,int edgeto);///<Adds an recombination edge to the tree
    int addRecEdge_FMA(double tfrom,double tto,unsigned int gstart,unsigned int gend,int edgefrom,int edgeto);///FMA_CHANGES: <Adds an recombination edge to the tree
    int addRecEdge(unsigned int gstart,unsigned int gend,int edgefrom,int edgeto);///< Adds a recombination edge with random start and end times between given edges
    int addRecEdge(std::string res,int sitesoffset=0);///<Add recedge from our output format
    void remRecEdge(int which);///<Removes a recombinant edge
    void setStart(int edge,int start);///<Sets the starting point of the given recedge
    void setEnd(int edge,int end);///<Sets the ending point of the given recedge
    void changeAge(int which,double dist);///<Changes the "which" node age, maintaining edge times
    inline void changeEdgeNodeFrom(int e,int newnode)
    {
        edge[e]->setTimeFrom(getEdgeTimeAbsFrom(e)-getNode(newnode)->getAge(),newnode);
    }///<changes the node to without changing the absolute time of the edge
    inline void changeEdgeNodeTo(int e,int newnode)
    {
        edge[e]->setTimeTo(getEdgeTimeAbsTo(e)-getNode(newnode)->getAge(),newnode);
    }///<changes the node to without changing the absolute time of the edge
    void chooseNodePath(int which, double dist, vector<int> *listLR, vector<int> *affedges);///< Sets up a specified order for which nodes to swap when they are disordered
    int moveNodeTime(int which, int *whichto,double dist,int oldwhich, vector<int> *listLR,  vector<int> *affedges);///<Updates the node list when the node which is moved a distance dist number of recedges moved to younger nodes - number moved to older nodes. edits affedges to make it the a vector of affected edges

    void moveEdges(int e1,int e2,double t1=0.0,double t2=-1.0);///<Moves all edge events on e1 between absolute times t1..t2 onto edge e2, maintaining their absolute times (t2=-1 means no upper limit)
// Information functions
    std::vector<std::vector<double> > pairwiseDistanceMatrix(long site);///< returns a pairwise distance matrix for a given site

};

} // end namespace weakarg
#endif


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

/* Author:  G. Jungman */

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

#include "error.h"

#include "chebyshev.h"
#include "cheb_eval.c"

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

/* Chebyshev expansion for log(1 + x(t))/x(t)
 *
 * x(t) = (4t-1)/(2(4-t))
 * t(x) = (8x+1)/(2(x+2))
 * -1/2 < x < 1/2
 * -1 < t < 1
 */
static double lopx_data[21] = {
  2.16647910664395270521272590407,
 -0.28565398551049742084877469679,
  0.01517767255690553732382488171,
 -0.00200215904941415466274422081,
  0.00019211375164056698287947962,
 -0.00002553258886105542567601400,
  2.9004512660400621301999384544e-06,
 -3.8873813517057343800270917900e-07,
  4.7743678729400456026672697926e-08,
 -6.4501969776090319441714445454e-09,
  8.2751976628812389601561347296e-10,
 -1.1260499376492049411710290413e-10,
  1.4844576692270934446023686322e-11,
 -2.0328515972462118942821556033e-12,
  2.7291231220549214896095654769e-13,
 -3.7581977830387938294437434651e-14,
  5.1107345870861673561462339876e-15,
 -7.0722150011433276578323272272e-16,
  9.7089758328248469219003866867e-17,
 -1.3492637457521938883731579510e-17,
  1.8657327910677296608121390705e-18
};
static cheb_series lopx_cs = {
  lopx_data,
  20,
  -1, 1,
  10
};

/* Chebyshev expansion for (log(1 + x(t)) - x(t))/x(t)^2
 *
 * x(t) = (4t-1)/(2(4-t))
 * t(x) = (8x+1)/(2(x+2))
 * -1/2 < x < 1/2
 * -1 < t < 1
 */
static double lopxmx_data[20] = {
 -1.12100231323744103373737274541,
  0.19553462773379386241549597019,
 -0.01467470453808083971825344956,
  0.00166678250474365477643629067,
 -0.00018543356147700369785746902,
  0.00002280154021771635036301071,
 -2.8031253116633521699214134172e-06,
  3.5936568872522162983669541401e-07,
 -4.6241857041062060284381167925e-08,
  6.0822637459403991012451054971e-09,
 -8.0339824424815790302621320732e-10,
  1.0751718277499375044851551587e-10,
 -1.4445310914224613448759230882e-11,
  1.9573912180610336168921438426e-12,
 -2.6614436796793061741564104510e-13,
  3.6402634315269586532158344584e-14,
 -4.9937495922755006545809120531e-15,
  6.8802890218846809524646902703e-16,
 -9.5034129794804273611403251480e-17,
  1.3170135013050997157326965813e-17
};
static cheb_series lopxmx_cs = {
  lopxmx_data,
  19,
  -1, 1,
  9
};

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

#ifndef HIDE_INLINE_STATIC
int
gsl_sf_log_e(const double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(x <= 0.0) {
    DOMAIN_ERROR(result);
  }
  else {
    result->val = log(x);
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
}

int
gsl_sf_log_abs_e(const double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(x == 0.0) {
    DOMAIN_ERROR(result);
  }
  else {
    result->val = log(fabs(x));
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
}
#endif

int
gsl_sf_complex_log_e(const double zr, const double zi, gsl_sf_result * lnr, gsl_sf_result * theta)
{
  /* CHECK_POINTER(lnr) */
  /* CHECK_POINTER(theta) */

  if(zr != 0.0 || zi != 0.0) {
    const double ax = fabs(zr);
    const double ay = fabs(zi);
    const double min = GSL_MIN(ax, ay);
    const double max = GSL_MAX(ax, ay);
    lnr->val = log(max) + 0.5 * log(1.0 + (min/max)*(min/max));
    lnr->err = 2.0 * GSL_DBL_EPSILON * fabs(lnr->val);
    theta->val = atan2(zi, zr);
    theta->err = GSL_DBL_EPSILON * fabs(lnr->val);
    return GSL_SUCCESS;
  }
  else {
    DOMAIN_ERROR_2(lnr, theta);
  }
}

int
gsl_sf_log_1plusx_e(const double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(x <= -1.0) {
    DOMAIN_ERROR(result);
  }
  else if(fabs(x) < GSL_ROOT6_DBL_EPSILON) {
    const double c1 = -0.5;
    const double c2 =  1.0/3.0;
    const double c3 = -1.0/4.0;
    const double c4 =  1.0/5.0;
    const double c5 = -1.0/6.0;
    const double c6 =  1.0/7.0;
    const double c7 = -1.0/8.0;
    const double c8 =  1.0/9.0;
    const double c9 = -1.0/10.0;
    const double t  =  c5 + x*(c6 + x*(c7 + x*(c8 + x*c9)));
    result->val = x * (1.0 + x*(c1 + x*(c2 + x*(c3 + x*(c4 + x*t)))));
    result->err = GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(fabs(x) < 0.5) {
    double t = 0.5*(8.0*x + 1.0)/(x+2.0);
    gsl_sf_result c;
    cheb_eval_e(&lopx_cs, t, &c);
    result->val = x * c.val;
    result->err = fabs(x * c.err);
    return GSL_SUCCESS;
  }
  else {
    result->val = log(1.0 + x);
    result->err = GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
}

int
gsl_sf_log_1plusx_mx_e(const double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(x <= -1.0) {
    DOMAIN_ERROR(result);
  }
  else if(fabs(x) < GSL_ROOT5_DBL_EPSILON) {
    const double c1 = -0.5;
    const double c2 =  1.0/3.0;
    const double c3 = -1.0/4.0;
    const double c4 =  1.0/5.0;
    const double c5 = -1.0/6.0;
    const double c6 =  1.0/7.0;
    const double c7 = -1.0/8.0;
    const double c8 =  1.0/9.0;
    const double c9 = -1.0/10.0;
    const double t  =  c5 + x*(c6 + x*(c7 + x*(c8 + x*c9)));
    result->val = x*x * (c1 + x*(c2 + x*(c3 + x*(c4 + x*t))));
    result->err = GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(fabs(x) < 0.5) {
    double t = 0.5*(8.0*x + 1.0)/(x+2.0);
    gsl_sf_result c;
    cheb_eval_e(&lopxmx_cs, t, &c);
    result->val = x*x * c.val;
    result->err = x*x * c.err;
    return GSL_SUCCESS;
  }
  else {
    const double lterm = log(1.0 + x);
    result->val = lterm - x;
    result->err = GSL_DBL_EPSILON * (fabs(lterm) + fabs(x));
    return GSL_SUCCESS;
  }
}

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

#include "eval.h"

double gsl_sf_log(const double x)
{
  EVAL_RESULT(gsl_sf_log_e(x, &result));
}

double gsl_sf_log_abs(const double x)
{
  EVAL_RESULT(gsl_sf_log_abs_e(x, &result));
}

double gsl_sf_log_1plusx(const double x)
{
  EVAL_RESULT(gsl_sf_log_1plusx_e(x, &result));
}

double gsl_sf_log_1plusx_mx(const double x)
{
  EVAL_RESULT(gsl_sf_log_1plusx_mx_e(x, &result));
}


text:
static char help[] =
"Coupled reaction-diffusion equations (Pearson 1993).  Option prefix -ptn_.\n"
"Demonstrates form  F(t,Y,dot Y) = G(t,Y)  where F() is IFunction and G() is\n"
"RHSFunction().  Implements IJacobian().  Defaults to ARKIMEX (= adaptive\n"
"Runge-Kutta implicit-explicit) TS type.\n\n";

// compare runs with
//    -dm_mat_type aij|baij|sbaij
//    -dm_mat_type sbaij -ksp_type cg -pc_type icc  # 20% speed up?

#include <petsc.h>

typedef struct {
  double u, v;
} Field;

typedef struct {
  double    L,     // domain side length
            Du,    // diffusion coefficient: u equation
            Dv,    //                        v equation
            phi,   // "dimensionless feed rate" (F in Pearson 1993)
            kappa; // "dimensionless rate constant" (k in Pearson 1993)
} PatternCtx;

extern PetscErrorCode InitialState(DM, Vec, double, PatternCtx*);
extern PetscErrorCode FormRHSFunctionLocal(DMDALocalInfo*, double, Field**,
                                           Field**, PatternCtx*);
extern PetscErrorCode FormIFunctionLocal(DMDALocalInfo*, double, Field**, Field**,
                                         Field **, PatternCtx*);
extern PetscErrorCode FormIJacobianLocal(DMDALocalInfo*, double, Field**, Field**,
                                         double, Mat, Mat, PatternCtx*);

int main(int argc,char **argv)
{
  PetscErrorCode ierr;
  PatternCtx     user;
  TS             ts;
  Vec            x;
  DM             da;
  DMDALocalInfo  info;
  double         noiselevel = -1.0;  // negative value means no initial noise

  PetscInitialize(&argc,&argv,(char*)0,help);

  // parameter values from pages 21-22 in Hundsdorfer & Verwer (2003)
  user.L      = 2.5;
  user.Du     = 8.0e-5;
  user.Dv     = 4.0e-5;
  user.phi    = 0.024;
  user.kappa  = 0.06;
  ierr = PetscOptionsBegin(PETSC_COMM_WORLD, "ptn_", "options for patterns", ""); CHKERRQ(ierr);
  ierr = PetscOptionsReal("-noisy_init",
           "initialize u,v with this much random noise (e.g. 0.2) on top of usual initial values",
           "pattern.c",noiselevel,&noiselevel,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsReal("-L","square domain side length; recommend L >= 0.5",
           "pattern.c",user.L,&user.L,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsReal("-Du","diffusion coefficient of first equation",
           "pattern.c",user.Du,&user.Du,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsReal("-Dv","diffusion coefficient of second equation",
           "pattern.c",user.Dv,&user.Dv,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsReal("-phi","dimensionless feed rate (=F in (Pearson, 1993))",
           "pattern.c",user.phi,&user.phi,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsReal("-kappa","dimensionless rate constant (=k in (Pearson, 1993))",
           "pattern.c",user.kappa,&user.kappa,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsEnd(); CHKERRQ(ierr);

  ierr = DMDACreate2d(PETSC_COMM_WORLD,
               DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC,
               DMDA_STENCIL_BOX,  // for 9-point stencil
               3,3,PETSC_DECIDE,PETSC_DECIDE,
               2, 1,              // degrees of freedom, stencil width
               NULL,NULL,&da); CHKERRQ(ierr);
  ierr = DMSetFromOptions(da); CHKERRQ(ierr);
  ierr = DMSetUp(da); CHKERRQ(ierr);
  ierr = DMDASetFieldName(da,0,"u"); CHKERRQ(ierr);
  ierr = DMDASetFieldName(da,1,"v"); CHKERRQ(ierr);
  ierr = DMDAGetLocalInfo(da,&info); CHKERRQ(ierr);
  if (info.mx != info.my) {
      SETERRQ(PETSC_COMM_WORLD,1,"pattern.c requires mx == my");
  }
  ierr = DMDASetUniformCoordinates(da, 0.0, user.L, 0.0, user.L, -1.0, -1.0); CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD,
           "running on %d x %d grid with square cells of side h = %.6f ...\n",
           info.mx,info.my,user.L/(double)(info.mx)); CHKERRQ(ierr);

//STARTTSSETUP
  ierr = TSCreate(PETSC_COMM_WORLD,&ts); CHKERRQ(ierr);
  ierr = TSSetProblemType(ts,TS_NONLINEAR); CHKERRQ(ierr);
  ierr = TSSetDM(ts,da); CHKERRQ(ierr);
  ierr = TSSetApplicationContext(ts,&user); CHKERRQ(ierr);
  ierr = DMDATSSetRHSFunctionLocal(da,INSERT_VALUES,
           (DMDATSRHSFunctionLocal)FormRHSFunctionLocal,&user); CHKERRQ(ierr);
  ierr = DMDATSSetIFunctionLocal(da,INSERT_VALUES,
           (DMDATSIFunctionLocal)FormIFunctionLocal,&user); CHKERRQ(ierr);
  ierr = DMDATSSetIJacobianLocal(da,
           (DMDATSIJacobianLocal)FormIJacobianLocal,&user); CHKERRQ(ierr);
  ierr = TSSetType(ts,TSARKIMEX); CHKERRQ(ierr);
  ierr = TSSetTime(ts,0.0); CHKERRQ(ierr);
  ierr = TSSetMaxTime(ts,200.0); CHKERRQ(ierr);
  ierr = TSSetTimeStep(ts,5.0); CHKERRQ(ierr);
  ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP); CHKERRQ(ierr);
  ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
//ENDTSSETUP

  ierr = DMCreateGlobalVector(da,&x); CHKERRQ(ierr);
  ierr = InitialState(da,x,noiselevel,&user); CHKERRQ(ierr);
  ierr = TSSolve(ts,x); CHKERRQ(ierr);

  VecDestroy(&x);  TSDestroy(&ts);  DMDestroy(&da);
  return PetscFinalize();
}

// Formulas from page 22 of Hundsdorfer & Verwer (2003).  Interpretation here is
// to always generate 0.5 x 0.5 non-trivial patch in (0,L) x (0,L) domain.
PetscErrorCode InitialState(DM da, Vec Y, double noiselevel, PatternCtx* user) {
  PetscErrorCode ierr;
  DMDALocalInfo  info;
  int            i,j;
  double         sx,sy;
  const double   ledge = (user->L - 0.5) / 2.0, // nontrivial initial values on
                 redge = user->L - ledge;       //   ledge < x,y < redge
  DMDACoor2d     **aC;
  Field          **aY;

  ierr = VecSet(Y,0.0); CHKERRQ(ierr);
  if (noiselevel > 0.0) {
      // noise added to usual initial condition is uniform on [0,noiselevel],
      //     independently for each location and component
      ierr = VecSetRandom(Y,NULL); CHKERRQ(ierr);
      ierr = VecScale(Y,noiselevel); CHKERRQ(ierr);
  }
  ierr = DMDAGetLocalInfo(da,&info); CHKERRQ(ierr);
  ierr = DMDAGetCoordinateArray(da,&aC); CHKERRQ(ierr);
  ierr = DMDAVecGetArray(da,Y,&aY); CHKERRQ(ierr);
  for (j = info.ys; j < info.ys+info.ym; j++) {
    for (i = info.xs; i < info.xs+info.xm; i++) {
      if ((aC[j][i].x >= ledge) && (aC[j][i].x <= redge)
              && (aC[j][i].y >= ledge) && (aC[j][i].y <= redge)) {
          sx = sin(4.0 * PETSC_PI * aC[j][i].x);
          sy = sin(4.0 * PETSC_PI * aC[j][i].y);
          aY[j][i].v += 0.5 * sx * sx * sy * sy;
      }
      aY[j][i].u += 1.0 - 2.0 * aY[j][i].v;
    }
  }
  ierr = DMDAVecRestoreArray(da,Y,&aY); CHKERRQ(ierr);
  ierr = DMDARestoreCoordinateArray(da,&aC); CHKERRQ(ierr);
  return 0;
}

// in system form  F(t,Y,dot Y) = G(t,Y),  compute G():
//     G^u(t,u,v) = - u v^2 + phi (1 - u)
//     G^v(t,u,v) = + u v^2 - (phi + kappa) v
//STARTRHSFUNCTION
PetscErrorCode FormRHSFunctionLocal(DMDALocalInfo *info,
                   double t, Field **aY, Field **aG, PatternCtx *user) {
  int            i, j;
  double         uv2;

  for (j = info->ys; j < info->ys + info->ym; j++) {
      for (i = info->xs; i < info->xs + info->xm; i++) {
          uv2 = aY[j][i].u * aY[j][i].v * aY[j][i].v;
          aG[j][i].u = - uv2 + user->phi * (1.0 - aY[j][i].u);
          aG[j][i].v = + uv2 - (user->phi + user->kappa) * aY[j][i].v;
      }
  }
  return 0;
}
//ENDRHSFUNCTION

// in system form  F(t,Y,dot Y) = G(t,Y),  compute F():
//     F^u(t,u,v,u_t,v_t) = u_t - D_u Laplacian u
//     F^v(t,u,v,u_t,v_t) = v_t - D_v Laplacian v
//STARTIFUNCTION
PetscErrorCode FormIFunctionLocal(DMDALocalInfo *info,
                   double t, Field **aY, Field **aYdot, Field **aF,
                   PatternCtx *user) {
  int            i, j;
  const double   h = user->L / (double)(info->mx),
                 Cu = user->Du / (6.0 * h * h),
                 Cv = user->Dv / (6.0 * h * h);
  double         u, v, lapu, lapv;

  for (j = info->ys; j < info->ys + info->ym; j++) {
      for (i = info->xs; i < info->xs + info->xm; i++) {
          u = aY[j][i].u;
          v = aY[j][i].v;
          lapu =     aY[j+1][i-1].u + 4.0*aY[j+1][i].u +   aY[j+1][i+1].u
                 + 4.0*aY[j][i-1].u -    20.0*u        + 4.0*aY[j][i+1].u
                 +   aY[j-1][i-1].u + 4.0*aY[j-1][i].u +   aY[j-1][i+1].u;
          lapv =     aY[j+1][i-1].v + 4.0*aY[j+1][i].v +   aY[j+1][i+1].v
                 + 4.0*aY[j][i-1].v -    20.0*v        + 4.0*aY[j][i+1].v
                 +   aY[j-1][i-1].v + 4.0*aY[j-1][i].v +   aY[j-1][i+1].v;
          aF[j][i].u = aYdot[j][i].u - Cu * lapu;
          aF[j][i].v = aYdot[j][i].v - Cv * lapv;
      }
  }
  return 0;
}
//ENDIFUNCTION

// in system form  F(t,Y,dot Y) = G(t,Y),  compute combined/shifted
// Jacobian of F():
//     J = (shift) dF/d(dot Y) + dF/dY
//STARTIJACOBIAN
PetscErrorCode FormIJacobianLocal(DMDALocalInfo *info,
                   double t, Field **aY, Field **aYdot, double shift,
                   Mat J, Mat P, PatternCtx *user) {
    PetscErrorCode ierr;
    int            i, j, s, c;
    const double   h = user->L / (double)(info->mx),
                   Cu = user->Du / (6.0 * h * h),
                   Cv = user->Dv / (6.0 * h * h);
    double         val[9], CC;
    MatStencil     col[9], row;

    for (j = info->ys; j < info->ys + info->ym; j++) {
        row.j = j;
        for (i = info->xs; i < info->xs + info->xm; i++) {
            row.i = i;
            for (c = 0; c < 2; c++) { // u,v equations are c=0,1
                row.c = c;
                CC = (c == 0) ? Cu : Cv;
                for (s = 0; s < 9; s++)
                    col[s].c = c;
                col[0].i = i;   col[0].j = j;
                val[0] = shift + 20.0 * CC;
                col[1].i = i-1; col[1].j = j;    val[1] = - 4.0 * CC;
                col[2].i = i+1; col[2].j = j;    val[2] = - 4.0 * CC;
                col[3].i = i;   col[3].j = j-1;  val[3] = - 4.0 * CC;
                col[4].i = i;   col[4].j = j+1;  val[4] = - 4.0 * CC;
                col[5].i = i-1; col[5].j = j-1;  val[5] = - CC;
                col[6].i = i-1; col[6].j = j+1;  val[6] = - CC;
                col[7].i = i+1; col[7].j = j-1;  val[7] = - CC;
                col[8].i = i+1; col[8].j = j+1;  val[8] = - CC;
                ierr = MatSetValuesStencil(P,1,&row,9,col,val,INSERT_VALUES); CHKERRQ(ierr);
            }
        }
    }

    ierr = MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    ierr = MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    if (J != P) {
        ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
        ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    }
    return 0;
}
//ENDIJACOBIAN


text:
#include <stdlib.h>
#include <gsl/block/gsl_block.h>

#define BASE_DOUBLE
#include <gsl/templates_on.h>
#include <gsl/block/init_source.c>
#include <gsl/templates_off.h>
#undef BASE_DOUBLE

max_stars_repo_name: MontyThibault/centre-of-mass-awareness

max_issues_repo_name: MontyThibault/centre-of-mass-awareness

max_forks_repo_name: MontyThibault/centre-of-mass-awareness


text:
/*
 * Copyright (c) 2020 Microsoft. All rights reserved.
 */

#pragma once

// a place for frequently used types, especially ones relied on for interfaces
// eg gsl::span<uint8_t> etc

#include "machine.h"

#include <vector>

#include <gsl/span>
#include <gsl/span_ext>

#include <gsl/util>

#include <type_traits> // for make_unsigned

template <typename T>
T RemoveUnsigned(typename std::make_unsigned<T>::type val) {
  return gsl::narrow_cast<T>(val);
}

template <typename T> typename std::make_unsigned<T>::type MakeUnsigned(T val) {
  return gsl::narrow_cast<typename std::make_unsigned<T>::type>(val);
}

// encapsulate some memory which may be written, use instead of a ptr, len
using MemRef = gsl::span<uint8_t>;
using MemRefSizeType = MemRef::size_type;

using SignedMemRef = gsl::span<char>;
using SignedMemRefSizeType = SignedMemRef::size_type;

using MemVector = std::vector<uint8_t>;

// encapsulate some read only memory
using CMemRef = gsl::span<const uint8_t>;
using CSignedMemRef = gsl::span<const char>;

// from any old class T make a span covering uint8_ts
// intended for the Params classes.
template <class T> MemRef MakeMemRef(T *p) noexcept {
  void *ptr = p;
  static_assert(sizeof(*p) != 1);
  return MemRef(static_cast<uint8_t *>(ptr), sizeof(*p));
}

template <class T> MemRef MakeMemRef(T *p, size_t len) noexcept {
  void *ptr = p;
  return MemRef(static_cast<uint8_t *>(ptr), len);
}

template <class T> CMemRef MakeCMemRef(const T *p, size_t len) noexcept {
  const void *ptr = p;
  return CMemRef(static_cast<const uint8_t *>(ptr), len);
}

template <class T> CMemRef MakeCMemRef(const T *p) noexcept {
  const void *ptr = p;
  return CMemRef(static_cast<const uint8_t *>(ptr), sizeof(*p));
}

template <class T, size_t TLength>
CMemRef MakeCMemRef2(const T (&p)[TLength]) noexcept {
  const void *ptr = p;
  return CMemRef(static_cast<const uint8_t *>(ptr), TLength);
}

template <class T> SignedMemRef MakeSignedMemRef(T *p) noexcept {
  void *ptr = p;
  static_assert(sizeof(*p) != 1);
  return SignedMemRef(static_cast<char *>(ptr), sizeof(*p));
}

template <class T> SignedMemRef MakeSignedMemRef(T *p, size_t len) noexcept {
  void *ptr = p;
  return SignedMemRef(static_cast<char *>(ptr), len);
}

template <class T> CSignedMemRef MakeCSignedMemRef(T *p) noexcept {
  const void *ptr = p;
  static_assert(sizeof(*p) != 1);
  return CSignedMemRef(static_cast<const char *>(ptr), sizeof(*p));
}

template <class T> CSignedMemRef MakeCSignedMemRef(T *p, size_t len) noexcept {
  const void *ptr = p;
  return CSignedMemRef(static_cast<const char *>(ptr), len);
}

// T will needs top be copyable
template <class T> T ReHydrate(CMemRef source) noexcept {
  const void *ptr = source.data();
  const T *tInPlace = static_cast<const T *>(ptr);
  T ret = *tInPlace;
  return ret;
}

// Work around C4686
using _FixedLengthString = std::array<char, 1024>;
class FixedLengthString : public _FixedLengthString {
  using _FixedLengthString::_FixedLengthString;
};

// helper for HexDump
void HexByte(uint8_t byte, char &c1, char &c2);
// hex dump a block of memory into another block of memory.
size_t HexDump(CMemRef src, SignedMemRef dst);

template <class T, int LEN> struct IsPODEtc {
  static constexpr bool value =
      sizeof(T) == LEN && std::is_pod<T>::value && std::is_trivial<T>::value &&
      std::is_trivially_copyable<T>::value && std::is_standard_layout<T>::value;
};

/*
    Compile time protection against objects being different between the
    host and HW.

    usage - put something like the following outside the class:

    static_assert(OkAsParam<MyClass, 666>::value);

*/

template <class T, int LEN> struct SizeIsOK {
  static constexpr bool value = sizeof(T) == LEN;
};

template <class T> struct LayoutIsOK {
  // These two checks fail for most of our types as they have user provided
  // default ctors or members with default initialisers.
  // static_assert(std::is_pod<T>::value);
  // static_assert(std::is_trivial<T>::value);

  static constexpr bool value =

      // but anything we can serialise by copying must be trivially copyable
      std::is_trivially_copyable<T>::value &&

      // and must be standard layout, ie if you can't pass it to a C program
      // what hope is there that you can pass it over the wire to the HW.
      std::is_standard_layout<T>::value;
};

template <class T, int LEN> struct OkAsParam {
  static constexpr bool value = SizeIsOK<T, LEN>::value && LayoutIsOK<T>::value;
};

// Thinking about labeling classes explicitly. TBC, ignore for now.

class CannotBeSerialsed {
public:
  constexpr bool canBeSerialised() const { return false; }
};

class CanBeSerialsed {
public:
  constexpr bool canBeSerialised() const { return true; }
};


text:
#ifndef QDM_TAU_H
#define QDM_TAU_H 1

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

typedef struct {
  bool use_table;

  double low;
  double high;

  size_t spline_df;
  const gsl_vector *knots;

  gsl_vector *reset;
  gsl_matrix *ispline_table;
  gsl_matrix *mspline_table;
} qdm_tau;

qdm_tau *
qdm_tau_alloc(
    int use_table,

    double low,
    double high,

    size_t spline_df,
    const gsl_vector *knots
);

void
qdm_tau_free(
    qdm_tau *t
);

double
qdm_tau_ispline_mmm(
    const qdm_tau *t,
    const double value,
    const gsl_vector *mmm
);

double
qdm_tau_mspline_mmm(
    const qdm_tau *t,
    const double value,
    const gsl_vector *mmm
);

double
qdm_tau_find(
    const qdm_tau *t,
    const double v,
    const gsl_vector *mmm
);

#endif /* QDM_TAU_H */

max_stars_repo_name: calebcase/qdm

max_issues_repo_name: calebcase/qdm

max_forks_repo_name: calebcase/qdm