text:
/* matrix/gsl_matrix_complex_double.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 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 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.
 */

#ifndef __GSL_MATRIX_COMPLEX_DOUBLE_H__
#define __GSL_MATRIX_COMPLEX_DOUBLE_H__

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_vector_complex_double.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;
  double * data;
  gsl_block_complex * block;
  int owner;
} gsl_matrix_complex ;

typedef struct
{
  gsl_matrix_complex matrix;
} _gsl_matrix_complex_view;

typedef _gsl_matrix_complex_view gsl_matrix_complex_view;

typedef struct
{
  gsl_matrix_complex matrix;
} _gsl_matrix_complex_const_view;

typedef const _gsl_matrix_complex_const_view gsl_matrix_complex_const_view;

/* Allocation */

GSL_EXPORT
gsl_matrix_complex *
gsl_matrix_complex_alloc (const size_t n1, const size_t n2);

GSL_EXPORT
gsl_matrix_complex *
gsl_matrix_complex_calloc (const size_t n1, const size_t n2);

GSL_EXPORT
gsl_matrix_complex *
gsl_matrix_complex_alloc_from_block (gsl_block_complex * b,
                                           const size_t offset,
                                           const size_t n1, const size_t n2, const size_t d2);

GSL_EXPORT
gsl_matrix_complex *
gsl_matrix_complex_alloc_from_matrix (gsl_matrix_complex * b,
                                            const size_t k1, const size_t k2,
                                            const size_t n1, const size_t n2);

GSL_EXPORT
gsl_vector_complex *
gsl_vector_complex_alloc_row_from_matrix (gsl_matrix_complex * m,
                                                const size_t i);

GSL_EXPORT
gsl_vector_complex *
gsl_vector_complex_alloc_col_from_matrix (gsl_matrix_complex * m,
                                                const size_t j);

GSL_EXPORT void gsl_matrix_complex_free (gsl_matrix_complex * m);

/* Views */

GSL_EXPORT
_gsl_matrix_complex_view
gsl_matrix_complex_submatrix (gsl_matrix_complex * m,
                            const size_t i, const size_t j,
                            const size_t n1, const size_t n2);

GSL_EXPORT
_gsl_vector_complex_view
gsl_matrix_complex_row (gsl_matrix_complex * m, const size_t i);

GSL_EXPORT
_gsl_vector_complex_view
gsl_matrix_complex_column (gsl_matrix_complex * m, const size_t j);

GSL_EXPORT
_gsl_vector_complex_view
gsl_matrix_complex_diagonal (gsl_matrix_complex * m);

GSL_EXPORT
_gsl_vector_complex_view
gsl_matrix_complex_subdiagonal (gsl_matrix_complex * m, const size_t k);

GSL_EXPORT
_gsl_vector_complex_view
gsl_matrix_complex_superdiagonal (gsl_matrix_complex * m, const size_t k);

GSL_EXPORT
_gsl_matrix_complex_view
gsl_matrix_complex_view_array (double * base,
                             const size_t n1,
                             const size_t n2);

GSL_EXPORT
_gsl_matrix_complex_view
gsl_matrix_complex_view_array_with_tda (double * base,
                                      const size_t n1,
                                      const size_t n2,
                                      const size_t tda);

GSL_EXPORT
_gsl_matrix_complex_view
gsl_matrix_complex_view_vector (gsl_vector_complex * v,
                              const size_t n1,
                              const size_t n2);

GSL_EXPORT
_gsl_matrix_complex_view
gsl_matrix_complex_view_vector_with_tda (gsl_vector_complex * v,
                                       const size_t n1,
                                       const size_t n2,
                                       const size_t tda);

GSL_EXPORT
_gsl_matrix_complex_const_view
gsl_matrix_complex_const_submatrix (const gsl_matrix_complex * m,
                                  const size_t i, const size_t j,
                                  const size_t n1, const size_t n2);

GSL_EXPORT
_gsl_vector_complex_const_view
gsl_matrix_complex_const_row (const gsl_matrix_complex * m,
                            const size_t i);

GSL_EXPORT
_gsl_vector_complex_const_view
gsl_matrix_complex_const_column (const gsl_matrix_complex * m,
                               const size_t j);

GSL_EXPORT
_gsl_vector_complex_const_view
gsl_matrix_complex_const_diagonal (const gsl_matrix_complex * m);

GSL_EXPORT
_gsl_vector_complex_const_view
gsl_matrix_complex_const_subdiagonal (const gsl_matrix_complex * m,
                                    const size_t k);

GSL_EXPORT
_gsl_vector_complex_const_view
gsl_matrix_complex_const_superdiagonal (const gsl_matrix_complex * m,
                                      const size_t k);

GSL_EXPORT
_gsl_matrix_complex_const_view
gsl_matrix_complex_const_view_array (const double * base,
                                   const size_t n1,
                                   const size_t n2);

GSL_EXPORT
_gsl_matrix_complex_const_view
gsl_matrix_complex_const_view_array_with_tda (const double * base,
                                            const size_t n1,
                                            const size_t n2,
                                            const size_t tda);

GSL_EXPORT
_gsl_matrix_complex_const_view
gsl_matrix_complex_const_view_vector (const gsl_vector_complex * v,
                                    const size_t n1,
                                    const size_t n2);

GSL_EXPORT
_gsl_matrix_complex_const_view
gsl_matrix_complex_const_view_vector_with_tda (const gsl_vector_complex * v,
                                             const size_t n1,
                                             const size_t n2,
                                             const size_t tda);

/* Operations */

GSL_EXPORT gsl_complex gsl_matrix_complex_get(const gsl_matrix_complex * m, const size_t i, const size_t j);
GSL_EXPORT void gsl_matrix_complex_set(gsl_matrix_complex * m, const size_t i, const size_t j, const gsl_complex x);

GSL_EXPORT gsl_complex * gsl_matrix_complex_ptr(gsl_matrix_complex * m, const size_t i, const size_t j);
GSL_EXPORT const gsl_complex * gsl_matrix_complex_const_ptr(const gsl_matrix_complex * m, const size_t i, const size_t j);

GSL_EXPORT void gsl_matrix_complex_set_zero (gsl_matrix_complex * m);
GSL_EXPORT void gsl_matrix_complex_set_identity (gsl_matrix_complex * m);
GSL_EXPORT void gsl_matrix_complex_set_all (gsl_matrix_complex * m, gsl_complex x);

GSL_EXPORT int gsl_matrix_complex_fread (FILE * stream, gsl_matrix_complex * m) ;
GSL_EXPORT int gsl_matrix_complex_fwrite (FILE * stream, const gsl_matrix_complex * m) ;
GSL_EXPORT int gsl_matrix_complex_fscanf (FILE * stream, gsl_matrix_complex * m);
GSL_EXPORT int gsl_matrix_complex_fprintf (FILE * stream, const gsl_matrix_complex * m, const char * format);

GSL_EXPORT int gsl_matrix_complex_memcpy(gsl_matrix_complex * dest, const gsl_matrix_complex * src);
GSL_EXPORT int gsl_matrix_complex_swap(gsl_matrix_complex * m1, gsl_matrix_complex * m2);

GSL_EXPORT int gsl_matrix_complex_swap_rows(gsl_matrix_complex * m, const size_t i, const size_t j);
GSL_EXPORT int gsl_matrix_complex_swap_columns(gsl_matrix_complex * m, const size_t i, const size_t j);
GSL_EXPORT int gsl_matrix_complex_swap_rowcol(gsl_matrix_complex * m, const size_t i, const size_t j);

GSL_EXPORT int gsl_matrix_complex_transpose (gsl_matrix_complex * m);
GSL_EXPORT int gsl_matrix_complex_transpose_memcpy (gsl_matrix_complex * dest, const gsl_matrix_complex * src);

GSL_EXPORT int gsl_matrix_complex_isnull (const gsl_matrix_complex * m);

GSL_EXPORT int gsl_matrix_complex_add (gsl_matrix_complex * a, const gsl_matrix_complex * b);
GSL_EXPORT int gsl_matrix_complex_sub (gsl_matrix_complex * a, const gsl_matrix_complex * b);
GSL_EXPORT int gsl_matrix_complex_mul_elements (gsl_matrix_complex * a, const gsl_matrix_complex * b);
GSL_EXPORT int gsl_matrix_complex_div_elements (gsl_matrix_complex * a, const gsl_matrix_complex * b);
GSL_EXPORT int gsl_matrix_complex_scale (gsl_matrix_complex * a, const gsl_complex x);
GSL_EXPORT int gsl_matrix_complex_add_constant (gsl_matrix_complex * a, const gsl_complex x);
GSL_EXPORT int gsl_matrix_complex_add_diagonal (gsl_matrix_complex * a, const gsl_complex x);

/***********************************************************************/
/* The functions below are obsolete                                    */
/***********************************************************************/
GSL_EXPORT int gsl_matrix_complex_get_row(gsl_vector_complex * v, const gsl_matrix_complex * m, const size_t i);
GSL_EXPORT int gsl_matrix_complex_get_col(gsl_vector_complex * v, const gsl_matrix_complex * m, const size_t j);
GSL_EXPORT int gsl_matrix_complex_set_row(gsl_matrix_complex * m, const size_t i, const gsl_vector_complex * v);
GSL_EXPORT int gsl_matrix_complex_set_col(gsl_matrix_complex * m, const size_t j, const gsl_vector_complex * v);

#ifdef HAVE_INLINE

extern inline
gsl_complex
gsl_matrix_complex_get(const gsl_matrix_complex * m,
                     const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  gsl_complex zero = {{0,0}};

  if (i >= m->size1)
    {
      GSL_ERROR_VAL("first index out of range", GSL_EINVAL, zero) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_VAL("second index out of range", GSL_EINVAL, zero) ;
    }
#endif
  return *(gsl_complex *)(m->data + 2*(i * m->tda + j)) ;
}

extern inline
void
gsl_matrix_complex_set(gsl_matrix_complex * m,
                     const size_t i, const size_t j, const gsl_complex x)
{
#if GSL_RANGE_CHECK
  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
  *(gsl_complex *)(m->data + 2*(i * m->tda + j)) = x ;
}

extern inline
gsl_complex *
gsl_matrix_complex_ptr(gsl_matrix_complex * m,
                             const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  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 (gsl_complex *)(m->data + 2*(i * m->tda + j)) ;
}

extern inline
const gsl_complex *
gsl_matrix_complex_const_ptr(const gsl_matrix_complex * m,
                                   const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  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 gsl_complex *)(m->data + 2*(i * m->tda + j)) ;
}

#endif /* HAVE_INLINE */

__END_DECLS

#endif /* __GSL_MATRIX_COMPLEX_DOUBLE_H__ */

max_stars_repo_name: dynaryu/vaws

max_issues_repo_name: dynaryu/vaws

max_forks_repo_name: dynaryu/vaws


text:
/*==========================================================
 * pca.c - Principal component analysis of diffraction data
 *
 * Calculates the two main components of the data
 * by Principal Component Analysis, as detailed in
 * Bernhardt et.al (2016), Biophys. J.
 *
 * The calling syntax is:
 *
 *		pca(data, mask, qx, qy, w, phi, m, n)
 *
 * Copyright 2017 Institute for X-ray Physics (University of Göttingen)
 *
 * Permission is hereby granted, free of charge, to any person obtaining
 * a copy of this software and associated documentation files (the "Software"),
 * to deal in the Software without restriction, including without limitation
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
 * and/or sell copies of the Software, and to permit persons to whom the
 * Software is furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included
 * in all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
 * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
 * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
 * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
 * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE
 * OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
 *
 *========================================================*/

#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_eigen.h>

#define PI 3.14159265

/* The computational routine */
void pca(double *x, double *mask, double *qx, double *qy, double *w, double *phi, int m, int n)
{
    int col;
    int row;
    double tr = 0;

    /* total sum of x */
    for (row=0; row<m; row++) {
        for (col=0; col<n; col++) {
            x[row*m + col] = x[row*m + col] * mask[row*m + col];
            tr += x[row*m + col];
        }
    }

    /* covariance matrix */
    double cov_qx = 0;
    double cov_mixed = 0;
    double cov_qy = 0;

    /* total sum of x */
    for (row=0; row<m; row++) {
        for (col=0; col<n; col++) {
            cov_qx += x[row*m + col] * pow(qx[row*m + col],2);
            cov_qy += x[row*m + col] * pow(qy[row*m + col],2);
            cov_mixed += x[row*m + col] * qx[row*m + col] * qy[row*m + col];
        }
    }

    double data[] = {cov_qx/tr   , cov_mixed/tr,
                     cov_mixed/tr, cov_qy/tr   };

    gsl_matrix_view mv
        = gsl_matrix_view_array (data, 2, 2);

    gsl_vector *eval = gsl_vector_alloc (2);
    gsl_matrix *evec = gsl_matrix_alloc (2, 2);

    gsl_eigen_symmv_workspace * ws
        = gsl_eigen_symmv_alloc (2);

    gsl_eigen_symmv (&mv.matrix, eval, evec, ws);

    gsl_eigen_symmv_free (ws);

    gsl_eigen_symmv_sort (eval, evec,
                        GSL_EIGEN_SORT_ABS_DESC);

    /* get eigenvalues and eigenvectors (assuming they are sorted) */
    double eval_i = gsl_vector_get (eval, 0);
    double eval_j = gsl_vector_get (eval, 1);

    gsl_vector_view evec_i
       = gsl_matrix_column (evec, 0);

    double xcomp = gsl_vector_get (&evec_i.vector, 0);
    double ycomp = gsl_vector_get (&evec_i.vector, 1);

    /* anisotropy */
    *w = (eval_i - eval_j) / (eval_i + eval_j);

    /* orientation */
    *phi = atan2(ycomp,xcomp) * 180.0 / PI;
    if (*phi < 0){
        *phi += 180;
    }

    gsl_vector_free (eval);
    gsl_matrix_free (evec);

    /* debugging information */
    /*
    printf ("%g\t%g\t%g\n", cov_qx/tr, cov_mixed/tr, cov_qy/tr );

    printf ("eigenvalue 1 = %g\n", eval_i);
    printf ("eigenvalue 2 = %g\n", eval_j);

    printf ("eigenvector = \n");
    gsl_vector_fprintf (stdout,
                        &evec_i.vector, "%g");

    printf ("x = %g\n", xcomp);
    printf ("y = %g\n", ycomp);
    printf ("w = %g\n", *w);
    printf ("phi = %g\n", *phi);
    */
}

ext: c

max_stars_repo_name: irpgoe/nanodiffraction

max_issues_repo_name: irpgoe/nanodiffraction

max_forks_repo_name: irpgoe/nanodiffraction


text:
/* Leaf bindings for the GNU Scientific Library

   Note that the GSL is GPL, not LGPL!.  If you include this code in
   the libprim build, the result will be covered under the GPL.
 */

#ifndef _LEAF_GSL_H_
#define _LEAF_GSL_H_

#include <gsl/gsl_math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_poly.h>
#include <leaf/grid.h>

/* Direct gsl_ -> grid_ adapters. */
int grid_linalg_SV_decomp(grid *A, grid *V, grid *S, grid *work);
int grid_linalg_SV_solve(grid *gU, grid *gV, grid *gS, grid *gb, grid *gx);

int grid_blas_dgemv(int transA, double alpha, grid *gA, grid *gx, double beta, grid *gy);
int grid_blas_dgemm(int transA, int transB, double alpha, grid *gA, grid *gB, double beta, grid *gC);
int grid_unfold(grid *gA, grid *gx, grid *gout);

#endif

ext: h

max_stars_repo_name: zwizwa/libprim

max_issues_repo_name: zwizwa/libprim

max_forks_repo_name: zwizwa/libprim


text:
#ifndef __RNG_H__
#define __RNG_H__
#include <gsl/gsl_rng.h>
extern const gsl_rng_type *rng_T;
extern gsl_rng *rng_R;
#endif

max_stars_repo_name: ai-ku/scode

max_issues_repo_name: ai-ku/scode

max_forks_repo_name: ai-ku/scode


text:
#include <stdio.h>
#include <gsl/gsl_fit.h>

int main (void)
{
  int n = 11;
  double x[11] = {10.0, 8.0, 13.0, 9.0,11.0,14.0,6.0,4.0 ,12.0,7.0,5.0};
  double y[11] = {  8.04,  6.95,7.68, 8.81, 8.33,9.96,7.24,4.26,10.84,4.82,5.68 };

  double c0, c1, cov00, cov01, cov11, sumsq;

  gsl_fit_linear (x, 1, y, 1, n,
                   &c0, &c1, &cov00, &cov01, &cov11,
                   &sumsq);

  printf ("best fit: Y = %g + %g X\n", c0, c1);
  printf ("covariance matrix:\n");
  printf ("[ %g, %g\n  %g, %g]\n",
          cov00, cov01, cov01, cov11);
  printf ("sumsq = %g\n", sumsq);

  printf ("\n");

  return 0;
}


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

#include <gsl/gsl_linalg.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_matrix_double.h>

int main() {
	int ret;
	gsl_matrix A;
	double data[9];
	int i, j;

	memset(&A, 0, sizeof(gsl_matrix));

	A.size1 = 3;
	A.size2 = 3;
	A.tda = 3;
	A.data = data;

	gsl_matrix_set(&A, 0, 0, 34.0);
	gsl_matrix_set(&A, 0, 1, 4.0);
	gsl_matrix_set(&A, 0, 2, 14.0);
	gsl_matrix_set(&A, 1, 0, 1.0);
	gsl_matrix_set(&A, 1, 1, 8.0);
	gsl_matrix_set(&A, 1, 2, 3.0);
	gsl_matrix_set(&A, 2, 0, 7.0);
	gsl_matrix_set(&A, 2, 1, 1.0);
	gsl_matrix_set(&A, 2, 2, 8.0);

	for (i=0; i<A.size1; i++) {
		printf((i==0) ? "A = (" : "    (");
		for (j=0; j<A.size2; j++) {
			printf(" %12.5g ", gsl_matrix_get(&A, i, j));
		}
		printf(")\n");
	}
	printf("\n");

	ret = gsl_linalg_cholesky_decomp(&A);

	for (i=0; i<A.size1; i++) {
		printf((i==0) ? "L = (" : "    (");
		for (j=0; j<A.size2; j++) {
			printf(" %12.5g ", (j <= i) ? gsl_matrix_get(&A, i, j) : 0.0);
		}
		printf(")\n");
	}
	printf("\n");

	return 0;
}

max_stars_repo_licenses: Net-SNMP Xnet

max_issues_repo_licenses: Net-SNMP Xnet

max_forks_repo_licenses: Net-SNMP Xnet


text:
#include <pygsl/solver.h>
#include <gsl/gsl_multiroots.h>

static PyObject * module = NULL;

/* accessor methods */
#define _SOLVER_METH(solver, name) \
static PyObject*  \
PyGSL_multiroot_ ## solver ## _ ## name(PyGSL_solver *self, PyObject *args) \
{ return PyGSL_solver_ret_vec(self, args, (ret_vec)gsl_multiroot_## solver ## _ ## name); }

#define FSOLVER(name)   _SOLVER_METH(fsolver, name)
#define FDFSOLVER(name) _SOLVER_METH(fdfsolver, name)
#define BSOLVER_METH(name) _SOLVER_METH(fsolver, name) _SOLVER_METH(fdfsolver, name)
BSOLVER_METH(root)
BSOLVER_METH(f)
BSOLVER_METH(dx)

static PyObject*
PyGSL_multiroot_fsolver_test_residual(PyGSL_solver *self, PyObject *args)
{
     int flag;
     double epsabs;
     gsl_multiroot_fsolver *s = self->solver;
     if(!PyArg_ParseTuple(args, "d", &epsabs))
	  return NULL;
     flag = gsl_multiroot_test_residual(s->f, epsabs);
     return PyGSL_ERROR_FLAG_TO_PYINT(flag);
}
static PyObject*
PyGSL_multiroot_fdfsolver_test_residual(PyGSL_solver *self, PyObject *args)
{
     int flag;
     double epsabs;
     gsl_multiroot_fdfsolver *s = self->solver;
     if(!PyArg_ParseTuple(args, "d", &epsabs))
	  return NULL;
     flag = gsl_multiroot_test_residual(s->f, epsabs);
     return PyGSL_ERROR_FLAG_TO_PYINT(flag);
}

static PyObject*
PyGSL_multiroot_fsolver_test_delta(PyGSL_solver *self, PyObject *args)
{
     int flag;
     double epsabs, epsrel;
     gsl_multiroot_fsolver *s = self->solver;
     if(!PyArg_ParseTuple(args, "dd", &epsabs, &epsrel))
	  return NULL;
     flag = gsl_multiroot_test_delta(s->dx, s->x, epsabs, epsrel);
     return PyGSL_ERROR_FLAG_TO_PYINT(flag);
}
static PyObject*
PyGSL_multiroot_fdfsolver_test_delta(PyGSL_solver *self, PyObject *args)
{
     int flag;
     double epsabs, epsrel;
     gsl_multiroot_fdfsolver *s = self->solver;
     if(!PyArg_ParseTuple(args, "dd", &epsabs, &epsrel))
	  return NULL;
     flag = gsl_multiroot_test_delta(s->dx, s->x, epsabs, epsrel);
     return PyGSL_ERROR_FLAG_TO_PYINT(flag);
}
int
PyGSL_multiroot_function_wrap(const gsl_vector *x, void *params, gsl_vector *f)
{
     int flag;
     PyGSL_solver *p = (PyGSL_solver *) params;
     FUNC_MESS_BEGIN();
     if(PyGSL_DEBUG_LEVEL() > 2){
       gsl_vector_fprintf(stderr, x, "x = %d");
     }
     flag =  PyGSL_function_wrap_Op_On(x, f, p->cbs[0], p->args, x->size,
				       x->size, __FUNCTION__);
     FUNC_MESS_END();
     return flag;
}

int
PyGSL_multiroot_function_wrap_df(const gsl_vector *x, void *params, gsl_matrix *J)
{
     int flag;
     PyGSL_solver *p = (PyGSL_solver *) params;
     FUNC_MESS_BEGIN();
     flag =  PyGSL_function_wrap_Op_Opn(x, J, p->cbs[1], p->args, x->size, x->size,
					__FUNCTION__);
     FUNC_MESS_END();
     return flag;
}

int
PyGSL_multiroot_function_wrap_fdf(const gsl_vector *x, void *params, gsl_vector *f, gsl_matrix *J)
{
     int flag;
     PyGSL_solver *p = (PyGSL_solver *) params;
     FUNC_MESS_BEGIN();
     flag = PyGSL_function_wrap_Op_On_Opn(x, f, J, p->cbs[2], p->args, x->size, x->size,
					  __FUNCTION__);
     FUNC_MESS_END();
     return flag;
}

PyObject *
PyGSL_multiroot_fsolver_set(PyGSL_solver *self, PyObject *pyargs, PyObject *kw)
{

     gsl_multiroot_function * c_sys;
     struct pygsl_solver_n_set info = {0, NULL, (set_m_t) gsl_multiroot_fsolver_set};
     PyObject * tmp;

     FUNC_MESS_BEGIN();
     if(self->c_sys == NULL){
	  if((c_sys = calloc(1, sizeof(gsl_multiroot_function))) == NULL){
	       PyGSL_ERROR_NULL("Could not allocate the memory for the c_sys", GSL_ENOMEM);
	  }
	  c_sys->n =self->problem_dimensions[0];
	  c_sys->f =PyGSL_multiroot_function_wrap;
	  c_sys->params=(void*)self;
	  info.c_sys = c_sys;
     }else{
	  info.c_sys = self->c_sys;
     }
     tmp =  PyGSL_solver_n_set(self, pyargs, kw, &info);
     FUNC_MESS_END();
     return tmp;
}

PyObject *
PyGSL_multiroot_fdfsolver_set(PyGSL_solver *self, PyObject *pyargs, PyObject *kw)
{

     gsl_multiroot_function_fdf * c_sys;
     struct pygsl_solver_n_set info = {1, NULL, (set_m_t)gsl_multiroot_fdfsolver_set};
     PyObject * tmp;

     FUNC_MESS_BEGIN();
     if(self->c_sys == NULL){
	  if((c_sys = calloc(1, sizeof(gsl_multiroot_function_fdf))) == NULL){
	       PyGSL_ERROR_NULL("Could not allocate the memory for the c_sys", GSL_ENOMEM);
	  }
	  c_sys->n=self->problem_dimensions[0];
	  c_sys->f  = PyGSL_multiroot_function_wrap;
	  c_sys->df = PyGSL_multiroot_function_wrap_df;
	  c_sys->fdf = PyGSL_multiroot_function_wrap_fdf;
	  c_sys->params=(void*)self;
	  info.c_sys = c_sys;
     }else{
	  info.c_sys = self->c_sys;
     }
     tmp = PyGSL_solver_n_set(self, pyargs, kw, &info);
     FUNC_MESS_END();
     return tmp;
}

static PyMethodDef PyGSL_multiroot_fmethods[] = {
     {"dx",   (PyCFunction)PyGSL_multiroot_fsolver_dx,   METH_NOARGS, NULL},
     {"f",    (PyCFunction)PyGSL_multiroot_fsolver_f,    METH_NOARGS, NULL},
     {"root", (PyCFunction)PyGSL_multiroot_fsolver_root, METH_NOARGS, NULL},
     {"set",  (PyCFunction)PyGSL_multiroot_fsolver_set,  METH_VARARGS|METH_KEYWORDS, NULL},
     {"test_delta",    (PyCFunction)PyGSL_multiroot_fsolver_test_delta, METH_VARARGS, NULL},
     {"test_residual", (PyCFunction)PyGSL_multiroot_fsolver_test_residual, METH_VARARGS, NULL},
     {NULL, NULL, 0, NULL}
};

static PyMethodDef PyGSL_multiroot_fdfmethods[] = {
     {"dx",            (PyCFunction)PyGSL_multiroot_fdfsolver_dx,   METH_NOARGS, NULL},
     {"f",             (PyCFunction)PyGSL_multiroot_fdfsolver_f,    METH_NOARGS, NULL},
     {"root",          (PyCFunction)PyGSL_multiroot_fdfsolver_root, METH_NOARGS, NULL},
     {"set",           (PyCFunction)PyGSL_multiroot_fdfsolver_set,  METH_VARARGS|METH_KEYWORDS, NULL},
     {"test_delta",    (PyCFunction)PyGSL_multiroot_fdfsolver_test_delta, METH_VARARGS, NULL},
     {"test_residual", (PyCFunction)PyGSL_multiroot_fdfsolver_test_residual, METH_VARARGS, NULL},
     {NULL, NULL, 0, NULL}
};

static const char multiroot_f_type_name[] = "F-MultiRootSolver";
static const char multiroot_fdf_type_name[] = "FdF-MultiRootSolver";

const struct _SolverStatic
multiroot_solver_f   = {{(void_m_t) gsl_multiroot_fsolver_free,
			/* gsl_multiroot_fsolver_restart */  (void_m_t) NULL,
			(name_m_t) gsl_multiroot_fsolver_name,
			 (int_m_t) gsl_multiroot_fsolver_iterate},
			1, PyGSL_multiroot_fmethods, multiroot_f_type_name},
multiroot_solver_fdf = {{(void_m_t) gsl_multiroot_fdfsolver_free,
			 /* gsl_multiroot_fdfsolver_restart (void_m_t) */ NULL,
			 (name_m_t) gsl_multiroot_fdfsolver_name,
			 (int_m_t)  gsl_multiroot_fdfsolver_iterate},
			3, PyGSL_multiroot_fdfmethods, multiroot_fdf_type_name};

static PyObject*
PyGSL_multiroot_f_init(PyObject *self, PyObject *args,
		      const gsl_multiroot_fsolver_type * type)
{

     PyObject *tmp=NULL;
     solver_alloc_struct s = {type, (void_an_t) gsl_multiroot_fsolver_alloc,
			      &multiroot_solver_f};
     FUNC_MESS_BEGIN();
     tmp = PyGSL_solver_dn_init(self, args, &s, 1);
     FUNC_MESS_END();
     return tmp;
}

static PyObject*
PyGSL_multiroot_fdf_init(PyObject *self, PyObject *args,
		      const gsl_multiroot_fdfsolver_type * type)
{

     PyObject *tmp=NULL;
     solver_alloc_struct s = {type, (void_an_t) gsl_multiroot_fdfsolver_alloc,
			      &multiroot_solver_fdf};
     FUNC_MESS_BEGIN();
     tmp = PyGSL_solver_dn_init(self, args, &s, 1);
     FUNC_MESS_END();
     return tmp;
}

#define AMROOTS_F(name)                                                  \
static PyObject* PyGSL_multiroot_init_ ## name (PyObject *self, PyObject *args)\
{                                                                             \
     PyObject *tmp = NULL;                                                    \
     FUNC_MESS_BEGIN();                                                       \
     tmp = PyGSL_multiroot_f_init(self, args,  gsl_multiroot_fsolver_ ## name); \
     if (tmp == NULL){                                                        \
	  PyGSL_add_traceback(module, __FILE__, __FUNCTION__, __LINE__); \
     }                                                                        \
     FUNC_MESS_END();                                                         \
     return tmp;                                                              \
}

#define AMROOTS_FDF(name)                                                  \
static PyObject* PyGSL_multiroot_init_ ## name (PyObject *self, PyObject *args)\
{                                                                             \
     PyObject *tmp = NULL;                                                    \
     FUNC_MESS_BEGIN();                                                       \
     tmp = PyGSL_multiroot_fdf_init(self, args,  gsl_multiroot_fdfsolver_ ## name); \
     if (tmp == NULL){                                                        \
	  PyGSL_add_traceback(module, __FILE__, __FUNCTION__, __LINE__); \
     }                                                                        \
     FUNC_MESS_END();                                                         \
     return tmp;                                                              \
}

AMROOTS_F(dnewton)
AMROOTS_F(broyden)
AMROOTS_F(hybrid)
AMROOTS_F(hybrids)

AMROOTS_FDF(newton)
AMROOTS_FDF(gnewton)
AMROOTS_FDF(hybridj)
AMROOTS_FDF(hybridsj)

static PyObject *
PyGSL_multiroot_test_delta(PyObject * self, PyObject * args)
{
     return PyGSL_solver_vvdd_i(self, args, gsl_multiroot_test_delta);
}

static PyObject *
PyGSL_multiroot_test_residual(PyObject * self, PyObject * args)
{
     return PyGSL_solver_vd_i(self, args, gsl_multiroot_test_residual);
}

static PyMethodDef mMethods[] = {
     /* solvers */
     {"dnewton" ,      PyGSL_multiroot_init_dnewton,  METH_VARARGS, NULL},
     {"broyden" ,      PyGSL_multiroot_init_broyden,  METH_VARARGS, NULL},
     {"hybrid"  ,      PyGSL_multiroot_init_hybrid ,  METH_VARARGS, NULL},
     {"hybrids" ,      PyGSL_multiroot_init_hybrids,  METH_VARARGS, NULL},
     {"newton"  ,      PyGSL_multiroot_init_newton ,  METH_VARARGS, NULL},
     {"gnewton" ,      PyGSL_multiroot_init_gnewton,  METH_VARARGS, NULL},
     {"hybridj" ,      PyGSL_multiroot_init_hybridj,  METH_VARARGS, NULL},
     {"hybridsj",      PyGSL_multiroot_init_hybridsj, METH_VARARGS, NULL},
     /* mutliroot funcs */
     {"test_delta",    PyGSL_multiroot_test_delta,     METH_VARARGS, NULL},
     {"test_residual", PyGSL_multiroot_test_residual,  METH_VARARGS, NULL},
     {NULL, NULL, 0, NULL}
};

static const char PyGSL_multiroot_module_doc[] = "XXX Missing\n";

void
initmultiroot(void)
{
     PyObject* m, *dict, *item;
     FUNC_MESS_BEGIN();

     m=Py_InitModule("multiroot", mMethods);
     module = m;
     assert(m);

     init_pygsl()
     import_pygsl_solver();
     assert(PyGSL_API);

     dict = PyModule_GetDict(m);
     if(!dict)
	  goto fail;

     if (!(item = PyString_FromString((char*)PyGSL_multiroot_module_doc))){
	  PyErr_SetString(PyExc_ImportError,
			  "I could not generate module doc string!");
	  goto fail;
     }
     if (PyDict_SetItemString(dict, "__doc__", item) != 0){
	  PyErr_SetString(PyExc_ImportError,
			  "I could not init doc string!");
	  goto fail;
     }

     FUNC_MESS_END();

 fail:
     FUNC_MESS("FAIL");
     return;
}

max_stars_repo_name: juhnowski/FishingRod

max_issues_repo_name: juhnowski/FishingRod

max_forks_repo_name: juhnowski/FishingRod


text:
#ifndef ARRUS_CORE_COMMON_COLLECTIONS_H
#define ARRUS_CORE_COMMON_COLLECTIONS_H

#include <string>
#include <vector>
#include <numeric>
#include <unordered_set>
#include <type_traits>
#include <bitset>
#include <stdexcept>

#include <gsl/span>
#include <boost/range/combine.hpp>

namespace arrus {

/**
 * Returns an array of range [start, end).
 */
template<typename T>
inline std::vector<T> getRange(T start, T end, T step = 1) {
    std::vector<T> values;
    for(T i = start; i < end; i += step) {
        values.push_back(i);
    }
    return values;
}

template<typename Out, typename In>
inline std::vector<Out> castTo(std::vector<In> values) {
    std::vector<Out> result(values.size());
    std::transform(
        std::begin(values), std::end(values),
        std::begin(result),
        [](In &value) { return Out(value); }
    );
    return result;
}

template<typename Out, typename Iterator>
inline std::vector<Out> castTo(const Iterator begin, const Iterator end) {
    std::vector<Out> result;
    std::transform(begin, end, std::back_inserter(result),
                   [](auto &value) { return Out(value); });
    return result;
}

/**
 * Returns an array that holds given value n times.
 */
template<typename T>
inline std::vector<T> getNTimes(const T value, size_t n) {
    std::vector<T> values;
    for(size_t i = 0; i < n; ++i) {
        values.push_back(value);
    }
    return values;
}

template<typename T>
inline size_t countUnique(const std::vector<T> &values) {
    return std::unordered_set<T>(std::begin(values), std::end(values)).size();
}

template<typename T>
inline bool
setContains(const std::unordered_set<T> &set, const T &value) {
    return set.find(value) != set.end();
}

template<typename T, typename U>
inline std::vector<std::pair<T, U>>
zip(const std::vector<T> &a, const std::vector<U> &b) {
    if(a.size() != b.size()) {
        throw std::runtime_error("Zipped vectors should have the same size.");
    }
    std::vector<std::pair<T, U>> res;
    res.reserve(a.size());
    for(const auto &[x, y] : boost::combine(a, b)) {
        res.emplace_back(x, y);
    }
    return res;
}

template<typename R>
inline std::vector<R>
generate(size_t nElements, std::function<R(size_t)> transformation) {
    std::vector<R> result;
    result.reserve(nElements);
    for(size_t i = 0; i < nElements; ++i) {
        result.emplace_back(transformation(i));
    }
    return result;
}

template<typename T>
inline std::vector<T>
concat(const std::vector<T> &a, const std::vector<T> &b) {
    std::vector<T> result;
    result.reserve(a.size() + b.size());
    result.insert(std::begin(result), std::begin(a), std::end(a));
    result.insert(std::end(result), std::begin(b), std::end(b));
    return result;

}

template<typename T>
inline std::vector<T>
concat(const std::vector<std::vector<T>> &a) {
    std::vector<T> result;
    size_t totalSize = 0;
    for(const auto &v : a) {
        totalSize += v.size();
    }
    result.reserve(totalSize);
    for(const auto &vec : a) {
        result.insert(std::end(result), std::begin(vec), std::end(vec));
    }
    return result;
}

template<typename T>
inline std::vector<T>
permute(const std::vector<T> &input, const std::vector<unsigned short> &perm) {
    std::vector<T> output(perm.size());
    for(size_t i = 0; i < static_cast<size_t>(perm.size()); ++i) {
        output[perm[i]] = input[i];
    }
    return output;
}

template<int size>
inline std::bitset<size>
toBitset(const std::vector<bool> &in) {
    std::bitset<size> result;
    for(size_t i = 0; i < size; ++i) {
        result[i] = in[i];
    }
    return result;
}

template<typename Map, typename K>
inline bool containsKey(Map map, const K &key) {
    return map.find(key) != std::end(map);
}

template<typename T>
inline void setValuesInRange(std::vector<T> &container, size_t start, size_t end, const T &value) {
    for(size_t i = start; i < end; ++i) {
        container[i] = value;
    }
}

template<size_t N>
inline void setValuesInRange(std::bitset<N> &container, size_t start, size_t end, const bool &value) {
    for(size_t i = start; i < end; ++i) {
        container[i] = value;
    }
}

template<typename T>
inline void setValuesInRange(std::vector<T> &container, size_t start, size_t end,
                             const std::function<T(size_t)> &generator) {
    for(size_t i = start; i < end; ++i) {
        container[i] = generator(i);
    }
}

template<typename T>
inline void setValuesInRange(std::vector<T> &container, size_t start, size_t end,
                             const std::vector<T>& source) {
    if(source.size != end-start) {
        throw std::runtime_error("Source vector should have exactly "
                            + std::to_string(end-start) + " elements "
                          "when assigning it to the selected range.");
    }
    for(size_t i = start; i < end; ++i) {
        container[i] = source[i-start];
    }
}

template<class InputIt, class T, class BinaryOp>
inline T reduce(InputIt first, InputIt last, T init, BinaryOp binaryOp) {
    T result = init;
    for(auto it = first; it != last; ++it) {
        result = binaryOp(result, *it);
    }
    return result;
}

}

#endif //ARRUS_CORE_COMMON_COLLECTIONS_H


text:
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_siman.h>

/* set up parameters for this simulated annealing run */

/* how many points do we try before stepping */
#define N_TRIES 200

/* how many iterations for each T? */
#define ITERS_FIXED_T 1000

/* max step size in random walk */
#define STEP_SIZE 1.0

/* Boltzmann constant */
#define K 1.0

/* initial temperature */
#define T_INITIAL 0.008

/* damping factor for temperature */
#define MU_T 1.003
#define T_MIN 2.0e-6

gsl_siman_params_t params
  = {N_TRIES, ITERS_FIXED_T, STEP_SIZE,
     K, T_INITIAL, MU_T, T_MIN};

/* now some functions to test in one dimension */
double E1(void *xp)
{
  double x = * ((double *) xp);

  return exp(-pow((x-1.0),2.0))*sin(8*x);
}

double M1(void *xp, void *yp)
{
  double x = *((double *) xp);
  double y = *((double *) yp);

  return fabs(x - y);
}

void S1(const gsl_rng * r, void *xp, double step_size)
{
  double old_x = *((double *) xp);
  double new_x;

  double u = gsl_rng_uniform(r);
  new_x = u * 2 * step_size - step_size + old_x;

  memcpy(xp, &new_x, sizeof(new_x));
}

void P1(void *xp)
{
  printf ("%12g", *((double *) xp));
}

int
main(int argc, char *argv[])
{
  const gsl_rng_type * T;
  gsl_rng * r;

  double x_initial = 15.5;

  gsl_rng_env_setup();

  T = gsl_rng_default;
  r = gsl_rng_alloc(T);

  gsl_siman_solve(r, &x_initial, E1, S1, M1, P1,
                  NULL, NULL, NULL,
                  sizeof(double), params);

  gsl_rng_free (r);
  return 0;
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
/*
 * \file base.h
 * \brief Header file with base requirements for codes
 *
 *
 */

#ifndef BASE_H
#define BASE_H

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <getopt.h>
#include <unistd.h>
#include <sys/stat.h>
#include <sys/types.h>
#include <time.h>

#include <bzlib.h>
#include <hdf5.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_eigen.h>

#define NAME_LENGTH 256
#define LONG_LENGTH 3000

#endif    /*  BASE_H  */

max_stars_repo_name: rcanasv/astrolibc

max_issues_repo_name: rcanasv/astrolibc

max_forks_repo_name: rcanasv/astrolibc


text:
#pragma once

#ifndef _NOGSL
#include <gsl/gsl_vector.h>
#else
#include "CustomSolver.h"
#endif
#include "ValueGrad.h"
#include "BooleanNodes.h"
#include "BooleanDAG.h"
#include "NodeVisitor.h"
#include "VarStore.h"
#include <map>
#include <set>
#include "SymbolicEvaluator.h"
#include "Interface.h"

#include <iostream>

using namespace std;

class BoolAutoDiff: public NodeVisitor, public SymbolicEvaluator
{
	BooleanDAG& bdag;
	map<string, int>& floatCtrls; // Maps float ctrl names to indices within grad vector
	int nctrls; // number of float ctrls
	gsl_vector* ctrls; // ctrl values
	vector<ValueGrad*> values; // Keeps track of values along with gradients for each node
	vector<DistanceGrad*> distances; // Keeps track of distance metric for boolean nodes
    Interface* inputValues;

public:
    int DEFAULT_INP = -32;
    BoolAutoDiff(BooleanDAG& bdag_p, map<string, int>& floatCtrls_p);

	~BoolAutoDiff(void);

	virtual void visit( SRC_node& node );
	virtual void visit( DST_node& node );
	virtual void visit( CTRL_node& node );
	virtual void visit( PLUS_node& node );
	virtual void visit( TIMES_node& node );
	virtual void visit( ARRACC_node& node );
	virtual void visit( DIV_node& node );
	virtual void visit( MOD_node& node );
	virtual void visit( NEG_node& node );
	virtual void visit( CONST_node& node );
	virtual void visit( LT_node& node );
	virtual void visit( EQ_node& node );
	virtual void visit( AND_node& node );
	virtual void visit( OR_node& node );
	virtual void visit( NOT_node& node );
	virtual void visit( ARRASS_node& node );
	virtual void visit( UFUN_node& node );
	virtual void visit( TUPLE_R_node& node );
	virtual void visit( ASSERT_node& node );

    virtual void setInputs(Interface* inputValues_p);

	virtual void run(const gsl_vector* ctrls_p, const set<int>& nodesSubset);
    virtual void run(const gsl_vector* ctrls_p);

    virtual double getErrorOnConstraint(int nodeid, gsl_vector* grad);
    double getSqrtError(bool_node* node, gsl_vector* grad);
    double getAssertError(bool_node* node, gsl_vector* grad);
    double getBoolCtrlError(bool_node* node, gsl_vector* grad);
    double getBoolExprError(bool_node* node, gsl_vector* grad);

    virtual double getErrorOnConstraint(int nodeid);
    double getSqrtError(bool_node* node);
    double getAssertError(bool_node* node);
    double getBoolCtrlError(bool_node* node);
    double getBoolExprError(bool_node* node);

    virtual double getErrorForAsserts(const set<int>& assertIds, gsl_vector* grad);
	virtual double getErrorForAssert(int assertId, gsl_vector* grad);

    double computeSingleError(int nodeid, double errorSoFar, const gsl_vector* errorGradSoFar, gsl_vector* grad);

	void setvalue(bool_node& bn, ValueGrad* v) {
		values[bn.id] = v;
	}

	ValueGrad* v(bool_node& bn) {
		ValueGrad* val = values[bn.id];
		if (val == NULL) {
			gsl_vector* g = gsl_vector_alloc(nctrls);
			val = new ValueGrad(0, g);
			setvalue(bn, val);
		}
		return val;
	}

	ValueGrad* v(bool_node* bn) {
		return v(*bn);
	}

	void setdistance(bool_node& bn, DistanceGrad* d) {
		distances[bn.id] = d;
	}

	DistanceGrad* d(bool_node& bn) {
		DistanceGrad* dist = distances[bn.id];
		if (dist == NULL) {
			gsl_vector* g = gsl_vector_alloc(nctrls);
			dist = new DistanceGrad(0, g);
			setdistance(bn, dist);
		}
		return dist;
	}

	DistanceGrad* d(bool_node* bn) {
		return d(*bn);
	}

	virtual double dist(int nid) {
		DistanceGrad* dg = d(bdag[nid]);
		return dg->dist;
	}

	virtual double dist(int nid, gsl_vector* grad) {
		DistanceGrad* dg = d(bdag[nid]);
		gsl_vector_memcpy(grad, dg->grad);
		return dg->dist;
	}

	virtual void print() {
		/*for (int i = 0; i < bdag.size(); i++) {
			if (bdag[i]->type == bool_node::ASSERT) {
				DistanceGrad* dist = d(bdag[i]->mother);
				if (dist->set) {
				double gmag = gsl_blas_dnrm2(dist->grad);
				if (dist->dist < 0.1) {
					cout << bdag[i]->mother->lprint() << endl;
					cout << dist->printFull() << endl;
				}
				}
			}
		}
		return;*/
		for (int i = 0; i < bdag.size(); i++) {
			cout << bdag[i]->lprint() << " ";
			if (bdag[i]->getOtype() == OutType::FLOAT) {
				if (v(bdag[i])->set) {
					cout << v(bdag[i])->print() << endl;
				} else {
					cout << "UNSET" << endl;
				}
			} else {
				if (d(bdag[i])->set) {
					cout << d(bdag[i])->print() << endl;
				} else {
					cout << "UNSET" << endl;
				}
			}
		}
	}
	virtual void printFull() {
		for (int i = 0; i < bdag.size(); i++) {
			cout << bdag[i]->lprint() << endl;
			if (bdag[i]->getOtype() == OutType::FLOAT) {
				if (v(bdag[i])->set) {
					cout << v(bdag[i])->printFull() << endl;
				} else {
					cout << "UNSET" <<endl;
				}
			} else {
				if(d(bdag[i])->set) {
					cout << d(bdag[i])->printFull() << endl;
				} else {
					cout << "UNSET" << endl;
				}
			}
			if (bdag[i]->type == bool_node::PLUS) {
				cout << v(bdag[i]->mother())->getVal() << " + " << v(bdag[i]->father())->getVal() << endl;
			}
		}
	}

	bool isFloat(bool_node& bn) {
		return (bn.getOtype() == OutType::FLOAT);
	}

	bool isFloat(bool_node* bn) {
		return (bn->getOtype() == OutType::FLOAT);
	}

	int getInputValue(bool_node& bn) {
		if (inputValues->hasValue(bn.id)) {
			int val = inputValues->getValue(bn.id);
			return val;
		} else {
			return DEFAULT_INP;
		}
	}

	int getInputValue(bool_node* bn) {
		return getInputValue(*bn);
	}

};

max_stars_repo_name: natebragg/sketch-backend

max_issues_repo_name: natebragg/sketch-backend

max_forks_repo_name: natebragg/sketch-backend


text:
/* specfunc/gsl_sf_result.h
 *
 * 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.
 */

/* Author:  G. Jungman */

#ifndef __GSL_SF_RESULT_H__
#define __GSL_SF_RESULT_H__

#include <gsl/gsl_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

struct gsl_sf_result_struct {
  double val;
  double err;
};
typedef struct gsl_sf_result_struct gsl_sf_result;

#define GSL_SF_RESULT_SET(r,v,e) do { (r)->val=(v); (r)->err=(e); } while(0)

struct gsl_sf_result_e10_struct {
  double val;
  double err;
  int    e10;
};
typedef struct gsl_sf_result_e10_struct gsl_sf_result_e10;

GSL_EXPORT int gsl_sf_result_smash_e(const gsl_sf_result_e10 * re, gsl_sf_result * r);

__END_DECLS

#endif /* __GSL_SF_RESULT_H__ */


text:
#pragma once
#include <gsl/span.h>
#include <msgpack.hpp>
#include <type_traits>
#include <typeindex>
#include <unordered_map>
#include <cstddef> // somehow needed for size_t?

namespace meta {
using packer = msgpack::packer<std::ostream>;
using unpacker = msgpack::unpacker;

struct Attribute {
  std::type_index typeindex;
  const void *data;
};

struct Field {
  const char *name;
  const char *qualName;
  std::type_index typeindex;
  std::size_t offset;
  void *(*getPtr)(void *object);
  gsl::span<const Attribute> attributes;

  template <typename T> T *getAs(void *object) const {
    return reinterpret_cast<T *>(getPtr(object));
  }

  template <typename T> const T *getAttribute() const {
	  for (auto&& a : attributes) {
		  if (a.typeindex == std::type_index{ typeid(T) }) {
			  return static_cast<const T*>(a.data);
		  }
	  }
	  return nullptr;
  }
};

struct Enumerator {
  const char *name;
  uint64_t value;
  gsl::span<const Attribute> attributes;

  template <typename T> const T *getAttribute() const {
	  for (auto&& a : attributes) {
		  if (a.typeindex == std::type_index{ typeid(T) }) {
			  return static_cast<const T*>(a.data);
		  }
	  }
	  return nullptr;
  }
};

enum class TypeKind { Enum, Struct, Union };

struct Enum;
struct Record;

struct Type {
  const char *name;
  TypeKind kind;
  gsl::span<const Attribute> attributes;
  void (*serialize)(packer &, const void *);
  void (*deserialize)(unpacker &, void *);

  template <typename T> const T *as() const {
    if (kind == T::t_kind)
      return static_cast<const T *>(this);
    return nullptr;
  }
};

struct Record : public Type {
  static constexpr TypeKind t_kind = TypeKind::Struct;
  Record(const char *name, gsl::span<const Field> public_fields_,
           gsl::span<const std::type_index> bases_,
           gsl::span<const Attribute> attributes_,
           void (*serialize_)(packer &, const void *),
           void (*deserialize_)(unpacker &, void *))
      : Type{name, t_kind, attributes_, serialize_, deserialize_},
        publicFields{public_fields_}, bases{bases_} {}

  gsl::span<const Field> publicFields;
  gsl::span<const std::type_index> bases;
};

struct Enum : public Type {
  static constexpr TypeKind t_kind = TypeKind::Enum;
  Enum(const char *name, gsl::span<const Enumerator> enumerators_,
         uint64_t (*getValue_)(const void *e),
         void (*setValue_)(void *, uint64_t),
         gsl::span<const Attribute> attributes_,
         void (*serialize_)(packer &, const void *),
         void (*deserialize_)(unpacker &, void *))
      : Type{name, t_kind, attributes_, serialize_, deserialize_},
        enumerators{enumerators_}, getValue{getValue_},
        setValue{setValue_} {}

  gsl::span<const Enumerator> enumerators;

  int findEnumeratorIndex(const void *obj) const {
    auto v = getValue(obj);
    for (int i = 0; i < enumerators.size(); ++i) {
      if (enumerators[i].value == v)
        return i;
    }
    return -1;
  }

  const Enumerator *findEnumerator(const void *obj) const {
    int i = findEnumeratorIndex(obj);
    if (i < 0)
      return nullptr;
    return &enumerators[i];
  }

  // manipulation functions
  uint64_t (*getValue)(const void *obj);
  void (*setValue)(void *obj, uint64_t value);
};

template <typename T> const Type *typeOf();
const Type *typeOf(std::type_index ti);

template <typename T> void serialize(packer& p, const T& data)
{
	serialize(p, typeOf<T>(), &data);
}

template <typename T> void serialize_dynamic(packer& p, const T& data)
{
	// always use the dynamic type
	serialize(p, typeOf(typeid(data)), &data);
}

inline void serialize(packer& p, const Type* ty, const void* data)
{
	if (!ty) return;
	ty->serialize(p, data);
}

// enum value-to-string and string-to-value
template <typename T, typename = std::enable_if_t<std::is_enum<T>::value>>
const char *getEnumeratorName(T constant);
template <typename T, typename = std::enable_if_t<std::is_enum<T>::value>>
T getEnumeratorValue(const char *name);
}

max_stars_repo_name: ennis/autograph-pipelines

max_issues_repo_name: ennis/autograph-pipelines

max_forks_repo_name: ennis/autograph-pipelines


text:
// align.h
//
// prototypes for align.cc
//
// dw, 12/12/3

#include <gsl/gsl_matrix.h>

#include "pdb.h"
#include <map>

// find rigid body transformation z=Rx+t minimizing
// least squares error ||y-z||**2
// find rigid body transformation z=Rx+t minimizing
// least squares error ||y-z||**2
void align(gsl_matrix_const_view , gsl_matrix_const_view, gsl_matrix *, gsl_vector *, double *);

double align(vector<atom>&, vector<atom>&, map<int, int>&, int, vector<atom>&,
    vector<atom>&, map<int, int>&, int);

max_stars_repo_name: tecdatalab/legacy

max_issues_repo_name: tecdatalab/legacy

max_forks_repo_name: tecdatalab/biostructure


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

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

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

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

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

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

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_sort_float.h>
#include <gsl/gsl_statistics_float.h>
#include <gsl/gsl_integration.h>
#include "psrsalsa.h"
int filterPApoints(datafile_definition *datafile, verbose_definition verbose)
{
  int dPa_polnr;
  long i, j, nrpoints;
  float *olddata;
  if(datafile->poltype != POLTYPE_ILVPAdPA && datafile->poltype != POLTYPE_PAdPA && datafile->poltype != POLTYPE_ILVPAdPATEldEl) {
    printerror(verbose.debug, "ERROR filterPApoints: Data doesn't appear to have poltype ILVPAdPA, PAdPA or ILVPAdPATEldEl.");
    return 0;
  }
  if(datafile->poltype == POLTYPE_ILVPAdPA && datafile->NrPols != 5) {
    printerror(verbose.debug, "ERROR filterPApoints: 5 polarization channels were expected, but there are only %ld.", datafile->NrPols);
    return 0;
  }else if(datafile->poltype == POLTYPE_ILVPAdPATEldEl && datafile->NrPols != 8) {
    printerror(verbose.debug, "ERROR filterPApoints: 8 polarization channels were expected, but there are only %ld.", datafile->NrPols);
    return 0;
  }else if(datafile->poltype == POLTYPE_PAdPA && datafile->NrPols != 2) {
    printerror(verbose.debug, "ERROR filterPApoints: 2 polarization channels were expected, but there are only %ld.", datafile->NrPols);
    return 0;
  }
  if(datafile->NrSubints > 1 || datafile->NrFreqChan > 1) {
    printerror(verbose.debug, "ERROR filterPApoints: Can only do this opperation if there is one subint and one frequency channel.");
    return 0;
  }
  if(datafile->tsampMode != TSAMPMODE_LONGITUDELIST) {
    printerror(verbose.debug, "ERROR filterPApoints: Expected pulse longitudes to be defined.");
    return 0;
  }
  if(datafile->poltype == POLTYPE_ILVPAdPA || datafile->poltype == POLTYPE_ILVPAdPATEldEl) {
    dPa_polnr = 4;
  }else if(datafile->poltype == POLTYPE_PAdPA) {
    dPa_polnr = 1;
  }
  nrpoints = 0;
  for(i = 0; i < datafile->NrBins; i++) {
    if(datafile->data[i+dPa_polnr*datafile->NrBins] > 0) {
      nrpoints++;
    }
  }
  if(verbose.verbose)
    printf("Keeping %ld significant PA points\n", nrpoints);
  olddata = datafile->data;
  datafile->data = (float *)malloc(nrpoints*datafile->NrPols*sizeof(float));
  if(datafile->data == NULL) {
    printerror(verbose.debug, "ERROR filterPApoints: Memory allocation error.");
    return 0;
  }
  j = 0;
  for(i = 0; i < datafile->NrBins; i++) {
    if(olddata[i+dPa_polnr*datafile->NrBins] > 0) {
      if(datafile->poltype == POLTYPE_ILVPAdPA) {
 datafile->data[j+0*nrpoints] = olddata[i+0*datafile->NrBins];
 datafile->data[j+1*nrpoints] = olddata[i+1*datafile->NrBins];
 datafile->data[j+2*nrpoints] = olddata[i+2*datafile->NrBins];
 datafile->data[j+3*nrpoints] = olddata[i+3*datafile->NrBins];
 datafile->data[j+4*nrpoints] = olddata[i+4*datafile->NrBins];
      }else if(datafile->poltype == POLTYPE_ILVPAdPATEldEl) {
 datafile->data[j+0*nrpoints] = olddata[i+0*datafile->NrBins];
 datafile->data[j+1*nrpoints] = olddata[i+1*datafile->NrBins];
 datafile->data[j+2*nrpoints] = olddata[i+2*datafile->NrBins];
 datafile->data[j+3*nrpoints] = olddata[i+3*datafile->NrBins];
 datafile->data[j+4*nrpoints] = olddata[i+4*datafile->NrBins];
 datafile->data[j+5*nrpoints] = olddata[i+5*datafile->NrBins];
 datafile->data[j+6*nrpoints] = olddata[i+6*datafile->NrBins];
 datafile->data[j+7*nrpoints] = olddata[i+7*datafile->NrBins];
      }else {
 datafile->data[j+0*nrpoints] = olddata[i+0*datafile->NrBins];
 datafile->data[j+1*nrpoints] = olddata[i+1*datafile->NrBins];
      }
      datafile->tsamp_list[j] = datafile->tsamp_list[i];
      j++;
    }
  }
  free(olddata);
  datafile->NrBins = nrpoints;
  return datafile->NrBins;
}
int make_paswing_fromIQUV(datafile_definition *datafile, int extended, pulselongitude_regions_definition onpulse, int normalize, int correctLbias, float correctQV, float correctV, int nolongitudes, float loffset, float paoffset, datafile_definition *rms_file, float rebin_factor, verbose_definition verbose)
{
  int indent, rms_file_specified;
  long i, j, NrOffpulseBins, pulsenr, freqnr, output_nr_pols;
  float ymax, baseline_intensity, RMSQ, RMSU, *Loffpulse, *Poffpulse, medianL, medianP, *newdata, *newdata_rms;
  rms_file_specified = 1;
  if(rms_file == NULL) {
    rms_file = datafile;
    rms_file_specified = 0;
  }
  if(verbose.verbose) {
    for(indent = 0; indent < verbose.indent; indent++) printf(" ");
    printf("Constructing PA and degree of linear polarization");
    if(extended)
      printf(", total polarization and ellipticity");
    if(rms_file_specified)
      printf(" (using a seperate file to determine the off-pulse rms)");
    printf("\n");
    for(indent = 0; indent < verbose.indent; indent++) printf(" ");
    printf("  Reference frequency for PA is ");
    if(datafile->isDeFarad) {
      if((datafile->freq_ref > -1.1 && datafile->freq_ref < -0.9) || (datafile->freq_ref > 0.99e10 && datafile->freq_ref < 1.01e10))
 printf("infinity\n");
      else if(datafile->freq_ref < 0)
 printf("unknown\n");
      else
 printf("%f MHz\n", datafile->freq_ref);
    }else {
      if(datafile->NrFreqChan == 1)
 printf("%lf MHz\n", get_centre_frequency(*datafile, verbose));
      else
 printf("observing frequencies of individual frequency channels\n");
    }
    for(indent = 0; indent < verbose.indent; indent++) printf(" ");
    printf("  ");
    switch(correctLbias) {
    case -1: printf("No L de-bias applied"); break;
    case 0: printf("De-bias L using median noise subtraction"); break;
    case 1: printf("De-bias L using Wardle & Kronberg correction"); break;
    default: printerror(verbose.debug, "ERROR make_paswing_fromIQUV: Undefined L de-bias method specified."); return 0;
    }
    if(correctQV != 1 || correctV != 1)
      printf(", Q correction factor %f, V correction factor %f", 1.0/correctQV, 1.0/(correctQV*correctV));
    if(normalize)
      printf(", output is normalised");
    if(loffset != 0)
      printf(", pulse longitude shifted by %f deg\n", loffset);
    if(paoffset != 0)
      printf(", PA shifted by %f deg\n", paoffset);
    printf("\n");
    if(extended) {
      printwarning(verbose.debug, "WARNING make_paswing_fromIQUV: Total polarization is computed, and the median off-pulse value is subtracted (probably not a very good idea).");
    }
  }
  if(datafile->NrPols != 4 || rms_file->NrPols != 4) {
    fflush(stdout);
    printerror(verbose.debug, "ERROR make_paswing_fromIQUV: Expected 4 input polarizations.");
    return 0;
  }
  if(datafile->poltype != POLTYPE_STOKES) {
    if(datafile->poltype == POLTYPE_UNKNOWN) {
      printwarning(verbose.debug, "WARNING make_paswing_fromIQUV: Polarization state unknown, it is assumed the data are Stokes parameters.");
    }else {
      printerror(verbose.debug, "ERROR make_paswing_fromIQUV: Convert data into Stokes parameters first.");
      return 0;
    }
  }
  if(rms_file_specified) {
    if(rms_file->poltype != POLTYPE_STOKES) {
      if(rms_file->poltype == POLTYPE_UNKNOWN) {
 printwarning(verbose.debug, "WARNING make_paswing_fromIQUV: Polarization state of the data to be used to determine the off-pulse rms is unknown, it is assumed the data are Stokes parameters.");
      }else {
 printerror(verbose.debug, "ERROR make_paswing_fromIQUV: Convert data to be used to determine the off-pulse rms into Stokes parameters first.");
 return 0;
      }
    }
  }
  if(datafile->tsampMode != TSAMPMODE_FIXEDTSAMP) {
    fflush(stdout);
    printerror(verbose.debug, "ERROR make_paswing_fromIQUV: Expects input to have a regular sampling.");
    return 0;
  }
  if(correctQV == 0) {
    fflush(stdout);
    printerror(verbose.debug, "ERROR make_paswing_fromIQUV: correctQV is set to zero, you probably want this to be 1.");
    return 0;
  }
  if(correctV == 0) {
    fflush(stdout);
    printerror(verbose.debug, "ERROR make_paswing_fromIQUV: correctV is set to zero, you probably want this to be 1.");
    return 0;
  }
  if(datafile->isDebase == 0) {
    fflush(stdout);
    printerror(verbose.debug, "ERROR make_paswing_fromIQUV: Please remove baseline first, i.e. use pmod -debase.");
    return 0;
  }else if(datafile->isDebase != 1) {
    fflush(stdout);
    printwarning(verbose.debug, "WARNING make_paswing_fromIQUV: Unknown baseline state. It is assumed the baseline has already removed from the data.");
  }
  if(rms_file_specified) {
    if(rms_file->isDebase == 0) {
      fflush(stdout);
      printerror(verbose.debug, "ERROR make_paswing_fromIQUV: Please remove baseline first, i.e. use pmod -debase.");
      return 0;
    }else if(rms_file->isDebase != 1) {
      fflush(stdout);
      printwarning(verbose.debug, "WARNING make_paswing_fromIQUV: Unknown baseline state. It is assumed the baseline has already removed from the data.");
    }
  }
  if(rms_file_specified) {
    if(datafile->NrSubints != rms_file->NrSubints) {
      printerror(verbose.debug, "ERROR make_paswing_fromIQUV: Number of subintegrations is different in the data to be used to determine the off-pulse rms.");
 return 0;
    }
    if(datafile->NrFreqChan != rms_file->NrFreqChan) {
      printerror(verbose.debug, "ERROR make_paswing_fromIQUV: Number of frequency channels is different in the data to be used to determine the off-pulse rms (%ld != %ld).", rms_file->NrFreqChan, datafile->NrFreqChan);
 return 0;
    }
    if(correctLbias == 0) {
      printerror(verbose.debug, "ERROR make_paswing_fromIQUV: Subtracting the median of L is not supported when a separate file is used for the offpulse statistics.");
      return 0;
    }
    if(extended) {
      printerror(verbose.debug, "ERROR make_paswing_fromIQUV: Subtracting the median of P is not supported when a separate file is used for the offpulse statistics.");
      return 0;
    }
  }
  if(datafile->offpulse_rms != NULL) {
    free(datafile->offpulse_rms);
  }
  Loffpulse = (float *)malloc((rms_file->NrBins)*sizeof(float));
  Poffpulse = (float *)malloc((rms_file->NrBins)*sizeof(float));
  if(extended) {
    output_nr_pols = 8;
  }else {
    output_nr_pols = 5;
  }
  newdata = (float *)malloc(datafile->NrBins*datafile->NrSubints*datafile->NrFreqChan*output_nr_pols*sizeof(float));
  if(rms_file_specified) {
    newdata_rms = (float *)malloc(rms_file->NrBins*rms_file->NrSubints*rms_file->NrFreqChan*output_nr_pols*sizeof(float));
  }else {
    newdata_rms = newdata;
  }
  datafile->offpulse_rms = (float *)malloc(datafile->NrSubints*datafile->NrFreqChan*output_nr_pols*sizeof(float));
  if(Loffpulse == NULL || Poffpulse == NULL || newdata == NULL || datafile->offpulse_rms == NULL || newdata_rms == NULL) {
    fflush(stdout);
    printerror(verbose.debug, "ERROR make_paswing_fromIQUV: Memory allocation error.");
    return 0;
  }
  if(nolongitudes == 0) {
    datafile->tsamp_list = (double *)malloc(datafile->NrBins*sizeof(double));
    if(datafile->tsamp_list == NULL) {
      fflush(stdout);
      printerror(verbose.debug, "ERROR make_paswing_fromIQUV: Memory allocation error.");
      return 0;
    }
    for(j = 0; j < (datafile->NrBins); j++) {
      datafile->tsamp_list[j] = get_pulse_longitude(*datafile, 0, j, verbose);
      datafile->tsamp_list[j] += loffset;
    }
  }
  if(normalize && (datafile->NrSubints > 1 || datafile->NrFreqChan > 1)) {
    fflush(stdout);
    printwarning(verbose.debug, "WARNING make_paswing_fromIQUV: Normalization will cause all subintegrations/frequency channels to be normalised individually. This may not be desired.");
  }
  long sindex_I, sindex_Q, sindex_U, sindex_V, sindex_I_rms, sindex_Q_rms, sindex_U_rms, sindex_V_rms, newindex_I, newindex_L, newindex_V, newindex_Pa, newindex_dPa, newindex_T, newindex_Ell, newindex_dEll, newindex_L_rms, newindex_T_rms;
  for(pulsenr = 0; pulsenr < datafile->NrSubints; pulsenr++) {
    for(freqnr = 0; freqnr < datafile->NrFreqChan; freqnr++) {
      sindex_I = datafile->NrBins*(0+datafile->NrPols*(freqnr+pulsenr*datafile->NrFreqChan));
      sindex_Q = datafile->NrBins*(1+datafile->NrPols*(freqnr+pulsenr*datafile->NrFreqChan));
      sindex_U = datafile->NrBins*(2+datafile->NrPols*(freqnr+pulsenr*datafile->NrFreqChan));
      sindex_V = datafile->NrBins*(3+datafile->NrPols*(freqnr+pulsenr*datafile->NrFreqChan));
      newindex_I = datafile->NrBins*(0+output_nr_pols*(freqnr+datafile->NrFreqChan*pulsenr));
      newindex_L = datafile->NrBins*(1+output_nr_pols*(freqnr+datafile->NrFreqChan*pulsenr));
      newindex_V = datafile->NrBins*(2+output_nr_pols*(freqnr+datafile->NrFreqChan*pulsenr));
      newindex_Pa = datafile->NrBins*(3+output_nr_pols*(freqnr+datafile->NrFreqChan*pulsenr));
      newindex_dPa = datafile->NrBins*(4+output_nr_pols*(freqnr+datafile->NrFreqChan*pulsenr));
      if(extended) {
 newindex_T = datafile->NrBins*(5+output_nr_pols*(freqnr+datafile->NrFreqChan*pulsenr));
 newindex_Ell = datafile->NrBins*(6+output_nr_pols*(freqnr+datafile->NrFreqChan*pulsenr));
 newindex_dEll = datafile->NrBins*(7+output_nr_pols*(freqnr+datafile->NrFreqChan*pulsenr));
      }
      if(rms_file_specified) {
 sindex_I_rms = rms_file->NrBins*(0+rms_file->NrPols*(freqnr+pulsenr*rms_file->NrFreqChan));
 sindex_Q_rms = rms_file->NrBins*(1+rms_file->NrPols*(freqnr+pulsenr*rms_file->NrFreqChan));
 sindex_U_rms = rms_file->NrBins*(2+rms_file->NrPols*(freqnr+pulsenr*rms_file->NrFreqChan));
 sindex_V_rms = rms_file->NrBins*(3+rms_file->NrPols*(freqnr+pulsenr*rms_file->NrFreqChan));
 newindex_L_rms = rms_file->NrBins*(1+output_nr_pols*(freqnr+rms_file->NrFreqChan*pulsenr));
 if(extended) {
   newindex_T_rms = rms_file->NrBins*(5+output_nr_pols*(freqnr+rms_file->NrFreqChan*pulsenr));
 }
      }else {
 sindex_I_rms = sindex_I;
 sindex_Q_rms = sindex_Q;
 sindex_U_rms = sindex_U;
 sindex_V_rms = sindex_V;
 newindex_L_rms = newindex_L;
 if(extended) {
   newindex_T_rms = newindex_T;
 }
      }
      if(normalize == 0) {
 ymax = 1;
      }else {
 ymax = datafile->data[sindex_I];
 for(j = 1; j < (datafile->NrBins); j++) {
   if(datafile->data[sindex_I+j] > ymax)
     ymax = datafile->data[sindex_I+j];
 }
      }
      if(ymax == 0)
 ymax = 1;
      for(j = 0; j < (datafile->NrBins); j++) {
 datafile->data[sindex_I + j] /= ymax;
 datafile->data[sindex_Q + j] /= correctQV*ymax;
 datafile->data[sindex_U + j] /= ymax;
 datafile->data[sindex_V + j] /= correctV*correctQV*ymax;
 newdata[j+newindex_L] = sqrt((datafile->data[sindex_Q+j])*(datafile->data[sindex_Q+j])+(datafile->data[sindex_U+j])*(datafile->data[sindex_U+j]));
 if(extended) {
   newdata[j+newindex_T] = sqrt(newdata[j+newindex_L]*newdata[j+newindex_L]+datafile->data[sindex_V+j]*datafile->data[sindex_V+j]);
 }
      }
      if(rms_file_specified) {
 for(j = 0; j < (rms_file->NrBins); j++) {
   rms_file->data[sindex_I_rms + j] /= ymax;
   rms_file->data[sindex_Q_rms + j] /= correctQV*ymax;
   rms_file->data[sindex_U_rms + j] /= ymax;
   rms_file->data[sindex_V_rms + j] /= correctV*correctQV*ymax;
   newdata_rms[j+newindex_L_rms] = sqrt((rms_file->data[sindex_Q_rms+j])*(rms_file->data[sindex_Q_rms+j])+(rms_file->data[sindex_U_rms+j])*(rms_file->data[sindex_U_rms+j]));
   if(extended) {
     newdata_rms[j+newindex_T_rms] = sqrt(newdata[j+newindex_L_rms]*newdata[j+newindex_L_rms]+rms_file->data[sindex_V_rms+j]*rms_file->data[sindex_V_rms+j]);
   }
 }
      }
      datafile->offpulse_rms[0+output_nr_pols*(freqnr + datafile->NrFreqChan*pulsenr)] = 0;
      datafile->offpulse_rms[1+output_nr_pols*(freqnr + datafile->NrFreqChan*pulsenr)] = 0;
      datafile->offpulse_rms[2+output_nr_pols*(freqnr + datafile->NrFreqChan*pulsenr)] = 0;
      datafile->offpulse_rms[3+output_nr_pols*(freqnr + datafile->NrFreqChan*pulsenr)] = -1;
      datafile->offpulse_rms[4+output_nr_pols*(freqnr + datafile->NrFreqChan*pulsenr)] = -1;
      if(extended) {
 datafile->offpulse_rms[5+output_nr_pols*(freqnr + datafile->NrFreqChan*pulsenr)] = 0;
 datafile->offpulse_rms[6+output_nr_pols*(freqnr + datafile->NrFreqChan*pulsenr)] = -1;
 datafile->offpulse_rms[7+output_nr_pols*(freqnr + datafile->NrFreqChan*pulsenr)] = -1;
      }
      baseline_intensity = 0;
      RMSQ = 0;
      RMSU = 0;
      NrOffpulseBins = 0;
      for(i = 0; i < (rms_file->NrBins); i++) {
 Loffpulse[i] = 0;
 if(checkRegions(i, &onpulse, 0, verbose) == 0) {
   NrOffpulseBins++;
   baseline_intensity += rms_file->data[sindex_I_rms+i];
   RMSQ += (rms_file->data[sindex_Q_rms+i])*(rms_file->data[sindex_Q_rms+i]);
   RMSU += (rms_file->data[sindex_U_rms+i])*(rms_file->data[sindex_U_rms+i]);
   datafile->offpulse_rms[0+output_nr_pols*(freqnr + datafile->NrFreqChan*pulsenr)] += (rms_file->data[sindex_I_rms+i])*(rms_file->data[sindex_I_rms+i]);
   datafile->offpulse_rms[1+output_nr_pols*(freqnr + datafile->NrFreqChan*pulsenr)] += (newdata_rms[i+newindex_L_rms])*(newdata_rms[i+newindex_L_rms]);
   datafile->offpulse_rms[2+output_nr_pols*(freqnr + datafile->NrFreqChan*pulsenr)] += (rms_file->data[sindex_V_rms+i])*(rms_file->data[sindex_V_rms+i]);
   if(extended) {
     datafile->offpulse_rms[5+output_nr_pols*(freqnr + datafile->NrFreqChan*pulsenr)] += (newdata_rms[i+newindex_L_rms])*(newdata_rms[i+newindex_L_rms]) + (rms_file->data[sindex_V_rms+i])*(rms_file->data[sindex_V_rms+i]);
   }
   Loffpulse[NrOffpulseBins-1] = newdata_rms[i+newindex_L_rms];
   if(extended)
     Poffpulse[NrOffpulseBins-1] = newdata_rms[i+newindex_T_rms];
 }
      }
      baseline_intensity /= (float)NrOffpulseBins;
      RMSQ = sqrt(RMSQ/(float)NrOffpulseBins);
      RMSU = sqrt(RMSU/(float)NrOffpulseBins);
      datafile->offpulse_rms[0+output_nr_pols*(freqnr + datafile->NrFreqChan*pulsenr)] = sqrt(datafile->offpulse_rms[0+output_nr_pols*(freqnr + datafile->NrFreqChan*pulsenr)]/(float)NrOffpulseBins);
      datafile->offpulse_rms[1+output_nr_pols*(freqnr + datafile->NrFreqChan*pulsenr)] = sqrt(datafile->offpulse_rms[1+output_nr_pols*(freqnr + datafile->NrFreqChan*pulsenr)]/(float)NrOffpulseBins);
      datafile->offpulse_rms[2+output_nr_pols*(freqnr + datafile->NrFreqChan*pulsenr)] = sqrt(datafile->offpulse_rms[2+output_nr_pols*(freqnr + datafile->NrFreqChan*pulsenr)]/(float)NrOffpulseBins);
      if(extended) {
 datafile->offpulse_rms[5+output_nr_pols*(freqnr + datafile->NrFreqChan*pulsenr)] = sqrt(datafile->offpulse_rms[5+output_nr_pols*(freqnr + datafile->NrFreqChan*pulsenr)]/(float)NrOffpulseBins);
      }
      if(rms_file_specified) {
 float scale = 1.0/sqrt(rebin_factor);
 RMSQ *= scale;
 RMSU *= scale;
 datafile->offpulse_rms[0+output_nr_pols*(freqnr + datafile->NrFreqChan*pulsenr)] *= scale;
 datafile->offpulse_rms[1+output_nr_pols*(freqnr + datafile->NrFreqChan*pulsenr)] *= scale;
 datafile->offpulse_rms[2+output_nr_pols*(freqnr + datafile->NrFreqChan*pulsenr)] *= scale;
 if(extended) {
   datafile->offpulse_rms[5+output_nr_pols*(freqnr + datafile->NrFreqChan*pulsenr)] *= scale;
 }
      }
      if(verbose.verbose) {
 if((freqnr == 0 && pulsenr == 0) || verbose.debug) {
   for(indent = 0; indent < verbose.indent; indent++) printf(" ");
   fprintf(stdout, "  PA conversion output for subint %ld frequency channel %ld:\n", pulsenr, freqnr);
   for(indent = 0; indent < verbose.indent; indent++) printf(" ");
   fprintf(stdout, "    Average baseline Stokes I: %f\n", baseline_intensity);
   for(indent = 0; indent < verbose.indent; indent++) printf(" ");
   fprintf(stdout, "    RMS I:                  %f\n", datafile->offpulse_rms[0+output_nr_pols*(freqnr + datafile->NrFreqChan*pulsenr)]);
   for(indent = 0; indent < verbose.indent; indent++) printf(" ");
   fprintf(stdout, "    RMS Q:                  %f\n", RMSQ);
   for(indent = 0; indent < verbose.indent; indent++) printf(" ");
   fprintf(stdout, "    RMS U:                  %f\n", RMSU);
   for(indent = 0; indent < verbose.indent; indent++) printf(" ");
   fprintf(stdout, "    RMS V:                  %f\n", datafile->offpulse_rms[2+output_nr_pols*(freqnr + datafile->NrFreqChan*pulsenr)]);
   for(indent = 0; indent < verbose.indent; indent++) printf(" ");
   fprintf(stdout, "    RMS L (before de-bias): %f\n", datafile->offpulse_rms[1+output_nr_pols*(freqnr + datafile->NrFreqChan*pulsenr)]);
   if(extended) {
     for(indent = 0; indent < verbose.indent; indent++) printf(" ");
     fprintf(stdout, "    RMS sqrt(Q^2+U^2+V^2):  %f\n", datafile->offpulse_rms[5+output_nr_pols*(freqnr + datafile->NrFreqChan*pulsenr)]);
   }
 }
      }
      gsl_sort_float (Loffpulse, 1, NrOffpulseBins);
      medianL = gsl_stats_float_median_from_sorted_data(Loffpulse, 1, NrOffpulseBins);
      fflush(stdout);
      if((verbose.verbose && freqnr == 0 && pulsenr == 0) || verbose.debug) {
 for(indent = 0; indent < verbose.indent; indent++) printf(" ");
 fprintf(stdout, "    Median L: %f\n", medianL);
      }
      if(extended) {
 gsl_sort_float (Poffpulse, 1, NrOffpulseBins);
 medianP = gsl_stats_float_median_from_sorted_data(Poffpulse, 1, NrOffpulseBins);
 fflush(stdout);
 if((verbose.verbose && freqnr == 0 && pulsenr == 0) || verbose.debug) {
   for(indent = 0; indent < verbose.indent; indent++) printf(" ");
   fprintf(stdout, "    Median sqrt(Q^2+U^2+V^2): %f\n", medianP);
 }
      }
      for(j = 0; j < (datafile->NrBins); j++) {
 if(correctLbias == 1) {
   float junk = (0.5*(RMSQ+RMSU)/newdata[j+newindex_L]);
   if(junk < 1)
     newdata[j+newindex_L] *= sqrt(1.0-junk*junk);
   else
     newdata[j+newindex_L] = 0.0;
 }else if(correctLbias == 0) {
   newdata[j+newindex_L] -= medianL;
 }
 if(extended)
   newdata[j+newindex_T] -= medianP;
      }
      for(i = 0; i < (datafile->NrBins); i++) {
 newdata[i+newindex_Pa] = 90.0*atan2(datafile->data[sindex_U+i],datafile->data[sindex_Q+i])/M_PI;
 if(paoffset) {
   newdata[i+newindex_Pa] += paoffset;
   newdata[i+newindex_Pa] = derotate_180_small_double(newdata[i+newindex_Pa]);
 }
 if(datafile->data[i+sindex_Q] != 0 || datafile->data[i+sindex_U] != 0) {
   newdata[i+newindex_dPa] = sqrt((datafile->data[sindex_Q+i]*RMSU)*(datafile->data[sindex_Q+i]*RMSU) + (datafile->data[sindex_U+i]*RMSQ)*(datafile->data[sindex_U+i]*RMSQ));
   newdata[i+newindex_dPa] /= 2.0*(datafile->data[i+sindex_Q]*datafile->data[i+sindex_Q] + datafile->data[i+sindex_U]*datafile->data[i+sindex_U]);
   newdata[i+newindex_dPa] *= 180.0/M_PI;
 }else {
   newdata[i+newindex_dPa] = 0;
 }
 if(extended) {
   newdata[i+newindex_Ell] = 0;
     newdata[i+newindex_dEll] = 0;
 }
      }
      for(j = 0; j < (datafile->NrBins); j++) {
 newdata[j+newindex_I] = datafile->data[j+sindex_I];
 newdata[j+newindex_V] = datafile->data[j+sindex_V];
      }
    }
  }
  free(datafile->data);
  datafile->data = newdata;
  if(rms_file_specified) {
    free(newdata_rms);
  }
  if(nolongitudes == 0) {
    datafile->tsampMode = TSAMPMODE_LONGITUDELIST;
  }
  datafile->NrPols = output_nr_pols;
  if(extended) {
    datafile->poltype = POLTYPE_ILVPAdPATEldEl;
  }else {
    datafile->poltype = POLTYPE_ILVPAdPA;
  }
  free(Loffpulse);
  free(Poffpulse);
  return 1;
}
int writePPOLHeader(datafile_definition datafile, int argc, char **argv, verbose_definition verbose)
{
  char *txt;
  txt = malloc(10000);
  if(txt == NULL) {
    fflush(stdout);
    printerror(verbose.debug, "ERROR writePPOLHeader: Memory allocation error.");
    return 0;
  }
  constructCommandLineString(txt, 10000, argc, argv, verbose);
  fprintf(datafile.fptr_hdr, "#ppol file: %s\n", txt);
  free(txt);
  return 1;
}
int readPPOLHeader(datafile_definition *datafile, int extended, verbose_definition verbose)
{
  float dummy_float;
  int ret, maxlinelength, nrwords;
  char *txt, *ret_ptr, *word_ptr;
  maxlinelength = 2000;
  txt = malloc(maxlinelength);
  if(txt == NULL) {
    fflush(stdout);
    printerror(verbose.debug, "ERROR readPPOLHeader: Memory allocation error.");
    return 0;
  }
  datafile->isFolded = 1;
  datafile->foldMode = FOLDMODE_FIXEDPERIOD;
  datafile->fixedPeriod = 0;
  datafile->tsampMode = TSAMPMODE_LONGITUDELIST;
  datafile->fixedtsamp = 0;
  datafile->tsubMode = TSUBMODE_FIXEDTSUB;
  if(datafile->tsub_list != NULL)
    free(datafile->tsub_list);
  datafile->tsub_list = (double *)malloc(sizeof(double));
  if(datafile->tsub_list == NULL) {
    fflush(stdout);
    printerror(verbose.debug, "ERROR readPPOLHeader: Memory allocation error");
    return 0;
  }
  datafile->tsub_list[0] = 0;
  datafile->NrSubints = 1;
  datafile->NrFreqChan = 1;
  datafile->datastart = 0;
  rewind(datafile->fptr);
  ret = fread(txt, 1, 3, datafile->fptr);
  txt[3] = 0;
  if(ret != 3) {
    fflush(stdout);
    printerror(verbose.debug, "ERROR readPPOLHeader: cannot read from file.");
    free(txt);
    return 0;
  }
  if(strcmp(txt, "#pp") != 0
) {
    fflush(stdout);
    printwarning(verbose.debug, "WARNING readPPOLHeader: File does not appear to be in PPOL or PPOLSHORT format. I will try to load file, but this will probably fail. Did you run ppol first?");
  }
  skipallhashedlines(datafile);
  datafile->NrBins = 0;
  dummy_float = 0;
  do {
    ret_ptr = fgets(txt, maxlinelength, datafile->fptr);
    if(ret_ptr != NULL) {
      if(txt[0] != '#') {
 if(extended) {
   word_ptr = pickWordFromString(txt, 2, &nrwords, 1, ' ', verbose);
   if(nrwords != 10 && nrwords != 14) {
     fflush(stdout);
     printerror(verbose.debug, "ERROR readPPOLHeader: Line should have 10 or 14 words, got %d", nrwords);
     if(nrwords == 3)
       printerror(verbose.debug, "                             Maybe file is in format %s?", returnFileFormat_str(PPOL_SHORT_format));
     printerror(verbose.debug, "                             Line: '%s'.", txt);
     free(txt);
     return 0;
   }
   if(nrwords == 10) {
     datafile->poltype = POLTYPE_ILVPAdPA;
     datafile->NrPols = 5;
   }else {
     datafile->poltype = POLTYPE_ILVPAdPATEldEl;
     datafile->NrPols = 8;
   }
 }else {
   word_ptr = pickWordFromString(txt, 1, &nrwords, 1, ' ', verbose);
   if(nrwords != 3) {
     fflush(stdout);
     printerror(verbose.debug, "ERROR readPPOLHeader: Line should have 3 words, got %d", nrwords);
     if(nrwords == 10)
       printerror(verbose.debug, "                             Maybe file is in format %s?", returnFileFormat_str(PPOL_format));
     printerror(verbose.debug, "                             Line: '%s'.", txt);
     free(txt);
     return 0;
   }
 }
 ret = sscanf(word_ptr, "%f", &dummy_float);
 if(ret != 1) {
   fflush(stdout);
   printerror(verbose.debug, "ERROR readPPOLHeader: Cannot interpret as a float: '%s'.", txt);
   free(txt);
   return 0;
 }
 if(dummy_float >= 360) {
   fflush(stdout);
   printwarning(verbose.debug, "WARNING: IGNORING POINTS AT PULSE LONGITUDES > 360 deg.");
 }else {
   (datafile->NrBins)++;
 }
      }
    }
  }while(ret_ptr != NULL && dummy_float < 360);
  if(extended == 0) {
    datafile->poltype = POLTYPE_PAdPA;
    datafile->NrPols = 2;
  }
  fflush(stdout);
  if(verbose.verbose) fprintf(stdout, "Going to load %ld points from %s\n", datafile->NrBins, datafile->filename);
  if(datafile->NrBins == 0) {
    fflush(stdout);
    printerror(verbose.debug, "ERROR readPPOLHeader: No data in %s", datafile->filename);
    free(txt);
    return 0;
  }
  fseek(datafile->fptr, datafile->datastart, SEEK_SET);
  free(txt);
  if(datafile->offpulse_rms != NULL) {
    free(datafile->offpulse_rms);
    datafile->offpulse_rms = NULL;
  }
  if(extended) {
    datafile->offpulse_rms = (float *)malloc(datafile->NrSubints*datafile->NrFreqChan*datafile->NrPols*sizeof(float));
    if(datafile->offpulse_rms == NULL) {
      fflush(stdout);
      printerror(verbose.debug, "ERROR readPPOLHeader: Memory allocation error");
      return 0;
    }
  }
  datafile->tsamp_list = (double *)malloc(datafile->NrBins*sizeof(double));
  if(datafile->tsamp_list == NULL) {
    fflush(stdout);
    printerror(verbose.debug, "ERROR readPPOLHeader: Memory allocation error");
    return 0;
  }
  return 1;
}
int readPPOLfile(datafile_definition *datafile, float *data, int extended, float add_longitude_shift, verbose_definition verbose)
{
  int maxlinelength;
  long i, k, dummy_long;
  char *txt, *ret_ptr;
  if(datafile->NrBins == 0) {
    fflush(stdout);
    printerror(verbose.debug, "ERROR readPPOLfile: No data in %s", datafile->filename);
    return 0;
  }
  maxlinelength = 2000;
  txt = malloc(maxlinelength);
  if(txt == NULL) {
    fflush(stdout);
    printerror(verbose.debug, "ERROR readPPOLfile: Memory allocation error.");
    return 0;
  }
  fseek(datafile->fptr, datafile->datastart, SEEK_SET);
  k = 0;
  if(extended) {
    datafile->offpulse_rms[3] = -1;
    datafile->offpulse_rms[4] = -1;
    if(datafile->NrPols == 8) {
      datafile->offpulse_rms[6] = -1;
      datafile->offpulse_rms[7] = -1;
    }
  }
  for(i = 0; i < datafile->NrBins; i++) {
    ret_ptr = fgets(txt, maxlinelength, datafile->fptr);
    if(ret_ptr == NULL) {
      fflush(stdout);
      printerror(verbose.debug, "ERROR readPPOLfile: Cannot read next line, should not happen after successfully reading in header");
      free(txt);
      return 0;
    }
    if(txt[0] != '#') {
      if(extended == 0) {
 sscanf(txt, "%lf %f %f", &(datafile->tsamp_list[k]), &(data[k]), &(data[k+datafile->NrBins]));
      }else {
 if(datafile->NrPols == 8) {
   sscanf(txt, "%ld %lf %f %f %f %f %f %f %f %f %f %f %f %f", &dummy_long, &(datafile->tsamp_list[k]), &(data[k]), &(datafile->offpulse_rms[0]), &(data[k+datafile->NrBins]), &(datafile->offpulse_rms[1]), &(data[k+2*datafile->NrBins]), &(datafile->offpulse_rms[2]), &(data[k+3*datafile->NrBins]), &(data[k+4*datafile->NrBins]), &(data[k+5*datafile->NrBins]), &(datafile->offpulse_rms[5]), &(data[k+6*datafile->NrBins]), &(data[k+7*datafile->NrBins]));
 }else {
   sscanf(txt, "%ld %lf %f %f %f %f %f %f %f %f", &dummy_long, &(datafile->tsamp_list[k]), &(data[k]), &(datafile->offpulse_rms[0]), &(data[k+datafile->NrBins]), &(datafile->offpulse_rms[1]), &(data[k+2*datafile->NrBins]), &(datafile->offpulse_rms[2]), &(data[k+3*datafile->NrBins]), &(data[k+4*datafile->NrBins]));
 }
      }
      datafile->tsamp_list[k] += add_longitude_shift;
      if(datafile->tsamp_list[k] >= 0 && datafile->tsamp_list[k] < 360) {
   k++;
      }else {
 fflush(stdout);
 printwarning(verbose.debug, "WARNING readPPOLfile: IGNORING POINTS AT PULSE LONGITUDES > 360 deg.");
      }
    }
  }
  if(k != datafile->NrBins) {
    fflush(stdout);
    printerror(verbose.debug, "WARNING readPPOLfile: The nr of bins read in is different as determined from header. Something is wrong.");
    return 0;
  }
  fflush(stdout);
  if(verbose.verbose) fprintf(stdout, "readPPOLfile: Accepted %ld points\n", datafile->NrBins);
     free(txt);
  return 1;
}
int writePPOLfile(datafile_definition datafile, float *data, int extended, int onlysignificantPA, int twoprofiles, float PAoffset, verbose_definition verbose)
{
  long j;
  if(datafile.poltype != POLTYPE_ILVPAdPA && datafile.poltype != POLTYPE_PAdPA && datafile.poltype != POLTYPE_ILVPAdPATEldEl) {
    printerror(verbose.debug, "ERROR writePPOLfile: Data doesn't appear to have poltype ILVPAdPA, PAdPA or ILVPAdPATEldEl (it is %d).", datafile.poltype);
    return 0;
  }
  if(datafile.poltype == POLTYPE_ILVPAdPA && datafile.NrPols != 5) {
    printerror(verbose.debug, "ERROR writePPOLfile: 5 polarization channels were expected, but there are %ld.", datafile.NrPols);
    return 0;
  }else if(datafile.poltype == POLTYPE_PAdPA && datafile.NrPols != 2) {
    printerror(verbose.debug, "ERROR writePPOLfile: 2 polarization channels were expected, but there are %ld.", datafile.NrPols);
    return 0;
  }else if(datafile.poltype == POLTYPE_ILVPAdPATEldEl && datafile.NrPols != 8) {
    printerror(verbose.debug, "ERROR writePPOLfile: 8 polarization channels were expected, but there are %ld.", datafile.NrPols);
    return 0;
  }
  if(datafile.NrSubints > 1 || datafile.NrFreqChan > 1) {
    printerror(verbose.debug, "ERROR writePPOLfile: Can only do this opperation if there is one subint and one frequency channel.");
    return 0;
  }
  if(datafile.tsampMode != TSAMPMODE_LONGITUDELIST) {
    printerror(verbose.debug, "ERROR writePPOLfile: Expected pulse longitudes to be defined.");
    return 0;
  }
  int pa_offset, dpa_offset;
  if(datafile.poltype == POLTYPE_ILVPAdPA || datafile.poltype == POLTYPE_ILVPAdPATEldEl) {
    pa_offset = 3;
    dpa_offset = 4;
  }else if(datafile.poltype == POLTYPE_PAdPA) {
    pa_offset = 0;
    dpa_offset = 1;
  }
  for(j = 0; j < datafile.NrBins; j++) {
    if(data[j+dpa_offset*datafile.NrBins] > 0 || onlysignificantPA == 0) {
      if(extended) {
 fprintf(datafile.fptr, "%ld %e %e %e %e %e %e %e %e %e", j, datafile.tsamp_list[j], data[j], datafile.offpulse_rms[0], data[j+datafile.NrBins], datafile.offpulse_rms[1], data[j+2*datafile.NrBins], datafile.offpulse_rms[2], data[j+pa_offset*datafile.NrBins]+PAoffset, data[j+dpa_offset*datafile.NrBins]);
 if(datafile.poltype == POLTYPE_ILVPAdPATEldEl)
   fprintf(datafile.fptr, " %e %e %e %e", data[j+5*datafile.NrBins], datafile.offpulse_rms[5], data[j+6*datafile.NrBins], data[j+7*datafile.NrBins]);
 fprintf(datafile.fptr, "\n");
      }else {
 fprintf(datafile.fptr, "%e %e %e\n", datafile.tsamp_list[j], data[j+pa_offset*datafile.NrBins]+PAoffset, data[j+dpa_offset*datafile.NrBins]);
      }
    }
  }
  if(twoprofiles) {
    for(j = 0; j < datafile.NrBins; j++) {
      if(data[j+dpa_offset*datafile.NrBins] > 0 || onlysignificantPA == 0) {
 if(extended) {
   fprintf(datafile.fptr, "%ld %e %e %e %e %e %e %e %e %e", j, datafile.tsamp_list[j]+360, data[j], datafile.offpulse_rms[0], data[j+datafile.NrBins], datafile.offpulse_rms[1], data[j+2*datafile.NrBins], datafile.offpulse_rms[2], data[j+pa_offset*datafile.NrBins]+PAoffset, data[j+dpa_offset*datafile.NrBins]);
   if(datafile.poltype == POLTYPE_ILVPAdPATEldEl)
     fprintf(datafile.fptr, " %e %e %e %e", data[j+5*datafile.NrBins], datafile.offpulse_rms[5], data[j+6*datafile.NrBins], data[j+7*datafile.NrBins]);
   fprintf(datafile.fptr, "\n");
 }else {
   fprintf(datafile.fptr, "%e %e %e\n", datafile.tsamp_list[j]+360, data[j]+PAoffset, data[j+dpa_offset*datafile.NrBins]);
 }
      }
    }
  }
  return 1;
}

max_stars_repo_name: David-McKenna/psrsalsa

max_issues_repo_name: David-McKenna/psrsalsa

max_forks_repo_name: David-McKenna/psrsalsa


text:
#include <getopt.h>
#include <stdio.h>
#include <list>
#include <gsl/gsl_rng.h>
//#include "utils.h"
#include "data.h"
#include "model.h"

// helper function to write out per-user info
void log_user(FILE* file, Data *data, int user, int heldout, double rmse, double mae,
    double rank, int first, double crr, double ncrr, double ndcg, bool stats);

void log_item(FILE* file, Data *data, int item, int heldout, double rmse, double mae,
    double rank, int first, double crr, double ncrr, double ndcg, bool stats);

// helper function to sort predictions properly
bool prediction_compare(const pair<double,int>& itemA,
    const pair<double, int>& itemB);

// take a prediction function as an argument
void eval(Model* model, double (Model::*prediction)(int,int), string outdir, Data* data, bool stats, unsigned long int seed, bool verbose, string label, bool write_rankings, bool mapped_ids);

max_stars_repo_name: ajbc/spf

max_issues_repo_name: ajbc/spf

max_forks_repo_name: ajbc/spf