text:
#include <complex.h>
#include <math.h>
//#include <gsl/gsl_sf.h>
#include <stdio.h>
#include <stdlib.h>
#include "specfunc.h"

/*
 * to compile as a library (that can be used by specfunc.py):
 *
 *      cc -fPIC -c -O2 specfunc.c -o specfunc.o
 *      cc --shared specfunc.o -o libspecfunc.so
 *
 * if the gamma-function stuff an dthe parabolic cylinder function is commented
 * in, you need
 *
 *      cc --shared specfunc.o -o libspecfunc.so -lgslcblas -lgsl
 *
 * instead
 *
 * hyp2f1 and hyp1f1 seem to work well, the parabolic cylinder functions pcfd
 * are problematic for anything but small arguments
 *
 * felix, june 2014
 */

#define PREC_WARN_LIMIT 1e99 // warn if individual terms in the series are larger than this
// XXX 4. 12. 14: set to high value because contrary to what I thought I had
// understood, there seem to be many perfectly fine cases where individual
// terms in the sum are higher than 10^15

// hypergeometric function for |z| < 1
double complex hyp2f1(const double complex a, const double complex b,
        const double complex c, const double complex z, sf_prms_and_info_t* p)
{
    double frac = p->tol;
    double complex n = 0;
    double complex summand = a*b/c * z;
    double complex sum = 1. + summand;

    p->prec_warning = 0;
    for (n = 2; (creal(n) < p->max_iter) && (frac >= p->tol); n++)
    {
        summand *= (a+n-1)*(b+n-1)/(c+n-1) * z/n;
        p->prec_warning = p->prec_warning | (creal(summand) > PREC_WARN_LIMIT)
                                          | (cimag(summand) > PREC_WARN_LIMIT);
        sum += summand;
        frac = fabs(creal(summand/sum)) + fabs(cimag((summand/sum)));
    }
    p->iters_needed = (int)creal(n);
    p->tol_achieved = frac;
    return sum;
}

// confluent hypergeometric function
double complex hyp1f1(const double complex a, const double complex b,
        const double complex z, sf_prms_and_info_t* p)
{
    double frac = p->tol;
    double complex n = 0;
    double complex summand = a/b * z;
    double complex sum = 1. + summand;
    p->prec_warning = 0;
    for (n = 2; (creal(n) < p->max_iter) && (frac >= p->tol) ; n++)
    {
        summand = summand/(b+n-1) * (a+n-1) * (z/n);
        p->prec_warning = p->prec_warning | (creal(summand) > PREC_WARN_LIMIT)
                                          | (cimag(summand) > PREC_WARN_LIMIT);
        sum += summand;
        frac = fabs(creal(summand/sum)) + fabs(cimag((summand/sum)));
    }
    p->iters_needed = (int)creal(n);
    p->tol_achieved = frac;
    return sum;
}

void hyp1f1_a_arr(const double complex* a, const double complex b,
        const double complex z, double complex* out,
        size_t n_arr, sf_prms_and_info_t* p)
{
    double complex f = 0.;
    double frac = 0.;
    double maxfrac = p->tol;
    double complex n = 0;
    double complex* summand =
        (double complex*)malloc(sizeof(double complex) * n_arr);
    int i = 0;

    for (i = 0; i < n_arr; i++)
    {
        summand[i] = a[i]/b * z;
        out[i] = 1. + summand[i];
    }

    p->prec_warning = 0;
    for (n = 2; (creal(n) < p->max_iter) && (maxfrac >= p->tol) ; n++)
    {
        maxfrac = 0;
        f = 1./(b+n-1) * z / n; // same for all a
        for (i = 0; i < n_arr; i++)
        {
            summand[i] *= f * (a[i]+n-1);
            p->prec_warning = p->prec_warning | (creal(summand[i]) > PREC_WARN_LIMIT)
                                              | (cimag(summand[i]) > PREC_WARN_LIMIT);
            out[i] += summand[i];
            frac = fabs(creal(summand[i]/out[i]))
                + fabs(cimag((summand[i]/out[i])));
            maxfrac = frac > maxfrac ? frac : maxfrac;
        }
    }
    p->iters_needed = (int)creal(n);
    p->tol_achieved = maxfrac;
    free(summand);
}

void hyp1f1_b_arr(const double complex a, const double complex* b,
        const double complex z, double complex* out,
        size_t n_arr, sf_prms_and_info_t* p)
{
    double complex f = 0.;
    double frac = 0.;
    double maxfrac = p->tol;
    double complex n = 0;
    double complex* summand =
        (double complex*)malloc(sizeof(double complex) * n_arr);
    int i = 0;

    for (i = 0; i < n_arr; i++)
    {
        summand[i] = a/b[i] * z;
        out[i] = 1. + summand[i];
    }

    p->prec_warning = 0;
    for (n = 2; (creal(n) < p->max_iter) && (maxfrac >= p->tol) ; n++)
    {
        maxfrac = 0;
        f = (a+n-1) * z / n; // same for all b
        for (i = 0; i < n_arr; i++)
        {
            summand[i] *= f / (b[i]+n-1);
            p->prec_warning = p->prec_warning | (creal(summand[i]) > PREC_WARN_LIMIT)
                                              | (cimag(summand[i]) > PREC_WARN_LIMIT);
            out[i] += summand[i];
            frac = fabs(creal(summand[i]/out[i]))
                + fabs(cimag((summand[i]/out[i])));
            maxfrac = frac > maxfrac ? frac : maxfrac;
        }
    }
    p->iters_needed = (int)creal(n);
    p->tol_achieved = maxfrac;
    free(summand);
}

void hyp1f1_z_arr(const double complex a, const double complex b,
        const double complex* z, double complex* out,
        size_t n_arr, sf_prms_and_info_t* p)
{
    double complex f = 0.;
    double frac = 0.;
    double maxfrac = p->tol;
    double complex n = 0;
    double complex* summand =
        (double complex*)malloc(sizeof(double complex) * n_arr);
    int i = 0;

    for (i = 0; i < n_arr; i++)
    {
        summand[i] = a/b * z[i];
        out[i] = 1. + summand[i];
    }

    p->prec_warning = 0;
    for (n = 2; (creal(n) < p->max_iter) && (maxfrac >= p->tol) ; n++)
    {
        maxfrac = 0;
        f = (a+n-1)/(b+n-1) / n; // same for all z
        for (i = 0; i < n_arr; i++)
        {
            summand[i] *= f * z[i];
            p->prec_warning = p->prec_warning | (creal(summand[i]) > PREC_WARN_LIMIT)
                                              | (cimag(summand[i]) > PREC_WARN_LIMIT);
            out[i] += summand[i];
            frac = fabs(creal(summand[i]/out[i]))
                + fabs(cimag((summand[i]/out[i])));
            maxfrac = frac > maxfrac ? frac : maxfrac;
        }
    }
    p->iters_needed = (int)creal(n);
    p->tol_achieved = maxfrac;
    free(summand);
}

void hyp1f1_all_arr(const double complex* a, const double complex* b,
        const double complex* z, double complex* out,
        size_t n_arr, sf_prms_and_info_t* p)
{
    int i = 0;
    for (i = 0; i < n_arr; i++)
        out[i] = hyp1f1(a[i], b[i], z[i], p);
}

/*
double complex cgamma(const double complex z)
{
    gsl_sf_result lnabs = {0, 0};
    gsl_sf_result arg = {0, 0};
    gsl_sf_lngamma_complex_e(creal(z), cimag(z), &lnabs, &arg);
    return cexp(lnabs.val+1j*arg.val);
}

// parabolic cylinder function -- this is only ok for small parameters
double complex pcfd(const double complex nu, const double complex z,
        sf_prms_and_info_t* p)
{
    return cpow(2, nu/2)*cexp(-z*z/4) * csqrt(M_PI) *
        (1./cgamma((1-nu)/2) * hyp1f1(-nu/2, 1./2, z*z/2, p)
            - csqrt(2.)*z/cgamma(-nu/2)
            * hyp1f1((1-nu)/2, 3./2, z*z/2, p));
}

void pcfd_nu_arr(const double complex* nu, const double complex z,
        double complex* out, size_t n_arr, sf_prms_and_info_t* p)
{
    int i = 0;
    for (i = 0; i < n_arr; i++)
    {
        out[i] = pcfd(nu[i], z, p);
    }
}

void pcfd_z_arr(const double complex nu, const double complex* z,
        double complex* out, size_t n_arr, sf_prms_and_info_t* p)
{
    int i = 0;
    for (i = 0; i < n_arr; i++)
    {
        out[i] = pcfd(nu, z[i], p);
    }
}
*/

void hyp2f1_a_arr(const double complex* a, const double complex b,
        const double complex c, const double complex z, double complex* out,
        size_t n_arr, sf_prms_and_info_t* p)
{
    double complex f = 0.;
    double frac = 0.;
    double maxfrac = p->tol;
    double complex n = 0;
    double complex* summand =
        (double complex*)malloc(sizeof(double complex) * n_arr);
    int i = 0;

    for (i = 0; i < n_arr; i++)
    {
        summand[i] = a[i]*b/c * z;
        out[i] = 1. + summand[i];
    }

    p->prec_warning = 0;
    for (n = 2; (creal(n) < p->max_iter) && (maxfrac >= p->tol) ; n++)
    {
        maxfrac = 0;
        f = (b+n-1)/(c+n-1) * z / n; // same for all a
        for (i = 0; i < n_arr; i++)
        {
            summand[i] *= f * (a[i]+n-1);
            p->prec_warning = p->prec_warning | (creal(summand[i]) > PREC_WARN_LIMIT)
                                              | (cimag(summand[i]) > PREC_WARN_LIMIT);
            out[i] += summand[i];
            frac = fabs(creal(summand[i]/out[i]))
                + fabs(cimag((summand[i]/out[i])));
            maxfrac = frac > maxfrac ? frac : maxfrac;
        }
    }
    p->iters_needed = (int)creal(n);
    p->tol_achieved = maxfrac;
    free(summand);
}

void hyp2f1_b_arr(const double complex a, const double complex* b,
        const double complex c, const double complex z, double complex* out,
        size_t n_arr, sf_prms_and_info_t* p)
{
    hyp2f1_a_arr(b, a, c, z, out, n_arr, p);
}

void hyp2f1_c_arr(const double complex a, const double complex b,
        const double complex* c, const double complex z, double complex* out,
        size_t n_arr, sf_prms_and_info_t* p)
{
    double complex f = 0.;
    double frac = 0.;
    double maxfrac = p->tol;
    double complex n = 0;
    double complex* summand =
        (double complex*)malloc(sizeof(double complex) * n_arr);
    int i = 0;

    for (i = 0; i < n_arr; i++)
    {
        summand[i] = a*b/c[i] * z;
        out[i] = 1. + summand[i];
    }

    p->prec_warning = 0;
    for (n = 2; (creal(n) < p->max_iter) && (maxfrac >= p->tol) ; n++)
    {
        maxfrac = 0;
        f = (a+n-1)*(b+n-1) * z / n; // same for all a
        for (i = 0; i < n_arr; i++)
        {
            summand[i] *= f / (c[i]+n-1);
            p->prec_warning = p->prec_warning | (creal(summand[i]) > PREC_WARN_LIMIT)
                                              | (cimag(summand[i]) > PREC_WARN_LIMIT);
            out[i] += summand[i];
            frac = fabs(creal(summand[i]/out[i]))
                + fabs(cimag((summand[i]/out[i])));
            maxfrac = frac > maxfrac ? frac : maxfrac;
        }
    }
    p->iters_needed = (int)creal(n);
    p->tol_achieved = maxfrac;
    free(summand);
}

// hyp2d1 for an array of z values
void hyp2f1_z_arr(const double complex a, const double complex b,
        const double complex c, const double complex* z, double complex* out,
        size_t n_arr, sf_prms_and_info_t* p)
{
    double complex f = 0.;
    double frac = 0.;
    double maxfrac = p->tol;
    double complex n = 0;
    double complex* summand =
        (double complex*)malloc(sizeof(double complex) * n_arr);
    int i = 0;

    for (i = 0; i < n_arr; i++)
    {
        summand[i] = a*b/c * z[i];
        out[i] = 1. + summand[i];
    }

    p->prec_warning = 0;
    for (n = 2; (creal(n) < p->max_iter) && (maxfrac >= p->tol) ; n++)
    {
        maxfrac = 0;
        f = (a+n-1)*(b+n-1)/(c+n-1) / n; // same for all zs
        for (i = 0; i < n_arr; i++)
        {
            summand[i] *= f * z[i];
            p->prec_warning = p->prec_warning | (creal(summand[i]) > PREC_WARN_LIMIT)
                                              | (cimag(summand[i]) > PREC_WARN_LIMIT);
            out[i] += summand[i];
            frac = fabs(creal(summand[i]/out[i]))
                + fabs(cimag((summand[i]/out[i])));
            maxfrac = frac > maxfrac ? frac : maxfrac;
        }
    }
    p->iters_needed = (int)creal(n);
    p->tol_achieved = maxfrac;
    free(summand);
}

void hyp2f1_all_arr(const double complex* a, const double complex* b,
        const double complex* c, const double complex* z, double complex* out,
        size_t n_arr, sf_prms_and_info_t* p)
{
    int i = 0;
    for (i = 0; i < n_arr; i++)
        out[i] = hyp2f1(a[i], b[i], c[i], z[i], p);
}


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

#ifndef __GSL_MATRIX_CHAR_H__
#define __GSL_MATRIX_CHAR_H__

#if !defined( GSL_FUN )
#  if !defined( GSL_DLL )
#    define GSL_FUN extern
#  elif defined( BUILD_GSL_DLL )
#    define GSL_FUN extern __declspec(dllexport)
#  else
#    define GSL_FUN extern __declspec(dllimport)
#  endif
#endif

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_inline.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_vector_char.h>
#include <gsl/gsl_blas_types.h>

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

__BEGIN_DECLS

typedef struct
{
  size_t size1;
  size_t size2;
  size_t tda;
  char * data;
  gsl_block_char * block;
  int owner;
} gsl_matrix_char;

typedef struct
{
  gsl_matrix_char matrix;
} _gsl_matrix_char_view;

typedef _gsl_matrix_char_view gsl_matrix_char_view;

typedef struct
{
  gsl_matrix_char matrix;
} _gsl_matrix_char_const_view;

typedef const _gsl_matrix_char_const_view gsl_matrix_char_const_view;

/* Allocation */

GSL_FUN gsl_matrix_char *
gsl_matrix_char_alloc (const size_t n1, const size_t n2);

GSL_FUN gsl_matrix_char *
gsl_matrix_char_calloc (const size_t n1, const size_t n2);

GSL_FUN gsl_matrix_char *
gsl_matrix_char_alloc_from_block (gsl_block_char * b,
                                   const size_t offset,
                                   const size_t n1,
                                   const size_t n2,
                                   const size_t d2);

GSL_FUN gsl_matrix_char *
gsl_matrix_char_alloc_from_matrix (gsl_matrix_char * m,
                                    const size_t k1,
                                    const size_t k2,
                                    const size_t n1,
                                    const size_t n2);

GSL_FUN gsl_vector_char *
gsl_vector_char_alloc_row_from_matrix (gsl_matrix_char * m,
                                        const size_t i);

GSL_FUN gsl_vector_char *
gsl_vector_char_alloc_col_from_matrix (gsl_matrix_char * m,
                                        const size_t j);

GSL_FUN void gsl_matrix_char_free (gsl_matrix_char * m);

/* Views */

GSL_FUN _gsl_matrix_char_view
gsl_matrix_char_submatrix (gsl_matrix_char * m,
                            const size_t i, const size_t j,
                            const size_t n1, const size_t n2);

GSL_FUN _gsl_vector_char_view
gsl_matrix_char_row (gsl_matrix_char * m, const size_t i);

GSL_FUN _gsl_vector_char_view
gsl_matrix_char_column (gsl_matrix_char * m, const size_t j);

GSL_FUN _gsl_vector_char_view
gsl_matrix_char_diagonal (gsl_matrix_char * m);

GSL_FUN _gsl_vector_char_view
gsl_matrix_char_subdiagonal (gsl_matrix_char * m, const size_t k);

GSL_FUN _gsl_vector_char_view
gsl_matrix_char_superdiagonal (gsl_matrix_char * m, const size_t k);

GSL_FUN _gsl_vector_char_view
gsl_matrix_char_subrow (gsl_matrix_char * m, const size_t i,
                         const size_t offset, const size_t n);

GSL_FUN _gsl_vector_char_view
gsl_matrix_char_subcolumn (gsl_matrix_char * m, const size_t j,
                            const size_t offset, const size_t n);

GSL_FUN _gsl_matrix_char_view
gsl_matrix_char_view_array (char * base,
                             const size_t n1,
                             const size_t n2);

GSL_FUN _gsl_matrix_char_view
gsl_matrix_char_view_array_with_tda (char * base,
                                      const size_t n1,
                                      const size_t n2,
                                      const size_t tda);

GSL_FUN _gsl_matrix_char_view
gsl_matrix_char_view_vector (gsl_vector_char * v,
                              const size_t n1,
                              const size_t n2);

GSL_FUN _gsl_matrix_char_view
gsl_matrix_char_view_vector_with_tda (gsl_vector_char * v,
                                       const size_t n1,
                                       const size_t n2,
                                       const size_t tda);

GSL_FUN _gsl_matrix_char_const_view
gsl_matrix_char_const_submatrix (const gsl_matrix_char * m,
                                  const size_t i, const size_t j,
                                  const size_t n1, const size_t n2);

GSL_FUN _gsl_vector_char_const_view
gsl_matrix_char_const_row (const gsl_matrix_char * m,
                            const size_t i);

GSL_FUN _gsl_vector_char_const_view
gsl_matrix_char_const_column (const gsl_matrix_char * m,
                               const size_t j);

GSL_FUN _gsl_vector_char_const_view
gsl_matrix_char_const_diagonal (const gsl_matrix_char * m);

GSL_FUN _gsl_vector_char_const_view
gsl_matrix_char_const_subdiagonal (const gsl_matrix_char * m,
                                    const size_t k);

GSL_FUN _gsl_vector_char_const_view
gsl_matrix_char_const_superdiagonal (const gsl_matrix_char * m,
                                      const size_t k);

GSL_FUN _gsl_vector_char_const_view
gsl_matrix_char_const_subrow (const gsl_matrix_char * m, const size_t i,
                               const size_t offset, const size_t n);

GSL_FUN _gsl_vector_char_const_view
gsl_matrix_char_const_subcolumn (const gsl_matrix_char * m, const size_t j,
                                  const size_t offset, const size_t n);

GSL_FUN _gsl_matrix_char_const_view
gsl_matrix_char_const_view_array (const char * base,
                                   const size_t n1,
                                   const size_t n2);

GSL_FUN _gsl_matrix_char_const_view
gsl_matrix_char_const_view_array_with_tda (const char * base,
                                            const size_t n1,
                                            const size_t n2,
                                            const size_t tda);

GSL_FUN _gsl_matrix_char_const_view
gsl_matrix_char_const_view_vector (const gsl_vector_char * v,
                                    const size_t n1,
                                    const size_t n2);

GSL_FUN _gsl_matrix_char_const_view
gsl_matrix_char_const_view_vector_with_tda (const gsl_vector_char * v,
                                             const size_t n1,
                                             const size_t n2,
                                             const size_t tda);

/* Operations */

GSL_FUN void gsl_matrix_char_set_zero (gsl_matrix_char * m);
GSL_FUN void gsl_matrix_char_set_identity (gsl_matrix_char * m);
GSL_FUN void gsl_matrix_char_set_all (gsl_matrix_char * m, char x);

GSL_FUN int gsl_matrix_char_fread (FILE * stream, gsl_matrix_char * m) ;
GSL_FUN int gsl_matrix_char_fwrite (FILE * stream, const gsl_matrix_char * m) ;
GSL_FUN int gsl_matrix_char_fscanf (FILE * stream, gsl_matrix_char * m);
GSL_FUN int gsl_matrix_char_fprintf (FILE * stream, const gsl_matrix_char * m, const char * format);

GSL_FUN int gsl_matrix_char_memcpy(gsl_matrix_char * dest, const gsl_matrix_char * src);
GSL_FUN int gsl_matrix_char_swap(gsl_matrix_char * m1, gsl_matrix_char * m2);
GSL_FUN int gsl_matrix_char_tricpy(CBLAS_UPLO_t Uplo, CBLAS_DIAG_t Diag, gsl_matrix_char * dest, const gsl_matrix_char * src);

GSL_FUN int gsl_matrix_char_swap_rows(gsl_matrix_char * m, const size_t i, const size_t j);
GSL_FUN int gsl_matrix_char_swap_columns(gsl_matrix_char * m, const size_t i, const size_t j);
GSL_FUN int gsl_matrix_char_swap_rowcol(gsl_matrix_char * m, const size_t i, const size_t j);
GSL_FUN int gsl_matrix_char_transpose (gsl_matrix_char * m);
GSL_FUN int gsl_matrix_char_transpose_memcpy (gsl_matrix_char * dest, const gsl_matrix_char * src);
GSL_FUN int gsl_matrix_char_transpose_tricpy (CBLAS_UPLO_t Uplo_src, CBLAS_DIAG_t Diag, gsl_matrix_char * dest, const gsl_matrix_char * src);

GSL_FUN char gsl_matrix_char_max (const gsl_matrix_char * m);
GSL_FUN char gsl_matrix_char_min (const gsl_matrix_char * m);
GSL_FUN void gsl_matrix_char_minmax (const gsl_matrix_char * m, char * min_out, char * max_out);

GSL_FUN void gsl_matrix_char_max_index (const gsl_matrix_char * m, size_t * imax, size_t *jmax);
GSL_FUN void gsl_matrix_char_min_index (const gsl_matrix_char * m, size_t * imin, size_t *jmin);
GSL_FUN void gsl_matrix_char_minmax_index (const gsl_matrix_char * m, size_t * imin, size_t * jmin, size_t * imax, size_t * jmax);

GSL_FUN int gsl_matrix_char_equal (const gsl_matrix_char * a, const gsl_matrix_char * b);

GSL_FUN int gsl_matrix_char_isnull (const gsl_matrix_char * m);
GSL_FUN int gsl_matrix_char_ispos (const gsl_matrix_char * m);
GSL_FUN int gsl_matrix_char_isneg (const gsl_matrix_char * m);
GSL_FUN int gsl_matrix_char_isnonneg (const gsl_matrix_char * m);

GSL_FUN int gsl_matrix_char_add (gsl_matrix_char * a, const gsl_matrix_char * b);
GSL_FUN int gsl_matrix_char_sub (gsl_matrix_char * a, const gsl_matrix_char * b);
GSL_FUN int gsl_matrix_char_mul_elements (gsl_matrix_char * a, const gsl_matrix_char * b);
GSL_FUN int gsl_matrix_char_div_elements (gsl_matrix_char * a, const gsl_matrix_char * b);
GSL_FUN int gsl_matrix_char_scale (gsl_matrix_char * a, const double x);
GSL_FUN int gsl_matrix_char_scale_rows (gsl_matrix_char * a, const gsl_vector_char * x);
GSL_FUN int gsl_matrix_char_scale_columns (gsl_matrix_char * a, const gsl_vector_char * x);
GSL_FUN int gsl_matrix_char_add_constant (gsl_matrix_char * a, const double x);
GSL_FUN int gsl_matrix_char_add_diagonal (gsl_matrix_char * a, const double x);

/***********************************************************************/
/* The functions below are obsolete                                    */
/***********************************************************************/
GSL_FUN int gsl_matrix_char_get_row(gsl_vector_char * v, const gsl_matrix_char * m, const size_t i);
GSL_FUN int gsl_matrix_char_get_col(gsl_vector_char * v, const gsl_matrix_char * m, const size_t j);
GSL_FUN int gsl_matrix_char_set_row(gsl_matrix_char * m, const size_t i, const gsl_vector_char * v);
GSL_FUN int gsl_matrix_char_set_col(gsl_matrix_char * m, const size_t j, const gsl_vector_char * v);
/***********************************************************************/

/* inline functions if you are using GCC */

GSL_FUN INLINE_DECL char   gsl_matrix_char_get(const gsl_matrix_char * m, const size_t i, const size_t j);
GSL_FUN INLINE_DECL void    gsl_matrix_char_set(gsl_matrix_char * m, const size_t i, const size_t j, const char x);
GSL_FUN INLINE_DECL char * gsl_matrix_char_ptr(gsl_matrix_char * m, const size_t i, const size_t j);
GSL_FUN INLINE_DECL const char * gsl_matrix_char_const_ptr(const gsl_matrix_char * m, const size_t i, const size_t j);

#ifdef HAVE_INLINE
INLINE_FUN
char
gsl_matrix_char_get(const gsl_matrix_char * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_VAL("first index out of range", GSL_EINVAL, 0) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_VAL("second index out of range", GSL_EINVAL, 0) ;
        }
    }
#endif
  return m->data[i * m->tda + j] ;
}

INLINE_FUN
void
gsl_matrix_char_set(gsl_matrix_char * m, const size_t i, const size_t j, const char x)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_VOID("first index out of range", GSL_EINVAL) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_VOID("second index out of range", GSL_EINVAL) ;
        }
    }
#endif
  m->data[i * m->tda + j] = x ;
}

INLINE_FUN
char *
gsl_matrix_char_ptr(gsl_matrix_char * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
        }
    }
#endif
  return (char *) (m->data + (i * m->tda + j)) ;
}

INLINE_FUN
const char *
gsl_matrix_char_const_ptr(const gsl_matrix_char * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
        }
    }
#endif
  return (const char *) (m->data + (i * m->tda + j)) ;
}

#endif

__END_DECLS

#endif /* __GSL_MATRIX_CHAR_H__ */

ext: h


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

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

#include "demof.h"
#include "demof.c"

int
main ()
{
  int status;
  int iterations = 0, max_iterations = 100;
  gsl_root_fdfsolver *s;
  double x0, x = 5.0, r_expected = sqrt (5.0);
  gsl_function_fdf FDF;
  struct quadratic_params params = {1.0, 0.0, -5.0};

  FDF.f = &quadratic;
  FDF.df = &quadratic_deriv;
  FDF.fdf = &quadratic_fdf;
  FDF.params = &params;

  s = gsl_root_fdfsolver_alloc (gsl_root_fdfsolver_newton);
  gsl_root_fdfsolver_set (s, &FDF, x);

  printf ("using %s method\n", gsl_root_fdfsolver_name (s));

  printf ("%-5s %10s %10s %10s %10s\n",
          "iter", "root", "actual", "err", "err(est)");
  do
    {
      iterations++;
      status = gsl_root_fdfsolver_iterate (s);
      x0 = x;
      x = gsl_root_fdfsolver_root (s);
      status = gsl_root_test_delta (x, x0, 0, 0.001);

      if (status == GSL_SUCCESS)
        printf ("Converged:\n");

      printf ("%5d %10.7f %10.7f %+10.7f %10.7f\n",
              iterations, x, r_expected, x - r_expected, x - x0);
    }
  while (status == GSL_CONTINUE && iterations < max_iterations);

}

ext: c

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim


text:
/* rng/taus.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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

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

/* This is a maximally equidistributed combined Tausworthe
   generator. The sequence is,

   x_n = (s1_n ^ s2_n ^ s3_n)

   s1_{n+1} = (((s1_n & 4294967294) <<12) ^ (((s1_n <<13) ^ s1_n) >>19))
   s2_{n+1} = (((s2_n & 4294967288) << 4) ^ (((s2_n << 2) ^ s2_n) >>25))
   s3_{n+1} = (((s3_n & 4294967280) <<17) ^ (((s3_n << 3) ^ s3_n) >>11))

   computed modulo 2^32. In the three formulas above '^' means
   exclusive-or (C-notation), not exponentiation. Note that the
   algorithm relies on the properties of 32-bit unsigned integers (it
   is formally defined on bit-vectors of length 32). I have added a
   bitmask to make it work on 64 bit machines.

   We initialize the generator with s1_1 .. s3_1 = s_n MOD m, where
   s_n = (69069 * s_{n-1}) mod 2^32, and s_0 = s is the user-supplied
   seed.

   The theoretical value of x_{10007} is 2733957125. The subscript
   10007 means (1) seed the generator with s=1 (2) do six warm-up
   iterations, (3) then do 10000 actual iterations.

   The period of this generator is about 2^88.

   From: P. L'Ecuyer, "Maximally Equidistributed Combined Tausworthe
   Generators", Mathematics of Computation, 65, 213 (1996), 203--213.

   This is available on the net from L'Ecuyer's home page,

   http://www.iro.umontreal.ca/~lecuyer/myftp/papers/tausme.ps
   ftp://ftp.iro.umontreal.ca/pub/simulation/lecuyer/papers/tausme.ps

   Update: April 2002

   There is an erratum in the paper "Tables of Maximally
   Equidistributed Combined LFSR Generators", Mathematics of
   Computation, 68, 225 (1999), 261--269:
   http://www.iro.umontreal.ca/~lecuyer/myftp/papers/tausme2.ps

        ... the k_j most significant bits of z_j must be non-
        zero, for each j. (Note: this restriction also applies to the
        computer code given in [4], but was mistakenly not mentioned in
        that paper.)

   This affects the seeding procedure by imposing the requirement
   s1 > 1, s2 > 7, s3 > 15.

   The generator taus2 has been added to satisfy this requirement.
   The original taus generator is unchanged.

   Update: November 2002

   There was a bug in the correction to the seeding procedure for s2.
   It affected the following seeds 254679140 1264751179 1519430319
   2274823218 2529502358 3284895257 3539574397 (s2 < 8).

*/

static inline unsigned long int taus_get (void *vstate);
static double taus_get_double (void *vstate);
static void taus_set (void *state, unsigned long int s);

typedef struct
  {
    unsigned long int s1, s2, s3;
  }
taus_state_t;

static inline unsigned long
taus_get (void *vstate)
{
  taus_state_t *state = (taus_state_t *) vstate;

#define MASK 0xffffffffUL
#define TAUSWORTHE(s,a,b,c,d) (((s &c) <<d) &MASK) ^ ((((s <<a) &MASK)^s) >>b)

  state->s1 = TAUSWORTHE (state->s1, 13, 19, 4294967294UL, 12);
  state->s2 = TAUSWORTHE (state->s2, 2, 25, 4294967288UL, 4);
  state->s3 = TAUSWORTHE (state->s3, 3, 11, 4294967280UL, 17);

  return (state->s1 ^ state->s2 ^ state->s3);
}

static double
taus_get_double (void *vstate)
{
  return taus_get (vstate) / 4294967296.0 ;
}

static void
taus_set (void *vstate, unsigned long int s)
{
  taus_state_t *state = (taus_state_t *) vstate;

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

#define LCG(n) ((69069 * n) & 0xffffffffUL)
  state->s1 = LCG (s);
  state->s2 = LCG (state->s1);
  state->s3 = LCG (state->s2);

  /* "warm it up" */
  taus_get (state);
  taus_get (state);
  taus_get (state);
  taus_get (state);
  taus_get (state);
  taus_get (state);
  return;
}

static void
taus2_set (void *vstate, unsigned long int s)
{
  taus_state_t *state = (taus_state_t *) vstate;

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

#define LCG(n) ((69069 * n) & 0xffffffffUL)
  state->s1 = LCG (s);
  if (state->s1 < 2) state->s1 += 2UL;
  state->s2 = LCG (state->s1);
  if (state->s2 < 8) state->s2 += 8UL;
  state->s3 = LCG (state->s2);
  if (state->s3 < 16) state->s3 += 16UL;

  /* "warm it up" */
  taus_get (state);
  taus_get (state);
  taus_get (state);
  taus_get (state);
  taus_get (state);
  taus_get (state);
  return;
}

static const gsl_rng_type taus_type =
{"taus",                        /* name */
 0xffffffffUL,                  /* RAND_MAX */
 0,                             /* RAND_MIN */
 sizeof (taus_state_t),
 &taus_set,
 &taus_get,
 &taus_get_double};

const gsl_rng_type *gsl_rng_taus = &taus_type;

static const gsl_rng_type taus2_type =
{"taus2",                       /* name */
 0xffffffffUL,                  /* RAND_MAX */
 0,                             /* RAND_MIN */
 sizeof (taus_state_t),
 &taus2_set,
 &taus_get,
 &taus_get_double};

const gsl_rng_type *gsl_rng_taus2 = &taus2_type;

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
#ifndef _H_MAIN
#define _H_MAIN

#define CPP 0

#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <math.h>
#include <gsl/gsl_linalg.h>
#include <GL/glew.h>
#include <GLFW/glfw3.h>

#if CPP
#include <glm/glm.hpp>
#include <glm/gtx/transform.hpp>
#include <glm/gtc/matrix_transform.hpp>
#include <glm/gtc/type_ptr.hpp>
#endif

#define WW 640
#define WH 480
#define RAD(d) (d) * (M_PI/180)
#define ZEROCHK(i) ((i) == 0)?0:i

extern GLFWwindow *window;
extern GLuint tr_vb;
extern GLuint tr_c_vb;
extern GLuint va;
extern GLuint vs;
extern GLuint fs;
extern GLuint shader_bin;
uint64_t timer[10];

extern void cookTriangle(void);
extern void cookCube(void);
extern uint64_t getCycles(void);
extern char * readFile(char *);
extern void debugSC(GLuint);
extern uint64_t getCycles(void);

extern gsl_matrix * m_new_diag(uint32_t);
extern gsl_matrix * m_new(uint32_t, uint32_t);
extern void m_setT(gsl_matrix *, double, double, double, uint8_t);
extern void m_setRy(gsl_matrix *, double, uint8_t);
extern void m_mul(gsl_matrix *, gsl_matrix *, gsl_matrix *);
extern void m_print(gsl_matrix *, uint8_t, uint8_t);
extern void m_array(gsl_matrix *, uint8_t, uint8_t, float *);
extern void glmPerspective(double, double, double, double, gsl_matrix *);
extern void normalize(double *, uint8_t, double *);
extern double getVectorLength(double *, uint8_t);
extern void getCrossV3(double *, double *, double *);
extern double scalar(double *, double *, uint8_t);
extern void subVec(double *, double *, uint8_t, double *);
extern void glmLookAt(double *, double *, double *, gsl_matrix *);
extern void m_set_all(gsl_matrix *, double);

#include "lib.c"
#include "mat.c"

#endif


text:
/* ode-initval/rk2.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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

/* Runge-Kutta 2(3), Euler-Cauchy */

/* Author:  G. Jungman
 */

/* Reference: Abramowitz & Stegun, section 25.5. Runge-Kutta 2nd (25.5.7)
   and 3rd (25.5.8) order methods */

#include <config.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_odeiv.h>

#include "odeiv_util.h"

typedef struct
{
  double *k1;
  double *k2;
  double *k3;
  double *ytmp;
}
rk2_state_t;

static void *
rk2_alloc (size_t dim)
{
  rk2_state_t *state = (rk2_state_t *) malloc (sizeof (rk2_state_t));

  if (state == 0)
    {
      GSL_ERROR_NULL ("failed to allocate space for rk2_state", GSL_ENOMEM);
    }

  state->k1 = (double *) malloc (dim * sizeof (double));

  if (state->k1 == 0)
    {
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for k1", GSL_ENOMEM);
    }

  state->k2 = (double *) malloc (dim * sizeof (double));

  if (state->k2 == 0)
    {
      free (state->k1);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for k2", GSL_ENOMEM);
    }

  state->k3 = (double *) malloc (dim * sizeof (double));

  if (state->k3 == 0)
    {
      free (state->k2);
      free (state->k1);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for k2", GSL_ENOMEM);
    }

  state->ytmp = (double *) malloc (dim * sizeof (double));

  if (state->ytmp == 0)
    {
      free (state->k3);
      free (state->k2);
      free (state->k1);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for k2", GSL_ENOMEM);
    }

  return state;
}

static int
rk2_apply (void *vstate,
           size_t dim,
           double t,
           double h,
           double y[],
           double yerr[],
           const double dydt_in[],
           double dydt_out[],
           const gsl_odeiv_system * sys)
{
  rk2_state_t *state = (rk2_state_t *) vstate;

  size_t i;

  double *const k1 = state->k1;
  double *const k2 = state->k2;
  double *const k3 = state->k3;
  double *const ytmp = state->ytmp;

  /* k1 step */
  /* k1 = f(t,y) */

  if (dydt_in != NULL)
    {
      DBL_MEMCPY (k1, dydt_in, dim);
    }
  else
    {
      int s = GSL_ODEIV_FN_EVAL (sys, t, y, k1);

      if (s != GSL_SUCCESS)
	{
	  return s;
	}
    }

  /* k2 step */
  /* k2 = f(t + 0.5*h, y + 0.5*k1) */

  for (i = 0; i < dim; i++)
    {
      ytmp[i] = y[i] + 0.5 * h * k1[i];
    }

  {
    int s = GSL_ODEIV_FN_EVAL (sys, t + 0.5 * h, ytmp, k2);

    if (s != GSL_SUCCESS)
      {
	return s;
      }
  }

  /* k3 step */
  /* for 3rd order estimates, is used for error estimation
     k3 = f(t + h, y - k1 + 2*k2) */

  for (i = 0; i < dim; i++)
    {
      ytmp[i] = y[i] + h * (-k1[i] + 2.0 * k2[i]);
    }

  {
    int s = GSL_ODEIV_FN_EVAL (sys, t + h, ytmp, k3);

    if (s != GSL_SUCCESS)
      {
	return s;
      }
  }

  /* final sum */

  for (i = 0; i < dim; i++)
    {
      /* Save original values if derivative evaluation below fails */
      ytmp[i] = y[i];

      {
	const double ksum3 = (k1[i] + 4.0 * k2[i] + k3[i]) / 6.0;
	y[i] += h * ksum3;
      }
    }

  /* Derivatives at output */

  if (dydt_out != NULL)
    {
      int s = GSL_ODEIV_FN_EVAL (sys, t + h, y, dydt_out);

      if (s != GSL_SUCCESS)
	{
	  /* Restore original values */
	  DBL_MEMCPY (y, ytmp, dim);

	  return s;
	}
    }

  /* Error estimation */

  for (i = 0; i < dim; i++)
    {
      const double ksum3 = (k1[i] + 4.0 * k2[i] + k3[i]) / 6.0;
      yerr[i] = h * (k2[i] - ksum3);
    }

  return GSL_SUCCESS;
}

static int
rk2_reset (void *vstate, size_t dim)
{
  rk2_state_t *state = (rk2_state_t *) vstate;

  DBL_ZERO_MEMSET (state->k1, dim);
  DBL_ZERO_MEMSET (state->k2, dim);
  DBL_ZERO_MEMSET (state->k3, dim);
  DBL_ZERO_MEMSET (state->ytmp, dim);

  return GSL_SUCCESS;
}

static unsigned int
rk2_order (void *vstate)
{
  rk2_state_t *state = (rk2_state_t *) vstate;
  state = 0; /* prevent warnings about unused parameters */
  return 2;
}

static void
rk2_free (void *vstate)
{
  rk2_state_t *state = (rk2_state_t *) vstate;
  free (state->k1);
  free (state->k2);
  free (state->k3);
  free (state->ytmp);
  free (state);
}

static const gsl_odeiv_step_type rk2_type = { "rk2",    /* name */
  1,                           /* can use dydt_in */
  1,                           /* gives exact dydt_out */
  &rk2_alloc,
  &rk2_apply,
  &rk2_reset,
  &rk2_order,
  &rk2_free
};

const gsl_odeiv_step_type *gsl_odeiv_step_rk2 = &rk2_type;


text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <float.h>
#include "compearth.h"
#ifdef COMPEARTH_USE_MKL
#ifdef __clang__
#pragma clang diagnostic push
#pragma clang diagnostic ignored "-Wreserved-id-macro"
#pragma clang diagnostic ignored "-Wstrict-prototypes"
#endif
#include <mkl_cblas.h>
#ifdef __clang__
#pragma clang diagnostic pop
#endif
#else
#include <cblas.h>
#endif
#define deg 180.0/M_PI;

/*!
 * @brief Returns the angle between two vectors, in degrees.
 *
 * @param[in] n    Length of vectors.
 * @param[in] va   Initial vector.  This is an array of dimension [n].
 * @param[in] vb   Rotate vector.  This is an array of dimension [n].
 *
 * @result Rotation angle in degrees between va and vb.
 *
 * @author Ben Baker
 *
 * @copyright MIT
 *
 */
double compearth_matlab_fangle(const int n,
                               const double *__restrict__ va,
                               const double *__restrict__ vb)
{
    double angle, magVa, magVb, xden, xnum;
    xnum = cblas_ddot(n, va, 1, vb, 1);
    magVa = cblas_dnrm2(n, va, 1);
    magVb = cblas_dnrm2(n, vb, 1);
    xden = magVa*magVb;
    angle = (double) NAN;
    if (xden > 0.0){angle = acos(xnum/xden)*deg;}
    return angle;
}

max_stars_repo_name: OUCyf/mtbeach

max_issues_repo_name: carltape/mtbeach

max_forks_repo_name: carltape/mtbeach


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

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

/* Gaussian function */
double tgauss(double x,double sigma)
{
  double y,z,a=2.506628273;
  z = x / sigma;
  y = exp((double)-z*z*0.5)/(sigma * a);
  return y;
}

/* create a Gaussian kernel for convolution */
gsl_vector *CGaussKernel(double sigma)
{
  int i,dim,n;
  double x,u,step;
  gsl_vector *kernel;

  /* no temporal smoothing */
  if (sigma < 0.001) return NULL;

  dim  = 4.0 * sigma + 1;
  n    = 2*dim+1;
  step = 1;

  kernel = gsl_vector_alloc(n);

  x = -(float)dim;
  for (i=0; i<n; i++) {
    u = tgauss(x,sigma);
    gsl_vector_set(kernel,i,u);
    x += step;
  }
  return kernel;
}

/* fill gaussian matrix */
void GaussMatrix(gsl_matrix *S,gsl_vector *kernel)
{
  int dim = kernel->size/2;
  int i,j,k;
  gsl_matrix_set_zero(S);
  for (i=0; i<S->size1; i++) {
    double sum = 0;
    for (j=0; j<S->size2; j++) {
      k = i-j+dim;
      if (k < 0 || k >= kernel->size) continue;
      double x = gsl_vector_get(kernel,k);
      sum += x;
      gsl_matrix_set(S,i,j,x);
    }
    /* normalize */
    for (j=0; j<S->size2; j++) {
      double x = gsl_matrix_get(S,i,j);
      gsl_matrix_set(S,i,j,x/sum);
    }
  }
}

void VectorConvolve(const double *src,gsl_vector *dst,gsl_vector *kernel)
{
  int i,j,k,len,n,dim;
  double sum;

  dim = kernel->size/2;
  len = dst->size;
  n   = len - dim;

  for (i=dim; i<n; i++) {
    sum = 0;
    k=0;
    for (j=i-dim; j<=i+dim; j++) {
      sum += src[j] * kernel->data[k];
      k++;
    }
    dst->data[i] = sum;
  }

  /* Randbehandlung */
  for (i=0; i<dim; i++) {
    dst->data[i] = src[i];
  }
  for (i=n; i<len; i++) {
    dst->data[i] = src[i];
  }
}

max_stars_repo_name: zrajna/lipsia

max_issues_repo_name: zrajna/lipsia

max_forks_repo_name: zrajna/lipsia


text:
#include <stdio.h>
#include <getopt.h>
#include <omp.h>
#include <math.h>
#include <string.h>

#ifdef __APPLE__

#include <cblas.h>

#define set_num_threads(x) openblas_set_num_threads(x)
#define get_num_threads() openblas_get_num_threads()
#else
#include <mkl.h>
#define set_num_threads(x) mkl_set_num_threads(x)
#define get_num_threads() mkl_get_num_threads()
#endif

double norm_inf_mat_Y11(double *X, int n, int extra);

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

    int n, i, j;
    char *filename = NULL;
    double tolerance = 1e-8;
    int option = 0;
    int threads_num = 1;
    while ((option = getopt(argc, argv, "t:f:")) != -1) {
        switch (option) {
            case 'f':
                filename = optarg;
                break;
            case 't':
                threads_num = atoi(optarg);
                break;
            default:
                printf("Usage: mexp -f string \n");
                return 0;
        }
    }
    if (filename == NULL) {
        printf("Usage: mexp -t num_threads -f string \n");
        return 0;
    }

    double *A = NULL, *R = NULL, *B = NULL, *W = NULL;
    double *x, *r, *p, *w;
    FILE *file = fopen(filename, "rb");
    if (!file) {
        printf("Unable to open file!");
        return 0;
    }

    fread(&n, sizeof(int), 1, file);
    A = (double *) malloc((n + 1) * (n + 1) * sizeof(double));
    R = (double *) malloc((n + 1) * (n + 1) * sizeof(double));
    B = (double *) malloc((n + 1) * sizeof(double));
    W = (double *) malloc(n * sizeof(double));
    x = (double *) malloc(n * sizeof(double));
    r = (double *) malloc(n * sizeof(double));
    p = (double *) malloc(n * sizeof(double));
    w = (double *) malloc(n * sizeof(double));

    printf("set_threads_num:%d\t", threads_num++);

    fseek(file, 0, SEEK_SET);
    fread(&n, sizeof(int), 1, file);

    //fseek(file, sizeof(int), SEEK_SET);
    memset(A, 0, sizeof(double) * (n + 1) * (n + 1));
    memset(R, 0, sizeof(double) * (n + 1) * (n + 1));
    for (i = 0; i < n; ++i) {
        fread(A + i * (n + 1), sizeof(double), (size_t) n, file);
    }
    memset(B, 0, sizeof(double) * (n + 1));
    fread(B, sizeof(double), (size_t) n, file);

    set_num_threads(threads_num);
    omp_set_num_threads(threads_num);

    cblas_dcopy(n, B, 1, A + n, n + 1);

    double omp_time = omp_get_wtime();
    double lambda = norm_inf_mat_Y11(A, n, 1);

    //init Y matrix
    for (i = 0; i < n; ++i) {
        double *A_T = A + i * (n + 1);
#pragma omp parallel for
        for (j = 0; j < n; ++j) {
            A_T[j] = -A_T[j] / lambda;
        }
    }
#pragma omp parallel for
    for (i = 0; i < (n + 1) * (n + 1); i += n + 2) {
        A[i] += 1;
    }
#pragma omp parallel for
    for (i = n; i < (n + 1) * n; i += n + 1) {
        A[i] /= lambda;
    }
#pragma omp parallel for
    for (i = 0; i < n; ++i) {
        double *A_T = A + n * (n + 1);
        A_T[i] = 0.0;
    }
    A[(n + 1) * (n + 1)] = 1.0;

    //loop
    double a, b, *T;
    size_t index;
    int itern = 0;
    while (1) {
        a = norm_inf_mat_Y11(A, n, 1);
        index = cblas_idamax(n, A + n, n + 1);
        b = A[(index + 1) * (n + 1) - 1];
        if (a / b < tolerance)
            break;
        else {
            //R = A*A
            cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, n + 1, n + 1,
                        n + 1, 1.0, A, n + 1, A, n + 1, 0.0, R, n + 1);
            T = A;
            A = R;
            R = T;
            itern++;
        }
    }

    omp_time = omp_get_wtime() - omp_time;

    cblas_dcopy(n, A + n, n + 1, B, 1);

    /*
    * check result
    */
    for (i = 0; i < n; ++i) {
        if (fabs(B[i] - 1) > 1e-3) {
            printf("MEXP ERR %f\t", B[i]);
            break;
        }
    }
    printf("mexp_iter:%d\t", itern);
    printf("omp_time:%f\t", omp_time);

    fseek(file, 0, SEEK_SET);
    fread(&n, sizeof(int), 1, file);

    fread(A, sizeof(double), (size_t) (n * n), file);
    memset(B, 0, sizeof(double) * (n + 1));
    fread(B, sizeof(double), (size_t) n, file);
#pragma omp parallel for
    for (i = 0; i < n; ++i) {
        W[i] = 1.0;
        //B[i] /= A[i*(n+1)];
        //cblas_dscal(n,1/A[i*(n+1)],A+i*n,1);
    }

    omp_time = omp_get_wtime();
    int iteration = 0;

    double alpha, beta, rho_new = 0.0, rho_old;
    memset(x, 0, sizeof(double) * n);
    cblas_dcopy(n, B, 1, r, 1); //r <- B

    rho_new = cblas_ddot(n, r, 1, r, 1);
    double w_norm2 = sqrt(rho_new);
    while (1) {
        if (sqrt(rho_new) / w_norm2 < tolerance) {
            break;
        }
        iteration++;
        if (iteration == 1) {
            cblas_dcopy(n, r, 1, p, 1); //p <- r
        } else {
            beta = rho_new / rho_old;
            cblas_dscal(n, beta, p, 1); // p <- beta * p
            cblas_daxpy(n, 1.0, r, 1, p, 1); // p <- r + pcd
        }
        cblas_dsymv(CblasRowMajor, CblasUpper, n, 1.0, A, n, p, 1, 0.0, w, 1);
        alpha = rho_new / cblas_ddot(n, p, 1, w, 1);
        cblas_daxpy(n, alpha, p, 1, x, 1);
        cblas_daxpy(n, -alpha, w, 1, r, 1);  // r = r - alpha* w
        rho_old = rho_new;
        rho_new = cblas_ddot(n, r, 1, r, 1); //rho_new = r' dot r;
    }

    printf("cg_it: %d\t", iteration);
    printf("cg_omp_time: %f", omp_get_wtime() - omp_time);

    /*
    * check result
    */
    for (i = 0; i < n; ++i) {
        if (fabs(x[i] - 1) > 1e-3) {
            printf(" CG ERR %f\t", x[i]);
            break;
        }
    }
    printf("\n");

    free(A);
    free(B);
    free(W);
    free(R);
    free(x);
    free(r);
    free(p);
    free(w);
    fclose(file);
}

double norm_inf_mat_Y11(double *X, int n, int extra) {
    double lambda = 0.0, lambda_t;
    int i;
    for (i = 0; i < n; ++i) {
        lambda_t = cblas_dasum(n, X + (i * (n + extra)), 1);
        if (lambda_t > lambda)
            lambda = lambda_t;
    }
    return lambda;
}

max_stars_repo_name: baishuai/MEXP

max_issues_repo_name: baishuai/MEXP

max_forks_repo_name: baishuai/MEXP


text:
#ifndef LOGGER_INCLUDE
#define LOGGER_INCLUDE

#include <gsl/string_span>

#include "log/log.h"
BOOST_LOG_GLOBAL_LOGGER(exec_helper_test_logger, execHelper::log::LoggerType);

static const gsl::czstring<> LOG_CHANNEL = "test";
#define LOG(x)                                                                 \
    BOOST_LOG_STREAM_CHANNEL_SEV(exec_helper_test_logger::get(), LOG_CHANNEL,  \
                                 execHelper::log::x)                           \
        << boost::log::add_value(fileLog, __FILE__)                            \
        << boost::log::add_value(lineLog, __LINE__)

#endif /* LOGGER_INCLUDE */

max_stars_repo_name: exec-helper/source

max_issues_repo_name: exec-helper/source

max_forks_repo_name: exec-helper/source


text:
#ifndef BASIS_PURSUIT_A_H
#define BASIS_PURSUIT_A_H

#include "Matrix.h"
#include "APPROX2.h"
#include "ALM_I_APPROX.h"
#include <string>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <stdio.h>      /* printf */
#include <time.h>
#include <fstream>
#include <algorithm>
#include <iomanip>
#include <ctime>
#include <math.h>

/*This class solves problem of the form f(x)+g(x) under the constraint Mx=c by IPALM_APPROX;
where f(x)=0 and g(x)=frac{lambda2}{2}|x|_2+lambda1 |x|_1.
*/

template<typename L, typename D>
class Basis_pursuit_a: public ALM_I_APPROX<L, D>
{
private:

  D lambda1;

  D lambda2;

  Matrix<L,D> my_M;

  Matrix<L,D> my_A;

  D val_lambda_f;

protected:

public:

  Basis_pursuit_a(const char* Matrix_file,D val_lambda1, D val_lambda2)
  :ALM_I_APPROX<L,D>(), my_M(Matrix_file)
  {
    lambda1=val_lambda1;
    lambda2=val_lambda2;
    val_lambda_f=0;
    this->val_mu_f=0;
    this->val_mu_g=lambda2;
  }

  L get_n(){return my_M.nfeatures;}
  L get_m(){return my_M.nsamples;}

  inline D gradient_of_phi_j(D x1, L i){
    return 0;
  }

  inline D value_of_g_i(D x, L i){
    return lambda2*x*x/2+lambda1*fabs(x);
  }

  inline D value_of_phi_j(D x1, L i){return 0;}

  inline D prox_of_g_i(D x1,D x2,D x3, L i){
    return compute_one_step(1.0/x2, x1, x3)-x3;
  }

  D compute_one_step(D tau, D u, D x){
    D new_x;
    if(x>tau*(lambda1+u))
    new_x=(x-tau*(lambda1+u))/(1+lambda2*tau);
    else if(x<tau*(u-lambda1))
    new_x=(x-tau*(u-lambda1))/(1+lambda2*tau);
    else
    new_x=0;
    return new_x;
  }

  inline D feasible_dual(vector<D> & x, vector<D> & y){

    if(lambda2+this->beta_s>0)
    {
      return 1;
    }
    else
    {
      D scal=1;
      L l=x.size();
      for(L i=0;i<l;i++)
      if(fabs(x[i]+y[i])>lambda1)
      scal=min(lambda1/fabs(x[i]+y[i]),scal);
      return scal;
    }
  }

  inline D rescale(vector<D> & x, vector<D> & y){
    D scal=1;
    L l=x.size();
    for(L i=0;i<l;i++){
      if(fabs(x[i]+y[i])>lambda1){
        scal=min(lambda1/fabs(x[i]+y[i]),scal);
      }
    }
    return scal;
  }

  inline D value_of_phistar_i(D x,L i) {return 0;}

  inline D value_of_g_tilde_star(D scal,vector<D> & x, vector<D> & y){

    if(lambda2+this->tau_s*this->beta_s>0)
    {
      D res=0;
      D ip=0;
      L l=x.size();
      for(L i=0;i<l;i++)
      {
        ip=max(fabs(x[i]+y[i]+this->tau_s*this->beta_s*this->x_s[i])-lambda1,0.0);
        res+=ip*ip;
      }
      D xsnorm=0;
      for(L i=0;i<l;i++)
      xsnorm+=this->tau_s*this->beta_s/2*this->x_s[i]*this->x_s[i];
      return res/2/(lambda2+this->tau_s*this->beta_s)-xsnorm;
    }
    else
    {
      return 0;
    }
  }

  inline D value_of_g_star(D scal,vector<D> & x, vector<D> & y){
    return 0;
  }

  inline void set_matrix_M(){
    this->data_M=my_M;
  }

  inline void set_matrix_A(){
    this->data_A.nsamples=0;
    this->data_A.nfeatures=this->data_M.nfeatures;
    this->data_A.nnz=0;
    this->data_A.ptr.resize(1,0);
    this->data_A.ptr_t.resize(this->data_M.nfeatures+1,0);
  }

  inline D distance_to_subgradient_of_g(){
    D res=0;
    D tmp;
    D xv;
    for(L i=0;i<this->n;i++){
      tmp=-this->gradient_of_f[i]-this->M_tlambda_s[i];
      xv=this->x_s[i];
      if(xv>0) res+=(tmp-lambda2*xv-lambda1)*(tmp-lambda2*xv-lambda1);
      else if(xv<0) res+=(tmp-lambda2*xv+lambda1)*(tmp-lambda2*xv+lambda1);
      else if(tmp-lambda2*xv>lambda1) res+=(tmp-lambda2*xv-lambda1)*(tmp-lambda2*xv-lambda1);
      else if(tmp-lambda2*xv<-lambda1) res+=(tmp-lambda2*xv+lambda1)*(tmp-lambda2*xv+lambda1);
      else res+=0;
    }
    return sqrt(res);
  }

  inline D subgradient_of_Tx(){
    D res=0;
    D tmp;
    D xv;
    for(L i=0;i<this->n;i++){
      xv=this->Tx[i];
      tmp=this->gradient_of_Tx[i]+ lambda2*xv+this->tau_s*this->beta_s*(xv- this->x_s[i]);
      if(xv>0) res+=(tmp+ lambda1)*(tmp+ lambda1);
      else if(xv<0) res+=(tmp-lambda1)*(tmp-lambda1);
      else if(tmp>lambda1) res+=(tmp-lambda1)*(tmp-lambda1);
      else if(tmp<-lambda1) res+=(tmp+lambda1)*(tmp+lambda1);
      else res+=0;
    }
    return sqrt(res);
  }

  void ALM_I_APPROX_solver(D beta_0, D epsilon_0,  D eta, D rho,vector<D> & x0,vector<D> & y0,L val_tau, L max_nb_outer, L p_N_1, L p_N_2,string filename1, string filename2, D time){
    this->ALM_I_APPROX_solve_with_APPROX(beta_0, epsilon_0,  eta, rho,x0,y0, val_tau, max_nb_outer,  p_N_1,  p_N_2,  val_lambda_f, filename1,  filename2, time);
  }

};

#endif


text:
/*
  Copyright (c) 2016-2017 Hong Xu

  This file is part of WCSPLift.

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

  WCSPLift 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 WCSPLift.  If not, see <http://www.gnu.org/licenses/>.
*/

/** \file WCSPInstance.h
 *
 * Define the definition of WCSP instances and the constraints in it. Algorithms that can be applied
 * directly on WCSP instances, i.e., the min-sum message passing algorithm \cite xkk17 and integer
 * linear programming \cite xkk17a, are also included in this file.
 */

#ifndef WCSPINSTANCE_H_
#define WCSPINSTANCE_H_

#include <array>
#include <cstdint>
#include <istream>
#include <map>
#include <numeric>
#include <set>
#include <stdexcept>
#include <utility>
#include <vector>

#include <boost/dynamic_bitset.hpp>
#include <boost/property_tree/json_parser.hpp>
#include <boost/property_tree/ptree.hpp>
#include <cblas.h>

#include "global.h"
#include "RunningTime.h"
#include "LinearProgramSolver.h"

/** \brief Weighted constraint.
 *
 * This class describes a weighted constraint. It consists of the IDs of the variables as well as
 * how the weight corresponding to each assignment.
 *
 * \tparam VarIdType The type of variable IDs. It must be an integer type and defaults to \p
 * intmax_t.
 *
 * \tparam WeightType The type of weights. It must be a numeric type and defaults to \p double.
 *
 * \tparam ValueType The type of values of all variables in the constraint. It should usually be a
 * bitset type and defaults to \p boost::dynamic_bitset<>.
 */
template < class VarIdType = intmax_t,
           class WeightType = double,
           class ValueType = boost::dynamic_bitset<>,  // The type of the values to represent each constraint
           class NonBooleanValueType = size_t>
class Constraint
{
public:
    /** Alias of \p VarIDType. */
    typedef VarIdType variable_id_t;

    /** Alias of \p ValueType. */
    typedef ValueType value_t;

    /** Alias of \p WeightType. */
    typedef WeightType weight_t;

    /** Alias of \p NonBooleanValueType. */
    typedef NonBooleanValueType non_boolean_value_t;

private:
    // This class is only used for Polynomial key comparison. This function ensures keys with
    // a smaller number of variables must precede larger number of variables.
    class PolynomialKeyComparison
    {
    public:
        bool operator () (const std::set<variable_id_t>& a, const std::set<variable_id_t>& b) const
        {
            if (a.size() > b.size())
                return true;

            if (a.size() < b.size())
                return false;

            return a < b;
        }
    };
public:
    /** The polynomial form of constraints consisting of the coefficient of each term (each term is
     * represented by an assignment of values to variables). */
    typedef std::map<std::set<variable_id_t>, weight_t, PolynomialKeyComparison> Polynomial;

private:
    std::vector<variable_id_t> nonBooleanVariables;
    std::vector<variable_id_t> variables;
    std::map<value_t, weight_t> weights;
    std::map<std::vector<non_boolean_value_t>, weight_t> nonBooleanWeights;

public:
    /** Alias of \p Constraint<VarIdType, WeightType, ValueType, NonBooleanValueType>, i.e., the
     * class type itself. */
    typedef Constraint<VarIdType, WeightType, ValueType, NonBooleanValueType> self_t;

    /** \brief Get the list of the IDs of non-Boolean variables in the constraint.
     *
     * \return A list of variable IDs.
     */
    inline const std::vector<variable_id_t>& getNonBooleanVariables() const noexcept
    {
        return nonBooleanVariables;
    }

    /** \brief Get the list of the IDs of variables in the constraint.
     *
     * \return A list of variable IDs.
     */
    inline const std::vector<variable_id_t>& getVariables() const noexcept
    {
        return variables;
    }

    /** \brief Get the ID of the <tt>i</tt>'th variable.
     *
     * \param[in] i The index of the variable of which to get the ID.
     *
     * \return The ID of the <tt>i</tt>'th variable.
     */
    inline variable_id_t getVariable(size_t i) const noexcept
    {
        return variables.at(i);
    }

    /** \brief Set the IDs of the variables in the constraint from a \p std::vector<variable_id_t>
     * object.
     *
     * \param[in] variables: A list of variable IDs to set.
     */
    inline void setVariables(const std::vector<variable_id_t>& variables) noexcept
    {
        this->variables = variables;
        this->weights.clear();
    }

    /** \brief Set the IDs of the non-Boolean variables in the constraint from a \p
     * std::vector<variable_id_t> object.
     *
     * \param[in] variables: A list of variable IDs to set.
     */
    inline void setNonBooleanVariables(const std::vector<variable_id_t>& variables) noexcept
    {
        this->nonBooleanVariables = variables;
        this->weights.clear();
    }

    /** \brief Set the IDs of the variables in the constraint from a range of iterators.
     *
     * \param[in] b the beginning iterator.
     *
     * \param[in] e the ending iterator.
     */
    template<class Iter>
    inline void setVariables(Iter b, Iter e) noexcept
    {
        variables.clear();
        std::copy(b, e, std::back_inserter(variables));
    }

    /** \brief Set the weight of a given assignment of values to variables specified by \p v to \p
     * w.
     *
     * \param[in] v An assignment of values to variables.
     *
     * \param[in] w A weight.
     *
     * \return A reference to the class object itself.
     */
    inline self_t& setWeight(const value_t& v, weight_t w)
    {
        weights[v] = w;
        return *this;
    }

    /** \brief Set the weight of a given assignment of values to non-Boolean variables specified by
     * \p v to \p w.
     *
     * \param[in] v An assignment of values to variables.
     *
     * \param[in] w A weight.
     *
     * \return A reference to the class object itself.
     */
    inline self_t& setWeight(const std::vector<non_boolean_value_t>& v, weight_t w)
    {
        nonBooleanWeights[v] = w;
        return *this;
    }

    /** \brief Get the weight of a given assignment of values to variables specified by \p v.
     *
     * \param[in] v The assignment of values to variables.
     *
     * \return The weight of a given assignment of values to variables specified by \p v.
     */
    inline weight_t getWeight(const value_t& v) const noexcept
    {
        auto it = weights.find(v);
        if (it == weights.end()) // not found, return 0
            return weight_t();

        return it->second;
    }

    /** \brief Get the weight of a given assignment of values to non-Boolean variables specified by
     * \p v.
     *
     * \param[in] v The assignment of values to variables.
     *
     * \return The weight of a given assignment of values to variables specified by \p v.
     */
    inline weight_t getWeight(const std::vector<non_boolean_value_t>& v) const noexcept
    {
        auto it = nonBooleanWeights.find(v);
        if (it == nonBooleanWeights.end()) // not found, return 0
            return weight_t();

        return it->second;
    }

    /** \brief Get a map object that maps assignments of values to variables to weights.
     *
     * \return A map object that maps assignments of values to variables to weights.
     */
    inline const std::map<value_t, weight_t> getWeights() const noexcept
    {
        return weights;
    }

    /** \brief Represent the constraint using a \p boost::property_tree::ptree object.
     *
     * \param[out] t A \p boost::property_tree::ptree object that represents the constraint.
     */
    void toPropertyTree(boost::property_tree::ptree& t) const
    {
        using namespace boost::property_tree;

        t.clear();

        // put the variables
        ptree vs;
        for (auto& v : variables)
            vs.push_back(
                std::make_pair("", ptree(std::to_string(v))));

        t.put_child("variables", vs);

        // put the weights
        ptree ws;
        for (auto& w : weights)
        {
            ptree we;

            for (size_t i = 0; i < getVariables().size(); ++ i)
            {
                if (w.first[i])
                    we.push_back(std::make_pair("", ptree("1")));
                else
                    we.push_back(std::make_pair("", ptree("0")));
            }

            we.push_back(std::make_pair("", ptree(std::to_string(w.second))));

            ws.push_back(std::make_pair("", we));
        }

        t.put_child("weights", ws);
    }

    /** \brief Compute the coefficients of the polynomial converted from a constraint according to
     * \cite k08.
     *
     * \param[out] p A Constraint::Polynomial object corresponding to the polynomial representation
     * of the constraint.
     */
    void toPolynomial(Polynomial& p) const noexcept
    {
        size_t s = getVariables().size();

        // we basically solve A * x = b

        size_t max_int = ((size_t) 1) << s;

        // matrix a
        double * a = new double[max_int * max_int];
        memset(a, 0, max_int * max_int * sizeof(double));

/// \cond
#define ENTRY(a, x, y) a[(y) * max_int + (x)]
/// \endcond
        // set the first column to 1
        for (size_t i = 0; i < max_int; ++ i)
            ENTRY(a, i, 0) = 1.0;
        // Set the lower triangle. The row number corresponds to the assignments of variables, and
        // the col number corresponds to terms of the polynomial in the following order:
        // X_1, X_2, X_1 X_2, X_3, X_1 X_3, X_2 X_3, X_1 X_2 X_3...
        for (size_t i = 1; i < max_int; ++ i)
            for (size_t j = 1; j <= i; ++ j)
                // (i|j)==i : the 1-bits integer i includes all the 1-bits of j
                ENTRY(a, i, j) = ((i | j) == i) ? 1.0 : 0.0;
#undef ENTRY

        // assign the vector b
        double * x = new double[max_int];
        for (size_t i = 0; i < max_int; ++ i)
            x[i] = getWeight(value_t(s, i));

        cblas_dtrsv(CblasColMajor, CblasLower, CblasNoTrans, CblasUnit, max_int, a, max_int, x, 1);
        delete[] a;

        // Insert all the coefficients into the polynomial.
        for (size_t i = 0; i < max_int; ++ i)
        {
            typename Polynomial::key_type k;

            for (size_t j = 0; j < s; ++ j)
                if ((((size_t) 1u) << j) & i)
                    k.insert(getVariable(j));

            p[k] += x[i];
        }

        delete[] x;
    }
};

/** \brief A WCSP instance.
 *
 * This class describes a WCSP instance. It consists a set of constraints.
 *
 * \tparam VarIdType The type of variable IDs. It must be an integer type and defaults to \p
 * intmax_t.
 *
 * \tparam WeightType The type of weights. It must be a numeric type and defaults to \p double.
 *
 * \tparam ConstraintValueType The type of values of all variables in the constraints in the WCSP
 * instance. It should usually be a bitset type and defaults to \p boost::dynamic_bitset<>.
 */
template <class VarIdType = intmax_t, class WeightType = double,
          class ConstraintValueType = boost::dynamic_bitset<>,
          class NonBooleanValueType = size_t>
class WCSPInstance
{
public:
    /** Alias of \p VarIdType. */
    typedef VarIdType variable_id_t;
    /** Alias of \p WeightType. */
    typedef WeightType weight_t;
    /** Alias of \p ConstraintValueType. */
    typedef ConstraintValueType constraint_value_t;
    /** Alias of \p NonBooleanValueType. */
    typedef NonBooleanValueType non_boolean_value_t;

public:
    /** Alias of the Constraint type in the WCSP instance, i.e., \p
     * Constraint<variable_id_t,weight_t,constraint_value_t>.
     */
    typedef Constraint<variable_id_t, weight_t, constraint_value_t> constraint_t;

private:
    std::vector<constraint_t> constraints;

    // map from non-Boolean variables to their Boolean variable representations
    std::vector<std::vector<variable_id_t> > nonBooleanVariables;

public:

    /** \brief Get the mapping from original variables to Boolean variables.
     *
     * \return The mapping.
     */
    inline const std::vector<std::vector<variable_id_t> >& getNonBooleanVariables()
        const noexcept
    {
        return nonBooleanVariables;
    }

    /** \brief Display the mapping from original variables to Boolean variables.
     */
    inline void displayNonBooleanVariableMapping() const noexcept
    {
        std:: cout << "--- Non-Boolean Variable Mapping BEGINS ---" << std::endl;
        for (auto i = 0; i < nonBooleanVariables.size(); ++ i)
        {
            std::cout << i << '\t';
            for (auto j : nonBooleanVariables[i])
            {
                std::cout << j << ' ';
            }
            std::cout << std::endl;
        }
        std:: cout << "--- Non-Boolean Variable Mapping ENDS ---" << std::endl;
    }

    /** \brief Get the list of constraints.
     *
     * \return A list of \p constraint_t objects.
     */
    inline const std::vector<constraint_t>& getConstraints() const noexcept
    {
        return constraints;
    }

    /** \brief Get the <tt>i</tt>'th constraint.
     *
     * \param[in] i The index of the Constraint object to get.
     *
     * \return The <tt>i</tt>'th constraint.
     */
    inline const constraint_t& getConstraint(size_t i) const noexcept
    {
        return constraints.at(i);
    }

    /** \brief Compute the total weight corresponding to an assignment of values to all variables.
     *
     * \param[in] assignments The assignment of values to all variables.
     *
     * \return The computed total weight.
     */
    inline weight_t computeTotalWeight(
        const std::map<variable_id_t, non_boolean_value_t>& assignments) const noexcept
    {
        weight_t tw = 0;

        for (const auto& c : constraints)
        {
            std::vector<non_boolean_value_t> val(c.getNonBooleanVariables().size());
            for (size_t i = 0; i < c.getNonBooleanVariables().size(); ++ i)
            {
                try
                {
                    val[i] = assignments.at(c.getNonBooleanVariables().at(i));
                } catch (std::out_of_range& e)
                {
                    // Here, it means the variable assignments does not hold
                    // c.getNonBooleanVariables().at(i) -- the CCG does not contain that variable.
                    // This may well be that the variable itself is a "dummy" variable -- the value
                    // of it does not affect the total weight. In this case, we always assign zero
                    // (or anything else) to it.
                    val[i] = 0;
                }
            }

            tw += c.getWeight(val);
        }

        return tw;
    }

public:
    /** \brief Load a problem instance from a stream.
     *
     * \param[in] f The input stream.
     */
    void load(std::istream& f);

    /** \brief Load the problem in DIMACS format. The format specification is at
     * http://graphmod.ics.uci.edu/group/WCSP_file_format
     *
     * \param[in] f The input stream.
     */
    void loadDimacs(std::istream& f);

    /** \brief Load the problem in UAI format. The format specification is at
     * http://www.hlt.utdallas.edu/~vgogate/uai14-competition/modelformat.html
     *
     * \param[in] f The input stream.
     */
    void loadUAI(std::istream& f);

    /** List of supported input file formats.
     */
    enum class Format
    {
        DIMACS,
        UAI
    };

    /** \brief Construct a WCSP instance from an input stream in a given format.
     *
     * \param[in] f The input stream.
     *
     * \param[in] format The format of \p f.
     */
    WCSPInstance(std::istream& f, Format format)
    {
        switch (format)
        {
        case Format::DIMACS:
            loadDimacs(f);
            break;
        case Format::UAI:
            loadUAI(f);
            break;
        }
    }

    /** \brief Convert this WCSP instance to a human-readable string.
     *
     * \return The human-readable string.
     */
    std::string toString() const noexcept
    {
        std::stringstream ss;

        ss << '{' << std::endl;
        for (auto& c : constraints)
        {
            ss << c;
        }

        return ss.str();
    }

    /** \brief Solve the WCSP instance using the min-sum message passing algorithm.
     *
     * \param[in] delta The threshold for convergence determination. That is, the algorithm
     * terminates iff all messages change within \p delta compared with the previous iteration.
     *
     * \return A solution to the WCSP instance.
     */
    std::map<variable_id_t, bool> solveUsingMessagePassing(weight_t delta) const noexcept
    {
        std::map<variable_id_t, std::set<const constraint_t*> > constraint_list_for_v;

        // Initialize all messages.
        std::map<std::pair<variable_id_t, const constraint_t*>, std::array<weight_t, 2> > msgs_v_to_c;
        std::map<std::pair<const constraint_t*, variable_id_t>, std::array<weight_t, 2> > msgs_c_to_v;
        for (const auto& c : constraints)
        {
            for (auto v : c.getVariables())
            {
                constraint_list_for_v[v].insert(&c);
                msgs_v_to_c[std::make_pair(v, &c)] = {0, 0};
                msgs_c_to_v[std::make_pair(&c, v)] = {0, 0};
            }
        }

        bool converged = false;
        uintmax_t num_iterations = 0;
        do
        {
            if (RunningTime::GetInstance().isTimeOut())  // time's up
                break;

            ++ num_iterations;

            converged = true;

            auto msgs_v_to_c0 = msgs_v_to_c;
            auto msgs_c_to_v0 = msgs_c_to_v;

            // Update all the messages from variables to constraints.
            for (auto& it : msgs_v_to_c)
            {
                auto v = it.first.first;
                auto c = it.first.second;

                std::array<weight_t, 2> m = {0, 0};

                for (auto nc : constraint_list_for_v.at(v))
                {
                    if (nc != c)
                    {
                        const auto m_to_v = msgs_c_to_v0.find(std::make_pair(nc, v))->second;

                        m[0] += m_to_v.at(0);
                        m[1] += m_to_v.at(1);
                    }
                }

                it.second = m;

                // is it now convergent?
                if (converged)
                {
                    auto old_m = msgs_v_to_c0.at(it.first);

                    if (std::abs(old_m.at(0) - m.at(0)) > delta ||
                        std::abs(old_m.at(1) - m.at(1)) > delta)
                        converged = false;
                }
            }

            // Update all the messages from constraints to variables.
            for (auto& it : msgs_c_to_v)
            {
                auto c = it.first.first;
                auto v = it.first.second;

                std::vector<std::array<weight_t, 2> > m_to_c;
                m_to_c.reserve(c->getVariables().size());
                std::vector<variable_id_t> v_ids;
                v_ids.reserve(c->getVariables().size());
                size_t self_index;   // index of the variable itself
                for (size_t i = 0; i < c->getVariables().size(); ++ i)
                {
                    auto nv = c->getVariable(i);

                    m_to_c.push_back(msgs_v_to_c0.at(std::make_pair(nv, c)));
                    v_ids.push_back(nv);

                    if (nv == v)
                        self_index = i;
                }

                std::array<weight_t, 2> m = {
                    std::numeric_limits<weight_t>::max(), std::numeric_limits<weight_t>::max()
                };

                // TODO: Change the iteration which does not have an upper limit of number of bits.
                for (unsigned long i = 0; i < (1ul << c->getVariables().size()); ++ i)
                {
                    constraint_value_t values(c->getVariables().size(), i);

                    weight_t sum = c->getWeight(values);

                    for (size_t j = 0; j < c->getVariables().size(); ++ j)
                        if (c->getVariable(j) != v)
                            sum += msgs_v_to_c0.at(
                                std::make_pair(c->getVariable(j), c))[values[j] ? 1 : 0];

                    size_t index_update = values[self_index] ? 1 : 0;

                    if (m[index_update] > sum)
                        m[index_update] = sum;
                }

                auto m_min = m[0] < m[1] ? m[0] : m[1];
                m[0] -= m_min;
                m[1] -= m_min;

                it.second = m;

                // is it now convergent?
                if (converged)
                {
                    auto old_m = msgs_c_to_v0.at(it.first);

                    if (std::abs(old_m.at(0) - m.at(0)) > delta ||
                        std::abs(old_m.at(1) - m.at(1)) > delta)
                        converged = false;
                }
            }
        } while (!converged);

        std::cout << "Number of iterations: " << num_iterations << std::endl;

        // compute the variable values
        std::map<variable_id_t, std::array<weight_t, 2> > assignments_w;

        converged = true;
        for (const auto& c : constraints)
            for (auto v : c.getVariables())
            {
                assignments_w[v][0] += msgs_c_to_v[std::make_pair(&c, v)][0];
                assignments_w[v][1] += msgs_c_to_v[std::make_pair(&c, v)][1];
                if (isinf(assignments_w[v][0]) || isinf(assignments_w[v][1]))
                    converged = false;
            }

        std::map<variable_id_t, bool> assignments;
        for (const auto& a : assignments_w)
            assignments[a.first] = a.second[0] > a.second[1] ? a.second[1] : a.second[0];

        if (!converged)
            std::cout << "*** Message passing not converged! ***" << std::endl;

        return assignments;
    }

    /** \brief Solve the WCSP instance using linear programming.
     *
     * \param[in] lps The linear program solver to be used.
     *
     * \return A solution to the WCSP instance.
     */
    std::map<variable_id_t, non_boolean_value_t> solveUsingLinearProgramming(
        LinearProgramSolver& lps) const noexcept
    {
        lps.reset();
        lps.setTimeLimit(RunningTime::GetInstance().getTimeLimit().count());

        lps.setObjectiveType(LinearProgramSolver::ObjectiveType::MIN);

        auto constraints = this->getConstraints();

        // Every variable needs a unary constraint. This is the indices of all unary constraints.
        std::vector<size_t> unary_indices(this->getNonBooleanVariables().size(),
                                                 std::numeric_limits<size_t>::max());
        size_t num_unary = 0;
        for (size_t i = 0; i < constraints.size(); ++ i)
        {
            auto& vs = constraints.at(i).getVariables();
            if (vs.size() == 1)  // unary constraint
            {
                unary_indices[vs.at(0)] = i;
                ++ num_unary;
            }
        }
        constraints.reserve(constraints.size() + unary_indices.size() - num_unary);
        for (variable_id_t i = 0; i < unary_indices.size(); ++ i)
        {
            // no unary constraint for this var
            if (unary_indices.at(i) == std::numeric_limits<size_t>::max())
            {
                constraint_t cons;
                cons.setNonBooleanVariables(std::vector<variable_id_t>{i});
                constraints.push_back(std::move(cons));
                unary_indices[i] = constraints.size() - 1;
            }
        }

        std::vector<std::vector<LinearProgramSolver::variable_id_t> > lp_variables(
            constraints.size());

        // add LP variables and constraints
        for (size_t i = 0; i < constraints.size(); ++ i)
        {
            size_t num_values = 1;
            // Cartesian product of all values
            for (const auto& nbv : constraints.at(i).getNonBooleanVariables())
                num_values *= nonBooleanVariables[nbv].size() + 1;
            lp_variables[i].reserve(num_values);
            for (size_t j = 0; j < num_values; ++ j)  // add all variables
            {
                std::vector<non_boolean_value_t> val;
                val.reserve(constraints.at(i).getNonBooleanVariables().size());

                size_t j0 = j;
                for (const auto& nbv : constraints.at(i).getNonBooleanVariables())
                {
                    size_t d = nonBooleanVariables[nbv].size() + 1;
                    val.push_back(j0 % d);
                    j0 /= d;
                }

                lp_variables[i].push_back(lps.addVariable(constraints.at(i).getWeight(val)));
            }

            // Only one value in a WCSP constraint can take effect.
            lps.addConstraint(lp_variables[i], std::vector<double>(num_values, 1.0), 1.0,
                              LinearProgramSolver::ConstraintType::EQUAL);
        }

        // Add constraints that make sure LP variables that represent overlapped WCSP variables are
        // consistent.
        for (size_t i = 0; i < lp_variables.size(); ++ i)
        {
            for (size_t j : unary_indices)
            {
                if (j <= i)  // Don't do a pair twice.
                    continue;

                // Get the overlap map.
                std::vector<std::pair<size_t, size_t> > overlapped_variables;

                // Whether a variable is overlapping?
                std::vector<bool> vs_overlap_p[2];
                vs_overlap_p[0].resize(constraints.at(i).getNonBooleanVariables().size());
                vs_overlap_p[1].resize(constraints.at(j).getNonBooleanVariables().size());
                std::array<std::vector<size_t>, 2> domain_sizes;
                domain_sizes[0].reserve(constraints.at(i).getNonBooleanVariables().size());
                domain_sizes[1].reserve(constraints.at(j).getNonBooleanVariables().size());
                std::vector<size_t> overlapped_domain_sizes;

                for (size_t z = 0; z < 2; ++ z)
                {
                    size_t c = z ? j : i;
                    for (size_t k = 0; k < constraints.at(c).getNonBooleanVariables().size();
                         ++ k)
                        domain_sizes[z].push_back(
                            nonBooleanVariables[
                                constraints.at(c).getNonBooleanVariables()[k]].size() + 1);
                }

                for (size_t k0 = 0; k0 < constraints.at(i).getNonBooleanVariables().size(); ++ k0)
                    for (size_t k1 = 0; k1 < constraints.at(j).getNonBooleanVariables().size(); ++ k1)
                    {
                        if (constraints.at(i).getNonBooleanVariables()[k0] ==
                            constraints.at(j).getNonBooleanVariables()[k1])
                        {
                            overlapped_variables.push_back(std::make_pair(k0, k1));
                            overlapped_domain_sizes.push_back(
                                nonBooleanVariables[
                                    constraints.at(i).getNonBooleanVariables()[k0]].size() + 1
                                );
                            vs_overlap_p[0][k0] = true;
                            vs_overlap_p[1][k1] = true;
                            break;
                        }
                    }

                std::vector<non_boolean_value_t> vs[2];
                vs[0].resize(constraints.at(i).getNonBooleanVariables().size());
                vs[1].resize(constraints.at(j).getNonBooleanVariables().size());

                if (overlapped_domain_sizes.empty())  // no overlapping
                    continue;
                // Cartesian product of all values of overlapped variables
                size_t num_values = std::accumulate(overlapped_domain_sizes.begin(),
                                                    overlapped_domain_sizes.end(),
                                                    1, std::multiplies<size_t>());

                for (size_t k = 0; k < num_values; ++ k)
                {
                    size_t k0 = k;
                    for (const auto& ov : overlapped_variables)
                    {
                        size_t d = nonBooleanVariables[
                            constraints.at(i).getNonBooleanVariables()[ov.first]].size() + 1;
                        non_boolean_value_t val = k0 % d;
                        k0 /= d;
                        vs[0][ov.first] = val;
                        vs[1][ov.second] = val;
                    }

                    // non-overlapped variables
                    std::array<std::vector<non_boolean_value_t>, 2> vs_non_overlapped;
                    vs_non_overlapped[0].resize((vs[0].size() - overlapped_variables.size()));
                    vs_non_overlapped[1].resize((vs[1].size() - overlapped_variables.size()));

                    std::vector<LinearProgramSolver::variable_id_t> lp_vs;  // LP variables
                    std::vector<double> coefs;

                    for (size_t z = 0; z < 2; ++ z)
                    {
                        size_t l;
                        do {
                            for (l = 0; l < vs[z].size(); ++ l)
                            {
                                if (vs_overlap_p[z][l])
                                    continue;

                                ++ vs[z][l];
                                if (vs[z][l] >= domain_sizes[z][l])
                                    vs[z][l] = 0;
                                else
                                    break;
                            }

                            size_t index = 0;  // index corresponding to the variable values vs[z]
                            size_t cur_base = 1;
                            for (size_t m = 0; m < vs[z].size(); ++ m)
                            {
                                index += cur_base * vs[z][m];
                                cur_base *= domain_sizes[z][m];
                            }

                            if (z == 0)
                            {
                                lp_vs.push_back(lp_variables[i][index]);
                                coefs.push_back(1.0);
                            }
                            else
                            {
                                lp_vs.push_back(lp_variables[j][index]);
                                coefs.push_back(-1.0);
                            }
                        } while (l != vs[z].size());
                    }

                    lps.addConstraint(lp_vs, coefs, 0.0,
                                      LinearProgramSolver::ConstraintType::EQUAL);
                }
            }
        }

        std::vector<double> assignments;
        lps.solve(assignments);

        // Compute the solutions from assignments
        std::map<variable_id_t, non_boolean_value_t> solution;
        for (size_t i : unary_indices)
            for (size_t j = 0; j < lp_variables[i].size(); ++ j)
            {
                if (assignments[lp_variables[i][j]] > 0.99)
                {
                    solution[constraints.at(i).getNonBooleanVariables().at(0)] = j;
                    break;
                }
            }

        return solution;
    }
};

template <class VarIdType, class WeightType, class ConstraintValueType, class NonBooleanValueType>
void WCSPInstance<VarIdType, WeightType, ConstraintValueType, NonBooleanValueType>::loadDimacs(
    std::istream& f)
{
    std::string line;
    std::string tmp;
    size_t max_domain_size;
    size_t nv, nc; // number of variables, number of constraints

    // first line: problem_name number_of_variables max_domain_size number_of_constraints
    // global_upper_bound
    std::getline(f, line);
    std::istringstream ss(line);
    ss >> tmp; // name
    ss >> nv;
    ss >> max_domain_size; // max domain size
    ss >> nc;
    ss >> tmp; // ignore global upper bound

    // domain size of all variables
    std::getline(f, line);
    ss.clear();
    ss.str(line);
    variable_id_t cur_var_id = 0;
    nonBooleanVariables.clear();
    nonBooleanVariables.reserve(nv);
    for (size_t i = 0; i < nv; ++ i)
    {
        size_t domain_size;
        ss >> domain_size;
        // Boolean variables corresponding to the current variable
        std::vector<variable_id_t> vs(domain_size-1);
        std::iota(vs.begin(), vs.end(), cur_var_id);
        cur_var_id += vs.size();
        nonBooleanVariables.push_back(std::move(vs));
    }

    constraints.clear();
    constraints.resize(nc);

    // iterate over constraints
    for (size_t i = 0; i < nc; ++ i)
    {
        size_t arity;
        weight_t default_cost;
        size_t ntuples; // number of tuples not having the default cost
        std::getline(f, line);
        ss.clear();
        ss.str(line);

        ss >> arity;

        // load the variables in the constraint
        std::vector<variable_id_t> non_bool_variables;
        non_bool_variables.reserve(arity);
        // load the corresponding Boolean variables in the constraint
        std::vector<variable_id_t> variables;
        variables.reserve(arity*max_domain_size);
        for (size_t j = 0; j < arity; ++ j)
        {
            variable_id_t vid;
            ss >> vid;
            non_bool_variables.push_back(vid);
            const auto& nbv = nonBooleanVariables[vid];
            variables.insert(variables.end(), nbv.begin(), nbv.end());
        }
        constraints[i].setVariables(std::move(variables));
        constraints[i].setNonBooleanVariables(std::move(non_bool_variables));
        ss >> default_cost;
        // TODO: More efficient handling
        if (default_cost > std::abs(1e-6))   // non-zero default cost
        {
            constraint_value_t values;
            values.resize(constraints[i].getVariables().size());
            unsigned long te = (1u << constraints[i].getVariables().size());
            for (unsigned long l = 0; l < te; ++ l)
            {
                for (unsigned long k = 0; k < values.size(); ++ k)
                    values.set(k, (l & (1u << k)) ? 1 : 0);
                constraints[i].setWeight(values, default_cost);
            }
        }
        ss >> ntuples;

        // load the entries in the constraints
        for (size_t j = 0; j < ntuples; ++ j)
        {
            weight_t cost;
            constraint_value_t values;
            values.resize(constraints[i].getVariables().size());
            values.set();
            std::vector<non_boolean_value_t> non_boolean_values(arity);
            std::getline(f, line);
            ss.clear();
            ss.str(line);
            for (size_t k = 0, cur_bit = 0; k < arity; ++ k)
            {
                int val;
                ss >> val;
                non_boolean_values[k] = val;
                size_t bvs = nonBooleanVariables[constraints[i].getNonBooleanVariables()[k]].size();
                if (bvs == 1)
                    values.set(cur_bit ++, val ? true : false);
                else
                {
                    if (val != 0)
                        values.set(cur_bit + val - 1, false);
                    cur_bit += bvs;
                }
            }
            ss >> cost;
            constraints[i].setWeight(std::move(non_boolean_values), cost);
            constraints[i].setWeight(std::move(values), cost);
        }
    }
}

template <class VarIdType, class WeightType, class ConstraintValueType, class NonBooleanValueType>
void WCSPInstance<VarIdType, WeightType, ConstraintValueType, NonBooleanValueType>::loadUAI(
    std::istream& f)
{
    std::string line;
    std::string tmp;
    size_t nv, nc; // number of variables, number of constraints

    // first line: MARKOV, just ignore it
    std::getline(f, tmp);

    // second line: number of variables
    {
        std::getline(f, line);
        std::istringstream ss(line);
        ss >> nv; // number of variables
    }

    // third line: arity
    {
        std::getline(f, line);
        std::istringstream ss(line);
        variable_id_t cur_var_id = 0;
        nonBooleanVariables.clear();
        nonBooleanVariables.reserve(nv);
        for (size_t i = 0; i < nv; ++ i)
        {
            size_t domain_size;
            ss >> domain_size;
            // Boolean varaibles corresponding to the current variable
            std::vector<variable_id_t> vs(domain_size-1);
            std::iota(vs.begin(), vs.end(), cur_var_id);
            cur_var_id += vs.size();
            nonBooleanVariables.push_back(std::move(vs));
        }
    }

    // fourth line: number of constraints
    {
        std::getline(f, line);
        std::istringstream ss(line);
        ss >> nc;
        constraints.clear();
        constraints.resize(nc);
    }

    // iterate over constraints
    for (size_t i = 0; i < nc; ++ i)
    {
        size_t arity;
        std::getline(f, line);
        std::istringstream ss(line);

        ss >> arity;

        // Load the variables in the constraint. We load it reversely because UAI format arrange
        // their constraint value reversely as we do.
        std::vector<variable_id_t> variables;
        std::vector<variable_id_t> non_bool_variables(arity);

        for (size_t j = 0; j < arity; ++ j)
        {
            variable_id_t vid;
            ss >> vid;
            non_bool_variables[arity-j-1] = vid;
            variables.insert(variables.end(),
                             nonBooleanVariables[vid].rbegin(),
                             nonBooleanVariables[vid].rend());
        }
        std::reverse(variables.begin(), variables.end());

        constraints[i].setVariables(std::move(variables));
        constraints[i].setNonBooleanVariables(std::move(non_bool_variables));
    }

    // iterate over constraints
    for (size_t i = 0; i < nc; ++ i)
    {
        size_t ntuples;
        f >> ntuples;  // number of entries

        // load the entries in the constraints
        std::vector<weight_t> costs(ntuples);
        for (size_t j = 0; j < ntuples; ++ j)
            f >> costs[j];

        // normalization constant
        weight_t sum(std::accumulate(costs.begin(), costs.end(), weight_t()));

        // set the value of each tuple
        for (size_t j = 0; j < ntuples; ++ j)
        {
            costs[j] = -std::log(costs[j] / sum);
            if (!std::isfinite(costs[j]))
                costs[j] = 1e6;
            constraint_value_t values;
            std::vector<non_boolean_value_t> non_boolean_values;
            non_boolean_values.reserve(constraints[i].getNonBooleanVariables().size());
            size_t j0 = j;
            for (const auto& nbv : constraints[i].getNonBooleanVariables())
            {
                auto d = nonBooleanVariables[nbv].size() + 1;
                size_t cur_val = j0 % d;
                non_boolean_values.push_back(cur_val);

                if (d == 2)
                    values.push_back(cur_val ? true : false);
                else
                {
                    for (size_t k = 1; k < d; ++ k)
                    {
                        if (cur_val == k)
                            values.push_back(false);
                        else
                            values.push_back(true);
                    }
                }
                j0 /= d;
            }
            constraints[i].setWeight(std::move(values), costs[j]);
            constraints[i].setWeight(std::move(non_boolean_values), costs[j]);
        }
    }
}

/** \brief Write a Constraint object to a stream in a human-readable form.
 *
 * \param[in] o The stream to write to.
 *
 * \param[in] c The constraint object to be write to \p o.
 *
 * \return The stream \p o.
 */
template <class ...T>
std::ostream& operator << (std::ostream& o, const Constraint<T...>& c)
{
    boost::property_tree::ptree t;
    c.toPropertyTree(t);
    boost::property_tree::write_json(o, t, true);
    return o;
}

#endif /* WCSPINSTANCE_H_ */


text:
/*
Performing approximate Bayesian computation sequential Monte Carlo (Toni et al.
2009) on a beta-binomial model.

Synthetic data is generated by `ground_truth_and_analysis.ipynb`.

This script writes particles.csv, where each row corresponds to a particle and
each column corresponds to a round of SMC.

Defining #DEBUG_MODE will silence all writing to stdout. One may then add
printf statements in the code, and perhaps write the output to file as:
`./run.sh > output.txt`

TODO: Export all functions which are specific to the beta-binomial model into
smc.h, and perhaps rename smc.h as beta_binomial.h

Author: Juvid Aryaman
*/

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

#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_sort_double.h>
#include <gsl/gsl_statistics.h>

#define N_DATA 50
#define N_TRUTH 10

#define PRIOR_ALPHA 0.5
#define PRIOR_BETA 0.5

#define N_PARTICLES 5000
#define N_ROUNDS_SMC 50
#define KERNEL_SD 0.05
#define QUANTILE_ACCEPT_DISTANCE 0.8

#define RND gsl_rng_uniform(r)
#define SEED 1
#define DISTANCE_THRESHOLD_INIT 10

#define OUTFILE_NAME "particles.csv"

//#define DEBUG_MODE

#include "smc.h"

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

/* set up GSL RNG */
gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937);
/* end of GSL setup */

gsl_rng_set(r, SEED);

/////////////////////////
/*Read data*/
/////////////////////////

FILE *data_pointer;

data_pointer = fopen("binom_data.csv", "r");

int data[N_DATA];
int i, j, read_error_status;
for (i=0; i < N_DATA; i++){
	read_error_status = fscanf(data_pointer, "%d\n", &data[i]);
}
if (read_error_status != 1){printf("Error reading data\n"); return 0;}

/////////////////////////
/*Initialise variables*/
/////////////////////////

/*Make a (N_ROUNDS_SMC X N_PARTICLES) array to store all particles at all rounds
of SMC*/
double** theta_particle;
theta_particle = (double**) malloc(N_ROUNDS_SMC * sizeof(double*));
for (i = 0; i < N_ROUNDS_SMC; i++){
	theta_particle[i] = (double*) malloc(N_PARTICLES * sizeof(double));
}

double distance_threshold = DISTANCE_THRESHOLD_INIT;
int *simulated_data = (int*) malloc(N_DATA * sizeof(int));
double *distance = malloc(N_PARTICLES * sizeof(double));
double *weight = malloc(N_PARTICLES * sizeof(double));
double weight_normalizer = 0.0;

int time_smc=0; // an index of each round of SMC
int param_index_chosen;

/////////////////////////
/*Perform ABC SMC*/
/////////////////////////

/*For every round of SMC*/
for (time_smc = 0; time_smc < N_ROUNDS_SMC; time_smc++) {
	#ifndef DEBUG_MODE
		printf("Round %d of SMC\n", time_smc);
	#endif

	/*Draw or perturb a particle and compute distance*/
	for (i = 0; i < N_PARTICLES; i++) {
		distance[i] = distance_threshold + 1.0; // reset distance of particle i
		while (distance[i] > distance_threshold) {
			if (time_smc == 0) {
				// Sample from the prior
				theta_particle[0][i] = gsl_ran_beta(r, PRIOR_ALPHA, PRIOR_BETA);
			}
			else{
				/*Sample from the old weights and perturb*/
				param_index_chosen = weighted_choice(r, weight);
				if ((param_index_chosen < 0)||(param_index_chosen >= N_PARTICLES)) {
					printf("Error in param_index_chosen\n"); return -1;
			}
				theta_particle[time_smc][i] =
					theta_particle[time_smc-1][param_index_chosen] +
					gsl_ran_gaussian(r, KERNEL_SD);

				// check if bounds of prior exceeded
				if((theta_particle[time_smc][i]<0) || (theta_particle[time_smc][i] > 1)) continue;
			}

			// Simulate a candidate dataset
			for (j = 0; j < N_DATA; j++) {
				simulated_data[j] = gsl_ran_binomial(r, theta_particle[time_smc][i], N_TRUTH);
			}

			// Compute distance between data and simulation
			distance[i] = distance_metric(data, simulated_data);

		}
	}
	#ifndef DEBUG_MODE
		printf("Particles sampled.\n");
	#endif

	/*Compute weights*/
	if (time_smc==0){ for (i = 0; i < N_PARTICLES; i++) weight[i] = 1.0;}
	else{
		weight_normalizer = 0.0;
		for (i = 0; i < N_PARTICLES; i++) {
			weight_normalizer += weight[i]*kernel_pdf(theta_particle[time_smc-1][i],
				theta_particle[time_smc][i]);
		}

		for (i = 0; i < N_PARTICLES; i++) {
			weight[i] = prior_pdf(theta_particle[time_smc][i])/weight_normalizer;
		}

	}

	/*Normalise weights*/
	weight_normalizer = 0.0;
	for (i = 0; i < N_PARTICLES; i++)	weight_normalizer += weight[i];
	for (i = 0; i < N_PARTICLES; i++)	weight[i] = weight[i]/weight_normalizer;

	/* Resample weights*/
	distance_threshold = update_distance_threshold(distance);

}

#ifndef DEBUG_MODE
	printf("Writing particles to file\n");
#endif
	write_particles_to_csv(theta_particle);
#ifndef DEBUG_MODE
	printf("Done!\n");
#endif

return 0; //return from main
} //close main