text:
/**
 * \author Sylvain Marsat, University of Maryland - NASA GSFC
 *
 * \brief C header for functions windowing and computing FFT/IFFT of time/frequency series.
 *
 */

#ifndef __GENERATELLVTD_H__
#define __GENERATELLVTD_H__ 1

#ifdef __GNUC__
#define UNUSED __attribute__ ((unused))
#else
#define UNUSED
#endif

#define _XOPEN_SOURCE 500

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <complex.h>
#include <time.h>
#include <unistd.h>
#include <getopt.h>
#include <stdbool.h>
#include <string.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_bspline.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_min.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_complex.h>

#include "constants.h"
#include "struct.h"
#include "waveform.h"
#include "LISAgeometry.h"

/********************************** Structures ******************************************/

typedef enum GenWavetag {
  hphcTD,
  hphcFD,
  hlm
} GenWavetag;

/* Parameters for the generation of a LISA waveform (in the form of a list of modes) */
typedef struct tagGenLLVTDparams {
  double lambda;             /* first angle for the position in the sky (rad, default 0) */
  double beta;               /* second angle for the position in the sky (rad, default 0) */
  double polarization;       /* polarization angle (rad, default 0) */
  int tagnetwork;            /* Tag selecting the desired detector network */
  int nsamplesinfile;        /* Number of lines of input file */
  int binaryin;              /* Tag for loading the data in gsl binary form instead of text (default false) */
  int binaryout;             /* Tag for outputting the data in gsl binary form instead of text (default false) */
  char indir[256];           /* Path for the input directory */
  char infile[256];          /* Path for the input file */
  char outdir[256];          /* Path for the output directory */
  char outfile[256];         /* Path for the output file */
} GenLLVTDparams;

#if 0
{ /* so that editors will match succeeding brace */
#elif defined(__cplusplus)
}
#endif

#endif /* _GENERATELLVTD_H */

max_stars_repo_name: titodalcanton/flare

max_issues_repo_name: titodalcanton/flare

max_forks_repo_name: titodalcanton/flare


text:
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

int main(void)
{
  int i; double z; gsl_rng *r;
  r = gsl_rng_alloc(gsl_rng_mt19937);
  gsl_rng_set(r,0);
  for (i=0;i<10;i++) {
    z = gsl_ran_gaussian(r,1.0);
    printf("z(%d) = %f\n",i,z);
  }
  exit(EXIT_SUCCESS);
}

ext: c


text:
/* linalg/choleskyc.c
 *
 * Copyright (C) 2007, 2019 Patrick Alken
 * Copyright (C) 2010 Huan Wu
 *
 * 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 <config.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_errno.h>

#include "recurse.h"

/*
 * This module contains routines related to the Cholesky decomposition
 * of a complex Hermitian positive definite matrix.
 */

static void cholesky_complex_conj_vector(gsl_vector_complex *v);
static int complex_cholesky_decomp_L2(gsl_matrix_complex * A);
static int complex_cholesky_decomp_L3(gsl_matrix_complex * A);

/*
gsl_linalg_complex_cholesky_decomp()
  Perform the Cholesky decomposition on a Hermitian positive definite
matrix using Level 3 BLAS.

Inputs: A - (input/output) complex positive definite matrix

Return: success or error

Notes:
1) The lower triangle of A is overwritten with the Cholesky decomposition
*/

int
gsl_linalg_complex_cholesky_decomp(gsl_matrix_complex *A)
{
  const size_t N = A->size1;

  if (N != A->size2)
    {
      GSL_ERROR("Cholesky decomposition requires square matrix", GSL_ENOTSQR);
    }
  else
    {
      return complex_cholesky_decomp_L3(A);
    }
}

/*
gsl_linalg_complex_cholesky_solve()
  Solve A x = b where A is in cholesky form
*/

int
gsl_linalg_complex_cholesky_solve (const gsl_matrix_complex * cholesky,
                                   const gsl_vector_complex * b,
                                   gsl_vector_complex * x)
{
  if (cholesky->size1 != cholesky->size2)
    {
      GSL_ERROR ("cholesky matrix must be square", GSL_ENOTSQR);
    }
  else if (cholesky->size1 != b->size)
    {
      GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
    }
  else if (cholesky->size2 != x->size)
    {
      GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
    }
  else
    {
      gsl_vector_complex_memcpy (x, b);
      return gsl_linalg_complex_cholesky_svx(cholesky, x);
    }
}

/*
gsl_linalg_complex_cholesky_svx()
  Solve A x = b in place where A is in cholesky form
*/

int
gsl_linalg_complex_cholesky_svx (const gsl_matrix_complex * cholesky,
                                 gsl_vector_complex * x)
{
  if (cholesky->size1 != cholesky->size2)
    {
      GSL_ERROR ("cholesky matrix must be square", GSL_ENOTSQR);
    }
  else if (cholesky->size2 != x->size)
    {
      GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
    }
  else
    {
      /* solve for y using forward-substitution, L y = b */
      gsl_blas_ztrsv (CblasLower, CblasNoTrans, CblasNonUnit, cholesky, x);

      /* perform back-substitution, L^H x = y */
      gsl_blas_ztrsv (CblasLower, CblasConjTrans, CblasNonUnit, cholesky, x);

      return GSL_SUCCESS;
    }
}

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

gsl_linalg_complex_cholesky_invert()
  Compute the inverse of an Hermitian positive definite matrix in
  Cholesky form.

Inputs: LLT - matrix in cholesky form on input
              A^{-1} = L^{-H} L^{-1} on output

Return: success or error
******************************************************************************/

int
gsl_linalg_complex_cholesky_invert(gsl_matrix_complex * LLT)
{
  if (LLT->size1 != LLT->size2)
    {
      GSL_ERROR ("cholesky matrix must be square", GSL_ENOTSQR);
    }
  else
    {
      int status;
      size_t N = LLT->size1;
      size_t i, j;

      /* invert the lower triangle of LLT */
      status = gsl_linalg_complex_tri_invert(CblasLower, CblasNonUnit, LLT);
      if (status)
        return status;

      /* compute A^{-1} = L^{-H} L^{-1} */
      status = gsl_linalg_complex_tri_LHL(LLT);
      if (status)
        return status;

      /* copy the Hermitian lower triangle to the upper triangle */
      for (i = 1; i < N; ++i)
        {
          for (j = 0; j < i; ++j)
            {
              gsl_complex z = gsl_matrix_complex_get(LLT, i, j);
              gsl_matrix_complex_set(LLT, j, i, gsl_complex_conjugate(z));
            }
        }

      return GSL_SUCCESS;
    }
}

/********************************************
 *           INTERNAL ROUTINES              *
 ********************************************/

static void
cholesky_complex_conj_vector(gsl_vector_complex *v)
{
  size_t i;

  for (i = 0; i < v->size; ++i)
    {
      gsl_complex * vi = gsl_vector_complex_ptr(v, i);
      GSL_IMAG(*vi) = -GSL_IMAG(*vi);
    }
}

/*
complex_cholesky_decomp_L2()
  Perform the Cholesky decomposition on a Hermitian positive definite
matrix using Level 2 BLAS. See Golub & Van Loan, "Matrix Computations" (3rd ed),
algorithm 4.2.2.

Inputs: A - (input/output) complex postive definite matrix

Return: success or error

Notes:
1) The lower triangle of A is overwritten with the Cholesky decomposition
*/

static int
complex_cholesky_decomp_L2(gsl_matrix_complex * A)
{
  const size_t N = A->size1;

  if (N != A->size2)
    {
      GSL_ERROR("Cholesky decomposition requires square matrix", GSL_ENOTSQR);
    }
  else
    {
      size_t j;
      gsl_complex z;
      double ajj;

      for (j = 0; j < N; ++j)
        {
          z = gsl_matrix_complex_get(A, j, j);
          ajj = GSL_REAL(z);

          if (j > 0)
            {
              gsl_vector_complex_const_view aj =
                gsl_matrix_complex_const_subrow(A, j, 0, j);

              gsl_blas_zdotc(&aj.vector, &aj.vector, &z);
              ajj -= GSL_REAL(z);
            }

          if (ajj <= 0.0)
            {
              GSL_ERROR("matrix is not positive definite", GSL_EDOM);
            }

          ajj = sqrt(ajj);
          GSL_SET_COMPLEX(&z, ajj, 0.0);
          gsl_matrix_complex_set(A, j, j, z);

          if (j < N - 1)
            {
              gsl_vector_complex_view av =
                gsl_matrix_complex_subcolumn(A, j, j + 1, N - j - 1);

              if (j > 0)
                {
                  gsl_vector_complex_view aj =
                    gsl_matrix_complex_subrow(A, j, 0, j);
                  gsl_matrix_complex_view am =
                    gsl_matrix_complex_submatrix(A, j + 1, 0, N - j - 1, j);

                  cholesky_complex_conj_vector(&aj.vector);

                  gsl_blas_zgemv(CblasNoTrans,
                                 GSL_COMPLEX_NEGONE,
                                 &am.matrix,
                                 &aj.vector,
                                 GSL_COMPLEX_ONE,
                                 &av.vector);

                  cholesky_complex_conj_vector(&aj.vector);
                }

              gsl_blas_zdscal(1.0 / ajj, &av.vector);
            }
        }

      return GSL_SUCCESS;
    }
}

/*
complex_cholesky_decomp_L3()
  Perform the Cholesky decomposition on a Hermitian positive definite
matrix using Level 3 BLAS.

Inputs: A - (input/output) complex postive definite matrix

Return: success or error

Notes:
1) The lower triangle of A is overwritten with the Cholesky decomposition

2) Based on ReLAPACK recursive variant with Level 3 BLAS
*/

static int
complex_cholesky_decomp_L3(gsl_matrix_complex * A)
{
  const size_t N = A->size1;

  if (N != A->size2)
    {
      GSL_ERROR("Cholesky decomposition requires square matrix", GSL_ENOTSQR);
    }
  else if (N <= CROSSOVER_CHOLESKY)
    {
      /* use unblocked Level 2 algorithm */
      return complex_cholesky_decomp_L2(A);
    }
  else
    {
      /*
       * partition matrix:
       *
       * A11 A12
       * A21 A22
       *
       * where A11 is N1-by-N1
       */
      int status;
      const size_t N1 = GSL_LINALG_SPLIT_COMPLEX(N);
      const size_t N2 = N - N1;
      gsl_matrix_complex_view A11 = gsl_matrix_complex_submatrix(A, 0, 0, N1, N1);
      gsl_matrix_complex_view A21 = gsl_matrix_complex_submatrix(A, N1, 0, N2, N1);
      gsl_matrix_complex_view A22 = gsl_matrix_complex_submatrix(A, N1, N1, N2, N2);

      /* recursion on A11 */
      status = complex_cholesky_decomp_L3(&A11.matrix);
      if (status)
        return status;

      /* A21 = A21 * A11^{-1} */
      gsl_blas_ztrsm(CblasRight, CblasLower, CblasConjTrans, CblasNonUnit, GSL_COMPLEX_ONE, &A11.matrix, &A21.matrix);

      /* A22 -= A21 A21^H */
      gsl_blas_zherk(CblasLower, CblasNoTrans, -1.0, &A21.matrix, 1.0, &A22.matrix);

      /* recursion on A22 */
      status = complex_cholesky_decomp_L3(&A22.matrix);
      if (status)
        return status;

      return GSL_SUCCESS;
    }
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_eigen.h>

#include "egsl.h"

val egsl_sub(val v1,val v2){
	return egsl_sum(v1, egsl_scale(-1.0,v2));
}

val egsl_compose_col(val v1, val v2){
	gsl_matrix *m1 = egsl_gslm(v1);
	gsl_matrix *m2 = egsl_gslm(v2);
	egsl_expect_size(v2, 0, m1->size2);
	val v3 = egsl_alloc(m1->size1+m2->size1,m1->size2);
	gsl_matrix *m3 = egsl_gslm(v3);
	size_t i,j;
	for(j=0;j<m1->size2;j++) {
		for(i=0;i<m1->size1;i++)
			gsl_matrix_set(m3, i, j, gsl_matrix_get(m1,i,j));

		for(i=0;i<m2->size1;i++)
			gsl_matrix_set(m3, m1->size1+i, j, gsl_matrix_get(m2,i,j));
	}
	return v3;
}

val egsl_compose_row(val v1, val v2){
	gsl_matrix *m1 = egsl_gslm(v1);
	gsl_matrix *m2 = egsl_gslm(v2);
	egsl_expect_size(v2, m1->size1, 0);
	val v3 = egsl_alloc(m1->size1, m1->size2 + m2->size2);
	gsl_matrix *m3 = egsl_gslm(v3);
	size_t i,j;
	for(i=0;i<m1->size1;i++) {
		for(j=0;j<m1->size2;j++)
			gsl_matrix_set(m3, i, j, gsl_matrix_get(m1,i,j));

		for(j=0;j<m2->size2;j++)
			gsl_matrix_set(m3, i, m1->size2+j, gsl_matrix_get(m2,i,j));
	}
	return v3;
}

void egsl_add_to(val v1, val v2) {
	gsl_matrix * m1 = egsl_gslm(v1);
	gsl_matrix * m2 = egsl_gslm(v2);
	gsl_matrix_add(m1,m2);
}

void egsl_add_to_col(val v1, size_t j, val v2) {
/*	egsl_print("m1",v1);
	egsl_print("m2",v2); */
	gsl_matrix * m1 = egsl_gslm(v1);
	gsl_matrix * m2 = egsl_gslm(v2);

/*	printf("m1 size = %d,%d j = %d\n",m1->size1,m1->size2,j); */
	egsl_expect_size(v2, m1->size1, 1);
	size_t i;
	for(i=0;i<m1->size1;i++) {
		*gsl_matrix_ptr(m1, i, j) += gsl_matrix_get(m2,i,0);
	}
}

val egsl_copy_val(val v1) {
	gsl_matrix * m1 = egsl_gslm(v1);
	val v2 = egsl_alloc(m1->size1,m1->size2);
	gsl_matrix * m2 = egsl_gslm(v2);
	gsl_matrix_memcpy(m2,m1);
	return v2;
}

val egsl_scale(double s, val v1){
	val v2 = egsl_copy_val(v1);
	gsl_matrix * m2 = egsl_gslm(v2);
	gsl_matrix_scale(m2, s);
	return v2;
}

val egsl_sum(val v1, val v2){
	gsl_matrix * m1 = egsl_gslm(v1);
	gsl_matrix * m2 = egsl_gslm(v2);
	val v3 = egsl_alloc(m1->size1,m1->size2);
	gsl_matrix * m3 = egsl_gslm(v3);
	gsl_matrix_memcpy(m3,m1);
	gsl_matrix_add(m3,m2);
	return v3;
}

val egsl_sum3(val v1, val v2, val v3){
	return egsl_sum(v1, egsl_sum(v2,v3));
}

val egsl_mult(val v1, val v2){
	gsl_matrix * a = egsl_gslm(v1);
	gsl_matrix * b = egsl_gslm(v2);
	val v = egsl_alloc(a->size1,b->size2);
	gsl_matrix * ab = egsl_gslm(v);
	gsl_blas_dgemm(CblasNoTrans,CblasNoTrans,1.0,a,b,0.0,ab);
	return v;
}

val egsl_transpose(val v1){
	gsl_matrix * m1 = egsl_gslm(v1);
	val v2 = egsl_alloc(m1->size2,m1->size1);
	gsl_matrix * m2 = egsl_gslm(v2);
	gsl_matrix_transpose_memcpy(m2,m1);
	return v2;
}

val egsl_inverse(val v1){
	gsl_matrix*A = egsl_gslm(v1);
	val v2 = egsl_alloc(A->size1,A->size1);
	gsl_matrix*invA = egsl_gslm(v2);
	size_t n = A->size1;
	gsl_matrix * m = gsl_matrix_alloc(n,n);
	gsl_matrix_memcpy(m,A);
	gsl_permutation * perm = gsl_permutation_alloc (n);
	/* Make LU decomposition of matrix m */
	int s;
	gsl_linalg_LU_decomp (m, perm, &s);
	/* Invert the matrix m */
	gsl_linalg_LU_invert (m, perm, invA);
	gsl_permutation_free(perm);
	gsl_matrix_free(m);
	return v2;
}

void egsl_symm_eig(val v, double* eigenvalues, val* eigenvectors) {
	gsl_matrix *m = egsl_gslm(v);
	size_t N = m->size1;
	/* Check for v to be square */

	gsl_matrix *A = gsl_matrix_alloc(N,N);
	gsl_matrix_memcpy(A, m);

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

	gsl_eigen_symmv_workspace * ws = gsl_eigen_symmv_alloc(N);
	gsl_eigen_symmv(A, eval, evec, ws);
	gsl_eigen_symmv_free(ws);

	gsl_eigen_symmv_sort(eval, evec, GSL_EIGEN_SORT_VAL_DESC);

	size_t j;
	for(j=0;j<N;j++) {
		eigenvalues[j] = gsl_vector_get(eval, j);
		eigenvectors[j] = egsl_alloc(N,1);
		size_t i;
		for(i=0;i<N;i++)
			*egsl_atmp(eigenvectors[j],i,0) = gsl_matrix_get(evec,i,j);
	}

	gsl_vector_free(eval);
	gsl_matrix_free(evec);
	gsl_matrix_free(A);
}


text:
/**
 * File: trfit_utils.c
 *
 */
#include "trfit_utils.h"
#include <stdio.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_sort_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_statistics.h>
#include "aXe_grism.h"
#include "aXe_utils.h"
#include "spce_sect.h"
#include "spce_fitting.h"
#include "spce_is_in.h"
#include "spc_back.h"
#include "spce_pathlength.h"
#include "lmmin.h"
#include "lm_eval.h"

#define MIN(x,y) (((x)<(y))?(x):(y))
#define MAX(x,y) (((x)>(y))?(x):(y))

int gauss_f(const gsl_vector * x, void *params, gsl_vector * f)
{
  int ndata      = ((struct function_data *)params)->n;
  double *x_data = ((struct function_data *)params)->x;
  double *y_data = ((struct function_data *)params)->y;
  //double *e_data = ((struct function_data *)params)->sig;

  double a  = gsl_vector_get (x, 0);
  double b  = gsl_vector_get (x, 1);
  double c  = gsl_vector_get (x, 2);
  double x0 = gsl_vector_get (x, 3);

  double Yi = 0.0;

  int i;

  for (i = 0; i < ndata; i++)
    {
      /* Model Yi = a * exp(-b*(Xi-X0) + c */
      Yi = a * exp (-b * (x_data[i]-x0)*(x_data[i]-x0)) + c;
      //      gsl_vector_set (f, i, (Yi - y_data[i])/e_data[i]);
      gsl_vector_set (f, i, Yi - y_data[i]);
      //      fprintf(stderr, "Yi: %e, f: %e\n", Yi, (Yi - y_data[i])/e_data[i]);
    }

  return GSL_SUCCESS;
}

int
gauss_df(const gsl_vector * x, void *params,
	 gsl_matrix * J)
{
  int ndata      = ((struct function_data *)params)->n;
  double *x_data = ((struct function_data *)params)->x;
  //double *y_data = ((struct function_data *)params)->y;
  //double *e_data = ((struct function_data *)params)->sig;

  double a  = gsl_vector_get (x, 0);
  double b  = gsl_vector_get (x, 1);
  //double c  = gsl_vector_get (x, 2);
  double x0 = gsl_vector_get (x, 3);

  int i;
  double value;
  double xdiff;

  for (i = 0; i < ndata; i++)
    {
      /* Jacobian matrix J(i,j) = dfi / dxj, */
      /* where fi = (Yi - yi)/sigma[i],      */
      /* Model Yi = a * exp(-b*(Xi-X0) + c   */
      /* and the xj are the parameters (a,b,c,x0) */

      xdiff = x_data[i]-x0;

      value = exp(-b*(xdiff)*(xdiff));

      gsl_matrix_set (J, i, 0, value);
      gsl_matrix_set (J, i, 1, -a*(xdiff)*(xdiff)*value);
      gsl_matrix_set (J, i, 2, 1.0);
      gsl_matrix_set (J, i, 3,a*b*value*2*xdiff);
    }
  return GSL_SUCCESS;
}

int
gauss_fdf (const gsl_vector * x, void *params, gsl_vector *f, gsl_matrix *J)
{
  gauss_f(x, params, f);
  gauss_df(x, params, J);

  return GSL_SUCCESS;
}

/**
 * Function: fit_beamtrace
 * The function determines a new trace solution for a given beam. A gaussian
 * is fitted to the pixel values around the nominal trace position.
 * A kappa-sigma-klipping algorithm eliminates fitted trace centers
 * which are off due to cosmics or similar.
 * Finally the a new trace solution is fitted to the set of trace centers
 * The new slope is directly written to the beam structure, the
 * new y-value of the reference point is returned.
 *
 * Parameters:
 * @param conf     - the aXe configuration structure
 * @param obs      - the current observation
 * @param beamID   - the ID of the current beam
 * @param act_beam - the beam to fit
 *
 * Returns:
 * @return yf - the y-value of the reference point for the fitted trace
 */

gsl_vector *
fit_beamtrace(const aperture_conf  *conf, observation *obs,
	      const int beamID, beam act_beam)
{
  //gsl_vector_int *yvec;
  gsl_vector     *fit_params;
  gsl_vector     *lin_fit;
  gsl_vector     *fit_result;
  trace_func     *tracefun;

  px_point        xborder;
  px_point        tpoint;

  fit_data       *f_data;

  double *x;
  double *y;
  double *w;

  double *trace;

  double yf, yf_err;

  int i;
  int index;

  // This is the return result
  fit_result = gsl_vector_alloc(5);
  gsl_vector_set(fit_result, 0, 0);
  gsl_vector_set(fit_result, 1, 0);
  gsl_vector_set(fit_result, 2, 0);
  gsl_vector_set(fit_result, 3, 0);
  gsl_vector_set(fit_result, 4, 0);

  // define the beam and the trace function
  tracefun = act_beam.spec_trace;
  trace = (double *)(tracefun->data);

  // If this beam's ignore flag is set to 1 then do nothing
  if (act_beam.ignore == 1)
    return fit_result;

  // determine the start and end point in x
  //xborder = get_fit_xrange(conf, obs, beamID, act_beam);
  xborder = get_fit_xrange(conf, obs, act_beam);

  // allocate arrays for the fitting
  x = (double*)malloc((xborder.y-xborder.x+1)*sizeof(double));
  y = (double*)malloc((xborder.y-xborder.x+1)*sizeof(double));
  w = (double*)malloc((xborder.y-xborder.x+1)*sizeof(double));

  // Loop over all columns
  index = 0;
  for (i = xborder.x; i < xborder.y; i++)
    //for (i=130; i < 150; i++)
    {
      // determine the pixel which is
      // closest to the nominal trace position
      tpoint.x = i;
      tpoint.y = (int)floor(tracefun->func((double)i-act_beam.refpoint.x,
					   tracefun->data)
			    + act_beam.refpoint.y+0.5);

      // get the data to be fitted
      // in the current column
      f_data = get_fitdata(obs, N_COLUMNS_FIT, tpoint);
      /*
      if (tpoint.x > 130 && tpoint.x < 150)
	{
	  fprintf(stdout, "x: %i\n", tpoint.x);
	  print_fit_data(f_data);
	  fprintf(stdout, "\n\n");
	}
      if (tpoint.x == 135)
      {*/
      //fprintf(stdout, "x: %i\n", tpoint.x);
      //print_fit_data(f_data);
      //fprintf(stdout, "\n\n");
	  // fit the gaussian to the data
	  fit_params = fit_wuttke(f_data);
	  /*	}
      else{
	fit_params = fit_wuttke_talk(f_data);
	}*/
      // check whether the fit was successful
      if (gsl_vector_get(fit_params, 4) == 1)
	{
	  // transfer the results of the fit
	  // to this column in the vector
	  // which contains all fit results
	  x[index] = (double)tpoint.x;
	  y[index] = gsl_vector_get(fit_params, 3);
	  w[index] = 1.0;
	  index++;
	}

      // freem memory
      free_fit_data(f_data);
      gsl_vector_free(fit_params);
    }

  //for (i=0; i < index; i++)
  //  fprintf(stdout, "%f %f %f\n", x[i], y[i], w[i]);

  // do a kappa-sigma klipping on
  // the array with the trace locations
  // to improve the robustness of the result
  comp_kappasigma_interp(x, y, w, index, 1, N_KAPPASIG_ITER,
			 N_KAPPASIG_SIG, obs, 0);

  // make a linear fit to the trace positions
  lin_fit = det_vector_linear(x, y, w, index, 0);
  //for (i=0; i < index; i++)
  //fprintf(stdout, "%f %f %f\n", x[i], y[i], w[i]);

  // compute the y-position of the beam reference position
  // according to the new trace solution
  gsl_fit_linear_est (act_beam.refpoint.x-gsl_vector_get(lin_fit, 0), gsl_vector_get(lin_fit, 1),
		      gsl_vector_get(lin_fit, 2), gsl_vector_get(lin_fit, 3),
		      gsl_vector_get(lin_fit, 4), gsl_vector_get(lin_fit, 5),
		      &yf, &yf_err);

  gsl_vector_set(fit_result, 0, yf);
  gsl_vector_set(fit_result, 1, yf_err);
  gsl_vector_set(fit_result, 2, gsl_vector_get(lin_fit, 2));
  gsl_vector_set(fit_result, 3, sqrt(gsl_vector_get(lin_fit, 5)));
  gsl_vector_set(fit_result, 4, gsl_vector_get(lin_fit, 6));

  fprintf(stdout, "\nFitting result XXX: xpos = %f , ypos = %f +- %f , slope = %f +- %f , chi^2 = %f, x_pivot = %f, y_abs = %f\n",
	  act_beam.refpoint.x, gsl_vector_get(fit_result, 0), gsl_vector_get(fit_result, 1),
	  gsl_vector_get(fit_result, 2), gsl_vector_get(fit_result, 3), gsl_vector_get(lin_fit, 6),
	  gsl_vector_get(lin_fit, 0), gsl_vector_get(lin_fit, 1));

  // copy the slope of the new trace solution
  // to the according position in the beam description
  trace[1] = gsl_vector_get(lin_fit, 2);
  trace[2] = sqrt(gsl_vector_get(lin_fit, 5));

  // free the memory
  gsl_vector_free(lin_fit);
  free(x);
  free(y);
  free(w);

  // return the new y-position of
  // the reference point
  //  return yf;
  return fit_result;
}

/**
 * Function: gagauss
 * The function computes and returns the value
 * of a gauss function y = a * exp(-b*(x-c)^2) - d
 * at a position x.
 *
 * Parameters:
 * @param x      - the x-poition to evaluate the function
 * @param params - the function parameters
 *
 * Returns:
 * @return - the function value
 */
double
gagauss(const double x, const double *params)
{
  // just compute and return the result
  return (params[0] * exp (-params[1] * (x-params[3])*(x-params[3])) + params[2]);
}

/**
 * Function: comp_shift
 *
 * Returns:
 * @return - the function value
 */
double
comp_shift(const double x, const double *params, const double *fpars)
{
  double dy;
  double tmp;

  // Parameters in the original
  // notation fo W. Freudling:
  // fpars[0] = 'contrast'
  // fpars[1] = 'exp'
  // fpars[2] = 'tan_ang'
  // gpars[0] = 'shift'

  // compute the shift for the
  // x-value
  dy = x * fpars[2] + params[0];

  // compute the sin()-expression
  tmp = (sin((2.0*dy-0.5)*M_PI) + fpars[1])/2.0;

  // compute and return the whole expression
  return tmp * tmp * fpars[0] + (1.0 - fpars[0]);
}

/**
 * Function: fit_wuttke
 * Fits a gaussian to the data given in the input structure.
 * Uses a library for nn-linear fits (made by Wuttke)
 *
 * Parameters:
 * @param f_data - the fit data
 *
 * Returns:
 * @return params - the fit results
 */
gsl_vector *
fit_wuttke(const fit_data *f_data)
{
  gsl_vector *params;

  int m_dat = f_data->n_data;
  int n_p =  4;

  double p[4];

  // auxiliary settings:
  lm_control_type control;
  lm_data_type data;

  // allocate the parameter vector
  params = gsl_vector_alloc(5);

  // give a first guess
  // for the gauss parameters
  p[0] = f_data->y_values[f_data->n_data / 2];
  p[1] = 1.0;
  p[2] = (f_data->y_values[f_data->n_data-1]-f_data->y_values[0])/2.0;
  p[3] = f_data->x_values[0] + (f_data->x_values[f_data->n_data-1]-f_data->x_values[0])/2.0;

  // initialize the fit
  lm_initialize_control( &control );

  // define the parameters for the fitting
  data.user_func = gagauss;
  data.user_t = f_data->x_values;
  data.user_y = f_data->y_values;

  // perform the fit:
  lm_minimize( m_dat, n_p, p, lm_evaluate_default, lm_print_nothing,
	       &data, &control );

  //fprintf(stdout, "Result: %i, amp: %f, width: %f, const:  %f, pos: %f\n", control.info, p[0], p[1], p[2], p[3]);

  // copy the fit results to the output structure
  gsl_vector_set(params, 0, p[0]);
  gsl_vector_set(params, 1, p[1]);
  gsl_vector_set(params, 2, p[2]);
  gsl_vector_set(params, 3, p[3]);
  gsl_vector_set(params, 4, (double)control.info);

  // return the output structure
  return params;
}
gsl_vector *
fit_wuttke_talk(const fit_data *f_data)
{
  gsl_vector *params;

  int m_dat = f_data->n_data;
  int n_p =  4;

  double p[4];

  // auxiliary settings:
  lm_control_type control;
  lm_data_type data;

  // allocate the parameter vector
  params = gsl_vector_alloc(5);

  // give a first guess
  // for the gauss parameters
  p[0] = 1.0;
  p[1] = 1.0;
  p[2] = 0.0;
  p[3] = f_data->x_values[0] + (f_data->x_values[f_data->n_data-1]-f_data->x_values[0])/2.0;

  // initialize the fit
  lm_initialize_control( &control );

  // define the parameters for the fitting
  data.user_func = gagauss;
  data.user_t = f_data->x_values;
  data.user_y = f_data->y_values;

  // perform the fit:
  lm_minimize( m_dat, n_p, p, lm_evaluate_default, lm_print_default,
	       &data, &control );

  // copy the fit results to the output structure
  gsl_vector_set(params, 0, p[0]);
  gsl_vector_set(params, 1, p[1]);
  gsl_vector_set(params, 2, p[2]);
  gsl_vector_set(params, 3, p[3]);
  gsl_vector_set(params, 4, (double)control.info);

  // return the output structure
  return params;
}

d_point
find_grav_center(const fit_data *f_data)
{
  d_point ret={0.0,0.0};

  int i=0;

  double sum    = 0.0;
  double w_sum    = 0.0;
  double weight = 0.0;

  for (i=0; i < f_data->n_data; i++)
    {
      weight = fabs(f_data->y_values[i]/f_data->e_values[i]);
      sum   += weight * f_data->x_values[i];
      w_sum += weight;
    }

  ret.x = sum/w_sum;
  ret.y = weight;

  return ret;
}

/*
 * Function: get_fit_xrange
 * The function determines and returns the mininmum
 * and the maximum x-values for a beam.
 *
 * Parameters:
 * @param conf     - the aXe configuration structure
 * @param obs      - the current observation
 * @param beamID   - the ID of the current beam
 * @param act_beam - the beam to fit
 *
 * Returns:
 * @return ret - x_min and x_max of the beam
 */
px_point
get_fit_xrange(const aperture_conf  *conf, const observation *obs,
	       beam act_beam)
  {
    px_point ret;

    // compute one extremum
    ret.x = (int)floor(act_beam.refpoint.x
		       + conf->beam[act_beam.ID].offset.dx0 +.5);

    // compute the secund extremum
    ret.y = (int)floor(act_beam.refpoint.x
		       + conf->beam[act_beam.ID].offset.dx1 +.5);

    // limit the start and end column
    // to the image size
    ret.x = MAX(ret.x, 0);
    ret.y = MIN(ret.y, (int)obs->grism->size1);

    // return the min and the max value
    return ret;
  }

gsl_vector *
fit_function(const fit_data *f_data)
{
  gsl_vector *fit_params = gsl_vector_alloc(4);

  const gsl_multifit_fdfsolver_type *T;
  gsl_multifit_fdfsolver *s;
  gsl_multifit_function_fdf f;
  gsl_vector *x_init = gsl_vector_alloc(4);

  int status;
  //size_t i,
  size_t iter = 0;

  const size_t n = f_data->n_data;
  const size_t p = 4;

  //gsl_matrix *covar = gsl_matrix_alloc (p, p);

  struct function_data d = {f_data->n_data, f_data->x_values,
			    f_data->y_values, f_data->e_values};

  fit_params = gsl_vector_alloc(p);

  gsl_vector_set(x_init, 0, 1.0);
  gsl_vector_set(x_init, 1, 1.0);
  gsl_vector_set(x_init, 2, 0.0);
  gsl_vector_set(x_init, 3,
		 f_data->x_values[0]
		 +(f_data->x_values[f_data->n_data-1]-f_data->x_values[0])/2.0);

  f.f      = gauss_f;
  f.df     = gauss_df;
  f.fdf    = gauss_fdf;
  f.n      = n;
  f.p      = p;
  f.params = &d;

  T = gsl_multifit_fdfsolver_lmsder;
  s = gsl_multifit_fdfsolver_alloc (T, n, p);
  gsl_multifit_fdfsolver_set (s, &f, x_init);

  print_fit_state (iter, s);

  do
    {
      iter++;
      status = gsl_multifit_fdfsolver_iterate (s);

      printf ("status %i = %s\n", status, gsl_strerror (status));

      print_fit_state (iter, s);

      if (status != GSL_CONTINUE)
       	break;

      status = gsl_multifit_test_delta (s->dx, s->x,
                                        1e-12, 1e-12);
      printf ("status  a = %s\n", gsl_strerror (status));
      fprintf(stdout, "stat: %i, GSL_CONT: %i\n", status, GSL_CONTINUE);
    }
  while (status == GSL_CONTINUE && iter < 500);

  gsl_multifit_fdfsolver_free (s);
  return fit_params;
}

/*
 * Function: get_fitdata
 * The function collects the pixel values of at most np pixels
 * on both sides of the nominal trace position and stores
 * these values in an apppropriate structure.
 *
 * Parameters:
 * @param obs      - the current observation
 * @param np       - the number of pixels to sample around the nominal position
 * @param tr_point - the nominal traxe pixel
 *
 * Returns:
 * @return f_data  - the pixel data around the nominal trace position
 *
 */
fit_data *
get_fitdata(observation *obs, int np, px_point tr_point)
{
  gsl_vector_int *tmp;

  fit_data *f_data;

  int np_act=0;

  //int y_low;
  //int y_upp;

  int l_space=1;
  int u_space=1;

  int l_act;
  int u_act;

  int l_np=0;
  int u_np=0;

  int ncols=obs->grism->size2;

  int ii;
  int index;

  // allocate the vector
  tmp = gsl_vector_int_alloc(2*np+2);
  gsl_vector_int_set_all(tmp, -1);

  // limit the starting ppoint of the search
  // to values within the image dimension
  tr_point.y = MAX(0,tr_point.y);
  tr_point.y = MIN((int)obs->grism->size2,tr_point.y);

  // initialize the row numbers
  // to search up- and downwards
  l_act = tr_point.y - 1;
  u_act = tr_point.y;

  // as long as interpolation points are missing
  // and one direction, either up or down,
  // is 'open', continue searching
  while (np_act < 2*np && (l_space || u_space))
    {

      // check whether the direction
      // downwards is still open
      if (l_space)
	{
	  // if you are at the end of the frame
	  if (l_act < 0)
	    {
	      // close the direction downwards
	      l_space=0;
	    }
	  else
	    {
	      if (!isnan(gsl_matrix_get(obs->grism, tr_point.x, l_act)))
		{
		  // else store the interpolation point,
		  // do the various increments
		  gsl_vector_int_set(tmp, np_act, l_act);
		  //		  fprintf(stderr, "put in %i\n", l_act);
		  np_act++;
		  l_np++;
		}
	      l_act--;
	    }
	}

      // check whether the direction
      // upwards is still open
      if (u_space)
	{
	  // if you are at the end of the frame
	  if (u_act >= ncols)
	    {
	      // close the direction upwards
	      u_space=0;
	    }
	  else
	    {
	      if (!isnan(gsl_matrix_get(obs->grism, tr_point.x, u_act)))
		{
		  // else store the interpolation point,
		  // do the various increments
		  gsl_vector_int_set(tmp, np_act, u_act);
		  np_act++;
		  u_np++;
		}
	      u_act++;
	    }
	}
    }

  // sort the pixel indices
  gsl_sort_vector_int(tmp);

  // allocate a data structure
  f_data = alloc_fit_data(np_act);

  // fill the pixel values into the data structure
  index = 0;
  for (ii=((int)tmp->size - np_act); ii < (int)tmp->size; ii++)
    {
      f_data->x_values[index] = (double)gsl_vector_int_get(tmp, ii);
      f_data->y_values[index] = gsl_matrix_get(obs->grism, tr_point.x,
					       gsl_vector_int_get(tmp, ii));
      f_data->e_values[index] = gsl_matrix_get(obs->pixerrs, tr_point.x,
					       gsl_vector_int_get(tmp, ii));
      index++;
    }

  // release memory
  gsl_vector_int_free(tmp);

  //  return the result
  return f_data;
}

/*
 * Function: get_ipc_coldata
 * The function checks for valid pixel data on both side of the trace.
 * The maximum distance away from the trace is given as a parameter.
 * Valid data is collected and returned in a vector.
 *
 * Parameters:
 * @param img_data - the current observation
 * @param np       - the number of pixels to sample around the nominal position
 * @param tr_point - the nominal traxe pixel
 *
 * Returns:
 * @return ipc_row  - the pixel data around the nominal trace position
 *
 */
gsl_vector *
get_ipc_coldata(gsl_matrix *img_data, int np, px_point tr_point)
{
  gsl_vector *tmp;

  gsl_vector *ipc_row = NULL;

  int np_act=0;

  int y_act;

  int ncols = img_data->size2;

  int ii;

  // check whether the central tracepoint
  // is in the image
  if (tr_point.y < 0 || tr_point.y > ncols-1)
    // give back the NULL array
    return ipc_row;

  // allocate the vector
  tmp = gsl_vector_alloc(2*np+1);
  gsl_vector_set_all(tmp, -1.0);

  // set the current y-position
  y_act = tr_point.y;

  // check whether there is good data
  // at the current position
  if (!isnan(gsl_matrix_get(img_data, tr_point.x, y_act)))
    {
      // store the data point and enhance the counter
      gsl_vector_set(tmp, np_act, gsl_matrix_get(img_data, tr_point.x, y_act));
      np_act++;
    }

  // check on both sides
  // of the central tracepoint
  for (ii=1; ii <= np; ii++)
    {
      // go down
      y_act = tr_point.y - ii;

      // check whether we are still inside the image and
      // whether we have valid data
      if (y_act > -1 && !isnan(gsl_matrix_get(img_data, tr_point.x, y_act)))
	{
	  // store the data point and enhance the counter
	  gsl_vector_set(tmp, np_act, gsl_matrix_get(img_data, tr_point.x, y_act));
	  np_act++;
	}

      // go up
      y_act = tr_point.y + ii;

      // check whether we are still inside the image and
      // whether we have valid data
      if (y_act < ncols && !isnan(gsl_matrix_get(img_data, tr_point.x, y_act)))
	{
	  // store the data point and enhance the counter
	  gsl_vector_set(tmp, np_act, gsl_matrix_get(img_data, tr_point.x, y_act));
	  np_act++;
	}
    }

  // check whether good data was found
  if (np_act > 0)
    {
      // allocate the return vector
      ipc_row = gsl_vector_alloc(np_act);

      // transfer the data to the
      // return vector
      for (ii=0; ii < np_act; ii++)
	gsl_vector_set(ipc_row, ii, gsl_vector_get(tmp, ii));
    }

  // release memory
  gsl_vector_free(tmp);

  //  return the vector
  return ipc_row;
}

/*
 * Function: alloc_fit_data
 * The function allocates memory for a fit-data structure
 *
 * Parameters:
 * @param n_data - the number of points to allocate
 *
 * Returns:
 * @return f_data - the allocated structure
 */
fit_data *
alloc_fit_data(const int n_data)
{
  fit_data *f_data;

  f_data = (fit_data *)malloc(sizeof(fit_data));

  f_data->x_values = (double *)malloc(n_data * sizeof(double));
  f_data->y_values = (double *)malloc(n_data * sizeof(double));
  f_data->e_values = (double *)malloc(n_data * sizeof(double));

  f_data->n_data = n_data;

  return f_data;
}

/*
 * Function: print_fit_data
 * The function prints the content of a fit-data structure
 * onto the screen (error-out).
 *
 * Parameters:
 * @param f_data - the allocated structure
 *
 * Returns:
 * @return -
 */
void
print_fit_data(fit_data *f_data)
{
  int i;

  for (i=0; i < f_data->n_data; i++)
    {
      fprintf(stdout, "%f %f %f\n", f_data->x_values[i],
	      f_data->y_values[i], f_data->e_values[i]);
      //      fprintf(stderr, "x: %f, y: %f, err: %f\n", f_data->x_values[i],
      //	      f_data->y_values[i], f_data->e_values[i]);
    }
}

/*
 * Function: print_fit_data
 * The function releases the memory allocated in
 * a fit-data structure
 *
 * Parameters:
 * @param f_data - the allocated structure
 *
 * Returns:
 * @return -
 */
void
free_fit_data(fit_data *f_data)
{
  free(f_data->x_values);
  free(f_data->y_values);
  free(f_data->e_values);

  free(f_data);

  f_data = NULL;
}

void
print_fit_state (size_t iter, gsl_multifit_fdfsolver * s)
{
  printf ("iter: %3zu x = %e %e %e %e "
          "|f(x)| = %g |dx| = %g\n\n", iter,
          gsl_vector_get (s->x, 0),
          gsl_vector_get (s->x, 1),
          gsl_vector_get (s->x, 2),
          gsl_vector_get (s->x, 3),
          gsl_blas_dnrm2 (s->f),
          gsl_blas_dnrm2 (s->dx));
}

/**
 * Function: comp_intrel_max
 *
 * Parameters:
 * @param f_data  - the input data
 *
 * Returns:
 * @return rel12 - "second brightes row pixel divided by brighest row pixel"
 */
double
comp_intrel_max(const fit_data *f_data)
{

  int i;
  double rel12=-100000.0;

  gsl_vector *val_vector;

  // check whether there is enough
  // data to make the relation
  if (f_data->n_data < 2)
    return rel12;

  // allocate a vector
  val_vector = gsl_vector_alloc(f_data->n_data);

  // fill in the data
  for (i=0; i < f_data->n_data; i++)
    gsl_vector_set(val_vector, i, f_data->y_values[i]);

  // sort the values
  gsl_sort_vector(val_vector);

  // prevent division by zero
  if (gsl_vector_get(val_vector, 0))
    // compute the quantity second brightes pixel
    // divided by brighest pixel
    rel12 = gsl_vector_get(val_vector, f_data->n_data-2) / gsl_vector_get(val_vector, f_data->n_data-1);

  // free the memory
  gsl_vector_free(val_vector);

  // return the result
  return rel12;
}

/**
 * Function: bcksub_observation
 *
 * Parameters:
 * @param obs           - the grism image
 * @param bck           - the background image
 *
 * Returns:
 * @return ret - background subtracted data
 */
gsl_matrix *
bcksub_observation(observation *obs, observation *bck)
{
  gsl_matrix *ret;

  int ii;
  int jj;

  // initialize the return
  ret = gsl_matrix_alloc(obs->grism->size1, obs->grism->size2);

  // set all values to NAN
  gsl_matrix_set_all(ret, GSL_NAN);

  // go over all cols
  for (ii=0; ii < (int)obs->grism->size1; ii++)
    // go over all rows
    for (jj=0; jj < (int)obs->grism->size2; jj++)
      {
	// check whether there is valid data in
	// object and background image
	if (!isnan(gsl_matrix_get(obs->grism, ii, jj))
	    && !isnan(gsl_matrix_get(bck->grism, ii, jj))
	    && gsl_matrix_get(bck->grism, ii, jj))

	  // compute and set the pixel
	  gsl_matrix_set(ret, ii, jj, gsl_matrix_get(obs->grism, ii, jj) - gsl_matrix_get(bck->grism, ii, jj));
      }

  // return matrix
  return ret;
}

/**
 * Function: evaluate_ipc_fit
 *
 * Parameters:
 * @param f_data  - the input for the independent values
 * @param gpars   - the fited parameters
 * @param fpars   - the fixed parameters
 *
 * Returns:
 * @return y_values - the function values
 */
double *
evaluate_ipc_fit(fit_data *f_data, double *gpars, double *fpars)
{
  int index;

  //double yvalue;

  double *y_values;

  // allocate memory for the function values
  y_values = (double *)malloc(f_data->n_data * sizeof(double));

  // go over all independent data
  for (index=0; index < f_data->n_data; index++)
      // compute and store a function value
      y_values[index] = comp_shift(f_data->x_values[index], gpars, fpars);

  // return the function values
  return y_values;
}

/**
 * Function: write_ipc_fit
 *
 * Parameters:
 * @param f_data   - the input data
 * @param y_values - the model values
 * @param ipc_file - name of the file to write to
 *
 * Returns:
 * @return -
 */
void
write_ipc_fit(fit_data *f_data, double *y_values, const char *ipc_file, const double qual, double *yfract)
{
  FILE *fout;
  char Buffer[MAXCHAR];

  int index;

  // open the file for the values
  fout = fopen (ipc_file, "w");

  // write norm (= quality of fit) to comment line
  sprintf (Buffer, "# norm of the residue vector: %e\n", qual);
  fputs (Buffer, fout);
  // write norm (= quality of fit) to comment line
  sprintf (Buffer, "#%e\n", qual);
  fputs (Buffer, fout);
  sprintf (Buffer, "# Wanderer kommst Du nach Sparta: %e\n", qual);
  fputs (Buffer, fout);
  sprintf (Buffer, "# berichtige dorten: %e\n", qual);
  fputs (Buffer, fout);
  sprintf (Buffer, "# Du habest uns hier liegen gesehen: %e\n", qual);
  fputs (Buffer, fout);
  sprintf (Buffer, "# wie es die Geschichte befahl. %e\n", qual);
  fputs (Buffer, fout);
  sprintf (Buffer, "#%e\n", qual);
  fputs (Buffer, fout);

  // go over the data
  for (index=0; index < f_data->n_data; index++)
    {
      // write the data point to the buffer
      sprintf (Buffer, "%e %e %e %e %e\n", f_data->x_values[index], f_data->y_values[index], y_values[index],  f_data->e_values[index], yfract[index]);
      // write the buffer to the file
      fputs (Buffer, fout);
    }

  // close the file
  fclose(fout);
}

/**
 * Function: fit_ipc_data
 *
 * Parameters:
 * @param f_data  - the input data
 * @param gpars   - the guess or fit parameters
 * @param n_gpars - number of guess or fit parameters
 * @param fpars   - the fixed parameters
 * @param n_fpars - number of fixed parameters
 *
 * Returns:
 * @return control.fnorm - goodness of the fit value
 */
double
fit_ipc_data(fit_data *f_data, double *gpars, int n_gpars,
	     double *fpars, int n_fpars)
{
  int m_dat = f_data->n_data;

  // auxiliary settings:
  lm_control_type control;
  lm_data_fpar_type data;

  // initialize the fit
  lm_initialize_control( &control );

  // define the parameters for the fitting
  data.user_func = comp_shift;
  data.user_t = f_data->x_values;
  data.user_y = f_data->y_values;
  data.fpars  = fpars;

  // print the initial
  // value
  fprintf(stdout, "input: %e, ", gpars[0]);

  // perform the fit:
  lm_minimize( m_dat, n_gpars, gpars, lm_evaluate_fpar, lm_print_nothing,
  	       &data, &control );

  // print the new value plus additional information
  fprintf(stdout, " output: %e, info: %i, niter: %i norm: %e\n",
	  gpars[0], control.info, control.nfev, control.fnorm);

  // return the variable
  // describing the goodness
  // of the fit
  return control.fnorm;
}

/**
 * Function: get_ipc_fdata
 *
 * Parameters:
 * @param conf          - the aXe configuratuion structure
 * @param obs           - the grism image
 * @param data_matrix   - background subtracted grism image
 * @param actbeam       - the beam to examine
 *
 * Returns:
 * @return f_final - the filled data structure
 */
fit_data *
get_ipc_fdata(const aperture_conf  *conf,  observation *obs,
	      gsl_matrix *data_vals, beam act_beam)
{
  trace_func     *tracefun;

  px_point        xborder;
  px_point        tpoint;

  gsl_vector      *ipc_data;

  double rel12;

  int i;
  int index;

  fit_data *fdata  = NULL;
  fit_data *f_final= NULL;

  // define the beam and the trace function
  tracefun = act_beam.spec_trace;

  // If this beam's ignore flag is set to 1 then do nothing
  if (act_beam.ignore == 1)
    return fdata;

  // determine the start and end point in x
  xborder = get_fit_xrange(conf, obs, act_beam);

  // an assertion to prevent
  // catastrophic failures
  if  (xborder.x > xborder.y)
    return fdata;
  // introduced after making release 1 reduction

  // allocate memory for the fitting structure
  fdata = alloc_fit_data(xborder.y-xborder.x);

  // Loop over all columns
  index = 0;
  for (i = xborder.x; i < xborder.y; i++)
    {
      // determine the pixel which is
      // closest to the nominal trace position
      tpoint.x = i;
      tpoint.y = (int)floor(tracefun->func((double)i-act_beam.refpoint.x,
					   tracefun->data)
			    + act_beam.refpoint.y+0.5);

      // collect pixel data
      ipc_data = get_ipc_coldata(data_vals, N_COLUMNS_FIT, tpoint);

      // go to the next
      // row if no data is available
      if (!ipc_data)
	continue;

      // sort the values
      gsl_sort_vector(ipc_data);

      // at least two values are needed;
      // prevent division by zero
      if (ipc_data->size > 1 && gsl_vector_get(ipc_data, ipc_data->size-1))
	{
	  // compute the relation
	  rel12 = gsl_vector_get(ipc_data, ipc_data->size-2) / gsl_vector_get(ipc_data, ipc_data->size-1);

	  fdata->x_values[index] = (double)i - act_beam.refpoint.x;
	  fdata->y_values[index] = rel12;
	  fdata->e_values[index] = 1.0;
	  index++;
	}

      // free the allocated memory
      gsl_vector_free(ipc_data);
    }

  // check whether the ipc data covers the minimum length
  if (fabs(fdata->x_values[fdata->n_data-1] - fdata->x_values[0]) > MIN_IPC_LENGTH)
    // reduce the data to what is necessary
    f_final = strip_fitdata(fdata, index);
  else
    fprintf(stdout, "beam tooo short\n");

  // release the original structure
  free_fit_data(fdata);

  // return the fit data
  return f_final;
}

/**
 * Function: get_wf_special
 *
 * Parameters:
 * @param conf          - the aXe configuratuion structure
 * @param obs           - the grism image
 * @param data_matrix   - background subtracted grism image
 * @param actbeam       - the beam to examine
 *
 * Returns:
 * @return f_final - the filled data structure
 */
double *
get_wf_special(const aperture_conf  *conf,  observation *obs,
	       gsl_matrix *data_vals, beam act_beam, double yoffs)
{
  trace_func     *tracefun;

  px_point        xborder;
  px_point        tpoint;

  gsl_vector      *ipc_data;

  double rel12;

  int i;
  int index;

  fit_data *fdata  = NULL;
  //fit_data *f_final= NULL;

  double yyy;

  double *yfract=NULL;
  double *yfinal = 0;

  // define the beam and the trace function
  tracefun = act_beam.spec_trace;

  // If this beam's ignore flag is set to 1 then do nothing
  if (act_beam.ignore == 1)
    return yfinal;

  // determine the start and end point in x
  xborder = get_fit_xrange(conf, obs, act_beam);

  // allocate memory for the fitting structure
  fdata = alloc_fit_data(xborder.y-xborder.x);

  // allocate memory
  yfract = (double *)malloc((xborder.y-xborder.x)*sizeof(double));

  // Loop over all columns
  index = 0;
  for (i = xborder.x; i < xborder.y; i++)
    {
      // determine the pixel which is
      // closest to the nominal trace position
      tpoint.x = i;
      tpoint.y = (int)floor(tracefun->func((double)i-act_beam.refpoint.x,
					   tracefun->data)
			    + act_beam.refpoint.y+0.5);

      yyy           = tracefun->func((double)i-act_beam.refpoint.x, tracefun->data) + act_beam.refpoint.y + yoffs;
      yfract[index] = yyy - floor(yyy);

      // collect pixel data
      ipc_data = get_ipc_coldata(data_vals, N_COLUMNS_FIT, tpoint);

      // go to the next
      // row if no data is available
      if (!ipc_data)
	continue;

      // sort the values
      gsl_sort_vector(ipc_data);

      // at least two values are needed;
      // prevent division by zero
      if (ipc_data->size > 1 && gsl_vector_get(ipc_data, ipc_data->size-1))
	{
	  // compute the relation
	  rel12 = gsl_vector_get(ipc_data, ipc_data->size-2) / gsl_vector_get(ipc_data, ipc_data->size-1);

	  fdata->x_values[index] = (double)i - act_beam.refpoint.x;
	  fdata->y_values[index] = rel12;
	  fdata->e_values[index] = 1.0;
	  index++;
	}

      // free the allocated memory
      gsl_vector_free(ipc_data);
    }

  yfinal = strip_wf_special(fdata, yfract, index);

  // release the original structure
  free_fit_data(fdata);
  if (yfract)
    free(yfract);

  // return the fit data
  return yfinal;
}

/**
 * Function: strip_wf_special
 *
 * Parameters:
 * @param old_fdata - the experimental data
 * @param i_max     - the maximum index to check for valid data
 *
 * Returns:
 * @return fdata - the new, completely filled data structure
 */
double *
strip_wf_special(fit_data *old_fdata, double *y_fract, int i_max)
{
  double *y_new=NULL;

  int ii;
  int i_count;

  // check that SOME values
  // are accepted
  if (i_max < 1)
    // return NULL if not
    return y_new;

  // count the number of
  // values with weight
  i_count=0;
  for (ii=0; ii < i_max; ii++)
    if (old_fdata->e_values[ii])
      i_count++;

  // check whether there
  // is data at all
  if (i_count < 1)
    // return NULL if not
    return y_new;

  // allocate a new structure
  y_new = (double *)malloc(i_count * sizeof(double));

  // set the counter
  i_count=0;

  // go over the old data
  // to the maximum index
  for (ii=0; ii < i_max; ii++)
    {
      // check whether the
      // data is valid
      if (old_fdata->e_values[ii])
	{
	  // transfer valid data
	  y_new[i_count] = y_fract[ii];

	  // enhance the counter
	  i_count++;
	}
    }

  // return the new structure
  return y_new;
}

/**
 * Function: strip_fitdata
 *
 * Parameters:
 * @param old_fdata - the experimental data
 * @param i_max     - the maximum index to check for valid data
 *
 * Returns:
 * @return fdata - the new, completely filled data structure
 */
fit_data *
strip_fitdata(fit_data *old_fdata, int i_max)
{
  fit_data *fdata = NULL;

  int ii;
  int i_count;

  // check that SOME values
  // are accepted
  if (i_max < 1)
    // return NULL if not
    return fdata;

  // count the number of
  // values with weight
  i_count=0;
  for (ii=0; ii < i_max; ii++)
    if (old_fdata->e_values[ii])
      i_count++;

  // check whether there
  // is data at all
  if (i_count < 1)
    // return NULL if not
    return fdata;

  // allocate a new structure
  fdata = alloc_fit_data(i_count);

  // set the counter
  i_count=0;

  // go over the old data
  // to the maximum index
  for (ii=0; ii < i_max; ii++)
    {
      // check whether the
      // data is valid
      if (old_fdata->e_values[ii])
	{
	  // transfer valid data
	  fdata->x_values[i_count] = old_fdata->x_values[ii];
	  fdata->y_values[i_count] = old_fdata->y_values[ii];
	  fdata->e_values[i_count] = old_fdata->e_values[ii];

	  // enhance the counter
	  i_count++;
	}
    }

  // return the new structure
  return fdata;
}

/**
 * Function: kappa_sigma_klipp_ipc
 *
 * Parameters:
 * @param conf          - the aXe configuratuion structure
 * @param obs           - the grism image
 * @param data_matrix   - background subtracted grism image
 * @param actbeam       - the beam to examine
 * @param ipc_file_path - file name for the phase data
 *
 * Returns:
 * @return ret - the change in y-shift and the goodness of the fit
 */
d_point
kappa_sigma_klipp_ipc(const aperture_conf  *conf, observation *obs,
		      gsl_matrix *data_matrix, beam act_beam,
		      const char *ipc_file)
{
  fit_data *f_data   = NULL;
  fit_data *f_work   = NULL;

  double   *y_values = NULL;
  //double   *test = NULL;
  double    sigma;

  d_point   ret;

  trace_func *tracefun  = act_beam.spec_trace;
  double     *tracedata =  (double *)(tracefun->data);

  double p[1];
  double f[3] = {0.6, 0.9, tracedata[2]};

  int index=0;
  int nold_data;

  // REMOVEME ASAP !!!!!
  double *yfract=NULL;

  // compute a first guess for the shift
  p[0] = act_beam.refpoint.y - floor(act_beam.refpoint.y + tracedata[1]);

  // remember the shift
  ret.x = p[0];

  // give infinity for
  // the goodness of the fit
  ret.y = 1.0e+20;

  // get the data for the fit
  f_data = get_ipc_fdata(conf, obs, data_matrix, act_beam);

  // check whether
  // there is data
  if (f_data)
    {

      // fix the number of initial data points
      nold_data= f_data->n_data;

      // make the frist fit
      ret.y = fit_ipc_data(f_data, p, 1, f, 3);

      // evaluate the fit; comute the sigma
      y_values = evaluate_ipc_fit(f_data, p, f);
      sigma = calculate_sigma(f_data, y_values);

      // start the iteration
      while (index < N_ITER_IPC)
	{

	  // reject points according to the sigma
	  // generate a new struct to work with
	  reject_kappasigma(f_data, y_values, KAPPA_IPC*sigma);
	  // free memory in the old work structure

	  if (f_work)
	    free_fit_data(f_work);
	  f_work = strip_fitdata(f_data, f_data->n_data);

	  // check whether there
	  // is still data
	  if (!f_work)
	    break;

	  // check whether data was rejected
	  // brek if 'yes' transfer the new number if 'no'
	  if (f_work->n_data == nold_data)
	    break;
	  else
	    nold_data = f_work->n_data;

	  // make a new fit
	  ret.y = fit_ipc_data(f_work, p, 1, f, 3);

	  // evaluate fit values; comute the sigma
	  if (y_values)
	    free(y_values);
	  y_values = evaluate_ipc_fit(f_data, p, f);
	  sigma = calculate_sigma(f_data, y_values);

	  // enhance
	  // the kappa counter
	  index++;
	}

      // evaluate the fit values
      if (y_values)
	free(y_values);
      y_values = evaluate_ipc_fit(f_data, p, f);

      // REMOVEME ASAP
      yfract = get_wf_special(conf, obs, data_matrix, act_beam, p[0] - ret.x);

      // CHANGEME back ASAP
      // write the phase shift file
      write_ipc_fit(f_data, y_values, ipc_file, ret.y, yfract);

      // release memory
      if (f_data)
	free_fit_data(f_data);
      if (f_work)
	free_fit_data(f_work);
      if (y_values)
	free(y_values);
      if (yfract)
	free(yfract);
    }

  // compute the
  // difference in y-shift
  ret.x = p[0] - ret.x;

  // return the result
  return ret;
}

/**
 * Function: reject_kappasigma
 *
 * Parameters:
 * @param f_data   - the experimental data
 * @param y_values - the model data
 * @param max_diff - the maximum allowed difference
 *                   between data and model
 *
 * Returns:
 * @return -
 */
void
reject_kappasigma(fit_data *f_data, double *y_values, double max_diff)
{
  int index;

  // go over all data
  for (index=0; index < f_data->n_data; index++)
    // check whether the data point is valid and OVER the maximum allowed difference
    if (f_data->e_values[index] && fabs(f_data->y_values[index] - y_values[index]) > max_diff)
      // make the data point INVALID
      f_data->e_values[index] = 0.0;
}

/**
 * Function: calculate_sigma
 *
 * Parameters:
 * @param f_data   - the experimental data
 * @param y_values - the model data
 *
 * Returns:
 * @return sigma - the standard deviation
 */
double
calculate_sigma(fit_data *f_data, double *y_values)
{
  double *tmp;

  double sigma;
  double mean=0.0;
  double diff;

  int n_data=0;
  int index;

  // count the number of valid data points
  for (index=0; index < f_data->n_data; index++)
    if (f_data->e_values[index])
      n_data++;

  // allocate a temporary vector for
  // the valid data
  tmp = (double *)malloc(n_data * sizeof(double));

  // re-set the counter;
  // go over all data points
  n_data=0;
  for (index=0; index < f_data->n_data; index++)
    {
      // check whether the data point
      // is valid
      if (f_data->e_values[index])
	{
	  // compute the difference between data and model
	  diff =  f_data->y_values[index] - y_values[index];

	  // store the difference
	  tmp[n_data] = diff;

	  // add to the mean
	  mean += diff;

	  // enhance the counter
	  n_data++;
	}
    }

  // comopute the mean
  mean /= (double)n_data;

  // compute standard deviation
  sigma = gsl_stats_sd_m(tmp, 1, n_data, mean);

  // relese the memory
  free(tmp);

  // return the sigma
  return sigma;
}

max_stars_repo_name: sosey/pyaxe

max_issues_repo_name: sosey/pyaxe

max_forks_repo_name: sosey/pyaxe


text:
/*
** Implementation of LISA algorithm
** for statistical inference of fMRI images
**
** 2nd level inference for grouped samples, two-sample ttest, paired t-test
**
** G.Lohmann, MPI-KYB, 2018
** Zalan Rajna, University of Oulu, 2020
*/
#include <viaio/Vlib.h>
#include <viaio/file.h>
#include <viaio/mu.h>
#include <viaio/option.h>
#include <viaio/os.h>
#include <viaio/VImage.h>
#include <via/via.h>

#include <gsl/gsl_cdf.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_histogram.h>
#include <gsl/gsl_histogram2d.h>

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

#ifdef _OPENMP
#include <omp.h>
#endif /*_OPENMP*/

#define ABS(x) ((x) > 0 ? (x) : -(x))
#define MIN(x,y) ((x) < (y) ? (x) : (y))
#define MAX(x,y) ((x) > (y) ? (x) : (y))

extern double ttest2(double *data1,double *data2,int n1,int n2);
extern double xtest2(double *data1,double *data2,int n);
extern double paired_ttest(double *data1,double *data2,int n);
extern double welchtest(double *data1,double *data2,int n1,int n2);
extern void   VIsolatedVoxels(VImage src,float threshold);
extern void   VHistogram(gsl_histogram *histogram,VString filename);
extern void   VCheckImage(VImage src);
extern void   FDR(VImage src,VImage dest,gsl_histogram *nullhist,gsl_histogram *realhist,double);
extern double ttest1(double *data1,int n);
extern void   ImageStats(VImage src,double *,double *,double *hmin,double *hmax);
extern void   VBilateralFilter(VImage src,VImage dest,int radius,double var1,double var2,int);
extern double VImageVar(VImage src);
extern void   VGetHistRange(VImage src,double *hmin,double *hmax);
extern void   VZScale(VImage src,float mode,float stddev);
extern float  VGetMode(VImage src);
extern void   HistoUpdate(VImage,gsl_histogram *);

/* generate permutation table */
int **genperm(long seed,int n1,int n2,int numperm,int testtype)
{
  int i,j,n=n1+n2;
  int **table = (int **) VCalloc(numperm,sizeof(int *));
  int *base = (int *) VCalloc(n,sizeof(int));
  for (j=0; j<n; j++) base[j] = j;

  gsl_rng_env_setup();
  const gsl_rng_type *T = gsl_rng_default;
  gsl_rng *rx = gsl_rng_alloc(T);
  gsl_rng_set(rx,(unsigned long int)seed);

  /* for each permutation:  */
  for (i = 0; i < numperm; i++) {

    /* unpaired test */
    if (testtype == 0 || testtype == 2) {
      table[i] = (int *) VCalloc(n,sizeof(int));
      gsl_ran_shuffle (rx,base,n,sizeof(int));
      for (j=0; j<n; j++) {
	table[i][j] = base[j];
      }
    }

    /* paired test */
    if (testtype == 1) {
      table[i] = (int *) VCalloc(n1,sizeof(int));
      for (j=0; j<n1; j++) {
	table[i][j] = 0;
	if (gsl_ran_bernoulli (rx,(double)0.5) == 1) table[i][j] = 1;
      }
    }
  }

  gsl_rng_free(rx);
  VFree(base);
  return table;
}

void TTest(VImage *src1,VImage *src2,int *permtable,VImage dest,int n1,int n2,int testtype)
{
  int i,b,r,c,nslices,nrows,ncols;
  double u=0,v=0,z=0;
  int k1=0,k2=0,k=0;
  int n = n1+n2;
  int min1 = n1-2;
  int min2 = n2-2;

  gsl_set_error_handler_off();

  nslices = VImageNBands(src1[0]);
  nrows   = VImageNRows(src1[0]);
  ncols   = VImageNColumns(src1[0]);
  VFillImage(dest,VAllBands,0);

  /* alloc */
  double *data1 = (double *) VCalloc(n1,sizeof(double));
  double *data2 = (double *) VCalloc(n2,sizeof(double));
  double *data = (double *) VCalloc(n,sizeof(double));

  /* for every voxel */
  for (b=0; b<nslices; b++) {
    for (r=0; r<nrows; r++) {
      for (c=0; c<ncols; c++) {

	/* unpaired test */
	if (testtype == 0 || testtype == 2) {

	  k = 0;
	  for (i=0; i<n1; i++) {
	    data[k++] = VPixel(src1[i],b,r,c,VFloat);
	  }
	  for (i=0; i<n2; i++) {
	    data[k++] = VPixel(src2[i],b,r,c,VFloat);
	  }

	  k1 = 0;
	  for (i=0; i<n1; i++) {
	    u = data[permtable[i]];
	    if (fabs(u) > 0) {
	      data1[k1++] = u;
	    }
	  }
	  if (k1 < min1) continue;

	  k2 = 0;
	  for (i=0; i<n2; i++) {
	    u = data[permtable[i+n1]];
	    if (fabs(u) > 0) {
	      data2[k2++] = u;
	    }
	  }
	  if (k2 < min2) continue;
	}

	/* paired test */
	if (testtype == 1) {
	  k=0;
	  for (i=0; i<n1; i++) {
	    if (permtable[i] == 0) {  /* no flip */
	      u = VPixel(src1[i],b,r,c,VFloat);
	      v = VPixel(src2[i],b,r,c,VFloat);
	    }
	    else {   /* flip */
	      u = VPixel(src2[i],b,r,c,VFloat);
	      v = VPixel(src1[i],b,r,c,VFloat);
	    }
	    if (fabs(u) > 0 && fabs(v) > 0) {
	      data1[k] = u;
	      data2[k] = v;
	      k++;
	    }
	  }
	  if (k < min1) continue;
	}

	/* t-test */
	z = 0;
	switch(testtype) {
	case 0 :
	  z = ttest2(data1,data2,k1,k2);     /* pooled variance */
	  break;
	case 1 :
	  z = paired_ttest(data1,data2,k);   /* paired */
	  break;
	case 2 :
	  z = welchtest(data1,data2,k1,k2);  /* unequal variance */
	  break;
	default:
	  VError(" unknown testtype");
	}
	VPixel(dest,b,r,c,VFloat) = z;
      }
    }
  }
  VFree(data);
  VFree(data1);
  VFree(data2);
}

VDictEntry TSTDict[] = {
  { "pooled", 0, 0,0,0,0 },
  { "paired", 1, 0,0,0,0 },
  { "welch", 2, 0,0,0,0  },
  { NULL, 0,0,0,0,0 }
};

int main (int argc, char *argv[])
{
  static VArgVector in_files1;
  static VArgVector in_files2;
  static VString  out_filename="";
  static VString  mask_filename="";
  static VFloat   alpha = 0.05;
  static VShort   testtype = 2;  /* default is welch */
  static VShort   radius = 2;
  static VFloat   rvar = 2.0;
  static VFloat   svar = 2.0;
  static VShort   numiter = 2;
  static VShort   numperm = 5000;
  static VLong    seed = 99402622;
  static VBoolean centering = FALSE;
  static VBoolean cleanup = TRUE;
  static VShort   nproc = 0;
  static VOptionDescRec options[] = {
    {"in1", VStringRepn, 0, & in_files1, VRequiredOpt, NULL,"Input files 1 (single file per group)" },
    {"in2", VStringRepn, 0, & in_files2, VRequiredOpt, NULL,"Input files 2 (single file per group)" },
    {"out", VStringRepn, 1, & out_filename, VRequiredOpt, NULL,"Output file" },
    {"alpha",VFloatRepn,1,(VPointer) &alpha,VOptionalOpt,NULL,"FDR significance level"},
    {"perm",VShortRepn,1,(VPointer) &numperm,VOptionalOpt,NULL,"Number of permutations"},
    {"mask", VStringRepn, 1, (VPointer) &mask_filename, VRequiredOpt, NULL, "Mask"},
    {"test",VShortRepn,1,(VPointer) &testtype,VOptionalOpt,TSTDict,"type of test"},
    {"seed",VLongRepn,1,(VPointer) &seed,VOptionalOpt,NULL,"Seed for random number generation"},
    {"radius",VShortRepn,1,(VPointer) &radius,VOptionalOpt,NULL,"Bilateral parameter (radius in voxels)"},
    {"rvar",VFloatRepn,1,(VPointer) &rvar,VOptionalOpt,NULL,"Bilateral parameter (radiometric)"},
    {"svar",VFloatRepn,1,(VPointer) &svar,VOptionalOpt,NULL,"Bilateral parameter (spatial)"},
    {"filteriterations",VShortRepn,1,(VPointer) &numiter,VOptionalOpt,NULL,"Bilateral parameter (number of iterations)"},
    {"cleanup",VBooleanRepn,1,(VPointer) &cleanup,VOptionalOpt,NULL,"Whether to apply cleanup"},
    {"j",VShortRepn,1,(VPointer) &nproc,VOptionalOpt,NULL,"Number of processors to use, '0' to use all"},
  };

  FILE *fp=NULL;
  VString in_filename,str1,str2;
  VAttrList list1=NULL,list2=NULL,out_list=NULL,geolist=NULL;
  int i,ngroups1=0,ngroups2=0,npix=0;
  char *prg_name=GetLipsiaName("vlisa_grouped_twosample");
  fprintf (stderr, "%s\n", prg_name);

  /* parse command line */
  if (! VParseCommand (VNumber (options), options, & argc, argv)) {
    VReportUsage (argv[0], VNumber (options), options, NULL);
    exit (EXIT_FAILURE);
  }
  if (argc > 1) {
    VReportBadArgs (argc, argv);
    exit (EXIT_FAILURE);
  }

  /* omp-stuff */
#ifdef _OPENMP
  int num_procs=omp_get_num_procs();
  if (nproc > 0 && nproc < num_procs) num_procs = nproc;
  fprintf(stderr," using %d cores\n",(int)num_procs);
  omp_set_num_threads(num_procs);
#endif /* _OPENMP */

  /* read mask */
  VImage mask = VReadImageFile(mask_filename);
  if (mask==NULL) VError("Error reading mask file %s",mask_filename);

  /* number of inputs */
  ngroups1 = in_files1.number;
  ngroups2 = in_files2.number;
  fprintf(stderr," ngroups= %d  %d\n",ngroups1,ngroups2);

  /* images 1  */
  VImage *src1 = (VImage *) VCalloc(ngroups1,sizeof(VImage));
  for (i = 0; i < ngroups1; i++) {
    in_filename = ((VString *) in_files1.vector)[i];
    list1   = VReadAttrList(in_filename,0L,TRUE,FALSE);
    VMaskMinval(list1,mask,0.0);
    src1[i] = VReadImage(list1);
    if (src1[i] == NULL) VError(" no input image found");
    if (VPixelRepn(src1[i]) != VFloatRepn) VError(" input pixel repn must be float");
    if (i == 0) npix = VImageNPixels(src1[i]);
    else if (npix != VImageNPixels(src1[i])) VError(" inconsistent image dimensions");

    /* use geometry info from 1st file */
    if (geolist == NULL) geolist = VGetGeoInfo(list1);
  }

  /* images 2  */
  VImage *src2 = (VImage *) VCalloc(ngroups2,sizeof(VImage));
  for (i = 0; i < ngroups2; i++) {
    in_filename = ((VString *) in_files2.vector)[i];
    list2   = VReadAttrList(in_filename,0L,TRUE,FALSE);
    VMaskMinval(list2,mask,0.0);
    src2[i] = VReadImage(list2);
    if (src2[i] == NULL) VError(" no input image found");
    if (VPixelRepn(src2[i]) != VFloatRepn) VError(" input pixel repn must be float");
    if (i == 0) npix = VImageNPixels(src2[i]);
    else if (npix != VImageNPixels(src2[i])) VError(" inconsistent image dimensions");
  }

  if (testtype == 1) {  /* print filenames to terminal */
    for (i = 0; i < ngroups1; i++) {
      str1 = ((VString *) in_files1.vector)[i];
      str2 = ((VString *) in_files2.vector)[i];
      fprintf(stderr," %3d:  %s   %s\n",i,str1,str2);
    }
  }
  else {
    fprintf(stderr," Group 1:\n");
    for (i = 0; i < ngroups1; i++) {
      str1 = ((VString *) in_files1.vector)[i];
      fprintf(stderr," %3d:  %s\n",i,str1);
    }
    fprintf(stderr,"\n Group 2:\n");
    for (i = 0; i < ngroups2; i++) {
      str2 = ((VString *) in_files2.vector)[i];
      fprintf(stderr," %3d:  %s\n",i,str2);
    }
  }

  /* random permutations */
  size_t n = ngroups1+ngroups2;
  int nperm=0;
  int **permtable = genperm((long)seed,(int)ngroups1,(int)ngroups2,(int)numperm,(int)testtype);
  int *nopermtable = (int *) VCalloc(n,sizeof(int));
  for (i=0; i<n; i++) nopermtable[i]=i;
  if (testtype == 1) for (i=0; i<n; i++) nopermtable[i]=0;

  /* estimate null variance based on first 30 permutations */
  double hmin=0,hmax=0;
  float stddev=1.0;
  if (numperm > 0) {
    int tstperm = 30;
    if (tstperm > numperm) tstperm = numperm;
    double varsum=0,nx=0;

#pragma omp parallel for shared(src1,permtable) schedule(dynamic)
    for (nperm = 0; nperm < tstperm; nperm++) {
      VImage zmap = VCreateImageLike(src1[0]);
      TTest(src1,src2,permtable[nperm],zmap,ngroups1,ngroups2,(int)testtype);
#pragma omp critical
      {
	varsum += VImageVar(zmap);
	nx++;
      }
      VDestroyImage(zmap);
    }
    double meanvar = varsum/nx;
    stddev = sqrt(meanvar);
  }

  /* no permutation */
  VImage dst1  = VCreateImageLike (src1[0]);
  VImage zmap1 = VCreateImageLike(src1[0]);
  VFillImage(zmap1,VAllBands,0);
  TTest(src1,src2,nopermtable,zmap1,ngroups1,ngroups2,(int)testtype);

  if (numperm == 0) {
    double z = VImageVar(zmap1);
    stddev = (float)(sqrt(z)); /* update stddev */
  }
  float mode=0;
  if (centering) mode = VGetMode(zmap1);
  if (numperm > 0) VZScale(zmap1,mode,stddev);
  VBilateralFilter(zmap1,dst1,(int)radius,(double)rvar,(double)svar,(int)numiter);

  /* ini histograms */
  VGetHistRange(dst1,&hmin,&hmax);
  size_t nbins = 20000;
  gsl_histogram *hist0 = gsl_histogram_alloc (nbins);
  gsl_histogram_set_ranges_uniform (hist0,hmin,hmax);
  gsl_histogram *histz = gsl_histogram_alloc (nbins);
  gsl_histogram_set_ranges_uniform (histz,hmin,hmax);
  HistoUpdate(dst1,histz);

#pragma omp parallel for shared(src1,src2,permtable) schedule(dynamic)
  for (nperm = 0; nperm < numperm; nperm++) {
    if (nperm%20 == 0) fprintf(stderr," perm  %4d  of  %d\r",nperm,(int)numperm);

    VImage zmap = VCreateImageLike(src1[0]);
    VImage dst  = VCreateImageLike (zmap);
    TTest(src1,src2,permtable[nperm],zmap,ngroups1,ngroups2,(int)testtype);

    float mode=0;
    if (centering) mode = VGetMode(zmap);
    VZScale(zmap,mode,stddev);
    VBilateralFilter(zmap,dst,(int)radius,(double)rvar,(double)svar,(int)numiter);

#pragma omp critical
    {
      HistoUpdate(dst,hist0);
    }
    VDestroyImage(dst);
    VDestroyImage(zmap);
  }

  /* apply fdr */
  VImage fdrimage = VCopyImage (dst1,NULL,VAllBands);
  if (numperm > 0) {
    FDR(dst1,fdrimage,hist0,histz,(double)alpha);
    if (cleanup && alpha < 1.0) {
      VIsolatedVoxels(fdrimage,(float)(1.0-alpha));
    }
  }

  /* write output to disk */
  out_list = VCreateAttrList ();
  VHistory(VNumber(options),options,prg_name,&list1,&out_list);
  VSetGeoInfo(geolist,out_list);
  VAppendAttr (out_list,"image",NULL,VImageRepn,fdrimage);
  fp = VOpenOutputFile (out_filename, TRUE);
  if (! VWriteFile (fp, out_list)) exit (1);
  fclose(fp);
  fprintf (stderr, "\n");
  fprintf (stderr, "%s: done.\n", argv[0]);
  exit(0);
}

max_stars_repo_name: zrajna/lipsia

max_issues_repo_name: zrajna/lipsia

max_forks_repo_name: zrajna/lipsia


text:
#ifndef __DESIGN_H
#define __DESIGN_H

#include <gsl/gsl_matrix.h>

//function prototypes
double distance(double *x,double *y,int D,double p);
double cost(gsl_matrix *data,int Npoints,int D,double p,double lambda);
double diagonalCost(int Npoints,double lambda);
double swap(gsl_matrix *data,int Npoints,int D,double p,double lambda,int i1,int i2, int d);
void swapBack(gsl_matrix *data,int i1,int i2,int d);

//main design sampler prototype
double sample(int Npoints,int D,double p,double lambda,int seed,int maxIterations,gsl_matrix *data,double *costValues);

#endif

ext: h

max_stars_repo_name: asabyr/LensTools

max_issues_repo_name: asabyr/LensTools

max_forks_repo_name: asabyr/LensTools


text:
#pragma once

#include "Core/Assert.h"
#include "Core/Delegate.h"
#include "Core/Pointers.h"
#include "Scene/Components/Component.h"

#include <gsl/span>

#include <string>
#include <utility>
#include <vector>

class Scene;

class Entity
{
public:
   using OnDestroyDelegate = MulticastDelegate<void, Entity*>;

   ~Entity();

   void destroy();

   void tick(float dt);

   template<typename T>
   T* createComponent();

   Component* createComponentByName(const std::string& className);
   bool destroyComponent(Component* componentToDestroy);

   template<typename T>
   T* getComponentByClass();

   template<typename T>
   const T* getComponentByClass() const;

   template<typename T>
   std::vector<T*> getComponentsByClass();

   template<typename T>
   std::vector<const T*> getComponentsByClass() const;

   DelegateHandle addOnDestroyDelegate(OnDestroyDelegate::FuncType&& function);
   void removeOnDestroyDelegate(const DelegateHandle& handle);

   Scene& getScene()
   {
      return scene;
   }

   const Scene& getScene() const
   {
      return scene;
   }

private:
   friend class Scene;

   template<typename... ComponentTypes>
   static UPtr<Entity> create(Scene& scene);

   static UPtr<Entity> create(gsl::span<std::string> componentClassNames, Scene& scene);

   Entity(Scene& owningScene)
      : scene(owningScene)
   {
   }

   template<typename First, typename... Rest>
   void constructComponentsHelper();

   template<typename... ComponentTypes>
   void constructComponents();

   void onInitialized();
   void onComponentCreated(Component* component);

   std::vector<UPtr<Component>> components;
   OnDestroyDelegate onDestroyDelegate;
   Scene& scene;
};

template<typename T>
T* Entity::createComponent()
{
   UPtr<T> newComponent = ComponentRegistry::instance().createComponent<T>(*this);
   ASSERT(newComponent, "Failed to create a typed component");

   T* newComponentRaw = newComponent.get();
   components.push_back(std::move(newComponent));
   onComponentCreated(newComponentRaw);
   return newComponentRaw;
}

template<typename T>
T* Entity::getComponentByClass()
{
   for (const UPtr<Component>& component : components)
   {
      if (T* castedComponent = dynamic_cast<T*>(component.get()))
      {
         return castedComponent;
      }
   }

   return nullptr;
}

template<typename T>
const T* Entity::getComponentByClass() const
{
   for (const UPtr<Component>& component : components)
   {
      if (const T* castedComponent = dynamic_cast<const T*>(component.get()))
      {
         return castedComponent;
      }
   }

   return nullptr;
}

template<typename T>
std::vector<T*> Entity::getComponentsByClass()
{
   std::vector<T*> componentsOfClass;

   for (const UPtr<Component>& component : components)
   {
      if (T* castedComponent = dynamic_cast<T*>(component.get()))
      {
         componentsOfClass.push_back(castedComponent);
      }
   }

   return componentsOfClass;
}

template<typename T>
std::vector<const T*> Entity::getComponentsByClass() const
{
   std::vector<const T*> componentsOfClass;

   for (const UPtr<Component>& component : components)
   {
      if (const T* castedComponent = dynamic_cast<const T*>(component.get()))
      {
         componentsOfClass.push_back(castedComponent);
      }
   }

   return componentsOfClass;
}

// static
template<typename... ComponentTypes>
UPtr<Entity> Entity::create(Scene& scene)
{
   UPtr<Entity> entity(new Entity(scene));

   entity->constructComponents<ComponentTypes...>();
   entity->onInitialized();

   return entity;
}

template<typename First, typename... Rest>
void Entity::constructComponentsHelper()
{
   components.push_back(ComponentRegistry::instance().createComponent<First>(*this));
   constructComponents<Rest...>();
}

template<typename... ComponentTypes>
void Entity::constructComponents()
{
   constructComponentsHelper<ComponentTypes...>();
}

template<>
inline void Entity::constructComponents()
{
}

max_stars_repo_name: aaronmjacobs/Swap

max_issues_repo_name: aaronmjacobs/Swap

max_forks_repo_name: aaronmjacobs/Swap


text:
// Copyright (c) Microsoft Corporation.
// Licensed under the MIT License.

#pragma once

#include <algorithm>
#include <array>
#include <vector>
#include <gsl/gsl>
#include <assert.h>
#include <boost/optional.hpp>
#include "Data/Intrinsics.h"
#include "TimeUnits.h"
#include "MageSettings.h"

namespace mage
{
    enum class TrackingState
    {
        INITIALIZING = 0,
        TRACKING = 1,
        RELOCALIZING = 2,
        SKIPPED = 3
    };

    enum class FuserMode
    {
        Invalid = 0 ,
        WaitForMageInit = 1,
        WaitForGravityConverge = 2,
        VisualTrackingLost = 3,
        VisualTrackingReacquired = 4,
        ScaleInit = 5,
        Tracking = 6
    };

    enum class PixelFormat
    {
        GRAYSCALE8,
        NV12
    };

    struct FrameId
    {
        std::uint64_t CorrelationId;    //this is the number used to match frames either from each stereo cameras, or from low-res to high-res in the mono case. Should be unique (and increasing) per camera.
        CameraIdentity Camera;

        FrameId()
            : CorrelationId{ 0 }, Camera{ mage::CameraIdentity::MONO }
        {}

        FrameId(const std::uint64_t& id, CameraIdentity cam)
            : CorrelationId{ id }, Camera{ cam }
        {}

        bool operator ==(const FrameId& other) const
        {
            return CorrelationId == other.CorrelationId && Camera == other.Camera;
        }

        bool operator !=(const FrameId& other) const
        {
            return !(*this == other);
        }

        bool operator <(const FrameId& other) const
        {
            if (CorrelationId < other.CorrelationId)
                return true;

            if (other.CorrelationId < CorrelationId)
                return false;

            return Camera < other.Camera;
        }
    };

    struct Size
    {
        size_t Width;
        size_t Height;
    };

    struct Rect
    {
        int X;
        int Y;
        size_t Width;
        size_t Height;
    };

    struct Matrix
    {
        float M11, M12, M13, M14;
        float M21, M22, M23, M24;
        float M31, M32, M33, M34;
        float M41, M42, M43, M44;
    };

    struct Position
    {
        float X;
        float Y;
        float Z;

        std::array<float, 3> AsArray() const
        {
            return{ X, Y, Z };
        }
    };

    struct ProjectedPoint
    {
        float X;
        float Y;
        float Depth;
        int32_t Id;

        ProjectedPoint(float x, float y, float depth, int32_t id)
            : X{ x }, Y{ y }, Depth{ depth }, Id{ id }
        {}

        ProjectedPoint() = default;
    };

    struct AxisAlignedVolume
    {
        Position Min;
        Position Max;
    };

    struct Direction
    {
        float X;
        float Y;
        float Z;
    };

    enum class DataFormat
    {
        BINARY,
        PNG,
        LIVE,
        BOB,
        MIDDLEBURY,
    };

    namespace calibration
    {
        enum class DistortionType
        {
            None,
            Poly3k,
            Rational6k
        };

        class CameraModel
        {
        public:

            CameraModel() = default;
            CameraModel(const Intrinsics& intrinsics, DistortionType distType) : m_intrinsics(intrinsics), m_distortionType(distType) {}
            const Intrinsics& GetIntrinsics() const { return m_intrinsics; }

            //not in opencv order
            virtual gsl::span<const float> GetDistortionCoefficients() const = 0;
            virtual DistortionType GetDistortionType() const { return m_distortionType; }

        private:
            Intrinsics m_intrinsics;
            DistortionType m_distortionType;
        };

        class PinholeCameraModel : public CameraModel
        {
        public:
            PinholeCameraModel() = delete;
            PinholeCameraModel(const Intrinsics& intrinsics) : CameraModel(intrinsics, DistortionType::None) {}

            gsl::span<const float> GetDistortionCoefficients() const { return {}; }
        };

        class Poly3KCameraModel : public CameraModel
        {
        public:
            Poly3KCameraModel() = delete;
            Poly3KCameraModel(const Intrinsics& intrinsics, std::vector<float> distortionCoefficients) :
                CameraModel(intrinsics, DistortionType::Poly3k),
                m_distortionCoefficients{ distortionCoefficients }
            {
                assert(m_distortionCoefficients.size()== 5 && "incorrect number of distortion coefficients. want k1, k2, k3, p1, p2");
            }

            //not in opencv order
            gsl::span<const float> GetDistortionCoefficients() const { return m_distortionCoefficients; }

            float GetK1() const { return m_distortionCoefficients[0];}
            float GetK2() const { return m_distortionCoefficients[1];}
            float GetK3() const { return m_distortionCoefficients[2];}
            float GetP1() const { return m_distortionCoefficients[3];}
            float GetP2() const { return m_distortionCoefficients[4];}

            void SetK1(float k1) { m_distortionCoefficients[0] = k1; }
            void SetK2(float k2) { m_distortionCoefficients[1] = k2; }
            void SetK3(float k3) { m_distortionCoefficients[2] = k3; }
            void SetP1(float p1) { m_distortionCoefficients[3] = p1; }
            void SetP2(float p2) { m_distortionCoefficients[4] = p2; }

        private:
            // k1, k2, k3, p1, p2
            std::vector<float> m_distortionCoefficients;
        };

        class Rational6KCameraModel : public CameraModel
        {
        public:
            Rational6KCameraModel() = delete;
            Rational6KCameraModel(const Intrinsics& intrinsics, std::vector<float> distortionCoefficients)
                : CameraModel(intrinsics, DistortionType::Rational6k),
                m_distortionCoefficients{ distortionCoefficients }
            {
                assert(m_distortionCoefficients.size() == 8 && "incorrect number of distortion coefficients. want k1, k2, k3, k4, k5, k6, p1, p2");
            }

            //not in opencv order
            gsl::span<const float> GetDistortionCoefficients() const { return m_distortionCoefficients; }

            float GetK1() const { return m_distortionCoefficients[0]; }
            float GetK2() const { return m_distortionCoefficients[1]; }
            float GetK3() const { return m_distortionCoefficients[2]; }
            float GetK4() const { return m_distortionCoefficients[3]; }
            float GetK5() const { return m_distortionCoefficients[4]; }
            float GetK6() const { return m_distortionCoefficients[5]; }
            float GetP1() const { return m_distortionCoefficients[6]; }
            float GetP2() const { return m_distortionCoefficients[7]; }

            void SetK1(float k1) { m_distortionCoefficients[0] = k1; }
            void SetK2(float k2) { m_distortionCoefficients[1] = k2; }
            void SetK3(float k3) { m_distortionCoefficients[2] = k3; }
            void SetK4(float k4) { m_distortionCoefficients[3] = k4; }
            void SetK5(float k5) { m_distortionCoefficients[4] = k5; }
            void SetK6(float k6) { m_distortionCoefficients[5] = k6; }
            void SetP1(float p1) { m_distortionCoefficients[6] = p1; }
            void SetP2(float p2) { m_distortionCoefficients[7] = p2; }

        private:
            // k1, k2, k3, k4, k5, k6, p1, p2
            std::vector<float> m_distortionCoefficients;
        };

        struct Line
        {
            float M;
            float B;
        };

        struct Bounds
        {
            float Lower;
            float Upper;
        };

        class LinearFocalLengthModel
        {
        public:

            LinearFocalLengthModel(Line fx, Line fy, float cx, float cy, Bounds focalbounds, Size calibrationSize,
                const std::vector<float>& distortionPoly3k = {},
                const std::vector<float>& distortionRational6k = {}) :
                m_fx(fx), m_fy(fy), m_cx(cx), m_cy(cy), m_focalBounds(focalbounds), m_calibrationSize(calibrationSize)
            {
                assert(distortionPoly3k.empty() || distortionRational6k.empty() && "can't pass two types currently. calibrations have different intrinsics depending on model");
                assert(cx >= 0 && cx <= 1 && "cx and cy need to be a ratio");
                assert(cy >= 0 && cy <= 1 && "cx and cy need to be a ratio");

                if (!distortionPoly3k.empty())
                {
                    assert(distortionPoly3k.size() == 5 && "expecting k0, k1, k2, p0, p1, for poly3k distortion");
                    std::copy(distortionPoly3k.begin(), distortionPoly3k.end(), std::back_inserter(m_distortionPoly3k));
                }

                if (!distortionRational6k.empty())
                {
                    assert(distortionRational6k.size() == 8 && "expecting k0, k1, k2, k3, k4, k5, p0, p1, for rational6k distortion");
                    std::copy(distortionRational6k.begin(), distortionRational6k.end(), std::back_inserter(m_distortionRational6k));
                }
            }

            LinearFocalLengthModel()
            {
                m_fx = { 0,0 };
                m_fy = { 0,0 };
                m_cx = 0;
                m_cy = 0;
                m_focalBounds = { 0,0 };
                m_calibrationSize = { 0, 0 };
            }

            Intrinsics CreateIntrinsics(boost::optional<uint32_t> lensPosition, size_t width, size_t height) const
            {
                assert(width / (float)height == m_calibrationSize.Width / (float)m_calibrationSize.Height && "aspect not equal need to crop to modify resolution");
                assert(lensPosition || (m_fx.M == 0 && m_fy.M == 0));

                float fx = m_fx.B;
                float fy = m_fy.B;

                if (lensPosition)
                {
                    fx += m_fx.M * lensPosition.value();
                    fy += m_fy.M * lensPosition.value();
                }

                float FxInPixelCoordinates = fx * width;
                float FyInPixelCoordinates = fy * height;
                float CxInPixelCoordinates = m_cx * width;
                float CyInPixelCoordinates = m_cy * height;

                return Intrinsics(CxInPixelCoordinates, CyInPixelCoordinates,
                    FxInPixelCoordinates, FyInPixelCoordinates,
                    gsl::narrow<uint32_t>(width), gsl::narrow<uint32_t>(height));
            }

            std::shared_ptr<const PinholeCameraModel> CreatePinholeCameraModel(const boost::optional<uint32_t>& lensPosition, size_t width, size_t height) const
            {
                return std::make_shared<const PinholeCameraModel>(CreateIntrinsics(lensPosition, width, height));
            }

            std::shared_ptr<const Poly3KCameraModel> CreatePoly3kCameraModel(const boost::optional<uint32_t>& lensPosition, size_t width, size_t height) const
            {
                return std::make_shared<const Poly3KCameraModel>(CreateIntrinsics(lensPosition, width, height), m_distortionPoly3k);
            }

            std::shared_ptr<const Rational6KCameraModel> CreateRational6kCameraModel(const boost::optional<uint32_t>& lensPosition, size_t width, size_t height) const
            {
                return std::make_shared<const Rational6KCameraModel>(CreateIntrinsics(lensPosition, width, height), m_distortionRational6k);
            }

            std::shared_ptr<const CameraModel> GetCameraModel(const boost::optional<uint32_t>& lensPosition, size_t width, size_t height)
            {
                if (HasPoly3kModel())
                {
                    return CreatePoly3kCameraModel(lensPosition, width, height);
                }
                else if (HasRational6kModel())
                {
                    return CreateRational6kCameraModel(lensPosition, width, height);
                }
                else
                {
                    return CreatePinholeCameraModel(lensPosition, width, height);
                }
            }

            Line GetFx() const { return m_fx; }
            Line GetFy() const { return m_fy; }
            float GetCx() const { return m_cx; }
            float GetCy() const { return m_cy; }
            Bounds GetFocalBounds() const { return m_focalBounds; }
            Size GetCalibrationSize() const { return m_calibrationSize; }

            bool HasPoly3kModel()
            {
                return m_distortionPoly3k.size() > 0;
            }

            bool HasRational6kModel()
            {
                return m_distortionRational6k.size() > 0;
            }

            std::vector<float> GetDistortionPoly3k() const { return m_distortionPoly3k; }
            std::vector<float> GetDistortionRational6k() const { return m_distortionRational6k; }

        private:
            //F = m*focusValue + b
            Line m_fx, m_fy;
            float m_cx, m_cy;
            Bounds m_focalBounds;
            Size m_calibrationSize;

            //distortion
            std::vector<float> m_distortionPoly3k;
            std::vector<float> m_distortionRational6k;
        };
    }

    struct Depth
    {
        static constexpr float INVALID_DEPTH = -1.0f;

        float NearPlaneDepth;
        float FarPlaneDepth;
        gsl::span<const ProjectedPoint> SparseDepth;

        Depth(float nearDepth, float farDepth, gsl::span<const ProjectedPoint> sparse)
            : NearPlaneDepth{ nearDepth }, FarPlaneDepth{ farDepth }, SparseDepth{ sparse }
        {}

        Depth()
            : NearPlaneDepth{ INVALID_DEPTH }, FarPlaneDepth{ INVALID_DEPTH }
        {}
    };
}

max_stars_repo_name: syntheticmagus/mageslam

max_issues_repo_name: syntheticmagus/mageslam

max_forks_repo_name: syntheticmagus/mageslam


text:
/* **************************************** *
 *		ODE integration functions			*
 *	Provides GSL ode integrator and			*
 *	a simple euler integrator				*
 * **************************************** */

#include "canonical.h"
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv.h>
#define DIM 2
#define NPTS 100.0 // number of points to print out, best as double notation

/* Many good algorithms don't need the Jacobian, so we'll just pass a null pointer */
double *jac = NULL;

void
gslode(void *params, double MAX_TIME, FILE *theory)
{
	PAR pars = *(PAR *)params;

	/* Create our ODE system */
	gsl_odeiv_system sys = {func, jac, DIM, params};

	/* Range of time to simulate*/
	double t = 0.0, t1 = MAX_TIME;
	/* Initial step size, will be modified as needed by adaptive alogorithm */
	double h = 1e-6;
	/* initial conditions */
	double y[DIM] = { pars.Xo, 0.0 };

	/* Define method as Embedded Runge-Kutta Prince-Dormand (8,9) method */
	const gsl_odeiv_step_type * T
	 = gsl_odeiv_step_rk4;
//	 = gsl_odeiv_step_rk8pd;
	/* allocate stepper for our method of correct dimension*/
	gsl_odeiv_step * s
	 = gsl_odeiv_step_alloc (T, DIM);
	/* control will maintain absolute and relative error */
	gsl_odeiv_control * c
	 = gsl_odeiv_control_y_new (1e-6, 0.0);
	/* allocate the evolver */
	gsl_odeiv_evolve * e
	 = gsl_odeiv_evolve_alloc (DIM);

	/*dummy to make easy to switch to regular printing */
	double ti = t1;
	int i;

	/* Uncomment the outer for loop to print *
	 * 100 sample pts at regular intervals   */
	for (i = 1; i <= NPTS; i++){   ti = i * t1 / NPTS;
		while (t < ti)
		{
			int status = gsl_odeiv_evolve_apply (e, c, s,
												&sys,
												&t, t1,
												&h, y);
			if (status != GSL_SUCCESS)
			   break;
//			fprintf(theory,"%.6e %.6e %.6e\n",t,y[0],y[1]); //adaptive mesh
		}
		fprintf(theory,"%.6e %.6e %.6e\n",t,y[0],y[1]);  }

	gsl_odeiv_evolve_free (e);
	gsl_odeiv_control_free (c);
	gsl_odeiv_step_free (s);
}

void euler(void *params, double MAX_TIME, FILE *theory){
	PAR pars = *(PAR *)params;

	double dt = .01;
	int sample;

	double y[DIM] = { pars.Xo, 0.0 };
	double f[DIM];
	double t=0.0;

	for(sample = 0; sample < MAX_TIME/dt; sample++){

		if( sample % (int) (1/dt) == 0)
		{
			fprintf
			(
				theory, "%.6e %.6e %.6e\n",
				dt*sample,
				y[0],
				y[1]
			);
		}
		int status = func(t, y, f, &pars);
		if(status != GSL_SUCCESS)
			break;
		y[0] = y[0] + dt*f[0];
		y[1] = y[1] + dt*f[1];
	}

}

max_stars_repo_name: cboettig/fluctuationDomains

max_issues_repo_name: cboettig/fluctuationDomains

max_forks_repo_name: cboettig/fluctuationDomains


text:
#ifndef PRIMAL_DUAL_LOOPLESS_KATYUSHA0_H
#define PRIMAL_DUAL_LOOPLESS_KATYUSHA0_H

#include "Primal_Dual_LOOPLESS.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 <sstream>

//This class implements the method loopless kATYUSHA with arbitrary sampling

/*
The optimization problem to solve is:

\sum_{i=1}^n \lambda_i\phi_i(A_i^{\top} w)+ g(w)
Assumption 1: For each i, \phi_i is 1-smooth
By default, lambda_i=1/n for all i.
*
* The dual problem is
*         g*(-sum_{i=1}^n \lambda_i \alpha_i A_i)+\sum_{i=1}^n \lambda_i \phi*_i(\alpha_i)
*/

template<typename L, typename D>
class Primal_Dual_LOOPLESS_Katyusha0: public Primal_Dual_LOOPLESS<L, D>
{

protected:
   D theta1;
   D theta2;
   D theta3;
   D eta; //the stepsize parameter alpha in the paper

   D x;
   D w;
   D y;
   D z;

   std::vector<D> primal_y;     // the y variable

   std::vector<D> primal_z;     // the z variable

   std::vector<D> primal_w;     // the w variable

public:

  virtual inline D gradient_of_phi_i(D, L){return D(NULL);}
  virtual inline D gradient_of_gstar_j(D, L){return D(NULL);}
  virtual inline D value_of_phi_i(D, L) {return D(NULL);}
  virtual inline D value_of_g_j(D, L){return D(NULL);}
  virtual inline D value_of_phistar_i(D,L) {return D(NULL);}
  virtual inline D value_of_gstar(vector<D> &){return D(NULL);}
  virtual inline D value_of_gstar_minus(vector<D> &, D){return D(NULL);}
  virtual inline D feasible_dual(vector<D> &){return D(NULL);}
  virtual inline D compute_delta_alpha(D,D,L){return D(NULL);}
  virtual inline void set_auxiliary_v(){}
  virtual inline D get_lambda1(){return D(NULL);}
  virtual inline D get_lambda2(){return D(NULL);}
  virtual inline void compute_just_in_time_prox_grad(D, D, D &, L,L, D, D &, L){}

  Primal_Dual_LOOPLESS_Katyusha0()
  : Primal_Dual_LOOPLESS<L,D>()
  {

  }

  Primal_Dual_LOOPLESS_Katyusha0(const char* matrix_file, const char* vector_file)
  : Primal_Dual_LOOPLESS<L,D>(matrix_file, vector_file)
  {

  }

  Primal_Dual_LOOPLESS_Katyusha0(const char* matrix_file)
  : Primal_Dual_LOOPLESS<L,D>(matrix_file)
  {

  }

  void set_theta(){
    theta2=this->L2/2/max(this->Lf,this->L2);
    if(this->Lf<=this->L2/this->p){
      D tmp=this->mu/this->L2/this->p;
      cout<<"this->L2="<<this->L2<<endl;
      cout<<"tmp="<<tmp<<endl;
      if(tmp>=1){
        theta1=theta2;
      }
      else{
        theta1=sqrt(tmp)*theta2;
      }
    }else
      theta1=min(sqrt(this->mu/this->Lf),this->p/2);
    theta3=1-theta1-theta2;
    eta=1./(theta1*(this->Lf+2*max(this->L2,this->Lf)));
    cout<<"this->Lf="<<this->Lf<<"; this->L2="<<this->L2<<"; theta1="<<theta1<<"; theta2="<<theta2<<"; theta3="<<theta3<<endl;
    cout<<"eta="<<eta<<endl;
  }

  inline D compute_current_xj_value(D ui, L j, L t0 , L t1)
  {
    w=primal_w[j];
    y=primal_y[j];
    z=primal_z[j];
    compute_just_in_time_prox_grad(eta, ui , z, t0, t1, w, y, j);
    x=theta1*z+theta2*w+theta3*y;
    return x;
  }

  inline void set_stepsize(){
    primal_y.clear();
    primal_z.clear();
    primal_w.clear();
    primal_y.resize(this->nfeatures,0);
    primal_z.resize(this->nfeatures,0);
    primal_w.resize(this->nfeatures,0);
    for(L j=0; j<this->nfeatures; j++){
      primal_y[j]=this->primal_x[j];
      primal_z[j]=this->primal_x[j];
      primal_w[j]=this->primal_x[j];
    }
    cout<<"w size="<<primal_w.size()<<endl;
    set_theta();
    cout<<"w 2="<<primal_w.size()<<endl;
  }

  inline void update_x(D xj,L j){
    this->primal_x[j]=xj;
    primal_z[j]=z;
    primal_y[j]=y;
  }

  inline void update_baralpha()
   {
     D yi=gsl_rng_uniform(this->rng);
     if(yi<=this->p){
         this->nb_iters+=this->nsamples/this->tau;
         D aitg=0;
         D deltaalphai=0;
         for(L i=0;i<this->nsamples;i++)
         {
           aitg=this->compute_AiTxk(i);
           deltaalphai=gradient_of_phi_i(aitg,i)-this->dual_alpha[i];
           this->dual_alpha[i]+=deltaalphai;
           for (L k = this->ptr[i]; k < this->ptr[i + 1];k++)
           {
             L j=this->row_idx[k];
             this->baralpha[j]+=this->lambda_f[i]*deltaalphai*this->A[k];
           }
         }

       for(L j=0;j<this->nfeatures;j++)
       primal_w[j]=this->primal_x[j];
     }
   }

  void loopless_Katyusha(vector<D> & x0, vector<D> & w0, string filename, vector<D> & L_phi, D val_mu, D val_epsilon, L max_nb, L nb_tau, L nb_c, L u, L p_mod, D scal_p)
  {
    filename="Katyusha"+filename;
    this->loopless( x0,  w0,  filename,  L_phi,  val_mu,  val_epsilon, max_nb, nb_tau, nb_c, u, p_mod,scal_p);
  }

  void loopless_Katyusha2(vector<D> & x0, vector<D> & w0, string filename, vector<D> & L_phi, D val_mu, D val_epsilon, L max_nb, L nb_tau, L nb_c, L u, L p_mod, D scal_p)
  {
    filename="Katyusha"+filename;
    this->loopless2( x0,  w0,  filename,  L_phi,  val_mu,  val_epsilon, max_nb, nb_tau, nb_c, u, p_mod,scal_p);
  }

  };

  #endif /* MIN_SMOOTH_CONVEX_H */


text:
#pragma once
#include "PreConfig.h"

#ifdef CRU_PLATFORM_WINDOWS
#ifdef CRU_BASE_EXPORT_API
#define CRU_BASE_API __declspec(dllexport)
#else
#define CRU_BASE_API __declspec(dllimport)
#endif
#else
#define CRU_BASE_API
#endif

#include <gsl/gsl>

#define CRU_UNUSED(entity) static_cast<void>(entity);

#define CRU__CONCAT(a, b) a##b
#define CRU_MAKE_UNICODE_LITERAL(str) CRU__CONCAT(u, #str)

#define CRU_DEFAULT_COPY(classname)      \
  classname(const classname&) = default; \
  classname& operator=(const classname&) = default;

#define CRU_DEFAULT_MOVE(classname) \
  classname(classname&&) = default; \
  classname& operator=(classname&&) = default;

#define CRU_DELETE_COPY(classname)      \
  classname(const classname&) = delete; \
  classname& operator=(const classname&) = delete;

#define CRU_DELETE_MOVE(classname) \
  classname(classname&&) = delete; \
  classname& operator=(classname&&) = delete;

#define CRU_DEFAULT_DESTRUCTOR(classname) ~classname() = default;

#define CRU_DEFAULT_CONSTRUCTOR_DESTRUCTOR(classname) \
  classname() = default;                              \
  ~classname() = default;

#define CRU_DEFINE_COMPARE_OPERATORS(classname)                           \
  inline bool operator==(const classname& left, const classname& right) { \
    return left.Compare(right) == 0;                                      \
  }                                                                       \
                                                                          \
  inline bool operator!=(const classname& left, const classname& right) { \
    return left.Compare(right) != 0;                                      \
  }                                                                       \
                                                                          \
  inline bool operator<(const classname& left, const classname& right) {  \
    return left.Compare(right) < 0;                                       \
  }                                                                       \
                                                                          \
  inline bool operator<=(const classname& left, const classname& right) { \
    return left.Compare(right) <= 0;                                      \
  }                                                                       \
                                                                          \
  inline bool operator>(const classname& left, const classname& right) {  \
    return left.Compare(right) > 0;                                       \
  }                                                                       \
                                                                          \
  inline bool operator>=(const classname& left, const classname& right) { \
    return left.Compare(right) >= 0;                                      \
  }

namespace cru {
class CRU_BASE_API Object {
 public:
  Object() = default;
  CRU_DEFAULT_COPY(Object)
  CRU_DEFAULT_MOVE(Object)
  virtual ~Object() = default;
};

struct CRU_BASE_API Interface {
  Interface() = default;
  CRU_DELETE_COPY(Interface)
  CRU_DELETE_MOVE(Interface)
  virtual ~Interface() = default;
};

[[noreturn]] void CRU_BASE_API UnreachableCode();

using Index = gsl::index;

// https://www.boost.org/doc/libs/1_54_0/doc/html/hash/reference.html#boost.hash_combine
template <class T>
inline void hash_combine(std::size_t& s, const T& v) {
  std::hash<T> h;
  s ^= h(v) + 0x9e3779b9 + (s << 6) + (s >> 2);
}

#define CRU_DEFINE_CLASS_LOG_TAG(tag) \
 private:                             \
  constexpr static const char16_t* kLogTag = tag;
}  // namespace cru

max_stars_repo_name: crupest/Cru

max_issues_repo_name: crupest/Cru

max_forks_repo_name: crupest/Cru