text:
/*****************************************************************************
Copyright (c) 2011-2014, The OpenBLAS Project
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 OpenBLAS project 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 OWNER 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 "openblas_utest.h"
#include <sys/types.h>
#include <sys/wait.h>
#include <cblas.h>

void* xmalloc(size_t n)
{
    void* tmp;
    tmp = malloc(n);
    if (tmp == NULL) {
        fprintf(stderr, "You are about to die\n");
        exit(1);
    } else {
        return tmp;
    }
}

void check_dgemm(double *a, double *b, double *result, double *expected, int n)
{
    int i;
    cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, n, n, n,
        1.0, a, n, b, n, 0.0, result, n);
    for(i = 0; i < n * n; ++i) {
        ASSERT_DBL_NEAR_TOL(expected[i], result[i], DOUBLE_EPS);
    }
}

CTEST(fork, safety)
{
    int n = 1000;
    int i;

    double *a, *b, *c, *d;
    size_t n_bytes;

    pid_t fork_pid;
    pid_t fork_pid_nested;

    n_bytes = sizeof(*a) * n * n;

    a = xmalloc(n_bytes);
    b = xmalloc(n_bytes);
    c = xmalloc(n_bytes);
    d = xmalloc(n_bytes);

    // Put ones in a and b
    for(i = 0; i < n * n; ++i) {
        a[i] = 1;
        b[i] = 1;
    }

    // Compute a DGEMM product in the parent process prior to forking to
    // ensure that the OpenBLAS thread pool is initialized.
    cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, n, n, n,
       1.0, a, n, b, n, 0.0, c, n);

    fork_pid = fork();
    if (fork_pid == -1) {
        CTEST_ERR("Failed to fork process.");
    } else if (fork_pid == 0) {
        // Compute a DGEMM product in the child process to check that the
        // thread pool as been properly been reinitialized after the fork.
        check_dgemm(a, b, d, c, n);

        // Nested fork to check that the pthread_atfork protection can work
        // recursively
        fork_pid_nested = fork();
        if (fork_pid_nested == -1) {
            CTEST_ERR("Failed to fork process.");
            exit(1);
        } else if (fork_pid_nested == 0) {
            check_dgemm(a, b, d, c, n);
            exit(0);
        } else {
            check_dgemm(a, b, d, c, n);
            int child_status = 0;
            pid_t wait_pid = wait(&child_status);
            ASSERT_EQUAL(wait_pid, fork_pid_nested);
            ASSERT_EQUAL(0, WEXITSTATUS (child_status));
            exit(0);
        }
    } else {
        check_dgemm(a, b, d, c, n);
        // Wait for the child to finish and check the exit code.
        int child_status = 0;
        pid_t wait_pid = wait(&child_status);
        ASSERT_EQUAL(wait_pid, fork_pid);
        ASSERT_EQUAL(0, WEXITSTATUS (child_status));
    }
}


text:
#include "kepler.h"

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

#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_roots.h>

#define OVERPHASE(x)  ( 2.0 * M_PI * floor ( 0.5 + 0.5 * x * M_1_PI ) )

#define DEBUG 0

/****************************************************************************/

struct delta_params { double ecosmu, esinmu; };

static double delta_f ( double delta, void *params ) {

  struct delta_params * op = ( struct delta_params * ) params;
  double ecosmu = op->ecosmu;
  double esinmu = op->esinmu;

  return ecosmu * sin(delta) + esinmu * cos(delta) - delta;

}

static double delta_df ( double delta, void *params ) {

  struct delta_params * op = ( struct delta_params * ) params;
  double ecosmu = op->ecosmu;
  double esinmu = op->esinmu;

  return ecosmu * cos(delta) - esinmu * sin(delta) - 1.0;

}

static void delta_fdf ( double delta, void *params, double * f, double * df ) {
  struct delta_params * op = ( struct delta_params * ) params;
  double ecosmu = op->ecosmu;
  double esinmu = op->esinmu;
  double cosdelta = cos(delta);
  double sindelta = sin(delta);
  *f  = ecosmu * sindelta + esinmu * cosdelta - delta;
  *df = ecosmu * cosdelta - esinmu * sindelta - 1.0;
}

double kepler_psiofmu ( double mu, double ecc ) {

  int status;
  int iter = 0, max_iter = 100;
  double x0, x, x_expected;
  gsl_function_fdf FDF;
  static gsl_root_fdfsolver * solver;
  static int solverexists = 0;
  struct delta_params params = { ecc * cos(mu), ecc * sin(mu) };

  if ( ! solverexists ) { solverexists = 1;
    solver = gsl_root_fdfsolver_alloc ( gsl_root_fdfsolver_newton );
  }

  if ( DEBUG ) printf( "kepler_psiofmu ( mu = %lg, ecc = %lg )\n", mu, ecc );

  /* early exit make this more tollerant */
  if ( 0.0 == ecc ) return mu;
  if ( 0.0 == mu - OVERPHASE(mu) ) return 0.0;

  FDF.f = &delta_f;
  FDF.df = &delta_df;
  FDF.fdf = &delta_fdf;
  FDF.params = &params;

  x_expected = x = (&params)->esinmu;

  gsl_root_fdfsolver_set (solver, &FDF, x);

  if ( DEBUG ) printf ("using %s method\n", gsl_root_fdfsolver_name (solver));
  if ( DEBUG ) printf ("%-5s %10s %10s\n", "iter", "root", "err");

  do {
    iter++;
    status = gsl_root_fdfsolver_iterate (solver);
    x0 = x;
    x = gsl_root_fdfsolver_root (solver);
    status = gsl_root_test_delta (x, x0, 0, 1e-3);
    if ( DEBUG && status == GSL_SUCCESS ) printf ("Converged:\n");
    if ( DEBUG ) printf ("%5d %10.7f %+10.7f\n", iter, x, x - x_expected);
  } while (status == GSL_CONTINUE && iter < max_iter);

  if ( DEBUG ) printf( "kepler_psiofmu ( mu = %lg, ecc = %lg ) = %lg\n",
    mu, ecc, mu+x );

  return  mu + x;

}

/****************************************************************************/

double kepler_thetaofmu ( double mu, double ecc ) {

  double overmu = OVERPHASE(mu), psi, theta;

  psi = kepler_psiofmu ( mu, ecc );
  theta = 2.0 * atan( sqrt((1.0+ecc)/(1.0-ecc)) * tan(0.5*psi) );

  return  overmu + theta;

}

/****************************************************************************/

double kepler_rv ( double mu, double ecc, double omega ) {

  double psi, theta;

  psi = kepler_psiofmu ( mu, ecc );
  theta = 2.0 * atan( sqrt((1.0+ecc)/(1.0-ecc)) * tan(0.5*psi) );

  return  cos(theta + omega) + ecc*cos(omega);
}

/****************************************************************************/

double kepler_lite ( double mu, double ecc, double omega ) {

  double psi, theta;

  psi = kepler_psiofmu ( mu, ecc );
  theta = 2.0 * atan( sqrt((1.0+ecc)/(1.0-ecc)) * tan(0.5*psi) );

  return  (1.0 - ecc*ecc) * sin(theta + omega) / (1.0 + ecc*cos(theta));
}

/****************************************************************************/

double kepler_sep ( double mu, double ecc, double omega, double inc ) {

  double psi, theta, tmpd;

  psi = kepler_psiofmu ( mu, ecc );
  theta = 2.0 * atan( sqrt((1.0+ecc)/(1.0-ecc)) * tan(0.5*psi) );
  tmpd = sin(theta + omega) * sin(inc);
  tmpd = sqrt(1.0 - tmpd*tmpd);

  return  tmpd * (1.0 - ecc*ecc) / (1.0 + ecc * cos(theta));
}

/****************************************************************************/

ext: c

max_stars_repo_name: mrawls/FDBinary-tools

max_issues_repo_name: mrawls/FDBinary-tools

max_forks_repo_name: mrawls/FDBinary-tools


text:
/*
* -----------------------------------------------------------------
*  pmsr_lib.h
*  Pairwise Mixing Stirred Reactor Library
*  Version: 2.0
*  Last Update: Oct 1, 2019
*
*  This is the header file of a computational library with
*  routines to implement a Pairwise Mixing Stirred Reactor
*  (PMSR) model.
* -----------------------------------------------------------------
*  Programmer: Americo Barbosa da Cunha Junior
*              americo.cunhajr@gmail.com
* -----------------------------------------------------------------
*  Copyright (c) 2019 by Americo Barbosa da Cunha Junior
* -----------------------------------------------------------------
*/

#ifndef __PMSR_LIB_H__
#define __PMSR_LIB_H__

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

/*
* ------------------------------------------------------------
*  Types: struct pmsr, pmsr_wrk
* ------------------------------------------------------------
* This structure contains fields of PMSR model.
*
* last update: Nov 5, 2009
* ------------------------------------------------------------
*/

typedef struct pmsr
{
    int *pmsr_idx;         /* particles index vector */
    gsl_vector **pmsr_ps;  /* particle system        */
} pmsr_wrk;
/*------------------------------------------------------------*/

/*
*------------------------------------------------------------
* function prototypes
*------------------------------------------------------------
*/

void pmsr_title();

int pmsr_input(unsigned long int *seed,
               unsigned int *Np,
               unsigned int *Ndt,
	       double *t0,
               double *delta_t,
               double *tau_res,
               double *tau_mix,
	       double *tau_pair,
               double *Tmax,
               double *Tmin,
               double *atol,
               double *rtol);

pmsr_wrk *pmsr_alloc();

int pmsr_init(int Np,
              int Neq,
              pmsr_wrk *pmsr);

void pmsr_free(int Np,
               void **pmsr_bl);

void pmsr_set_all(int Np,
                  gsl_vector *phi,
                  pmsr_wrk *pmsr);

int pmsr_Nexc(int Np,
              double delta_t,
              double tau_res);

int pmsr_Npair(int Np,
               double delta_t,
               double tau_pair);

void pmsr_meanvar(int Np,
                  int k,
                  gsl_vector **ps,
                  double *mean,
                  double *var);

void pmsr_mixture(int Np,
                  int Neq,
                  double delta_t,
                  double tau_mix,
                  int *idx,
                  gsl_vector **ps);

void pmsr_iops(int Np,
               int Nexc,
               int Npair,
               gsl_rng *r,
               gsl_vector *phi,
               pmsr_wrk *pmsr);

#endif /* __PMSR_LIB_H__ */

ext: h

max_stars_repo_name: americocunhajr/CRFlowLib

max_issues_repo_name: americocunhajr/CRFlowLib

max_forks_repo_name: americocunhajr/CRFlowLib


text:
/*
 Copyright (C) 2002 M. Marques, A. Castro, A. Rubio, G. Bertsch, M. Oliveira

 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, or (at your option)
 any later version.

 This program is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 GNU General Public License for more details.

 You should have received a copy of the GNU General Public License
 along with this program; if not, write to the Free Software
 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
 02110-1301, USA.

*/

#include <config.h>

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

#include <gsl/gsl_sf.h>
#include <gsl/gsl_combination.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

#include "string_f.h"

#include <fortran_types.h>

/* Numerical threshold for oct_bessel_k0 and oct_bessel_k1 */
#define  BESSEL_K_THRES  1.0e2

/* ---------------------- Interface to GSL functions ------------------------ */

/* Special Functions */
double FC_FUNC_(oct_gamma, OCT_GAMMA)
     (const double *x)
{
  return gsl_sf_gamma(*x);
}

double FC_FUNC_(oct_incomplete_gamma, OCT_INCOMPLETE_GAMMA)
     (const double *a, const double *x)
{
  return gsl_sf_gamma_inc_Q(*a, *x);
}

double FC_FUNC_(oct_sph_bessel, OCT_SPH_BESSEL)
     (const fint *l, const double*x)
{
  return gsl_sf_bessel_jl(*l, *x);
}

double FC_FUNC_(oct_bessel, OCT_BESSEL)
     (const fint *n, const double *x)
{
  return gsl_sf_bessel_Jn(*n, *x);
}

double FC_FUNC_(oct_bessel_in, OCT_BESSEL_IN)
     (const fint *n, const double *x)
{
  return gsl_sf_bessel_In(*n, *x);
}

double FC_FUNC_(oct_bessel_j0, OCT_BESSEL_J0)
     (const double *x)
{
  return gsl_sf_bessel_J0(*x);
}

double FC_FUNC_(oct_bessel_j1, OCT_BESSEL_J1)
     (const double *x)
{
  return gsl_sf_bessel_J1(*x);
}

double FC_FUNC_(oct_bessel_k0, OCT_BESSEL_K0)
     (const double *x)
{
  if( *x > BESSEL_K_THRES )
  {
    return 0.0e0;
  } else {
    return gsl_sf_bessel_K0(*x);
  }
}

double FC_FUNC_(oct_bessel_k1, OCT_BESSEL_K1)
     (const double *x)
{
  if( *x > BESSEL_K_THRES )
  {
    return 0.0e0;
  } else {
    return gsl_sf_bessel_K1(*x);
  }
}

/* the GSL headers specify double x, not const double x */
double FC_FUNC_(oct_erfc, OCT_ERFC)
     (const double *x)
{
  /* avoid floating invalids in the asymptotic limit */
  if(*x >  20.0) return  0.0;
  if(*x < -20.0) return  2.0;
  /* otherwise call gsl */
  return gsl_sf_erfc(*x);
}

/* the GSL headers specify double x, not const double x */
double FC_FUNC_(oct_erf, OCT_ERF)
     (const double *x)
{
  /* avoid floating invalids in the asymptotic limit */
  if(*x >  20.0) return  1.0;
  if(*x < -20.0) return -1.0;
  /* otherwise call gsl */
  return gsl_sf_erf(*x);
}

double FC_FUNC_(oct_legendre_sphplm, OCT_LEGENDRE_SPHPLM)
     (const fint *l, const int *m, const double *x)
{
  return gsl_sf_legendre_sphPlm(*l, *m, *x);
}

/* generalized Laguerre polynomials */
double FC_FUNC_(oct_sf_laguerre_n, OCT_SF_LAGUERRE_N)
     (const int *n, const double *a, const double *x)
{
  return gsl_sf_laguerre_n(*n, *a, *x);
}

/* Combinations */

void FC_FUNC_(oct_combination_init, OCT_COMBINATION_INIT)
     (gsl_combination **c, const fint *n, const fint *k)
{
  *c = gsl_combination_calloc (*n, *k);
}

void FC_FUNC_(oct_combination_next, OCT_COMBINATION_NEXT)
     (gsl_combination **c, fint *next)
{
  *next = gsl_combination_next (((gsl_combination *)(*c)));
}

void FC_FUNC_(oct_get_combination, OCT_GET_COMBINATION)
     (gsl_combination **c, fint *comb)
{
  int i;
  for (i=0;i< ((gsl_combination *)(*c))->k; i++) {
    comb[i]=(fint)((gsl_combination *)(*c))->data[i];
  }
}

void FC_FUNC_(oct_combination_end, OCT_COMBINATION_END)
     (gsl_combination **c)
{
  gsl_combination_free (((gsl_combination *)(*c)));
}

/* Random Number Generation */
void FC_FUNC_(oct_ran_init, OCT_RAN_INIT)
     (gsl_rng **r)
{
  gsl_rng_env_setup();
  *r = gsl_rng_alloc(gsl_rng_default);
}

void FC_FUNC_(oct_ran_end, OCT_RAN_END)
     (gsl_rng **r)
{
  gsl_rng_free(*r);
}

/* Random Number Distributions */
double FC_FUNC_(oct_ran_gaussian, OCT_RAN_GAUSSIAN)
     (const gsl_rng **r, const double *sigma)
{
  return gsl_ran_gaussian(*r, *sigma);
}

double FC_FUNC_(oct_ran_flat, OCT_RAN_FLAT)
     (const gsl_rng **r, const double *a, const double *b)
{
  return gsl_ran_flat(*r, *a, *b);
}

void FC_FUNC_(oct_strerror, OCT_STRERROR)
     (const fint *err, STR_F_TYPE res STR_ARG1)
{
  const char *c;

  c = gsl_strerror(*err);
  TO_F_STR1(c, res);
}

max_stars_repo_name: shunsuke-sato/octopus

max_issues_repo_name: shunsuke-sato/octopus

max_forks_repo_name: shunsuke-sato/octopus


text:
/* rng/ranlxs.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 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.
 */

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_rng.h>

/* This is an implementation of M. Luescher's second generation
   version of the RANLUX generator.

   Thanks to Martin Luescher for providing information on this
   generator.

 */

static unsigned long int ranlxs_get (void *vstate);
static inline double ranlxs_get_double (void *vstate);
static void ranlxs_set_lux (void *state, unsigned long int s, unsigned int luxury);
static void ranlxs0_set (void *state, unsigned long int s);
static void ranlxs1_set (void *state, unsigned long int s);
static void ranlxs2_set (void *state, unsigned long int s);

static const int next[12] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 0};
static const int snext[24] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,
                              14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 0};

static const double sbase = 16777216.0;		/* 2^24 */
static const double sone_bit = 1.0 / 16777216.0;	/* 1/2^24 */
static const double one_bit = 1.0 / 281474976710656.0;	/* 1/2^48 */

static const double shift = 268435456.0;	/* 2^28 */

#define RANLUX_STEP(x1,x2,i1,i2,i3)      \
          x1=xdbl[i1] - xdbl[i2];        \
          if (x2 < 0)                    \
          {                              \
            x1-=one_bit;                 \
            x2+=1;                       \
          }                              \
          xdbl[i3]=x2

typedef struct
  {
    double xdbl[12], ydbl[12];  /* doubles first so they are 8-byte aligned */
    double carry;
    float xflt[24];
    unsigned int ir;
    unsigned int jr;
    unsigned int is;
    unsigned int is_old;
    unsigned int pr;
  }
ranlxs_state_t;

static void increment_state (ranlxs_state_t * state);

static void
increment_state (ranlxs_state_t * state)
{
  int k, kmax, m;
  double x, y1, y2, y3;

  float *xflt = state->xflt;
  double *xdbl = state->xdbl;
  double *ydbl = state->ydbl;
  double carry = state->carry;
  unsigned int ir = state->ir;
  unsigned int jr = state->jr;

  for (k = 0; ir > 0; ++k)
    {
      y1 = xdbl[jr] - xdbl[ir];
      y2 = y1 - carry;
      if (y2 < 0)
	{
	  carry = one_bit;
	  y2 += 1;
	}
      else
	{
	  carry = 0;
	}
      xdbl[ir] = y2;
      ir = next[ir];
      jr = next[jr];
    }

  kmax = state->pr - 12;

  for (; k <= kmax; k += 12)
    {
      y1 = xdbl[7] - xdbl[0];
      y1 -= carry;

      RANLUX_STEP (y2, y1, 8, 1, 0);
      RANLUX_STEP (y3, y2, 9, 2, 1);
      RANLUX_STEP (y1, y3, 10, 3, 2);
      RANLUX_STEP (y2, y1, 11, 4, 3);
      RANLUX_STEP (y3, y2, 0, 5, 4);
      RANLUX_STEP (y1, y3, 1, 6, 5);
      RANLUX_STEP (y2, y1, 2, 7, 6);
      RANLUX_STEP (y3, y2, 3, 8, 7);
      RANLUX_STEP (y1, y3, 4, 9, 8);
      RANLUX_STEP (y2, y1, 5, 10, 9);
      RANLUX_STEP (y3, y2, 6, 11, 10);

      if (y3 < 0)
	{
	  carry = one_bit;
	  y3 += 1;
	}
      else
	{
	  carry = 0;
	}
      xdbl[11] = y3;
    }

  kmax = state->pr;

  for (; k < kmax; ++k)
    {
      y1 = xdbl[jr] - xdbl[ir];
      y2 = y1 - carry;
      if (y2 < 0)
	{
	  carry = one_bit;
	  y2 += 1;
	}
      else
	{
	  carry = 0;
	}
      xdbl[ir] = y2;
      ydbl[ir] = y2 + shift;
      ir = next[ir];
      jr = next[jr];
    }

  ydbl[ir] = xdbl[ir] + shift;

  for (k = next[ir]; k > 0;)
    {
      ydbl[k] = xdbl[k] + shift;
      k = next[k];
    }

  for (k = 0, m = 0; k < 12; ++k)
    {
      x = xdbl[k];
      y2 = ydbl[k] - shift;
      if (y2 > x)
	y2 -= sone_bit;
      y1 = (x - y2) * sbase;

      xflt[m++] = (float) y1;
      xflt[m++] = (float) y2;
    }

  state->ir = ir;
  state->is = 2 * ir;
  state->is_old = 2 * ir;
  state->jr = jr;
  state->carry = carry;
}

static inline double
ranlxs_get_double (void *vstate)
{
  ranlxs_state_t *state = (ranlxs_state_t *) vstate;

  const unsigned int is = snext[state->is];

  state->is = is;

  if (is == state->is_old)
    increment_state (state);

  return state->xflt[state->is];
}

static unsigned long int
ranlxs_get (void *vstate)
{
  return ranlxs_get_double (vstate) * 16777216.0;	/* 2^24 */
}

static void
ranlxs_set_lux (void *vstate, unsigned long int s, unsigned int luxury)
{
  ranlxs_state_t *state = (ranlxs_state_t *) vstate;

  int ibit, jbit, i, k, m, xbit[31];
  double x, y;

  long int seed;

  if (s == 0)
    s = 1;			/* default seed is 1 */

  seed = s;

  i = seed & 0xFFFFFFFFUL;

  for (k = 0; k < 31; ++k)
    {
      xbit[k] = i % 2;
      i /= 2;
    }

  ibit = 0;
  jbit = 18;

  for (k = 0; k < 12; ++k)
    {
      x = 0;

      for (m = 1; m <= 48; ++m)
	{
	  y = (double) xbit[ibit];
	  x += x + y;
	  xbit[ibit] = (xbit[ibit] + xbit[jbit]) % 2;
	  ibit = (ibit + 1) % 31;
	  jbit = (jbit + 1) % 31;
	}
      state->xdbl[k] = one_bit * x;
    }

  state->carry = 0;
  state->ir = 0;
  state->jr = 7;
  state->is = 23;
  state->is_old = 0;
  state->pr = luxury;
}

static void
ranlxs0_set (void *vstate, unsigned long int s)
{
  ranlxs_set_lux (vstate, s, 109);
}

void
ranlxs1_set (void *vstate, unsigned long int s)
{
  ranlxs_set_lux (vstate, s, 202);
}

static void
ranlxs2_set (void *vstate, unsigned long int s)
{
  ranlxs_set_lux (vstate, s, 397);
}

static const gsl_rng_type ranlxs0_type =
{"ranlxs0",			/* name */
 0x00ffffffUL,			/* RAND_MAX */
 0,				/* RAND_MIN */
 sizeof (ranlxs_state_t),
 &ranlxs0_set,
 &ranlxs_get,
 &ranlxs_get_double};

static const gsl_rng_type ranlxs1_type =
{"ranlxs1",			/* name */
 0x00ffffffUL,			/* RAND_MAX */
 0,				/* RAND_MIN */
 sizeof (ranlxs_state_t),
 &ranlxs1_set,
 &ranlxs_get,
 &ranlxs_get_double};

static const gsl_rng_type ranlxs2_type =
{"ranlxs2",			/* name */
 0x00ffffffUL,			/* RAND_MAX */
 0,				/* RAND_MIN */
 sizeof (ranlxs_state_t),
 &ranlxs2_set,
 &ranlxs_get,
 &ranlxs_get_double};

const gsl_rng_type *gsl_rng_ranlxs0 = &ranlxs0_type;
const gsl_rng_type *gsl_rng_ranlxs1 = &ranlxs1_type;
const gsl_rng_type *gsl_rng_ranlxs2 = &ranlxs2_type;


text:
/*
  code to do weak lensing power spectrum computations

  Matthew Becker, UofC 2011
*/

#include <gsl/gsl_spline.h>

#ifndef _WEAKLENS_
#define _WEAKLENS_

typedef struct {
  int wlNum;
  double lmin;
  double lmax;
  double tmin;
  double tmax;
} weaklensData;

extern weaklensData wlData;

typedef struct {
  int initFlag;
  int currCosmoNum;
  int currWLNum;
  double zs1;
  double zs2;
  double chis1;
  double chis2;
  double chiLim;
  double sn;
  double ell; //temp var for integrals
  gsl_spline *spline;
  gsl_interp_accel *accel;
} _lensPowerSpectra;

typedef _lensPowerSpectra *lensPowerSpectra;

typedef struct {
  int initFlag;
  int currCosmoNum;
  int currWLNum;
  lensPowerSpectra lps;
  double theta; //temp var for integrals
  gsl_spline *splineP;
  gsl_interp_accel *accelP;
  gsl_spline *splineM;
  gsl_interp_accel *accelM;
} _lensCorrFunc;

typedef _lensCorrFunc * lensCorrFunc;

//functions to compute lensing power and cross spectra
double nonlinear_powspec_for_lens(double k, double a);
double lens_power_spectrum(double ell, lensPowerSpectra lps);
lensPowerSpectra init_lens_power_spectrum(double zs1, double zs2);
void free_lens_power_spectrum(lensPowerSpectra lps);

//functions to compute lens corr func.
double lens_corr_func_minus(double theta, lensCorrFunc lcf);
double lens_corr_func_plus(double theta, lensCorrFunc lcf);
lensCorrFunc init_lens_corr_func(lensPowerSpectra lps);
void free_lens_corr_func(lensCorrFunc lcf);

#endif /* _WEAKLENS_ */

max_stars_repo_name: beckermr/cosmocalc

max_issues_repo_name: beckermr/cosmocalc

max_forks_repo_name: beckermr/cosmocalc


text:
/*
* -----------------------------------------------------------------
*  Utilities Library --- util_lib.h
*  Version: 1.6180
*  Date: Feb 19, 2010
* -----------------------------------------------------------------
*  Programmer: Americo Barbosa da Cunha Junior
*              americo.cunhajr@gmail.com
* -----------------------------------------------------------------
*  Copyright (c) 2010 by Americo Barbosa da Cunha Junior
*
*  This program is free software: you can redistribute it and/or
*  modify it under the terms of the GNU General Public License as
*  published by the Free Software Foundation, either version 3 of
*  the License, or (at your option) any later version.
*
*  This program is distributed in the hope that it will be useful,
*  but WITHOUT ANY WARRANTY; without even the implied warranty of
*  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
*  GNU General Public License for more details.
*
*  A copy of the GNU General Public License is available in
*  LICENSE.txt or http://www.gnu.org/licenses/.
* -----------------------------------------------------------------
*  This is the header file for a library with miscelaneous
*  functions witch have a lot of applications.
* -----------------------------------------------------------------
*/

#ifndef __UTIL_LIB_H__
#define __UTIL_LIB_H__

#include <time.h>

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

/*
*------------------------------------------------------------
* function prototypes
*------------------------------------------------------------
*/

int gsl_fprint_vector(FILE* file,
                      gsl_vector *v);

int gsl_print_vector(gsl_vector *v);

int gsl_print_matrix_line(unsigned int i,
                          gsl_matrix *M);

int gsl_print_matrix_column(unsigned int j,
                            gsl_matrix *M);

int gsl_print_matrix(gsl_matrix *M);

int gsl_rand_vector(gsl_rng *r,
                    gsl_vector *v);

int gsl_rand_matrix(gsl_rng *r,
                    gsl_matrix *M);

int gsl_diagonal_matrix(gsl_vector *v,
                        gsl_matrix *D);

void util_time_used(clock_t cpu_start,
                    clock_t cpu_end,
			        time_t wall_start,
                    time_t wall_end);

#endif /* __UTIL_LIB_H__ */


text:
/* specfunc/gegenbauer.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_gegenbauer.h"

#include "error.h"

/* See: [Thompson, Atlas for Computing Mathematical Functions] */

int
gsl_sf_gegenpoly_1_e(double lambda, double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

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

int
gsl_sf_gegenpoly_2_e(double lambda, double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(lambda == 0.0) {
    const double txx = 2.0*x*x;
    result->val  = -1.0 + txx;
    result->err  = 2.0 * GSL_DBL_EPSILON * fabs(txx);
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    result->val = lambda*(-1.0 + 2.0*(1.0+lambda)*x*x);
    result->err = GSL_DBL_EPSILON * (2.0 * fabs(result->val) + fabs(lambda));
    return GSL_SUCCESS;
  }
}

int
gsl_sf_gegenpoly_3_e(double lambda, double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(lambda == 0.0) {
    result->val = x*(-2.0 + 4.0/3.0*x*x);
    result->err = GSL_DBL_EPSILON * (2.0 * fabs(result->val) + fabs(x));
    return GSL_SUCCESS;
  }
  else {
    double c = 4.0 + lambda*(6.0 + 2.0*lambda);
    result->val = 2.0*lambda * x * ( -1.0 - lambda + c*x*x/3.0 );
    result->err = GSL_DBL_EPSILON * (2.0 * fabs(result->val) + fabs(lambda * x));
    return GSL_SUCCESS;
  }
}

int
gsl_sf_gegenpoly_n_e(int n, double lambda, double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(lambda <= -0.5 || n < 0) {
    DOMAIN_ERROR(result);
  }
  else if(n == 0) {
    result->val = 1.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else if(n == 1) {
    return gsl_sf_gegenpoly_1_e(lambda, x, result);
  }
  else if(n == 2) {
    return gsl_sf_gegenpoly_2_e(lambda, x, result);
  }
  else if(n == 3) {
    return gsl_sf_gegenpoly_3_e(lambda, x, result);
  }
  else {
    if(lambda == 0.0 && (x >= -1.0 || x <= 1.0)) {
      /* 2 T_n(x)/n */
      const double z = n * acos(x);
      result->val = 2.0 * cos(z) / n;
      result->err = 2.0 * GSL_DBL_EPSILON * fabs(z * result->val);
      return GSL_SUCCESS;
    }
    else {
      int k;
      gsl_sf_result g2;
      gsl_sf_result g3;
      int stat_g2 = gsl_sf_gegenpoly_2_e(lambda, x, &g2);
      int stat_g3 = gsl_sf_gegenpoly_3_e(lambda, x, &g3);
      int stat_g  = GSL_ERROR_SELECT_2(stat_g2, stat_g3);
      double gkm2 = g2.val;
      double gkm1 = g3.val;
      double gk = 0.0;
      for(k=4; k<=n; k++) {
        gk = (2.0*(k+lambda-1.0)*x*gkm1 - (k+2.0*lambda-2.0)*gkm2) / k;
	gkm2 = gkm1;
	gkm1 = gk;
      }
      result->val = gk;
      result->err = 2.0 * GSL_DBL_EPSILON * 0.5 * n * fabs(gk);
      return stat_g;
    }
  }
}

int
gsl_sf_gegenpoly_array(int nmax, double lambda, double x, double * result_array)
{
  int k;

  /* CHECK_POINTER(result_array) */

  if(lambda <= -0.5 || nmax < 0) {
    GSL_ERROR("domain error", GSL_EDOM);
  }

  /* n == 0 */
  result_array[0] = 1.0;
  if(nmax == 0) return GSL_SUCCESS;

  /* n == 1 */
  if(lambda == 0.0)
    result_array[1] = 2.0*x;
  else
    result_array[1] = 2.0*lambda*x;

  /* n <= nmax */
  for(k=2; k<=nmax; k++) {
    double term1 = 2.0*(k+lambda-1.0) * x * result_array[k-1];
    double term2 = (k+2.0*lambda-2.0)	  * result_array[k-2];
    result_array[k] = (term1 - term2) / k;
  }

  return GSL_SUCCESS;
}

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

#include "eval.h"

double gsl_sf_gegenpoly_1(double lambda, double x)
{
  EVAL_RESULT(gsl_sf_gegenpoly_1_e(lambda, x, &result));
}

double gsl_sf_gegenpoly_2(double lambda, double x)
{
  EVAL_RESULT(gsl_sf_gegenpoly_2_e(lambda, x, &result));
}

double gsl_sf_gegenpoly_3(double lambda, double x)
{
  EVAL_RESULT(gsl_sf_gegenpoly_3_e(lambda, x, &result));
}

double gsl_sf_gegenpoly_n(int n, double lambda, double x)
{
  EVAL_RESULT(gsl_sf_gegenpoly_n_e(n, lambda, x, &result));
}


text:
#include <gsl/gsl_vector.h>
#include <gsl/gsl_multiroots.h>

max_stars_repo_name: datamole-ai/gomez

max_issues_repo_name: datamole-ai/gomez

max_forks_repo_name: datamole-ai/gomez


text:
#ifndef SIGNAL_H_
#define SIGNAL_H_

#include <gsl/gsl_statistics_double.h>
#include <complex>
#include <fftw3.h>
#include <math.h>
#include <iostream>
#include "fourier_transform.h"
#include "array.h"
#include "rng.h"
#include <assert.h>

// this is slowly turning into sth that should maybe rather be called
// "stochastic processes" rather than signal

namespace neurophys {

class SpectrumDescriptor
{
public:
    virtual double get(const double f) const = 0;
};

class WhiteSpectrum: public SpectrumDescriptor
{
public:
    WhiteSpectrum(const double s0): s0_(s0) {}
    virtual double get(const double f) const { return s0_; }
private:
    const double s0_;
};

class BandLimitedFlatSpectrum: public SpectrumDescriptor
{
public:
    BandLimitedFlatSpectrum(const double s0, const double f0, const double f1):
        s0_(s0), f0_(f0), f1_(f1)
        {}
    virtual double get(const double f) const
    {
        return (f >= f0_ && f < f1_) ? s0_ : 0.;
    }
private:
    const double s0_;
    const double f0_;
    const double f1_;
};

class PowerLawSpectrum: public SpectrumDescriptor
{
public:
    PowerLawSpectrum(const double expo, const double f0, const double f1):
        expo_(expo), f0_(f0), f1_(f1), A_(0)
    {
        // normalization, so that variance = 1
        A_ = 1./2 * (1.-expo) / (pow(f1,1.-expo) - expo*pow(f0,1.-expo));
    }
    virtual double get(const double f) const
    {
        if (f < f0_) return A_ * pow(f0_, -expo_);
        else if (f < f1_) return A_ * pow(f, -expo_);
        else return 0.;
    }
private:
    const double expo_;
    const double f0_;
    const double f1_;
    double A_;
};

namespace signal {

Array<double> generate_gaussian(RNG& rng,
        const C2RFourierTransform& c2rft, const SpectrumDescriptor& spec,
        const size_t N, const double T);

};

}

#endif /* SIGNAL_H */


text:
/*
  @copyright Russell Standish 2000-2013
  @author Russell Standish
  This file is part of EcoLab

  Open source licensed under the MIT license. See LICENSE for details.
*/

#include <classdesc_access.h>
#include <pack_stl.h>
#include <analysis.h>

#if defined(GNUSL)
#include <gsl/gsl_errno.h>
#include <gsl/gsl_odeiv2.h>
#endif

using namespace ecolab;
typedef map<NautyRep,double> archetype;
using namespace classdesc;
using array_ns::array;

struct Gsl_matrix;
struct GaussianMarkovModel
{
  /// lags and replications used calculating causal density and ilk
  int lags;
  Gsl_matrix* runData;
  typedef ConcreteGraph<DiGraph> G;
  G& g;
  gaussrand gran;

  GaussianMarkovModel(G& g): lags(10), runData(NULL), g(g) {}
  ~GaussianMarkovModel();
  /// run the model for lags+repetitions steps to populate the matrices
  void run(TCL_args repetitions);
  /// output the rundata into a file: each node has separate row, each
  /// sample column
  void outputRunData(TCL_args);
  /// calculate the causal density for a Gaussian Markov process
  /// running on the network.
  double causalDensity();
  double bipartitionCausalDensity();
};

// RAII data structure for GSL Runge-Kutta algorithm
struct RKdata
{
#if defined(GNUSL)
  gsl_odeiv2_system sys;
  gsl_odeiv2_driver* driver;
  template <class DynSysGen>
  RKdata(const DynSysGen& dynSys) {
    sys.function=DynSysGen::RKfunction;
    sys.dimension=dynSys.dim;
    sys.params=(void*)&dynSys;
    const gsl_odeiv2_step_type* stepper=gsl_odeiv2_step_rkf45;
    driver = gsl_odeiv2_driver_alloc_y_new
        (&sys, stepper, dynSys.stepMax, dynSys.epsAbs,
         dynSys.epsRel);
      gsl_odeiv2_driver_set_hmax(driver, dynSys.stepMax);
      gsl_odeiv2_driver_set_hmin(driver, dynSys.stepMin);
    }
    ~RKdata() {gsl_odeiv2_driver_free(driver);}
#else
  int sys, driver; // dummies, to keep classdesc happy
#endif
  };

// dynamic system network generators
template <class DynSys>
struct DynSysGen: public DynSys, public NetworkFromTimeSeries
{
  static int RKfunction(double t, const double x[], double f[], void *params)
  {
#if defined(GNUSL)
    if (params==NULL) return GSL_EBADFUNC;
    try
      {
        ((const DynSys*)params)->f(t,x,f);
      }
    catch (std::exception& e)
      {
        Tcl_AppendResult(interp(),e.what(),NULL);
        Tcl_AppendResult(interp(),"\n",NULL);
        return GSL_EBADFUNC;
      }
    return GSL_SUCCESS;
#else
    return 0;
#endif
  }
  // RK parameters
  double stepMin; // minimum step size
  double stepMax; // maximum step size
  double epsAbs;  // absolute error
  double epsRel;  // relative error
  DynSysGen(): stepMin(1e-5), stepMax(1e-3), epsAbs(1e-3), epsRel(1e-2)
  {
    resize(DynSys::dim);
    n.resize(DynSys::dim,100);
  }

  void fillNet(int nsteps) {
#if defined(GNUSL)
    RKdata rkData(*this);
    vector<double> x(DynSys::dim);
    for (size_t i=0; i<x.size(); ++i)
      x[i]=uni.rand();
    double t=0;
    gsl_odeiv2_driver_set_nmax(rkData.driver, 1);
    for (int i=0; i<nsteps; ++i)
      {
        gsl_odeiv2_driver_apply(rkData.driver,&t,1,&x[0]);
        for (size_t j=0; j<size(); ++j)
          (*this)[j]<<=x[j];
      }
    constructNet();
#endif
  }
  double complexity() const {return ::complexity(net);}
  urand uni;
  void random_rewire() {::random_rewire(net,uni);}
};

// example dynamical systems defined here
struct Lorenz
{
  static const unsigned dim=3;
  double sigma, rho, beta;
  Lorenz(): sigma(10), rho(28), beta(8/3.0) {}
  void f(double t, const double x[3], double y[3]) const
  {
    y[0]=sigma*(x[1]-x[0]);
    y[1]=x[0]*(rho-x[2])-x[1];
    y[2]=x[0]*x[1]-beta*x[2];
  }
};

struct HenonHeiles
{
  static const unsigned dim=4;
  double lambda;
  HenonHeiles(): lambda(1) {}
  void f(double t, const double x[4], double y[4]) const
  {
    y[0]=x[1];
    y[1]=-x[0]-2*lambda*x[0]*x[2];
    y[2]=x[3];
    y[3]=-x[2]-lambda*(x[0]*x[0]-x[2]*x[2]);
  }
};

// small 1D cellular automata
struct Small1DCA
{
  unsigned rule; ///< rule in 1..255
  unsigned nCells; ///< no cells < 0..8*sizeof(state)
  unsigned long state;
  /// 3 bit neighbourhood of cell i
  unsigned nbhd(unsigned i)
  {
    if (i==0)
      return ((state&3)<<1) | ((state>>(nCells-1)) & 1);
    else if (i==nCells-1)
      return ((state>>(i-1))&3) | ((state&1)<<2);
    else
      return state>>(i-1)&7;
  }
  // ratio of 1 bits out of 8
  double lambda() {
    int sum=0;
    for (int i=0; i<8; ++i)
      sum+=(rule>>i)&1;
    return sum/8.0;
  }
  void step() {
    unsigned long newState=0;
    for (unsigned i=0; i<nCells; ++i)
      {
        unsigned n=nbhd(i);
        newState |= ((rule>>n) & 1) << i;
      }
    state=newState;
  }
};

struct Larger1DCA
{
  // number of states a cell can take, and neighbourhood size (always odd)
  unsigned nStates;
  int nbhdSz;
  vector<unsigned> rule,cells;
  urand uni;
  size_t nbhd(unsigned i)
  {
    size_t r=0;
    for (int j=-nbhdSz/2; j<=nbhdSz/2; ++j)
      {
        int ii=i+j;
        if (ii<0) ii+=cells.size();
        if (ii>int(cells.size())-1) ii-=cells.size();
        r*=nStates;
        r+=cells[ii];
      }
    return r;
  }
  // ratio of non-quiescent states in rule to total
  double lambda() {
    int sum=0;
    for (unsigned i=0; i<rule.size(); ++i)
      sum+=rule[i]>0;
    return double(sum)/rule.size();
  }
  // resize rule vector, given nStates and nbhdSz
  void resizeRule() {
    rule.resize(pow(nStates,nbhdSz));
  }
  // installs a random rule matching a given lambda
  void initRandomRule(double lambda) {
    resizeRule();
    for (size_t i=0; i<rule.size(); ++i)
      if (uni.rand()>lambda)
        rule[i]=0;
      else
        rule[i]=(nStates-1)*uni.rand()+1;
  }
  void initRandom() {
    for (size_t i=0; i<cells.size(); ++i)
      cells[i]=nStates*uni.rand();
  }
  void step() {
    vector<unsigned> newState(cells.size());
    for (unsigned i=0; i<cells.size(); ++i)
      newState[i] = rule[nbhd(i)];
    cells.swap(newState);
  }

};

template <class T>
struct NodeMap: map<T,unsigned>
{
  unsigned nextId;
  NodeMap(): nextId(0) {}
  unsigned operator[](const T& x) {
    typename map<T,unsigned>::iterator i=this->find(x);
    if (i!=this->end())
      return i->second;
    else
      return map<T,unsigned>::operator[](x)=nextId++;
  }
};

struct NetworkFromSmall1DCA:
  public Small1DCA, public NetworkFromTimeSeries
{
  void fillNet(int nsteps, int transient) {
    resize(1); back().clear();
    state=pow(2,nCells)*uni.rand();
    NodeMap<unsigned long> nodeMap;

    // throw away the intial transient
    for (int i=0; i<transient; ++i) step();
    for (int i=0; i<nsteps; ++i)
      {
        step();
        back()<<=nodeMap[state];
 //       if (nodeMap.nextId>10)
 //         cout << back() << " ";
      }
    //   cout<<endl;
    n.resize(1); n[0]=nodeMap.nextId;
    constructNet();
  }

  double complexity() const {return ::complexity(net);}
  urand uni;
  void random_rewire() {::random_rewire(net,uni);}

};

struct NetworkFromLarger1DCA:
  public Larger1DCA, public NetworkFromTimeSeries
{
  void fillNet(int nsteps, int transient) {
    resize(1); back().clear();
    initRandom();
    // throw away the intial transient
    for (int i=0; i<transient; ++i) step();
    // map of state to node tag
    NodeMap<vector<unsigned> > nodeMap;
    unsigned tagid=1;
    for (int i=0; i<nsteps; ++i)
      {
        step();
        back()<<=nodeMap[cells];
      }
    n.resize(1); n[0]=nodeMap.nextId;
    constructNet();
  }

  double complexity() const {return ::complexity(net);}
  urand uni;
  void random_rewire() {::random_rewire(net,uni);}

};

class netc_t
{
  archetype archc;
  Iterator<archetype> ic;
  bool fill_curr();
  CLASSDESC_ACCESS(netc_t);
public:
  netc_t(): weight_dist("one"), gProcess(g) {}
  unsigned currc;
  void reset_arch() {archc.clear();}
  bool first();
  bool next();
  bool last() {return ic==archc.end();}
  unsigned size() {return archc.size();}
  double ODcomplexity() {return ::ODcomplexity(g);}
  double MA() {return ::MA(g);}
  void exhaust_search(TCL_args);
  //  void gen_uniq_rand(TCL_args);

  ConcreteGraph<DiGraph> g;
  Degrees degrees;
  void update_degrees() {degrees=::degrees(g);}
  urand uni;
  double complexity() {return ::complexity(g);}
  double lnomega() {BitRep eCanon; return canonical(eCanon,g);}
  double nauty_complexity() {return ::complexity(NautyRep(g));}
  double nauty_lnomega() {return NautyRep(g).lnomega();}
  double saucy_lnomega() {return ::saucy_lnomega(g);}

  /** race algorithms against each other in separate processes. Best
      run on a multicore machine with as many or more cores than there
      are algorithms (currently 4).

      @param time (in seconds) limiting the race
      @return the number of te winning algorithm, or 0 if the race times out
  */
  int automorphism_race(TCL_args);
  /**
     races variants of the bliss algorithm with different values of
     the splitting heuristic
     @param time (in seconds) limiting the  race
     @return the winning splitting heuristic, or -1 if the race times out
  */
  int bliss_race(TCL_args);
  /// returns a string representation of the algorithm given by the
  /// integer argument
  string automorphism_algorithm_name(TCL_args);

  TCL_obj_ref<random_gen> weight_dist;
  array<float> weights() {
    array<float> r;
    for (ConcreteGraph<DiGraph>::const_iterator e=g.begin(); e!=g.end(); ++e)
      r<<=e->weight;
    return r;
  }
  void genRandom(TCL_args args) {
    ErdosRenyi_generate(g,args[0],args[1],uni,*weight_dist);}
  void prefAttach(TCL_args args) {
    TCL_obj_ref<random_gen> indegree_dist;
    indegree_dist.set(args[1]);
    PreferentialAttach_generate(g,args[0],uni,*indegree_dist,*weight_dist);
  }
  void random_rewire() {::random_rewire(g,uni);}
  /// generate a fully connected directed graph, with each node having
  /// same indegree and outdegree. Only works if graph degree is odd.
  void genFullDirected(TCL_args args);
  bool testAgainstNauty();

  GaussianMarkovModel gProcess;

  DynSysGen<Lorenz> lorenzGen;
  DynSysGen<HenonHeiles> henonHeilesGen;
  NetworkFromSmall1DCA small1DCA;
  NetworkFromLarger1DCA large1DCA;
};

max_stars_repo_name: digiperfect/ecolab

max_issues_repo_name: digiperfect/ecolab

max_forks_repo_name: digiperfect/ecolab


text:
#ifndef LUM_H
#define LUM_H

#include <LHAPDF/LHAPDF.h>
#include <cmath>
#include <iostream>
#include <gsl/gsl_chebyshev.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_sf_log.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_sf_zeta.h>
#include <gsl/gsl_sf_psi.h>
#include <gsl/gsl_sf_dilog.h>
#include <algorithm>
#include <gsl/gsl_integration.h>

class Luminosity {
	public:
		typedef std::complex<long double> dcomplex;
		typedef LHAPDF::PDF* PDF_ptr;

		Luminosity(PDF_ptr thePDF, double MUF, unsigned short NF = 5, std::size_t order=20);
		Luminosity(const Luminosity& Lumi);
		virtual ~Luminosity();

		void Cheb_Lum(double );
		void c_large();

		dcomplex CLum_N(dcomplex N, unsigned channel) const;
		long double CLum_x(long double z, unsigned channel) const;

		dcomplex CLum_gg_N(dcomplex N) const {	return CLum_N(N,0);	}
		long double CLum_gg_x(long double z) const { return CLum_x(z,0);  }

		dcomplex CLum_gq_N(dcomplex N) const {	return CLum_N(N,1);	}
		long double CLum_gq_x(long double z) const { return CLum_x(z,1);  }

		dcomplex CLum_qqbar_N(dcomplex N) const {	return CLum_N(N,2);	}
		long double CLum_qqbar_x(long double z) const { return CLum_x(z,2);  }

		dcomplex CLum_qq_N(dcomplex N) const {	return CLum_N(N,3);	}
		long double CLum_qq_x(long double z) const { return CLum_x(z,3);  }

		dcomplex CLum_qqbarprime_N(dcomplex N) const {	return CLum_N(N,4);	}
		long double CLum_qqbarprime_x(long double z) const {	return CLum_x(z,4);	}

		dcomplex CLum_qqprime_N(dcomplex N) const {	return CLum_N(N,5);	}
		long double CLum_qqprime_x(long double z) const { return CLum_x(z,5);  }

		double xfxQ(int channel, double x){
		  return (_PDF->xfxQ(channel,x,_Muf)/x);
		}

		inline PDF_ptr get_PDF(){
		  return _PDF;
		}
		inline double get_muf(){
		  return _Muf;
		}
		inline size_t get_n(){
		  return _n;
		}

		struct par_struct{
			double Muf;
			PDF_ptr PDF;
			double t;
			unsigned int Nf;
		};

		void setNf(unsigned short NF){
		  if ( (NF < 1) || (NF > 6) ){
		    std::cout << "ERROR: INCONGRUENT NUMBER OF FLAVOUR (1-6 ACCEPTED)" << std::endl;
		    NF=5.;
		  }
		  _Nf=NF;
		}
		//void cheb_approx(double , double);

	private:
		size_t _n;	//Order of Chebishev approximation
		std::array< gsl_cheb_series*,3>  Chebseries;
		std::array< std::vector<double> ,3>  C_larges;
//		std::array< gsl_cheb_series*,6>  Chebseries;
//		std::array< std::vector<double> ,6>  C_larges;

		std::vector< std::vector< double > > _T;
		LHAPDF::PDF* _PDF;

		double _Muf;
		unsigned short _Nf;

};

#endif

max_stars_repo_name: gvita/RHEHpt

max_issues_repo_name: gvita/RHEHpt

max_forks_repo_name: gvita/RHEHpt


text:
//c++ network  author : liqi
//Nangjing University of Posts and Telecommunications
//date 2017.5.21,20:27
#ifndef NETWORK_H
#define NETWORK_H
#include "opencv2/imgproc/imgproc.hpp"
#include "opencv2/highgui/highgui.hpp"
#include <algorithm>
#include <stdlib.h>
#include <memory.h>
#include <fstream>
#include <cstring>
#include <cblas.h>
#include <string>
#include <math.h>
#include "pBox.h"
using namespace cv;

void addbias(struct pBox *pbox, mydataFmt *pbias);
void image2Matrix(const Mat &image, const struct pBox *pbox);
void featurePad(const pBox *pbox, const pBox *outpBox, const int pad);
void feature2Matrix(const pBox *pbox, pBox *Matrix, const Weight *weight);
void maxPooling(const pBox *pbox, pBox *Matrix, int kernelSize, int stride);
void relu(struct pBox *pbox, mydataFmt *pbias);
void prelu(struct pBox *pbox, mydataFmt *pbias, mydataFmt *prelu_gmma);
void convolution(const Weight *weight, const pBox *pbox, pBox *outpBox, const struct pBox *matrix);
void fullconnect(const Weight *weight, const pBox *pbox, pBox *outpBox);
void readData(string filename, long dataNumber[], mydataFmt *pTeam[]);
long initConvAndFc(struct Weight *weight, int schannel, int lchannel, int kersize, int stride, int pad);
void initpRelu(struct pRelu *prelu, int width);
void softmax(const struct pBox *pbox);

void image2MatrixInit(Mat &image, struct pBox *pbox);
void featurePadInit(const pBox *pbox, pBox *outpBox, const int pad);
void maxPoolingInit(const pBox *pbox, pBox *Matrix, int kernelSize, int stride);
void feature2MatrixInit(const pBox *pbox, pBox *Matrix, const Weight *weight);
void convolutionInit(const Weight *weight, const pBox *pbox, pBox *outpBox, const struct pBox *matrix);
void fullconnectInit(const Weight *weight, pBox *outpBox);

bool cmpScore(struct orderScore lsh, struct orderScore rsh);
void nms(vector<struct Bbox> &boundingBox_, vector<struct orderScore> &bboxScore_, const float overlap_threshold, string modelname = "Union");
void refineAndSquareBbox(vector<struct Bbox> &vecBbox, const int &height, const int &width);

#endif


text:
#ifndef ADD_TO_CONFIG_INCLUDE
#define ADD_TO_CONFIG_INCLUDE

#include <optional>
#include <stdexcept>
#include <string>
#include <utility>
#include <vector>

#include <gsl/pointers>

#include "config/variablesMap.h"

#include "base-utils/nonEmptyString.h"

namespace execHelper::test {
inline void
addToConfig(const execHelper::config::SettingsKeys& key,
            const std::string& value,
            gsl::not_null<execHelper::config::VariablesMap*> config) {
    if(!config->add(key, value)) {
        throw std::runtime_error("Failed to add key " + key.back() +
                                 " with value '" + value + "' to config");
    }
}

inline void
addToConfig(const execHelper::config::SettingsKeys& key,
            const NonEmptyString& value,
            gsl::not_null<execHelper::config::VariablesMap*> config) {
    if(!config->add(key, *value)) {
        throw std::runtime_error("Failed to add key " + key.back() +
                                 " with value '" + *value + "' to config");
    }
}

template <typename T>
inline void
addToConfig(const execHelper::config::SettingsKeys& key,
            const std::vector<T>& value,
            gsl::not_null<execHelper::config::VariablesMap*> config) {
    if(!config->add(key, value)) {
        throw std::runtime_error("Failed to add key " + key.back() +
                                 " with first value '" + value.front() +
                                 "' to config");
    }
}

inline void
addToConfig(const execHelper::config::SettingsKeys& key, bool value,
            gsl::not_null<execHelper::config::VariablesMap*> config) {
    if(value) {
        if(!config->add(key, "yes")) {
            throw std::runtime_error("Failed to add key " + key.back() +
                                     " with value 'true' to config");
        }
    }
    if(!config->add(key, "no")) {
        throw std::runtime_error("Failed to add key " + key.back() +
                                 " with value 'false' to config");
    }
}

inline void
addToConfig(execHelper::config::SettingsKeys key,
            const std::pair<std::string, std::string>& value,
            gsl::not_null<execHelper::config::VariablesMap*> config) {
    key.push_back(value.first);
    if(!config->add(key, value.second)) {
        throw std::runtime_error("Failed to add key " + value.first +
                                 " with first value '" + value.second +
                                 "' to config");
    }
}

template <typename T>
inline void
addToConfig(const execHelper::config::SettingsKeys& key,
            const std::optional<T>& value,
            gsl::not_null<execHelper::config::VariablesMap*> config) {
    if(value) {
        addToConfig(key, *value, config);
    }
}

template <typename T>
inline void
addToConfig(const execHelper::config::SettingsKey& key, const T& value,
            gsl::not_null<execHelper::config::VariablesMap*> config) {
    addToConfig(execHelper::config::SettingsKeys({key}), value, config);
}
} // namespace execHelper::test

#endif /* ADD_TO_CONFIG_INCLUDE */

max_stars_repo_name: exec-helper/source

max_issues_repo_name: exec-helper/source

max_forks_repo_name: exec-helper/source


text:
/*
 *
 *
 * The functions in this file are directly taken from
 https://github.com/s5248/cppconvnet
 The license and the copying notice are replicated below.
----------
LICENSE
----------
Copyright (c) 2015, BiaoZhi Huang.
All rights reserved.

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

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

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

* Neither the name of cppconvnet 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.

---------------
COPYING
---------------

Copyright (c) 2014 The MatConvNet team.
Copyright (c) 2015 BiaoZhi Huang.
All rights reserved.

Redistribution and use in source and binary forms are permitted
provided that the above copyright notice and this paragraph are
duplicated in all such forms and that any documentation,
advertising materials, and other materials related to such
distribution and use acknowledge that the software was developed
by the <organization>. The name of the
<organization> may not be used to endorse or promote products derived
from this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 */

#include "common_types.h"
#include <stdbool.h>
#include <stdlib.h>
#include <stdio.h>
#include <cblas.h>
#include "utils.h"

static inline int floor_divide(int a, int b) {
  if (a >= 0) return a / b;
  else return (a - b + 1) / b;
}

static inline int ceil_divide(int a, int b) {
  if (a >= 0) return (a + b - 1) / b;
  else return a / b;
}

static inline int min (int a, int b) {
  return a < b ? a : b;
}

static void CppConvnetIm2Row(float* stacked, const float* data, int numPatchesX,
                      int numPatchesY, int numRows, const TensorDim input_dim,
                      const TensorDim filter_dim,
                      int stride, int pad)
{

  /*
  Fill a row of the stacked image at a time. Since patches are stored
  aint the columns, scanning a row menas visiting all patche once.
  Each row corresponds to a particular offset within each patch.

  In this manner, as we fill a row
  we tend to access spatially adiacent elements
  in the input image, particulary for small strides.
  */
  float *im2row = stacked;
  // numRows = KxKXC, numPatchesX = H, numPatchesY = W
  for (int row = 0; row < numRows; ++row) {
    /*
    Get the patch offset corresponding to this row of the stacked
    image.
    */
    int u = row;
    int v = u / filter_dim.w;
    int z = v / filter_dim.h; // filter channel no
    u %= filter_dim.w; // filter col no
    v %= filter_dim.h; // filter row no

    /*
    Filling this row amounts to visiting all the pixels in the input
    image that appear at a given offset in the outut patches. Accounting
    for the subsampling of the output patches and input padding,
    these pixels are given by

    x_data(x) = x * m_stride[1] + u - m_pad[2],  0 <= x < numPatchesX
    y_data(y) = y * m_stride[0] + v - m_pad[0],   0 <= y < numPatchesY
    z_data(z) = z.

    Here (x,y) are the spatial indexes of the output patches. Depending
    on the padding, some of these values will read pixels outised
    the input image, which should default to 0. In particular, x lands
    on a x_data(x) within the image if x0 <= x < x1 where:

    x_data(x) >= 0 <=> x >= (m_pad[2] - u) / stride
    <=> x >= ceil((m_pad[2] - u) / stride) = x0
    x_data(x) <= m_inputDims[1]-1 <=> x <= (m_inputDims[1]-1 + m_pad[2] - u) / stride
    <=> x <= floor((m_inputDims[1]-1 + m_pad[2] - u) / stride)
    <=> x <  floor((m_inputDims[1]-1 + m_pad[2] - u) / stride) + 1 = x1

    and the same for y. Note that, while usually x0 <= x1, there are
    special cases for which x1 < x0. This is accounted for in the loops
    below.
    */

    int x0 = min(numPatchesX, (int)ceil_divide(pad - u, stride));
    int y0 = min(numPatchesY, (int)ceil_divide(pad - v, stride));
    int x1 = min(numPatchesX, (int)floor_divide(input_dim.w - 1 + pad - u, stride) + 1);
    int y1 = min(numPatchesY, (int)floor_divide(input_dim.h - 1 + pad - v, stride) + 1);
    int x;
    int y;

    for (y = 0; y < y0; ++y) {
      for (x = 0; x < numPatchesX; ++x) {
        *stacked++ = 0;
      }
    }
    for (; y < y1; ++y) {
      for (x = 0; x < x0; ++x) {
        *stacked++ = 0;
      }
      int y_data = y * stride + v - pad;
      int x_data = x * stride + u - pad;
      float const * b = data + (z * input_dim.h + y_data) * input_dim.w + x_data;
      for (; x < x1; ++x) {
        *stacked++ = *b;
        b += stride;
      }
      for (; x < numPatchesX; ++x) {
        *stacked++ = 0;
      }
    }
    for (; y < numPatchesY; ++y) {
      for (x = 0; x < numPatchesX; ++x) {
        *stacked++ = 0;
      }
    }
  }

  PrintMat("im2row", im2row, 1, numRows*numPatchesX*numPatchesY, CblasRowMajor);

}

// Not working for group != 1
bool CppConvnetConvLayer(const float *in_data, const float *filters,
                         const float *bias, TensorDim in_dim,
                         TensorDim filt_dim, int stride, int pad, int group,
                         float *output)
{
  float dataMult=1.0;
  float outputMult=0.0;
  TensorDim out_dim;
  out_dim.w = (in_dim.w + (pad + pad) - filt_dim.w) / stride + 1;
  out_dim.h = (in_dim.h + (pad + pad) - filt_dim.h) / stride + 1;
  out_dim.c = filt_dim.n;
  out_dim.n = in_dim.n;
  int m_biasDims = out_dim.c;

  //int numGroups = in_dim.c / filt_dim.c;// assume = 1
  bool fullyConnectedMode = (out_dim.w == 1 &&
      out_dim.h == 1 &&
    stride == 1 &&
    pad == 0 &&
    group == 1);
  if (fullyConnectedMode){
    return false;
  }
  printf("---------Inputs------\n");
  PrintTensor(in_data, in_dim);
  printf("---------------------\n");
  printf("-------Filters-------\n");
  PrintTensor(filters, filt_dim);
  int numFiltersPerGroup = filt_dim.n / group; // M
  int numOutputPixels = out_dim.w * out_dim.h; // H*W
  int filtersVolume = filt_dim.w * filt_dim.h * filt_dim.c/group; // K*K*C
  int tempVolume = numOutputPixels * filtersVolume * group; // H*W*K*K*C --> looks to be size of im2row buffer
  float* tempMemory = malloc( sizeof(float) *tempVolume); // im2row buffer?
  float* tempOnes = malloc(sizeof(float)*numOutputPixels); // buffer of size H*W all set to 1
  for (int i = 0; i < numOutputPixels; i++) tempOnes[i] = 1.0;

  for (int image = 0; image < in_dim.n; ++image) {
    int dataOffset = (in_dim.w * in_dim.h * in_dim.c) * image;// pointer to batch of images
    int outputOffset = (out_dim.w * out_dim.h * out_dim.c) * image;// pointer to batch of output maps.
    //under sample
    CppConvnetIm2Row(tempMemory, in_data + dataOffset, out_dim.w, out_dim.h,
                     filtersVolume*group, in_dim, filt_dim, stride, pad);

    for (int g = 0; g < group; ++g) {
      int filterGrpOffset = filtersVolume * numFiltersPerGroup * g;//0
      int tempGrpOffset = numOutputPixels * filtersVolume * g;//0
      int outputGrpOffset = numOutputPixels * numFiltersPerGroup * g;//0
      float alpha = dataMult;
      float beta = outputMult;
      //convolution
      // A -> MxK B-> KxN
      PrintMat("im2row", tempMemory + tempGrpOffset, numOutputPixels, filtersVolume, CblasColMajor);
      PrintMat("filters", filters + filterGrpOffset, filtersVolume, numFiltersPerGroup, CblasColMajor);
      //PrintMat("im2row", tempMemory, 1, filtersVolume* numOutputPixels, CblasColMajor);
      cblas_sgemm(CblasColMajor, CblasNoTrans, CblasNoTrans,
                  numOutputPixels/*m = HxW*/,
                  numFiltersPerGroup/*n = M*/, filtersVolume/*k = KxKxC/group*/,
                  alpha, tempMemory + tempGrpOffset, /*A*/
                  numOutputPixels/*lda*/,
                  filters + filterGrpOffset, // B
                  filtersVolume, // ldb
                  beta,
                  output + outputOffset + outputGrpOffset, // C
                  numOutputPixels // ldc
                  );

    }
    if (bias!=NULL) {
      float alpha = 1;
      float beta = 1;
      int k = 1;
      //add bias
      cblas_sgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, numOutputPixels,
                  m_biasDims, k, alpha, tempOnes, numOutputPixels,
                  bias, k, beta, output + outputOffset, numOutputPixels);
    }
  }
  free(tempMemory);
  free(tempOnes);
  return true;
}

max_stars_repo_name: gplhegde/convolution-flavors

max_issues_repo_name: chayitw/convolution-flavors

max_forks_repo_name: chayitw/convolution-flavors


text:
/*
 # This file is part of the Astrometry.net suite.
 # Licensed under a 3-clause BSD style license - see LICENSE
 */

#include <stdio.h>
#include <assert.h>
#include <string.h>
#include <math.h>

#include <gsl/gsl_matrix_double.h>
#include <gsl/gsl_vector_double.h>

#include "sip-utils.h"
#include "tweak.h"
#include "healpix.h"
#include "dualtree_rangesearch.h"
#include "kdtree_fits_io.h"
#include "mathutil.h"
#include "log.h"
#include "permutedsort.h"
#include "gslutils.h"
#include "errors.h"
#include "fit-wcs.h"

// TODO:
//
//  1. Write the document which explains every step of tweak in detail with
//     comments relating exactly to the code.
//  2. Implement polynomial terms - use order parameter to zero out A matrix
//  3. Make it robust to outliers
//     - Make jitter evolve as fit evolves
//     - Sigma clipping
//  4. Need test image with non-trivial rotation to test CD transpose problem
//
//
//  - Put CD inverse into its own function
//
//  BUG? transpose of CD matrix is similar to CD matrix!
//  BUG? inverse when computing sx/sy (i.e. same transpose issue)
//  Ability to fit without re-doing correspondences
//  Split fit x/y (i.e. two fits one for x one for y)

#define KERNEL_SIZE 5
#define KERNEL_MARG ((KERNEL_SIZE-1)/2)

sip_t* tweak_just_do_it(const tan_t* wcs, const starxy_t* imagexy,
                        const double* starxyz,
                        const double* star_ra, const double* star_dec,
                        const double* star_radec,
                        int nstars, double jitter_arcsec,
                        int order, int inverse_order, int iterations,
                        anbool weighted, anbool skip_shift) {
    tweak_t* twee = NULL;
    sip_t* sip = NULL;

    twee = tweak_new();
    twee->jitter = jitter_arcsec;
    twee->sip->a_order  = twee->sip->b_order  = order;
    twee->sip->ap_order = twee->sip->bp_order = inverse_order;
    twee->weighted_fit = weighted;
    if (skip_shift)
        tweak_skip_shift(twee);

    tweak_push_image_xy(twee, imagexy);
    if (starxyz)
        tweak_push_ref_xyz(twee, starxyz, nstars);
    else if (star_ra && star_dec)
        tweak_push_ref_ad(twee, star_ra, star_dec, nstars);
    else if (star_radec)
        tweak_push_ref_ad_array(twee, star_radec, nstars);
    else {
        logerr("Need starxyz, (star_ra and star_dec), or star_radec");
        return NULL;
    }
    tweak_push_wcs_tan(twee, wcs);
    tweak_iterate_to_order(twee, order, iterations);

    // Steal the resulting SIP structure
    sip = twee->sip;
    twee->sip = NULL;

    tweak_free(twee);
    return sip;
}

static void get_dydx_range(double* ximg, double* yimg, int nimg,
                           double* xcat, double* ycat, int ncat,
                           double *mindx, double *mindy,
                           double *maxdx, double *maxdy) {
    int i, j;
    *maxdx = -1e100;
    *mindx = 1e100;
    *maxdy = -1e100;
    *mindy = 1e100;

    for (i = 0; i < nimg; i++) {
        for (j = 0; j < ncat; j++) {
            double dx = ximg[i] - xcat[j];
            double dy = yimg[i] - ycat[j];
            *maxdx = MAX(dx, *maxdx);
            *maxdy = MAX(dy, *maxdy);
            *mindx = MIN(dx, *mindx);
            *mindy = MIN(dy, *mindy);
        }
    }
}

static void get_shift(double* ximg, double* yimg, int nimg,
                      double* xcat, double* ycat, int ncat,
                      double mindx, double mindy, double maxdx, double maxdy,
                      double* xshift, double* yshift) {
    int i, j;
    int themax, themaxind, ys, xs;

    // hough transform
    int hsz = 1000; // hough histogram size (per side)
    int *hough = calloc(hsz * hsz, sizeof(int)); // allocate bins
    int kern[] = {0, 2, 3, 2, 0,  // approximate gaussian smoother
                  2, 7, 12, 7, 2,       // should be KERNEL_SIZE x KERNEL_SIZE
                  3, 12, 20, 12, 3,
                  2, 7, 12, 7, 2,
                  0, 2, 3, 2, 0};

    assert(sizeof(kern) == KERNEL_SIZE * KERNEL_SIZE * sizeof(int));

    for (i = 0; i < nimg; i++) {    // loop over all pairs of source-catalog objs
        for (j = 0; j < ncat; j++) {
            double dx = ximg[i] - xcat[j];
            double dy = yimg[i] - ycat[j];
            int hszi = hsz - 1;
            int iy = hszi * ( (dy - mindy) / (maxdy - mindy) ); // compute deltay using implicit floor
            int ix = hszi * ( (dx - mindx) / (maxdx - mindx) ); // compute deltax using implicit floor

            // check to make sure the point is in the box
            if (KERNEL_MARG <= iy && iy < hsz - KERNEL_MARG &&
                KERNEL_MARG <= ix && ix < hsz - KERNEL_MARG) {
                int kx, ky;
                for (ky = -2; ky <= 2; ky++)
                    for (kx = -2; kx <= 2; kx++)
                        hough[(iy - ky)*hsz + (ix - kx)] += kern[(ky + 2) * 5 + (kx + 2)];
            }
        }
    }

    // find argmax in hough
    themax = 0;
    themaxind = -1;
    for (i = 0; i < hsz*hsz; i++) {
        if (themax < hough[i]) {
            themaxind = i;
            themax = hough[i];
        }
    }
    // which hough bin is the max?
    ys = themaxind / hsz;
    xs = themaxind % hsz;

    *yshift = (ys / (double)hsz) * (maxdy - mindy) + mindy;
    *xshift = (xs / (double)hsz) * (maxdx - mindx) + mindx;
    debug("xs = %d, ys = %d\n", xs, ys);
    debug("get_shift: mindx=%g, maxdx=%g, mindy=%g, maxdy=%g\n", mindx, maxdx, mindy, maxdy);
    debug("get_shift: xs=%g, ys=%g\n", *xshift, *yshift);

    free(hough);
}

static sip_t* do_entire_shift_operation(tweak_t* t, double rho) {
    get_shift(t->x, t->y, t->n,
              t->x_ref, t->y_ref, t->n_ref,
              rho*t->mindx, rho*t->mindy, rho*t->maxdx, rho*t->maxdy,
              &t->xs, &t->ys);
    wcs_shift(&(t->sip->wcstan), t->xs, t->ys);
    return NULL;
}

/* This function is intended only for initializing newly allocated tweak
 * structures, NOT for operating on existing ones.*/
void tweak_init(tweak_t* t) {
    memset(t, 0, sizeof(tweak_t));
    t->sip = sip_create();
}

tweak_t* tweak_new() {
    tweak_t* t = malloc(sizeof(tweak_t));
    tweak_init(t);
    return t;
}

void tweak_iterate_to_order(tweak_t* t, int maxorder, int iterations) {
    int order;
    int k;

    for (order=1; order<=maxorder; order++) {
        logverb("\n");
        logverb("--------------------------------\n");
        logverb("Order %i\n", order);
        logverb("--------------------------------\n");

        t->sip->a_order  = t->sip->b_order  = order;
        //t->sip->ap_order = t->sip->bp_order = order;
        tweak_go_to(t, TWEAK_HAS_CORRESPONDENCES);

        for (k=0; k<iterations; k++) {
            logverb("\n");
            logverb("--------------------------------\n");
            logverb("Iterating tweak: order %i, step %i\n", order, k);
            t->state &= ~TWEAK_HAS_LINEAR_CD;
            tweak_go_to(t, TWEAK_HAS_LINEAR_CD);
            tweak_clear_correspondences(t);
        }
    }
}

#define CHECK_STATE(x) if (state & x) sl_append(s, #x)

static char* state_string(unsigned int state) {
    sl* s = sl_new(4);
    char* str;
    CHECK_STATE(TWEAK_HAS_SIP);
    CHECK_STATE(TWEAK_HAS_IMAGE_XY);
    CHECK_STATE(TWEAK_HAS_IMAGE_XYZ);
    CHECK_STATE(TWEAK_HAS_IMAGE_AD);
    CHECK_STATE(TWEAK_HAS_REF_XY);
    CHECK_STATE(TWEAK_HAS_REF_XYZ);
    CHECK_STATE(TWEAK_HAS_REF_AD);
    CHECK_STATE(TWEAK_HAS_CORRESPONDENCES);
    CHECK_STATE(TWEAK_HAS_COARSLY_SHIFTED);
    CHECK_STATE(TWEAK_HAS_FINELY_SHIFTED);
    CHECK_STATE(TWEAK_HAS_REALLY_FINELY_SHIFTED);
    CHECK_STATE(TWEAK_HAS_LINEAR_CD);
    str = sl_join(s, " ");
    sl_free2(s);
    return str;
}

char* tweak_get_state_string(const tweak_t* t) {
    return state_string(t->state);
}

#undef CHECK_STATE

void tweak_clear_correspondences(tweak_t* t) {
    if (t->state & TWEAK_HAS_CORRESPONDENCES) {
        // our correspondences are also now toast
        assert(t->image);
        assert(t->ref);
        assert(t->dist2);
        il_free(t->image);
        il_free(t->ref);
        dl_free(t->dist2);
        if (t->weight)
            dl_free(t->weight);
        t->image    = NULL;
        t->ref      = NULL;
        t->dist2    = NULL;
        t->weight   = NULL;
        t->state &= ~TWEAK_HAS_CORRESPONDENCES;
    }
    assert(!t->image);
    assert(!t->ref);
    assert(!t->dist2);
    assert(!t->weight);
}

void tweak_clear_on_sip_change(tweak_t* t) {
    // tweak_clear_correspondences(t);
    tweak_clear_image_ad(t);
    tweak_clear_ref_xy(t);
    tweak_clear_image_xyz(t);
}

// ref_xy are the catalog star positions in image coordinates
void tweak_clear_ref_xy(tweak_t* t) {
    if (t->state & TWEAK_HAS_REF_XY) {
        //assert(t->x_ref);
        free(t->x_ref);
        //assert(t->y_ref);
        t->x_ref = NULL;
        free(t->y_ref);
        t->y_ref = NULL;
        t->state &= ~TWEAK_HAS_REF_XY;
    }
    assert(!t->x_ref);
    assert(!t->y_ref);
}

// radec of catalog stars
void tweak_clear_ref_ad(tweak_t* t) {
    if (t->state & TWEAK_HAS_REF_AD) {
        assert(t->a_ref);
        free(t->a_ref);
        t->a_ref = NULL;
        assert(t->d_ref);
        free(t->d_ref);
        t->d_ref = NULL;
        t->n_ref = 0;
        tweak_clear_correspondences(t);
        tweak_clear_ref_xy(t);
        t->state &= ~TWEAK_HAS_REF_AD;
    }
    assert(!t->a_ref);
    assert(!t->d_ref);
}

// image objs in ra,dec according to current tweak
void tweak_clear_image_ad(tweak_t* t) {
    if (t->state & TWEAK_HAS_IMAGE_AD) {
        assert(t->a);
        free(t->a);
        t->a = NULL;
        assert(t->d);
        free(t->d);
        t->d = NULL;
        t->state &= ~TWEAK_HAS_IMAGE_AD;
    }
    assert(!t->a);
    assert(!t->d);
}

void tweak_clear_image_xyz(tweak_t* t) {
    if (t->state & TWEAK_HAS_IMAGE_XYZ) {
        assert(t->xyz);
        free(t->xyz);
        t->xyz = NULL;
        t->state &= ~TWEAK_HAS_IMAGE_XYZ;
    }
    assert(!t->xyz);
}

void tweak_clear_image_xy(tweak_t* t) {
    if (t->state & TWEAK_HAS_IMAGE_XY) {
        assert(t->x);
        free(t->x);
        t->x = NULL;
        assert(t->y);
        free(t->y);
        t->y = NULL;
        t->state &= ~TWEAK_HAS_IMAGE_XY;
    }
    assert(!t->x);
    assert(!t->y);
}

// tell us (from outside tweak) where the catalog stars are
void tweak_push_ref_ad(tweak_t* t, const double* a, const double *d, int n) {
    assert(a);
    assert(d);
    assert(n);
    tweak_clear_ref_ad(t);
    assert(!t->a_ref);
    assert(!t->d_ref);
    t->a_ref = malloc(sizeof(double) * n);
    t->d_ref = malloc(sizeof(double) * n);
    memcpy(t->a_ref, a, n*sizeof(double));
    memcpy(t->d_ref, d, n*sizeof(double));
    t->n_ref = n;
    t->state |= TWEAK_HAS_REF_AD;
}

void tweak_push_ref_ad_array(tweak_t* t, const double* ad, int n) {
    int i;
    assert(ad);
    assert(n);
    tweak_clear_ref_ad(t);
    assert(!t->a_ref);
    assert(!t->d_ref);
    t->a_ref = malloc(sizeof(double) * n);
    t->d_ref = malloc(sizeof(double) * n);
    for (i=0; i<n; i++) {
        t->a_ref[i] = ad[2*i + 0];
        t->d_ref[i] = ad[2*i + 1];
    }
    t->n_ref = n;
    t->state |= TWEAK_HAS_REF_AD;
}

static void ref_xyz_from_ad(tweak_t* t) {
    int i;
    assert(t->state & TWEAK_HAS_REF_AD);
    assert(!t->xyz_ref);
    t->xyz_ref = malloc(sizeof(double) * 3 * t->n_ref);
    assert(t->xyz_ref);
    for (i = 0; i < t->n_ref; i++)
        radecdeg2xyzarr(t->a_ref[i], t->d_ref[i], t->xyz_ref + 3 * i);
    t->state |= TWEAK_HAS_REF_XYZ;
}

static void ref_ad_from_xyz(tweak_t* t) {
    int i, n;
    assert(t->state & TWEAK_HAS_REF_XYZ);
    assert(!t->a_ref);
    assert(!t->d_ref);
    n = t->n_ref;
    t->a_ref = malloc(sizeof(double) * n);
    t->d_ref = malloc(sizeof(double) * n);
    assert(t->a_ref);
    assert(t->d_ref);
    for (i=0; i<n; i++)
        xyzarr2radecdeg(t->xyz_ref + 3*i, t->a_ref + i, t->d_ref + i);
    t->state |= TWEAK_HAS_REF_XYZ;
}

// tell us (from outside tweak) where the catalog stars are
void tweak_push_ref_xyz(tweak_t* t, const double* xyz, int n) {
    assert(xyz);
    assert(n);
    tweak_clear_ref_ad(t);
    assert(!t->xyz_ref);
    t->xyz_ref = malloc(sizeof(double) * 3 * n);
    assert(t->xyz_ref);
    memcpy(t->xyz_ref, xyz, 3*n*sizeof(double));
    t->n_ref = n;
    t->state |= TWEAK_HAS_REF_XYZ;
}

void tweak_push_image_xy(tweak_t* t, const starxy_t* xy) {
    tweak_clear_image_xy(t);
    t->x = starxy_copy_x(xy);
    t->y = starxy_copy_y(xy);
    t->n = starxy_n(xy);
    t->state |= TWEAK_HAS_IMAGE_XY;
}

void tweak_skip_shift(tweak_t* t) {
    t->state |= (TWEAK_HAS_COARSLY_SHIFTED | TWEAK_HAS_FINELY_SHIFTED |
                 TWEAK_HAS_REALLY_FINELY_SHIFTED);
}

// DualTree RangeSearch callback. We want to keep track of correspondences.
// Potentially the matching could be many-to-many; we allow this and hope the
// optimizer can take care of it.
static void dtrs_match_callback(void* extra, int image_ind, int ref_ind, double dist2) {
    tweak_t* t = extra;
    image_ind = kdtree_permute(t->kd_image, image_ind);
    ref_ind   = kdtree_permute(t->kd_ref,   ref_ind);
    il_append(t->image, image_ind);
    il_append(t->ref, ref_ind);
    dl_append(t->dist2, dist2);
    if (t->weight)
        dl_append(t->weight, exp(-dist2 / (2.0 * t->jitterd2)));
}

void tweak_push_correspondence_indices(tweak_t* t, il* image, il* ref, dl* distsq, dl* weight) {
    t->image = image;
    t->ref = ref;
    t->dist2 = distsq;
    t->weight = weight;
    t->state |= TWEAK_HAS_CORRESPONDENCES;
}

// The jitter is in radians
static void find_correspondences(tweak_t* t, double jitter) {
    double dist;
    double* data_image = malloc(sizeof(double) * t->n * 3);
    double* data_ref = malloc(sizeof(double) * t->n_ref * 3);

    assert(t->state & TWEAK_HAS_IMAGE_XYZ);
    assert(t->state & TWEAK_HAS_REF_XYZ);
    tweak_clear_correspondences(t);

    memcpy(data_image, t->xyz, 3*t->n*sizeof(double));
    memcpy(data_ref, t->xyz_ref, 3*t->n_ref*sizeof(double));

    t->kd_image = kdtree_build(NULL, data_image, t->n, 3, 4, KDTT_DOUBLE,
                               KD_BUILD_BBOX);

    t->kd_ref = kdtree_build(NULL, data_ref, t->n_ref, 3, 4, KDTT_DOUBLE,
                             KD_BUILD_BBOX);

    // Storage for correspondences
    t->image = il_new(600);
    t->ref = il_new(600);
    t->dist2 = dl_new(600);
    if (t->weighted_fit)
        t->weight = dl_new(600);

    dist = rad2dist(jitter);

    logverb("search radius = %g arcsec\n", rad2arcsec(jitter));

    // Find closest neighbours
    dualtree_rangesearch(t->kd_image, t->kd_ref,
                         RANGESEARCH_NO_LIMIT, dist, FALSE, NULL,
                         dtrs_match_callback, t,
                         NULL, NULL);

    kdtree_free(t->kd_image);
    kdtree_free(t->kd_ref);
    t->kd_image = NULL;
    t->kd_ref = NULL;
    free(data_image);
    free(data_ref);

    logverb("Number of correspondences: %zu\n", il_size(t->image));
}

static double correspondences_rms_arcsec(tweak_t* t, int weighted) {
    double err2 = 0.0;
    int i;
    double totalweight = 0.0;
    for (i=0; i<il_size(t->image); i++) {
        double imgxyz[3];
        double refxyz[3];
        double weight;
        int refi, imgi;
        if (weighted && t->weight)
            weight = dl_get(t->weight, i);
        else
            weight = 1.0;
        totalweight += weight;

        imgi = il_get(t->image, i);
        sip_pixelxy2xyzarr(t->sip, t->x[imgi], t->y[imgi], imgxyz);

        refi = il_get(t->ref, i);
        radecdeg2xyzarr(t->a_ref[refi], t->d_ref[refi], refxyz);

        err2 += weight * distsq(imgxyz, refxyz, 3);
    }
    return distsq2arcsec( err2 / totalweight );
}

#if 0
// in arcseconds^2 on the sky (chi-sq)
static double figure_of_merit(tweak_t* t, double *rmsX, double *rmsY) {
    double sqerr = 0.0;
    int i;
    for (i = 0; i < il_size(t->image); i++) {
        double a, d;
        double xyzpt[3];
        double xyzpt_ref[3];
        sip_pixelxy2radec(t->sip, t->x[il_get(t->image, i)],
                          t->y[il_get(t->image, i)], &a, &d);

        // xref and yref should be intermediate WC's not image x and y!
        radecdeg2xyzarr(a, d, xyzpt);
        radecdeg2xyzarr(t->a_ref[il_get(t->ref, i)],
                        t->d_ref[il_get(t->ref, i)], xyzpt_ref);
        sqerr += distsq(xyzpt, xyzpt_ref, 3);
    }
    return rad2arcsec(1)*rad2arcsec(1)*sqerr;
}

static double figure_of_merit2(tweak_t* t) {
    // find error in pixels^2
    double sqerr = 0.0;
    int i;
    for (i = 0; i < il_size(t->image); i++) {
        double x, y, dx, dy;
        Unused anbool ok;
        ok = sip_radec2pixelxy(t->sip, t->a_ref[il_get(t->ref, i)], t->d_ref[il_get(t->ref, i)], &x, &y);
        assert(ok);
        dx = t->x[il_get(t->image, i)] - x;
        dy = t->y[il_get(t->image, i)] - y;
        sqerr += dx * dx + dy * dy;
    }
    // convert to arcsec^2
    return sqerr * square(sip_pixel_scale(t->sip));
}
#endif

// FIXME: adapt this function to take as input the correspondences to use VVVVV
//    wic is World Intermediate Coordinates, either along ra or dec
//       i.e. canonical image coordinates
//    wic_corr is the list of corresponding indexes for wip
//    pix_corr is the list of corresponding indexes for pixels
//    pix is raw image pixels (either u or v)
//    siporder is the sip order up to MAXORDER (defined in sip.h)
//    the correspondences are passed so that we can stick RANSAC around the whole
//    thing for better estimation.

// Run a polynomial tweak
static void do_sip_tweak(tweak_t* t) {
    sip_t sipout;
    size_t i, M;

    // a_order and b_order should be the same!
    assert(t->sip->a_order == t->sip->b_order);

    debug("do_sip_tweak starting.\n");
    logverb("RMS error of correspondences: %g arcsec\n",
            correspondences_rms_arcsec(t, 0));
    if (t->weighted_fit)
        logverb("Weighted RMS error of correspondences: %g arcsec\n",
                correspondences_rms_arcsec(t, 1));

    M = il_size(t->image);
    double* starxyz = malloc(M * 3 * sizeof(double));
    double* fieldxy = malloc(M * 2 * sizeof(double));
    double* weights = NULL;
    if (t->weighted_fit)
        weights = malloc(M * sizeof(double));

    int result;
    for (i=0; i<M; i++) {
        int refi;
        int imi;
        imi = il_get(t->image, i);
        fieldxy[2*i + 0] = t->x[imi];
        fieldxy[2*i + 1] = t->y[imi];
        refi = il_get(t->ref, i);
        radecdeg2xyzarr(t->a_ref[refi], t->d_ref[refi], starxyz + i*3);
        if (t->weighted_fit)
            weights[i] = dl_get(t->weight, i);
    }

    int doshift = 1;
    result = fit_sip_wcs(starxyz, fieldxy, weights, M,
                         &(t->sip->wcstan), t->sip->a_order, t->sip->ap_order, doshift,
                         &sipout);
    free(starxyz);
    free(fieldxy);
    free(weights);
    if (result) {
        ERROR("fit_sip_wcs failed\n");
        return;
    }
    memcpy(t->sip, &sipout, sizeof(sip_t));

    // recalc using new SIP
    tweak_clear_on_sip_change(t);
    tweak_go_to(t, TWEAK_HAS_IMAGE_AD);
    tweak_go_to(t, TWEAK_HAS_REF_XY);

    logverb("RMS error of correspondences: %g arcsec\n",
            correspondences_rms_arcsec(t, 0));
    if (t->weighted_fit)
        logverb("Weighted RMS error of correspondences: %g arcsec\n",
                correspondences_rms_arcsec(t, 1));
}

// Really what we want is some sort of fancy dependency system... DTDS!
// Duct-tape dependencey system (DTDS)
#define done(x) t->state |= x; return x;
#define want(x)                                 \
    if (flag == x && t->state & x)              \
        return x;                               \
    else if (flag == x)
#define ensure(x)                               \
    if (!(t->state & x)) {                      \
        return tweak_advance_to(t, x);          \
    }

unsigned int tweak_advance_to(tweak_t* t, unsigned int flag) {
    want(TWEAK_HAS_IMAGE_AD) {
        int jj;
        ensure(TWEAK_HAS_SIP);
        ensure(TWEAK_HAS_IMAGE_XY);
        debug("Satisfying TWEAK_HAS_IMAGE_AD\n");
        // Convert to ra dec
        assert(!t->a);
        assert(!t->d);
        t->a = malloc(sizeof(double) * t->n);
        t->d = malloc(sizeof(double) * t->n);
        for (jj = 0; jj < t->n; jj++)
            sip_pixelxy2radec(t->sip, t->x[jj], t->y[jj], t->a + jj, t->d + jj);
        done(TWEAK_HAS_IMAGE_AD);
    }

    want(TWEAK_HAS_REF_AD) {
        if (!(t->a_ref && t->d_ref)) {
            ensure(TWEAK_HAS_REF_XYZ);
            debug("Satisfying TWEAK_HAS_REF_AD\n");
            ref_ad_from_xyz(t);
        }
        assert(t->a_ref && t->d_ref);
        done(TWEAK_HAS_REF_AD);
    }

    want(TWEAK_HAS_REF_XYZ) {
        if (!t->xyz_ref) {
            ensure(TWEAK_HAS_REF_AD);
            debug("Satisfying TWEAK_HAS_REF_XYZ\n");
            ref_xyz_from_ad(t);
        }
        assert(t->xyz_ref);
        done(TWEAK_HAS_REF_XYZ);
    }

    want(TWEAK_HAS_REF_XY) {
        int jj;
        ensure(TWEAK_HAS_REF_AD);
        debug("Satisfying TWEAK_HAS_REF_XY\n");
        assert(t->state & TWEAK_HAS_REF_AD);
        assert(t->n_ref);
        assert(!t->x_ref);
        assert(!t->y_ref);
        t->x_ref = malloc(sizeof(double) * t->n_ref);
        t->y_ref = malloc(sizeof(double) * t->n_ref);
        for (jj = 0; jj < t->n_ref; jj++) {
            Unused anbool ok;
            ok = sip_radec2pixelxy(t->sip, t->a_ref[jj], t->d_ref[jj],
                                   t->x_ref + jj, t->y_ref + jj);
            assert(ok);
        }
        done(TWEAK_HAS_REF_XY);
    }

    want(TWEAK_HAS_IMAGE_XYZ) {
        int i;
        ensure(TWEAK_HAS_IMAGE_AD);
        debug("Satisfying TWEAK_HAS_IMAGE_XYZ\n");
        assert(!t->xyz);
        t->xyz = malloc(3 * t->n * sizeof(double));
        for (i = 0; i < t->n; i++)
            radecdeg2xyzarr(t->a[i], t->d[i], t->xyz + 3*i);
        done(TWEAK_HAS_IMAGE_XYZ);
    }

    want(TWEAK_HAS_COARSLY_SHIFTED) {
        ensure(TWEAK_HAS_REF_XY);
        ensure(TWEAK_HAS_IMAGE_XY);
        debug("Satisfying TWEAK_HAS_COARSLY_SHIFTED\n");
        get_dydx_range(t->x, t->y, t->n, t->x_ref, t->y_ref, t->n_ref,
                       &t->mindx, &t->mindy, &t->maxdx, &t->maxdy);
        do_entire_shift_operation(t, 1.0);
        tweak_clear_image_ad(t);
        tweak_clear_ref_xy(t);
        done(TWEAK_HAS_COARSLY_SHIFTED);
    }

    want(TWEAK_HAS_FINELY_SHIFTED) {
        ensure(TWEAK_HAS_REF_XY);
        ensure(TWEAK_HAS_IMAGE_XY);
        ensure(TWEAK_HAS_COARSLY_SHIFTED);
        debug("Satisfying TWEAK_HAS_FINELY_SHIFTED\n");
        // Shrink size of hough box
        do_entire_shift_operation(t, 0.3);
        tweak_clear_image_ad(t);
        tweak_clear_ref_xy(t);
        done(TWEAK_HAS_FINELY_SHIFTED);
    }

    want(TWEAK_HAS_REALLY_FINELY_SHIFTED) {
        ensure(TWEAK_HAS_REF_XY);
        ensure(TWEAK_HAS_IMAGE_XY);
        ensure(TWEAK_HAS_FINELY_SHIFTED);
        debug("Satisfying TWEAK_HAS_REALLY_FINELY_SHIFTED\n");
        // Shrink size of hough box
        do_entire_shift_operation(t, 0.03);
        tweak_clear_image_ad(t);
        tweak_clear_ref_xy(t);
        done(TWEAK_HAS_REALLY_FINELY_SHIFTED);
    }

    want(TWEAK_HAS_CORRESPONDENCES) {
        ensure(TWEAK_HAS_REF_XYZ);
        ensure(TWEAK_HAS_IMAGE_XYZ);
        debug("Satisfying TWEAK_HAS_CORRESPONDENCES\n");
        t->jitterd2 = arcsec2distsq(t->jitter);
        find_correspondences(t, 6.0 * arcsec2rad(t->jitter));
        done(TWEAK_HAS_CORRESPONDENCES);
    }

    want(TWEAK_HAS_LINEAR_CD) {
        ensure(TWEAK_HAS_SIP);
        ensure(TWEAK_HAS_REALLY_FINELY_SHIFTED);
        ensure(TWEAK_HAS_REF_XY);
        ensure(TWEAK_HAS_REF_AD);
        ensure(TWEAK_HAS_IMAGE_XY);
        ensure(TWEAK_HAS_CORRESPONDENCES);
        debug("Satisfying TWEAK_HAS_LINEAR_CD\n");
        do_sip_tweak(t);
        tweak_clear_on_sip_change(t);
        done(TWEAK_HAS_LINEAR_CD);
    }

    // small memleak -- but it's a major bug if this happens, so suck it up.
    logerr("die for dependence: %s\n", state_string(flag));
    assert(0);
    return -1;
}

void tweak_push_wcs_tan(tweak_t* t, const tan_t* wcs) {
    memcpy(&(t->sip->wcstan), wcs, sizeof(tan_t));
    t->state |= TWEAK_HAS_SIP;
}

void tweak_go_to(tweak_t* t, unsigned int dest_state) {
    while (!(t->state & dest_state))
        tweak_advance_to(t, dest_state);
}

#define SAFE_FREE(xx) {free((xx)); xx = NULL;}
void tweak_clear(tweak_t* t) {
    if (!t)
        return ;
    SAFE_FREE(t->a);
    SAFE_FREE(t->d);
    SAFE_FREE(t->x);
    SAFE_FREE(t->y);
    SAFE_FREE(t->xyz);
    SAFE_FREE(t->a_ref);
    SAFE_FREE(t->d_ref);
    SAFE_FREE(t->x_ref);
    SAFE_FREE(t->y_ref);
    SAFE_FREE(t->xyz_ref);
    if (t->sip) {
        sip_free(t->sip);
        t->sip = NULL;
    }
    il_free(t->image);
    il_free(t->ref);
    dl_free(t->dist2);
    if (t->weight)
        dl_free(t->weight);
    t->image = NULL;
    t->ref = NULL;
    t->dist2 = NULL;
    t->weight = NULL;
    kdtree_free(t->kd_image);
    kdtree_free(t->kd_ref);
}

void tweak_free(tweak_t* t) {
    tweak_clear(t);
    free(t);
}

max_stars_repo_licenses: Net-SNMP Xnet

max_issues_repo_licenses: Net-SNMP Xnet

max_forks_repo_licenses: Net-SNMP Xnet


text:
/* Copyright (c) 2020 Felix Kutzner (github.com/fkutzner)

 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.

 Except as contained in this notice, the name(s) of the above copyright holders
 shall not be used in advertising or otherwise to promote the sale, use or
 other dealings in this Software without prior written authorization.

*/

#pragma once

#include <libincmonk/verifier/Clause.h>

#include <cassert>
#include <cstdint>
#include <functional>
#include <gsl/span>
#include <iterator>
#include <limits>
#include <optional>
#include <ostream>

//
// Implementation of inline function declared in Clause.h
//

namespace incmonk::verifier {

constexpr Var::Var(uint32_t id) : m_rawValue(id) {}

constexpr Var::Var() : m_rawValue{0} {}

constexpr auto Var::getRawValue() const -> uint32_t
{
  return m_rawValue;
}

constexpr auto Var::operator==(Var rhs) const -> bool
{
  return m_rawValue == rhs.m_rawValue;
}

constexpr auto Var::operator!=(Var rhs) const -> bool
{
  return m_rawValue != rhs.m_rawValue;
}

constexpr auto Var::operator<(Var rhs) const -> bool
{
  return m_rawValue < rhs.m_rawValue;
}

constexpr auto Var::operator<=(Var rhs) const -> bool
{
  return m_rawValue <= rhs.m_rawValue;
}

constexpr auto Var::operator>(Var rhs) const -> bool
{
  return m_rawValue > rhs.m_rawValue;
}

constexpr auto Var::operator>=(Var rhs) const -> bool
{
  return m_rawValue >= rhs.m_rawValue;
}

inline auto operator""_Var(unsigned long long cnfValue) -> Var
{
  return Var{static_cast<uint32_t>(cnfValue)};
}

constexpr auto Key<Var>::get(Var const& item) -> std::size_t
{
  return item.getRawValue();
}

constexpr Lit::Lit(Var v, bool positive) : m_rawValue{(v.getRawValue() << 1) + (positive ? 1 : 0)}
{
}

constexpr Lit::Lit() : m_rawValue{0} {}

constexpr auto Lit::getRawValue() const -> uint32_t
{
  return m_rawValue;
}

constexpr auto Lit::getVar() const -> Var
{
  return Var{m_rawValue >> 1};
}

constexpr auto Lit::isPositive() const -> bool
{
  return (m_rawValue & 1) == 1;
}

constexpr auto Lit::operator-() const -> Lit
{
  return Lit{getVar(), !isPositive()};
}

constexpr auto Lit::operator==(Lit rhs) const -> bool
{
  return m_rawValue == rhs.m_rawValue;
}

constexpr auto Lit::operator!=(Lit rhs) const -> bool
{
  return m_rawValue != rhs.m_rawValue;
}

constexpr auto Lit::operator<(Lit rhs) const -> bool
{
  return m_rawValue < rhs.m_rawValue;
}

constexpr auto Lit::operator<=(Lit rhs) const -> bool
{
  return m_rawValue <= rhs.m_rawValue;
}

constexpr auto Lit::operator>(Lit rhs) const -> bool
{
  return m_rawValue > rhs.m_rawValue;
}

constexpr auto Lit::operator>=(Lit rhs) const -> bool
{
  return m_rawValue >= rhs.m_rawValue;
}

inline auto operator""_Lit(unsigned long long cnfValue) -> Lit
{
  Var var{static_cast<uint32_t>(cnfValue)};
  return Lit{var, cnfValue > 0};
}

constexpr auto Key<Lit>::get(Lit const& item) -> std::size_t
{
  return item.getRawValue();
}

inline auto Clause::operator[](size_type idx) noexcept -> Lit&
{
  return *((&m_firstLit) + idx);
}

inline auto Clause::operator[](size_type idx) const noexcept -> Lit const&
{
  return *((&m_firstLit) + idx);
}

inline auto Clause::getLiterals() noexcept -> gsl::span<Lit>
{
  return gsl::span<Lit>{&m_firstLit, (&m_firstLit) + m_size};
}

inline auto Clause::getLiterals() const noexcept -> gsl::span<Lit const>
{
  return gsl::span<Lit const>{&m_firstLit, (&m_firstLit) + m_size};
}

inline auto Clause::size() const noexcept -> size_type
{
  return m_size;
}

inline auto Clause::empty() const noexcept -> bool
{
  return m_size == 0;
}

inline void Clause::setState(ClauseVerificationState state) noexcept
{
  uint8_t rawState = static_cast<uint8_t>(state);
  m_flags = (m_flags & 0xFFFFFFFC) | rawState;
}

inline auto Clause::getState() const noexcept -> ClauseVerificationState
{
  return static_cast<ClauseVerificationState>(m_flags & 3);
}

inline auto Clause::getAddIdx() const noexcept -> ProofSequenceIdx
{
  return m_pointOfAdd;
}

inline Clause::Clause(size_type size,
                      ClauseVerificationState initialState,
                      ProofSequenceIdx addIdx) noexcept
  : m_size{size}, m_pointOfAdd{addIdx}, m_firstLit{Var{0}, false}
{
  setState(initialState);
}

inline auto ClauseCollection::Ref::operator==(Ref rhs) const noexcept -> bool
{
  return m_offset == rhs.m_offset;
}

inline auto ClauseCollection::Ref::operator!=(Ref rhs) const noexcept -> bool
{
  return !(*this == rhs);
}

inline ClauseCollection::RefIterator::RefIterator(char const* m_allocatorMemory,
                                                  std::size_t highWaterMark) noexcept
  : m_clausePtr{m_allocatorMemory}, m_distanceToEnd{highWaterMark}, m_currentRef{}
{
  if (m_distanceToEnd == 0) {
    m_clausePtr = nullptr;
  }
  // Otherwise, the iterator is valid and currentRef is 0, referring to the first clause
}

inline ClauseCollection::RefIterator::RefIterator() noexcept
  : m_clausePtr{nullptr}, m_distanceToEnd{0}, m_currentRef{}
{
}

inline auto ClauseCollection::RefIterator::operator*() const noexcept -> Ref const&
{
  assert(m_clausePtr != nullptr && "Dereferencing a non-dereferencaable RefIterator");
  return m_currentRef;
}

inline auto ClauseCollection::RefIterator::operator->() const noexcept -> Ref const*
{
  assert(m_clausePtr != nullptr && "Dereferencing a non-dereferencaable RefIterator");
  return &m_currentRef;
}

inline auto ClauseCollection::RefIterator::operator++(int) noexcept -> RefIterator
{
  RefIterator old = *this;
  ++(*this);
  return old;
}

inline auto ClauseCollection::RefIterator::operator==(RefIterator const& rhs) const noexcept -> bool
{
  return (this == &rhs) || (m_clausePtr == nullptr && rhs.m_clausePtr == nullptr) ||
         (m_clausePtr == rhs.m_clausePtr && m_distanceToEnd == rhs.m_distanceToEnd &&
          m_currentRef == rhs.m_currentRef);
}

inline auto ClauseCollection::RefIterator::operator!=(RefIterator const& rhs) const noexcept -> bool
{
  return !(*this == rhs);
}

}

namespace std {
template <>
struct hash<incmonk::verifier::Lit> {
  auto operator()(incmonk::verifier::Lit lit) const noexcept -> std::size_t
  {
    return std::hash<uint32_t>{}(lit.getRawValue());
  }
};

template <>
struct hash<incmonk::verifier::Var> {
  auto operator()(incmonk::verifier::Var var) const noexcept -> std::size_t
  {
    return std::hash<uint32_t>{}(var.getRawValue());
  }
};

template <>
struct hash<incmonk::verifier::ClauseCollection::Ref> {
  auto operator()(incmonk::verifier::ClauseCollection::Ref ref) const noexcept -> std::size_t
  {
    return std::hash<std::size_t>{}(ref.m_offset);
  }
};
}

max_stars_repo_name: fkutzner/IncrementalMonkey

max_issues_repo_name: fkutzner/IncrementalMonkey

max_forks_repo_name: fkutzner/IncrementalMonkey


text:
/*---
  Flow*: A Verification Tool for Cyber-Physical Systems.
  Authors: Xin Chen, Sriram Sankaranarayanan, and Erika Abraham.
  Email: Xin Chen <chenxin415@gmail.com> if you have questions or comments.

  The code is released as is under the GNU General Public License (GPL).
---*/

#ifndef INCLUDE_H_
#define INCLUDE_H_

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sys/stat.h>
#include <sys/types.h>

#include <cmath>
#include <mpfr.h>
#include <vector>
#include <string>
#include <list>
#include <iostream>
#include <cassert>
#include <map>
#include <time.h>
#include <algorithm>
#include <gsl/gsl_poly.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <glpk.h>

const int normal_precision	=	53;
const int high_precision	=	256;

#define MAX_REFINEMENT_STEPS	49
#define MAX_DYN_REF				49

#define MAX_WIDTH				1e-12

#define THRESHOLD_HIGH			1e-12
#define THRESHOLD_LOW			1e-20

#define STOP_RATIO				0.99
#define ABS_STOP_RATIO			0.99

#define PN 						15 			// the number of digits printed
#define INVALID 				-1e10
#define UNBOUNDED				1e10

#define DC_THRESHOLD_SEARCH 	1e-6
#define DC_THRESHOLD_IMPROV		0.9

#define ID_PRE			0
#define QR_PRE			1

#define UNIFORM			0
#define MULTI			1

#define LAMBDA_DOWN		0.5
#define LAMBDA_UP		1.1

#define NAME_SIZE		100

#define UNSAFE			-1
#define SAFE			0
#define UNKNOWN			1

#define COMPLETED_UNSAFE		1
#define COMPLETED_SAFE			2
#define COMPLETED_UNKNOWN		3
#define UNCOMPLETED_SAFE		4
#define UNCOMPLETED_UNSAFE		5
#define UNCOMPLETED_UNKNOWN		6

#define PLOT_GNUPLOT	0
#define PLOT_MATLAB		1

#define PLOT_INTERVAL	0
#define PLOT_OCTAGON	1
#define PLOT_GRID		2

#define INTERVAL_AGGREG	0
#define PARA_AGGREG		1
#define PCA_AGGREG		2

#define MSG_SIZE		100
#define NUM_LENGTH		50

#define ONLY_PICARD				1
#define LOW_DEGREE				2
#define HIGH_DEGREE				3
#define NONPOLY_TAYLOR			4
#define LTI						5
#define LTV						6
#define ONLY_PICARD_SYMB		7
#define NONPOLY_TAYLOR_SYMB		8

#define POLY_DYN		1
#define NONPOLY_DYN		2
#define LTI_DYN			3
#define LTV_DYN			4

#define REFINEMENT_PREC	1e-5

#define RESET_COLOR		"\033[0m"
#define BLACK_COLOR		"\033[30m"
#define RED_COLOR		"\033[31m"
#define GREEN_COLOR		"\033[32m"
#define BLUE_COLOR		"\033[34m"
#define BOLD_FONT		"\e[1m"

const char str_pi_up[]	=	"3.14159265358979323846264338327950288419716939937511";
const char str_pi_lo[]	=	"3.14159265358979323846264338327950288419716939937510";

const char outputDir[] = "./outputs/";
const char imageDir[] = "./images/";
const char counterexampleDir[] = "./counterexamples/";
const char local_var_name[] = "local_var_";

const char str_prefix_taylor_picard[] = "taylor picard { ";
const char str_prefix_taylor_remainder[] = "taylor remainder { ";
const char str_prefix_taylor_polynomial[] = "taylor polynomial { ";
const char str_prefix_center[] = "nonpolynomial center { ";

const char str_prefix_combination_picard[] = "combination picard { ";
const char str_prefix_combination_remainder[] = "combination remainder { ";
const char str_prefix_combination_polynomial[] = "combination polynomial { ";

const char str_prefix_univariate_polynomial[] = "univariate polynomial { ";
const char str_prefix_multivariate_polynomial[] = "multivariate polynomial { ";

const char str_prefix_expression[] = "expression { ";

const char str_prefix_matrix[] = "matrix { ";

const char str_suffix[] = " }";

const char str_counterexample_dumping_name_suffix[] = ".counterexample";

extern int lineNum;

extern void parseODE();

#endif /* INCLUDE_H_ */


text:
#include "../EngineTest.h"

#include <Catch2>

#include <gsl/gsl_math.h>

TEST_CASE("Lag-1 Autocorrelate Function Evaluation Tests", "[lag1]") {

    // SECTION("Empty Test"){
    //     requireIsEqual("lag1()", "Insufficient Number of Arguments for Function: lag1");
    // }

    SECTION("`lag1` Test 1"){
        requireIsEqual("lag1(3.95, 6.81, -3.05, -1.64, -3.5, 4.95, -8.23, 2.92, 6.06)", -0.2164409953);
    }

    SECTION("`autocorr` Test 2"){
        requireIsEqual("autocorr(1.61, 7.86, -2.26, -5.62, 8.1, -4.55, 8.97)", -0.5295711223);
    }

    SECTION("`autocorr` Test 3"){
        requireIsEqual("autocorr(-2.79, -2.1)", -0.5);
    }

    SECTION("`autocorr` Test 4"){
        requireIsEqual("autocorr(7.27, -7.71, -6.63, -6.2, -0.98, 4.56)", 0.0020706076);
    }

    SECTION("`lag1` Test 5"){
        requireIsEqual("lag1(8.81, -2.71)", -0.5);
    }

    SECTION("`lag1` Test 6"){
        requireIsEqual("lag1(5.72, -3.33, 9.2, 2.92, -9.47, -8.05, -3.63)", 0.153525765);
    }

    SECTION("`autocorr` Test 7"){
        requireIsEqual("autocorr(1.85, 6.89, -9.39, 5.94, -7.5)", -0.657389584);
    }

    SECTION("`autocorr` Test 8"){
        requireIsEqual("autocorr(-5.14, -7.63, 5.89, 0.62, -1.38, 8.92, 7.13)", 0.1620685975);
    }

    SECTION("`autocorr` Test 9"){
        requireIsEqual("autocorr(-3.61, 8.1, -9.53, -6.09, -1.66, 2.79, 7.55)", -0.0886150351);
    }

    SECTION("`lag1` Test 10"){
        requireIsEqual("lag1(-0.76, 0.63)", -0.5);
    }

}

max_stars_repo_name: antoniojkim/CalcPlusPlus

max_issues_repo_name: antoniojkim/CalcPlusPlus

max_forks_repo_name: antoniojkim/CalcPlusPlus