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

#include <config.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_vector_complex.h>
#include <gsl/gsl_matrix.h>

/*
 * This module computes the eigenvalues of a real nonsymmetric
 * matrix, using the double shift Francis method.
 *
 * See the references in francis.c.
 *
 * This module gets the matrix ready by balancing it and
 * reducing it to Hessenberg form before passing it to the
 * francis module.
 */

/*
gsl_eigen_nonsymm_alloc()

Allocate a workspace for solving the nonsymmetric eigenvalue problem.
The size of this workspace is O(2n)

Inputs: n - size of matrix

Return: pointer to workspace
*/

gsl_eigen_nonsymm_workspace *
gsl_eigen_nonsymm_alloc(const size_t n)
{
  gsl_eigen_nonsymm_workspace *w;

  if (n == 0)
    {
      GSL_ERROR_NULL ("matrix dimension must be positive integer",
                      GSL_EINVAL);
    }

  w = (gsl_eigen_nonsymm_workspace *)
      malloc (sizeof (gsl_eigen_nonsymm_workspace));

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

  w->size = n;
  w->Z = NULL;
  w->do_balance = 0;

  w->diag = gsl_vector_alloc(n);

  if (w->diag == 0)
    {
      GSL_ERROR_NULL ("failed to allocate space for balancing vector", GSL_ENOMEM);
    }

  w->tau = gsl_vector_alloc(n);

  if (w->tau == 0)
    {
      GSL_ERROR_NULL ("failed to allocate space for hessenberg coefficients", GSL_ENOMEM);
    }

  w->francis_workspace_p = gsl_eigen_francis_alloc();

  if (w->francis_workspace_p == 0)
    {
      GSL_ERROR_NULL ("failed to allocate space for francis workspace", GSL_ENOMEM);
    }

  return (w);
} /* gsl_eigen_nonsymm_alloc() */

/*
gsl_eigen_nonsymm_free()
  Free workspace w
*/

void
gsl_eigen_nonsymm_free (gsl_eigen_nonsymm_workspace * w)
{
  gsl_vector_free(w->tau);

  gsl_vector_free(w->diag);

  gsl_eigen_francis_free(w->francis_workspace_p);

  free(w);
} /* gsl_eigen_nonsymm_free() */

/*
gsl_eigen_nonsymm_params()
  Set some parameters which define how we solve the eigenvalue
problem.

Inputs: compute_t - 1 if we want to compute T, 0 if not
        balance   - 1 if we want to balance the matrix, 0 if not
        w         - nonsymm workspace
*/

void
gsl_eigen_nonsymm_params (const int compute_t, const int balance,
                          gsl_eigen_nonsymm_workspace *w)
{
  gsl_eigen_francis_T(compute_t, w->francis_workspace_p);
  w->do_balance = balance;
} /* gsl_eigen_nonsymm_params() */

/*
gsl_eigen_nonsymm()

Solve the nonsymmetric eigenvalue problem

A x = \lambda x

for the eigenvalues \lambda using the Francis method.

Here we compute the real Schur form

T = Z^t A Z

with the diagonal blocks of T giving us the eigenvalues.
Z is a matrix of Schur vectors which is not computed by
this algorithm. See gsl_eigen_nonsymm_Z().

Inputs: A    - general real matrix
        eval - where to store eigenvalues
        w    - workspace

Return: success or error

Notes: If T is computed, it is stored in A on output. Otherwise
       the diagonal of A contains the 1-by-1 and 2-by-2 eigenvalue
       blocks.
*/

int
gsl_eigen_nonsymm (gsl_matrix * A, gsl_vector_complex * eval,
                   gsl_eigen_nonsymm_workspace * w)
{
  const size_t N = A->size1;

  /* check matrix and vector sizes */

  if (N != A->size2)
    {
      GSL_ERROR ("matrix must be square to compute eigenvalues", GSL_ENOTSQR);
    }
  else if (eval->size != N)
    {
      GSL_ERROR ("eigenvalue vector must match matrix size", GSL_EBADLEN);
    }
  else
    {
      int s;

      if (w->do_balance)
        {
          /* balance the matrix */
          gsl_linalg_balance_matrix(A, w->diag);
        }

      /* compute the Hessenberg reduction of A */
      gsl_linalg_hessenberg(A, w->tau);

      if (w->Z)
        {
          /*
           * initialize the matrix Z to U, which is the matrix used
           * to construct the Hessenberg reduction.
           */

          /* compute U and store it in Z */
          gsl_linalg_hessenberg_unpack(A, w->tau, w->Z);

          /* find the eigenvalues and Schur vectors */
          s = gsl_eigen_francis_Z(A, eval, w->Z, w->francis_workspace_p);

          if (w->do_balance)
            {
              /*
               * The Schur vectors in Z are the vectors for the balanced
               * matrix. We now must undo the balancing to get the
               * vectors for the original matrix A.
               */
              gsl_linalg_balance_accum(w->Z, w->diag);
            }
        }
      else
        {
          /* find the eigenvalues only */
          s = gsl_eigen_francis(A, eval, w->francis_workspace_p);
        }

      w->n_evals = w->francis_workspace_p->n_evals;

      return s;
    }
} /* gsl_eigen_nonsymm() */

/*
gsl_eigen_nonsymm_Z()

Solve the nonsymmetric eigenvalue problem

A x = \lambda x

for the eigenvalues \lambda.

Here we compute the real Schur form

T = Z^t A Z

with the diagonal blocks of T giving us the eigenvalues.
Z is the matrix of Schur vectors.

Inputs: A    - general real matrix
        eval - where to store eigenvalues
        Z    - where to store Schur vectors
        w    - workspace

Return: success or error

Notes: If T is computed, it is stored in A on output. Otherwise
       the diagonal of A contains the 1-by-1 and 2-by-2 eigenvalue
       blocks.
*/

int
gsl_eigen_nonsymm_Z (gsl_matrix * A, gsl_vector_complex * eval,
                     gsl_matrix * Z, gsl_eigen_nonsymm_workspace * w)
{
  /* check matrix and vector sizes */

  if (A->size1 != A->size2)
    {
      GSL_ERROR ("matrix must be square to compute eigenvalues", GSL_ENOTSQR);
    }
  else if (eval->size != A->size1)
    {
      GSL_ERROR ("eigenvalue vector must match matrix size", GSL_EBADLEN);
    }
  else if ((Z->size1 != Z->size2) || (Z->size1 != A->size1))
    {
      GSL_ERROR ("Z matrix has wrong dimensions", GSL_EBADLEN);
    }
  else
    {
      int s;

      w->Z = Z;

      s = gsl_eigen_nonsymm(A, eval, w);

      w->Z = NULL;

      return s;
    }
} /* gsl_eigen_nonsymm_Z() */


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

#include <config.h>
#include <math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

/* Logarithmic distribution

   prob(n) =   p^n / (n log(1/(1-p)) for n = 1, 2, 3, ...

   We use Kemp's second accelerated generator, from Luc Devroye's book
   on "Non-Uniform Random Variate Generation", Springer */

unsigned int
gsl_ran_logarithmic (const gsl_rng * r, const double p)
{
  double c = log (1-p) ;

  double v = gsl_rng_uniform_pos (r);

  if (v >= p)
    {
      return 1 ;
    }
  else
    {
      double u = gsl_rng_uniform_pos (r);
      double q = 1 - exp (c * u);

      if (v <= q*q)
	{
	  double x = 1 + log(v)/log(q) ;
	  return x ;
	}
      else if (v <= q)
	{
	  return 2;
	}
      else
	{
	  return 1 ;
	}
    }
}

double
gsl_ran_logarithmic_pdf (const unsigned int k, const double p)
{
  if (k == 0)
    {
      return 0 ;
    }
  else
    {
      double P = pow(p, (double)k) / (double) k / log(1/(1-p)) ;
      return P;
    }
}

ext: c


text:
#include <gsl/gsl_combination.h>
#include <gsl/gsl_errno.h>
#include <stdio.h>

int mgsl_combination_fwrite(const char *filename, const gsl_combination *p)
{
  FILE *fp;
  if((fp = fopen(filename, "w")) == NULL) return GSL_EFAILED;
  if(gsl_combination_fwrite(fp, p) != GSL_SUCCESS) return GSL_EFAILED;
  fclose(fp);
  return GSL_SUCCESS;
}

int mgsl_combination_fread(const char *filename, gsl_combination *p)
{
  FILE *fp;
  if((fp = fopen(filename, "r")) == NULL) return GSL_EFAILED;
  if(gsl_combination_fread(fp, p) != GSL_SUCCESS) return GSL_EFAILED;
  fclose(fp);
  return GSL_SUCCESS;
}

int mgsl_combination_fprintf(const char *filename, const gsl_combination *p, const char *format)
{
  FILE *fp;
  if((fp = fopen(filename, "w")) == NULL) return GSL_EFAILED;
  if(gsl_combination_fprintf(fp, p, format) != GSL_SUCCESS) return GSL_EFAILED;
  fclose(fp);
  return GSL_SUCCESS;
}

int mgsl_combination_fscanf(const char *filename, gsl_combination *p)
{
  FILE *fp;
  if((fp = fopen(filename, "r")) == NULL) return GSL_EFAILED;
  if(gsl_combination_fscanf(fp, p) != GSL_SUCCESS) return GSL_EFAILED;
  fclose(fp);
  return GSL_SUCCESS;
}

max_stars_repo_name: frithnanth/raku-Math-Libgsl-Combination

max_issues_repo_name: frithnanth/raku-Math-Libgsl-Combination

max_forks_repo_name: frithnanth/raku-Math-Libgsl-Combination


text:
/**
 * @file square_root.h
 * @brief Updating complex Cholesky factorizations.
 * @author John McDonough and Kenichi Kumatani
 */

#ifndef SQUARE_ROOT_H
#define SQUARE_ROOT_H

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

#include <gsl/gsl_block.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include "common/jexception.h"

void cholesky_backsub(const gsl_matrix* A, gsl_vector* x);

void vector_matrix_product(const gsl_vector_complex* vec,
			 const gsl_matrix_complex* mat, gsl_matrix* D);

void make_conjugate_symmetric(gsl_matrix_complex* mat);

void cholesky_forwardsub(const gsl_matrix* A, gsl_vector* x);

void cholesky_forwardsub_complex(const gsl_matrix_complex* lt,
			       const gsl_vector_complex* rhs,
			       gsl_vector_complex* lhs,
			       bool conjugate = false);

void cholesky_backsub_complex(const gsl_matrix_complex* lt,
                              const gsl_vector_complex* rhs,
                              gsl_vector_complex* lhs,
                              bool conjugate = false);

void rank_one_update_cholesky_factor(gsl_matrix_complex* A11,
                                     const double alpha_m,
                                     const gsl_vector_complex* c_m);

void propagate_covar_square_root_real(gsl_matrix* A11,
                                      gsl_matrix* A12,
                                      gsl_matrix* A21,
                                      gsl_matrix* A22, bool flag = false);

void sweep_lower_triangular(gsl_matrix* A, gsl_matrix* B);

void propagate_covar_square_root_step1(gsl_matrix_complex* A12,
                                       gsl_matrix_complex* A22);

void propagate_covar_square_root_step2a(gsl_matrix_complex* A11,
                                        gsl_matrix_complex* A12,
                                        gsl_matrix_complex* A21,
                                        gsl_matrix_complex* A22);

void propagate_covar_square_root_step2b(gsl_matrix_complex* A22);

void propagate_covar_square_root(gsl_matrix_complex* A11,
                                 gsl_matrix_complex* A12,
                                 gsl_matrix_complex* A21,
                                 gsl_matrix_complex* A22);

void propagate_info_square_root(gsl_matrix_complex* sqrt_Pm_inv,
                                gsl_matrix_complex* A12,
                                gsl_vector_complex* a_21,
                                gsl_vector_complex* a_22, bool rankOneA12 = true);

void propagate_info_square_root_step2_rls(gsl_matrix_complex* sqrt_Pm_inv,
                                          gsl_vector_complex* a_12,
                                          gsl_vector_complex* a_21,
                                          gsl_complex a_22);

void propagate_info_square_root_rls(gsl_matrix_complex* sqrt_Pm_inv,
                                    gsl_vector_complex* a_12,
                                    gsl_vector_complex* a_21,
                                    gsl_complex a_22);

void add_diagonal_loading(gsl_matrix_complex* sqrt_Pm_inv, int dim, double wght);

gsl_vector* cholesky_diagonal(gsl_vector* v, const gsl_matrix* m);
gsl_vector* square_diagonal(gsl_vector* v, const gsl_matrix* m);

#endif

ext: h


text:
#include <assert.h>
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <stdint.h>
#include <string.h>

#include <gsl/gsl_math.h>

#include "gaussianInt2.h"
#include "blit3.h"

#ifndef verbose
#define verbose 1
#endif

double rand_range(double A, double B)
  // A random number between A and B
{
  if(A==B)
    return A;
  if(B>A)
    return A + (double) rand() / RAND_MAX * (B-A);
  //if(A>B)
  return B + (double) rand() / RAND_MAX * (A-B);
}

int blit3g(double * T, uint32_t Tm, uint32_t Tn, uint32_t Tp,
    double * D, size_t nD, photonw pw, int one_indexing)
  // T:  target volume
  // nD: number of dots
  // D:  list of dots (x,y,z, nphot, sigmax, sigmay, sigmaz)
{

  const int nF = 7; // stride for D

  double sigma_max = -1;
  for(size_t kk = 0; kk < nD; kk++)
    {
      for(int ll = 0; ll < 3; ll++)
      {
      if( D[nF*kk+4 + ll] > sigma_max)
        sigma_max = D[nF*kk+4 + ll];
      }
    }

  const size_t w = ceil(sigma_max)*3*2 + 1;
  const int hw = (w-1)/2;
  // allocate temporary space for the Gaussian
  double * G = malloc(w*w*w*sizeof(double));

  double sigma[] = {1.6,1.6,2.1};
  double mu[] = {0,0,0};
  double nPhot;

  for(uint32_t pp = 0; pp<nD; pp++)  {

    // Round the location
    int64_t Pm =(int64_t) (D[nF*pp+0]-.5);
    int64_t Pn =(int64_t) (D[nF*pp+1]-.5);
    int64_t Pp =(int64_t) (D[nF*pp+2]-.5);

    //    printf("%lu %lu %lu\n", Pm,Pn,Pp);

    // Get the offset from the centre of the pixel
    mu[0] = -roundf(D[nF*pp+0])+D[nF*pp+0];
    mu[1] = -roundf(D[nF*pp+1])+D[nF*pp+1];
    mu[2] = -roundf(D[nF*pp+2])+D[nF*pp+2];

    //    printf("%f %f %f\n", mu[0], mu[1], mu[1]);

    //   printf("mu: %f %f %f\n", mu[0], mu[1], mu[2]);
    sigma[0] = D[nF*pp+4];
    sigma[1] = D[nF*pp+5];
    sigma[2] = D[nF*pp+6];

    nPhot = D[nF*pp+3];
    // Create the kernel
    if(!gaussianInt3(G, mu,  sigma, w))
    {
      printf("Failed to allocate space for Gaussian kernel\n");
      return 1;
    }

    if(pw == xyz_volume)
      for(size_t tt =0; tt<w*w*w; tt++)
        G[tt] = G[tt]*nPhot;

    if(pw == xy_plane)
    {
      double wgt = 0;
      for(size_t tt =0; tt<w*w; tt++)
        wgt += G[tt+w*w*(w-1)/2];

      for(size_t tt =0; tt<w*w*w; tt++)
        G[tt] = G[tt]*nPhot/wgt;
    }

    if(pw == mid_point)
    {
      const double wgt = G[(w*w*w-1)/2];
      for(size_t tt =0; tt<w*w*w; tt++)
        G[tt] = G[tt]*nPhot/wgt;
    }

    // blit it
    if(one_indexing){
      blit3(T, Tm, Tn, Tp,
          G, w, w, w,
          Pm-hw, Pn-hw, Pp-hw,
          0);
    } else {
blit3(T, Tm, Tn, Tp,
          G, w, w, w,
          Pm-hw+1, Pn-hw+1, Pp-hw+1,
          0);
    }

  }
  free(G);
  return 0;
}

void blit3(double * T, uint32_t Tm, uint32_t Tn, uint32_t Tp,
    double * S, uint32_t Sm, uint32_t Sn, uint32_t Sp,
    int64_t Pm, int64_t Pn, int64_t Pp,
    int8_t anchor)
  // T: target image
{
  // anchor, 0=corner, 1=center of S
  assert(anchor == 0);
  // Which ranges in T to be visited
  // in [Pma, Pmb]
  uint32_t Tma = GSL_MAX(0, Pm);
  uint32_t Tmb = GSL_MIN(Tm-1, Pm+Sm-1);
  uint32_t Tna = GSL_MAX(0, Pn);
  uint32_t Tnb = GSL_MIN(Tn-1, Pn+Sn-1);
  uint32_t Tpa = GSL_MAX(0, Pp);
  uint32_t Tpb = GSL_MIN(Tp-1, Pp+Sp-1);

  uint32_t Sma = 0;
  if(Pm<0)
    Sma = -Pm;

  uint32_t Sna = 0;
  if(Pn<0)
    Sna = -Pn;

  uint32_t Spa = 0;
  if(Pp<0)
    Spa = -Pp;

#if verbose > 0

  uint32_t Smb = GSL_MIN(Sm-1, Tm-Pm-1);
  uint32_t Snb = GSL_MIN(Sn-1, Tn-Pn-1);
  uint32_t Spb = GSL_MIN(((int64_t) Sp)-1, Tp-Pp-1);
  printf("Size of T: %u x %u x %u\n", Tm, Tn, Tp);
  printf("Size of S: %u x %u x %u\n", Sm, Sn, Sp);
  printf("P: %ld %ld %ld\n", Pm, Pn, Pp);
  printf("Will visit:\n");
  printf("T: [%u %u]x[%u %u]x[%u %u] (%u)\n", Tma, Tmb, Tna, Tnb, Tpa, Tpb, (Tmb-Tma+1)*(Tnb-Tna+1)*(Tpb-Tpa+1));
  printf("S: [%u %u]x[%u %u]x[%u %u] (%u)\n", Sma, Smb, Sna, Snb, Spa, Spb, (Smb-Sma+1)*(Snb-Sna+1)*(Spb-Spa+1));
#endif

  uint32_t ppos, npos; // For pre-caluclation of multiplications

  uint32_t smm=Sma, snn=Sna, spp=Spa;

  for(uint32_t pp = Tpa; pp<=Tpb; pp++, spp++) {
    snn = Sna;

    ppos = pp*Tm*Tn;
    for(uint32_t nn = Tna; nn<=Tnb; nn++, snn++) {
      npos = nn*Tm;
      smm = Sma;

      for(uint32_t mm = Tma; mm<=Tmb; mm++, smm++) {
#if verbose > 1
        printf("T(%u %u %u) := S(%u %u %u)\n", mm, nn, pp, smm, snn, spp);
#endif
        T[mm + npos + ppos] += S[smm + snn*Sm + spp*Sm*Sn];

      }
    }
  }

}

#ifdef standalone
int unit_tests()
{
  printf("Testing \n");
  uint32_t Tm = 10;
  uint32_t Tn = 11;
  uint32_t Tp = 12;

  double * T = malloc(Tm*Tn*Tp*sizeof(double));
  for(size_t kk = 0; kk<Tm*Tn*Tp; kk++) {
    T[kk] = 0;
  }

  uint32_t Sm = 13;
  uint32_t Sn = 15;
  uint32_t Sp = 17;

  double * S = malloc(Sm*Sn*Sp*sizeof(double));

  for(size_t kk = 0; kk<Sm*Sn*Sp; kk++) {
    S[kk] = 1;
  }

  int64_t Pm = -7;
  int64_t Pn = -8;
  int64_t Pp = -9;

  blit3(T, Tm, Tn, Tp,   S, Sm, Sn, Sp,  Pm, Pn, Pp,  0);

  if(Pm>-1 && Pn>-1 && Pp>-1)
    if(Pm<Tm && Pn<Tn && Pp<Tp)
      printf("T(%ld,%ld,%ld) = %f\n", Pm, Pn, Pp, T[Pm + Pn*Tm + Pp*Tm*Tn]);

  free(S);
  free(T);

printf("--> blit3g\n");
size_t Gm = 11;
size_t Gn = 11;
size_t Gp = 11;
double * G = calloc(Gm*Gn*Gp,sizeof(double));
double * Dot = calloc(7,sizeof(double));
Dot[0] = 7; Dot[1] = 7; Dot[2] = 7;
Dot[3] = 100;
Dot[4] = 1; Dot[5]=1; Dot[6]=1;
blit3g(G, Gm, Gn, Gp, Dot, 1, mid_point, 0);

double maxg = 0;
for(size_t kk=0; kk<Gm*Gn*Gp; kk++)
  maxg = GSL_MAX(maxg, G[kk]);
printf("max(G): %f\n", maxg);
double midg = G[7+7*Gm+7*Gm*Gn];
printf("G(7,7,7)=%f\n", midg);
assert(midg == maxg);

free(G);
free(Dot);

  return 0;
}

int main(int argc, char ** argv)
{
  printf("%s\n", argv[0]);
  if(argc == 1)
    return unit_tests();
}
#endif

max_stars_repo_name: elgw/dotter

max_issues_repo_name: elgw/dotter

max_forks_repo_name: elgw/dotter


text:
/* normal.c
 *
 * Copyright (C) 2015, 2016 Patrick Alken
 *
 * 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_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_multifit.h>
#include <gsl/gsl_multilarge.h>

typedef struct
{
  size_t p;              /* number of columns of LS matrix */
  gsl_matrix *ATA;       /* A^T A, p-by-p */
  gsl_vector *ATb;       /* A^T b, p-by-1 */
  double normb;          /* || b || */
  gsl_matrix *work_ATA;  /* workspace for chol(ATA), p-by-p */
  gsl_vector *workp;     /* workspace size p */
  gsl_vector *work3p;    /* workspace size 3*p */
  gsl_vector *D;         /* scale factors for ATA, size p */
  gsl_vector *c;         /* solution vector for L-curve */
  int eigen;             /* 1 if eigenvalues computed */
  double eval_min;       /* minimum eigenvalue */
  double eval_max;       /* maximum eigenvalue */
  gsl_eigen_symm_workspace *eigen_p;
} normal_state_t;

static void *normal_alloc(const size_t p);
static void normal_free(void *vstate);
static int normal_reset(void *vstate);
static int normal_accumulate(gsl_matrix * A, gsl_vector * b,
                             void * vstate);
static int normal_solve(const double lambda, gsl_vector * x,
                        double * rnorm, double * snorm,
                        void * vstate);
static int normal_rcond(double * rcond, void * vstate);
static int normal_lcurve(gsl_vector * reg_param, gsl_vector * rho,
                         gsl_vector * eta, void * vstate);
static const gsl_matrix * normal_ATA(const void * vstate);
static const gsl_vector * normal_ATb(const void * vstate);
static int normal_solve_system(const double lambda, gsl_vector * x,
                               normal_state_t *state);
static int normal_solve_cholesky(gsl_matrix * ATA, const gsl_vector * ATb,
                                 gsl_vector * x, normal_state_t *state);
static int normal_calc_norms(const gsl_vector *x, double *rnorm,
                             double *snorm, normal_state_t *state);
static int normal_eigen(normal_state_t *state);

/*
normal_alloc()
  Allocate workspace for solving large linear least squares
problems using the normal equations approach

Inputs: p    - number of columns of LS matrix

Return: pointer to workspace
*/

static void *
normal_alloc(const size_t p)
{
  normal_state_t *state;

  if (p == 0)
    {
      GSL_ERROR_NULL("p must be a positive integer",
                     GSL_EINVAL);
    }

  state = calloc(1, sizeof(normal_state_t));
  if (!state)
    {
      GSL_ERROR_NULL("failed to allocate normal state", GSL_ENOMEM);
    }

  state->p = p;

  state->ATA = gsl_matrix_alloc(p, p);
  if (state->ATA == NULL)
    {
      normal_free(state);
      GSL_ERROR_NULL("failed to allocate ATA matrix", GSL_ENOMEM);
    }

  state->work_ATA = gsl_matrix_alloc(p, p);
  if (state->work_ATA == NULL)
    {
      normal_free(state);
      GSL_ERROR_NULL("failed to allocate temporary ATA matrix", GSL_ENOMEM);
    }

  state->ATb = gsl_vector_alloc(p);
  if (state->ATb == NULL)
    {
      normal_free(state);
      GSL_ERROR_NULL("failed to allocate ATb vector", GSL_ENOMEM);
    }

  state->D = gsl_vector_alloc(p);
  if (state->D == NULL)
    {
      normal_free(state);
      GSL_ERROR_NULL("failed to allocate D vector", GSL_ENOMEM);
    }

  state->workp = gsl_vector_alloc(p);
  if (state->workp == NULL)
    {
      normal_free(state);
      GSL_ERROR_NULL("failed to allocate temporary ATb vector", GSL_ENOMEM);
    }

  state->work3p = gsl_vector_alloc(3 * p);
  if (state->work3p == NULL)
    {
      normal_free(state);
      GSL_ERROR_NULL("failed to allocate work3p", GSL_ENOMEM);
    }

  state->c = gsl_vector_alloc(p);
  if (state->c == NULL)
    {
      normal_free(state);
      GSL_ERROR_NULL("failed to allocate c vector", GSL_ENOMEM);
    }

  state->eigen_p = gsl_eigen_symm_alloc(p);
  if (state->eigen_p == NULL)
    {
      normal_free(state);
      GSL_ERROR_NULL("failed to allocate eigen workspace", GSL_ENOMEM);
    }

  normal_reset(state);

  return state;
}

static void
normal_free(void *vstate)
{
  normal_state_t *state = (normal_state_t *) vstate;

  if (state->ATA)
    gsl_matrix_free(state->ATA);

  if (state->work_ATA)
    gsl_matrix_free(state->work_ATA);

  if (state->ATb)
    gsl_vector_free(state->ATb);

  if (state->D)
    gsl_vector_free(state->D);

  if (state->workp)
    gsl_vector_free(state->workp);

  if (state->work3p)
    gsl_vector_free(state->work3p);

  if (state->c)
    gsl_vector_free(state->c);

  if (state->eigen_p)
    gsl_eigen_symm_free(state->eigen_p);

  free(state);
}

static int
normal_reset(void *vstate)
{
  normal_state_t *state = (normal_state_t *) vstate;

  gsl_matrix_set_zero(state->ATA);
  gsl_vector_set_zero(state->ATb);
  state->normb = 0.0;
  state->eigen = 0;
  state->eval_min = 0.0;
  state->eval_max = 0.0;

  return GSL_SUCCESS;
}

/*
normal_accumulate()
  Add a new block of rows to the normal equations system

Inputs: A      - new block of rows, n-by-p
        b      - new rhs vector n-by-1
        vstate - workspace

Return: success/error
*/

static int
normal_accumulate(gsl_matrix * A, gsl_vector * b, void * vstate)
{
  normal_state_t *state = (normal_state_t *) vstate;
  const size_t n = A->size1;

  if (A->size2 != state->p)
    {
      GSL_ERROR("columns of A do not match workspace", GSL_EBADLEN);
    }
  else if (n != b->size)
    {
      GSL_ERROR("A and b have different numbers of rows", GSL_EBADLEN);
    }
  else
    {
      int s;

      /* ATA += A^T A, using only the lower half of the matrix */
      s = gsl_blas_dsyrk(CblasLower, CblasTrans, 1.0, A, 1.0, state->ATA);
      if (s)
        return s;

      /* ATb += A^T b */
      s = gsl_blas_dgemv(CblasTrans, 1.0, A, b, 1.0, state->ATb);
      if (s)
        return s;

      /* update || b || */
      state->normb = gsl_hypot(state->normb, gsl_blas_dnrm2(b));

      return GSL_SUCCESS;
    }
}

/*
normal_solve()
  Solve normal equations system:

(A^T A + \lambda^2 I) x = A^T b

using Cholesky decomposition

Inputs: lambda - regularization parameter
        x      - (output) solution vector p-by-1
        rnorm  - (output) residual norm ||b - A x||
        snorm  - (output) solution norm ||x||
        vstate - workspace

Return: success/error
*/

static int
normal_solve(const double lambda, gsl_vector * x,
             double * rnorm, double * snorm,
             void * vstate)
{
  normal_state_t *state = (normal_state_t *) vstate;

  if (x->size != state->p)
    {
      GSL_ERROR("solution vector does not match workspace", GSL_EBADLEN);
    }
  else
    {
      int status;

      /* solve system (A^T A) x = A^T b */
      status = normal_solve_system(lambda, x, state);
      if (status)
        {
          GSL_ERROR("failed to solve normal equations", status);
        }

      /* compute residual norm ||y - X c|| and solution norm ||x|| */
      normal_calc_norms(x, rnorm, snorm, state);

      return GSL_SUCCESS;
    }
}

static int
normal_rcond(double * rcond, void * vstate)
{
  normal_state_t *state = (normal_state_t *) vstate;
  int status = GSL_SUCCESS;
  double rcond_ATA;

  status = gsl_linalg_cholesky_rcond(state->work_ATA, &rcond_ATA, state->work3p);
  if (status == GSL_SUCCESS)
    *rcond = sqrt(rcond_ATA);

  return status;
}

/*
normal_lcurve()
  Compute L-curve of least squares system

Inputs: reg_param - (output) vector of regularization parameters
        rho       - (output) vector of residual norms
        eta       - (output) vector of solution norms
        vstate    - workspace

Return: success/error
*/

static int
normal_lcurve(gsl_vector * reg_param, gsl_vector * rho,
              gsl_vector * eta, void * vstate)
{
  normal_state_t *state = (normal_state_t *) vstate;
  int status;
  double smin, smax; /* minimum/maximum singular values */
  size_t i;

  if (state->eigen == 0)
    {
      status = normal_eigen(state);
      if (status)
        return status;
    }

  if (state->eval_max < 0.0)
    {
      GSL_ERROR("matrix is not positive definite", GSL_EDOM);
    }

  /* compute singular values which are sqrts of eigenvalues */
  smax = sqrt(state->eval_max);
  if (state->eval_min > 0.0)
    smin = sqrt(state->eval_min);
  else
    smin = 0.0;

  /* compute vector of regularization parameters */
  gsl_multifit_linear_lreg(smin, smax, reg_param);

  /* solve normal equations for each regularization parameter */
  for (i = 0; i < reg_param->size; ++i)
    {
      double lambda = gsl_vector_get(reg_param, i);
      double rnorm, snorm;

      status = normal_solve_system(lambda, state->c, state);
      if (status)
        return status;

      /* compute ||y - X c|| and ||c|| */
      normal_calc_norms(state->c, &rnorm, &snorm, state);

      gsl_vector_set(rho, i, rnorm);
      gsl_vector_set(eta, i, snorm);
    }

  return GSL_SUCCESS;
}

static const gsl_matrix *
normal_ATA(const void * vstate)
{
  const normal_state_t *state = (const normal_state_t *) vstate;
  return state->ATA;
}

static const gsl_vector *
normal_ATb(const void * vstate)
{
  const normal_state_t *state = (const normal_state_t *) vstate;
  return state->ATb;
}

/*
normal_solve_system()
  Compute solution to normal equations:

(A^T A + lambda^2*I) x = A^T b

using LDL decomposition.

Inputs: x     - (output) solution vector
        state - workspace

Return: success/error
*/

static int
normal_solve_system(const double lambda, gsl_vector * x, normal_state_t *state)
{
  int status;
  const double lambda_sq = lambda * lambda;
  gsl_vector_view d = gsl_matrix_diagonal(state->work_ATA);

  /* copy ATA matrix to temporary workspace and regularize */
  gsl_matrix_tricpy(CblasLower, CblasNonUnit, state->work_ATA, state->ATA);
  gsl_vector_add_constant(&d.vector, lambda_sq);

  /* solve with Cholesky decomposition */
  status = normal_solve_cholesky(state->work_ATA, state->ATb, x, state);
  if (status)
    return status;

  return status;
}

static int
normal_solve_cholesky(gsl_matrix * ATA, const gsl_vector * ATb,
                      gsl_vector * x, normal_state_t *state)
{
  int status;

  status = gsl_linalg_cholesky_decomp2(ATA, state->D);
  if (status)
    return status;

  status = gsl_linalg_cholesky_solve2(ATA, state->D, ATb, x);
  if (status)
    return status;

  return GSL_SUCCESS;
}

/*
normal_calc_norms()
  Compute residual norm ||y - X c|| and solution
norm ||c||

Inputs: x     - solution vector
        rnorm - (output) residual norm ||y - X c||
        snorm - (output) solution norm ||c||
        state - workspace
*/

static int
normal_calc_norms(const gsl_vector *x, double *rnorm,
                  double *snorm, normal_state_t *state)
{
  double r2;

  /* compute solution norm ||x|| */
  *snorm = gsl_blas_dnrm2(x);

  /* compute residual norm ||b - Ax|| */

  /* compute: A^T A x - 2 A^T b */
  gsl_vector_memcpy(state->workp, state->ATb);
  gsl_blas_dsymv(CblasLower, 1.0, state->ATA, x, -2.0, state->workp);

  /* compute: x^T A^T A x - 2 x^T A^T b */
  gsl_blas_ddot(x, state->workp, &r2);

  /* add b^T b */
  r2 += state->normb * state->normb;

  *rnorm = sqrt(r2);

  return GSL_SUCCESS;
}

/*
normal_eigen()
  Compute eigenvalues of A^T A matrix, which
are stored in state->workp on output. Also,
state->eval_min and state->eval_max are set
to the minimum/maximum eigenvalues
*/

static int
normal_eigen(normal_state_t *state)
{
  int status;

  /* copy lower triangle of ATA to temporary workspace */
  gsl_matrix_tricpy(CblasLower, CblasNonUnit, state->work_ATA, state->ATA);

  /* compute eigenvalues of ATA */
  status = gsl_eigen_symm(state->work_ATA, state->workp, state->eigen_p);
  if (status)
    return status;

  gsl_vector_minmax(state->workp, &state->eval_min, &state->eval_max);

  state->eigen = 1;

  return GSL_SUCCESS;
}

static const gsl_multilarge_linear_type normal_type =
{
  "normal",
  normal_alloc,
  normal_reset,
  normal_accumulate,
  normal_solve,
  normal_rcond,
  normal_lcurve,
  normal_ATA,
  normal_ATb,
  normal_free
};

const gsl_multilarge_linear_type * gsl_multilarge_linear_normal =
  &normal_type;

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim


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

#ifndef __GSL_BST_RB_H__
#define __GSL_BST_RB_H__

#include <gsl/gsl_math.h>

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

__BEGIN_DECLS

#ifndef GSL_BST_RB_MAX_HEIGHT
#define GSL_BST_RB_MAX_HEIGHT 48
#endif

/* red-black node */
struct gsl_bst_rb_node
{
  struct gsl_bst_rb_node *rb_link[2]; /* subtrees */
  void *rb_data;                      /* pointer to data */
  unsigned char rb_color;             /* color */
};

/* red-black tree data structure */
typedef struct
{
  struct gsl_bst_rb_node *rb_root;   /* tree's root */
  gsl_bst_cmp_function *rb_compare;  /* comparison function */
  void *rb_param;                    /* extra argument to |rb_compare| */
  const gsl_bst_allocator *rb_alloc; /* memory allocator */
  size_t rb_count;                   /* number of items in tree */
  unsigned long rb_generation;       /* generation number */
} gsl_bst_rb_table;

/* red-black traverser structure */
typedef struct
{
  const gsl_bst_rb_table *rb_table;  /* tree being traversed */
  struct gsl_bst_rb_node *rb_node;   /* current node in tree */
  struct gsl_bst_rb_node *rb_stack[GSL_BST_RB_MAX_HEIGHT];
                                     /* all the nodes above |rb_node| */
  size_t rb_height;                  /* number of nodes in |rb_parent| */
  unsigned long rb_generation;       /* generation number */
} gsl_bst_rb_traverser;

__END_DECLS

#endif /* __GSL_BST_RB_H__ */

max_stars_repo_name: vinej/sml

max_issues_repo_name: vinej/sml

max_forks_repo_name: vinej/sml


text:
/* randrand - how random is rand(), mostly for comparison with the
 * random numbers signal/jitterbits.c produces. produces 256 different
 * random values which in an ideal world should be evenly distributed
 *
 *   ./randrand | r-fu equichisq
 *
 * even with gsl_rng_uniform_int() and the "maximally equidistributed"
 * Tausworthe generator (mark II), the p-values over 100 runs of this
 * version of the program on Mac OS X vary wildly, with a mean of 0.55
 * and standard deviation of 0.28 (and similar stats from OpenBSD) */

#include <err.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sysexits.h>
#include <unistd.h>

#include <gsl/gsl_rng.h>

#define TIMES 10000000

int Return_Value = EXIT_SUCCESS;

int main(void)
{
    unsigned long i;
    gsl_rng *gslrng;

#ifdef __OpenBSD__
    if (pledge("stdio", NULL) == -1)
        err(1, "pledge failed");
#endif

    if (!(gslrng = gsl_rng_alloc(gsl_rng_taus2)))
        err(EX_SOFTWARE, "call to gsl_rng_alloc() failed");
    gsl_rng_set(gslrng, arc4random());

    /* for instead rand() */
    //srandomdev();
    for (i = 0; i < TIMES; i++) {
        /* for instead rand() */
        //printf("%u\n", rand() % 256);

        /* for instead arc4random() */
        //printf("%u\n", arc4random() % 256);

        printf("%lu\n", gsl_rng_uniform_int(gslrng, 256));

        /* reseeding influences things? (no?) */
        //if (i % 100 == 0)
        //  gsl_rng_set(gslrng, arc4random());
    }

    exit(Return_Value);
}

max_stars_repo_name: thrig/scripts

max_issues_repo_name: thrig/scripts

max_forks_repo_name: thrig/scripts


text:
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>
#include <galpy_potentials.h>
// Reflex Motion potential - based on Moving Object Potential
// 3 arguments: amp, t0, tf
void constrain_range2(double * d) {
  // Constrains index to be within interpolation range
  if (*d < 0) *d = 0.0;
  if (*d > 1) *d = 1.0;
}
double ReflexMotionRforce(double R,double z, double phi,
				   double t,
				   struct potentialArg * potentialArgs){
  double amp,t0,tf,d_ind,x,y,obj_x,obj_y,obj_z,obj_R,obj_phi,RF;
  double * args= potentialArgs->args;
  //Get args
  amp= *args;
  t0= *(args+1);
  tf= *(args+2);
  d_ind= (t-t0)/(tf-t0);
  x= R*cos(phi);
  y= R*sin(phi);
  constrain_range2(&d_ind);
  // Interpolate x, y, z
  obj_x= gsl_spline_eval(*potentialArgs->spline1d,d_ind,*potentialArgs->acc1d);
  obj_y= gsl_spline_eval(*(potentialArgs->spline1d+1),d_ind,
			 *(potentialArgs->acc1d+1));
  obj_z= gsl_spline_eval(*(potentialArgs->spline1d+2),d_ind,
			 *(potentialArgs->acc1d+2));
  // Rdist= pow(pow(x-obj_x, 2)+pow(y-obj_y, 2), 0.5);
  obj_R = pow(pow(obj_x,2) + pow(obj_y,2),0.5);
  obj_phi = atan2(obj_y,obj_x);
  // Calculate R force
  RF= calcRforce(obj_R,obj_z,obj_phi,t,potentialArgs->nwrapped,
		 potentialArgs->wrappedPotentialArg);
  return -amp*RF*(cos(phi)*cos(obj_phi)+sin(phi)*sin(obj_phi));
}

double ReflexMotionzforce(double R,double z,double phi,
				      double t,
				      struct potentialArg * potentialArgs){
  double amp,t0,tf,d_ind,x,y,obj_x,obj_y,obj_z, obj_R, obj_phi;
  double * args= potentialArgs->args;
  //Get args
  amp= *args;
  t0= *(args+1);
  tf= *(args+2);
  d_ind= (t-t0)/(tf-t0);
  x= R*cos(phi);
  y= R*sin(phi);
  constrain_range2(&d_ind);
  // Interpolate x, y, z
  obj_x= gsl_spline_eval(*potentialArgs->spline1d,d_ind,*potentialArgs->acc1d);
  obj_y= gsl_spline_eval(*(potentialArgs->spline1d+1),d_ind,
			 *(potentialArgs->acc1d+1));
  obj_z= gsl_spline_eval(*(potentialArgs->spline1d+2),d_ind,
			 *(potentialArgs->acc1d+2));
  // Rdist= pow(pow(x-obj_x, 2)+pow(y-obj_y, 2), 0.5);
  obj_R = pow(pow(obj_x,2) + pow(obj_y,2),0.5);
  obj_phi = atan2(obj_y,obj_x);

  // Calculate z force
  return -amp * calczforce(obj_R,obj_z,obj_phi,t,potentialArgs->nwrapped,
			   potentialArgs->wrappedPotentialArg);
}

double ReflexMotionphiforce(double R,double z,double phi,
					double t,
					struct potentialArg * potentialArgs){
  double amp,t0,tf,d_ind,x,y,obj_x,obj_y,obj_z, obj_R, obj_phi, RF;
  double * args= potentialArgs->args;
  //Get args
  amp= *args;
  t0= *(args+1);
  tf= *(args+2);
  d_ind= (t-t0)/(tf-t0);
  x= R*cos(phi);
  y= R*sin(phi);
  constrain_range2(&d_ind);
  // Interpolate x, y, z
  obj_x= gsl_spline_eval(*potentialArgs->spline1d,d_ind,*potentialArgs->acc1d);
  obj_y= gsl_spline_eval(*(potentialArgs->spline1d+1),d_ind,
			 *(potentialArgs->acc1d+1));
  obj_z= gsl_spline_eval(*(potentialArgs->spline1d+2),d_ind,
			 *(potentialArgs->acc1d+2));
  // Rdist= pow(pow(x-obj_x, 2)+pow(y-obj_y, 2), 0.5);
  obj_R = pow(pow(obj_x,2) + pow(obj_y,2),0.5);
  obj_phi = atan2(obj_y,obj_x);

  // Calculate phiforce
  RF= calcRforce(obj_R,obj_z,obj_phi,t,potentialArgs->nwrapped,
		 potentialArgs->wrappedPotentialArg);
  return amp*RF*(sin(phi)*cos(obj_phi)-cos(phi)*sin(obj_phi));
}

// double MovingObjectPotentialPlanarRforce(double R, double phi,
// 				      double t,
// 				      struct potentialArg * potentialArgs){
//  double amp,t0,tf,d_ind,x,y,obj_x,obj_y,Rdist,RF;
//  double * args= potentialArgs->args;
  //Get args
//  amp= *args;
//  t0= *(args+1);
//  tf= *(args+2);
//  d_ind= (t-t0)/(tf-t0);
//  x= R*cos(phi);
//  y= R*sin(phi);
//  constrain_range(&d_ind);
  // Interpolate x, y
//  obj_x= gsl_spline_eval(*potentialArgs->spline1d,d_ind,*potentialArgs->acc1d);
//  obj_y= gsl_spline_eval(*(potentialArgs->spline1d+1),d_ind,
//			 *(potentialArgs->acc1d+1));
//  Rdist= pow(pow(x-obj_x, 2)+pow(y-obj_y, 2), 0.5);
  // Calculate R force
//  RF= calcPlanarRforce(Rdist, phi, t, potentialArgs->nwrapped,
//		       potentialArgs->wrappedPotentialArg);
//  return -amp*RF*(cos(phi)*(obj_x-x)+sin(phi)*(obj_y-y))/Rdist;
//}

// double MovingObjectPotentialPlanarphiforce(double R, double phi,
//					double t,
//					struct potentialArg * potentialArgs){
//  double amp,t0,tf,d_ind,x,y,obj_x,obj_y,Rdist,RF;
//  double * args= potentialArgs->args;
  // Get args
//  amp= *args;
//  t0= *(args+1);
//  tf= *(args+2);
//  d_ind= (t-t0)/(tf-t0);
//  x= R*cos(phi);
//  y= R*sin(phi);
//  constrain_range(&d_ind);
  // Interpolate x, y
//  obj_x= gsl_spline_eval(*potentialArgs->spline1d,d_ind,*potentialArgs->acc1d);
//  obj_y= gsl_spline_eval(*(potentialArgs->spline1d+1),d_ind,
//			 *(potentialArgs->acc1d+1));
//  Rdist= pow(pow(x-obj_x, 2)+pow(y-obj_y, 2), 0.5);
  // Calculate phiforce
//  RF= calcPlanarRforce(Rdist, phi, t, potentialArgs->nwrapped,
//		       potentialArgs->wrappedPotentialArg);
//  return -amp*RF*R*(cos(phi)*(obj_y-y)-sin(phi)*(obj_x-x))/Rdist;
//}

max_stars_repo_name: radsouza/galpy

max_issues_repo_name: radsouza/galpy

max_forks_repo_name: radsouza/galpy


text:
#include <math.h>
#include <time.h>
#include <string.h>
#include <stdarg.h>
#include <stdio.h>
#include <gsl/gsl_vector.h>
// #include "libiomp/omp.h"
#include <omp.h>
#include "../include/type.h"
#include "../include/util.h"
#include "../include/global.h"

extern clock_t dur;
extern char *optarg;
extern int optind;

char c_getopt(int argc, char *argv[], char *pattern) {
    static size_t i = 0;
    int j;
    char ret;
    ++optind;

    if (optind > (int)strlen(pattern)) {
        return -1;
    }

    ret = pattern[i++];

    for (j = 0; j < argc; ++j) {
        if (argv[j][0] == '-' && argv[j][1] == ret) {
            optarg = argv[j+1];
        }
    }

    return ret;
}

/**
 * @brief initiator of odiff.
 *
 * @param inc1
 * @param inc2
 * @param skp1
 * @param skp2
 * @param inclu_len
 * @param skip_len
 * @param flag
 * @param id
 *
 * @return
 */
odiff* diff_alloc(gsl_vector* inc1, gsl_vector* inc2,
                  gsl_vector* skp1, gsl_vector* skp2,
                  int inclu_len, int skip_len, int flag, char* id) {
    odiff *diff = (odiff*)malloc(sizeof(odiff));
    diff->inc1 = inc1, diff->skp1 = skp1, diff->inc2 = inc2, diff->skp2 = skp2;
    diff->inclu_len = inclu_len, diff->skip_len = skip_len;
    diff->flag = flag, diff->id = (char*)malloc(sizeof(char)*(strlen(id)+1));
    strcpy(diff->id, id);
    return diff;
}

/**
 * @brief performing unary operation on gsl_vector.
 * @param vec
 * @param fun
 */
void unary_operate_vec(gsl_vector *vec, double (*fun) (double)) {
    size_t i;
    for (i = 0; i < vec->size; ++i) {
        gsl_vector_set(vec, i, fun(gsl_vector_get(vec, i)));
    }
    return;
}

/**
 * @brief performing binary operation on gsl_vector.
 *
 * @param vec
 * @param pa
 * @param fun
 */
void binary_operate_vec(gsl_vector *vec, double pa, double (*fun) (double, const double)) {
    size_t i;
    for (i = 0; i < vec->size; ++i) {
        gsl_vector_set(vec, i, fun(gsl_vector_get(vec, i), pa));
    }
    return;
}

/**
 * @brief performing element-wise operation on gsl_vector.
 *
 * @param vec
 * @param fun
 * @param argc
 * @param ...
 */
void element_wise_vec(gsl_vector *vec, double (*fun) (double, va_list), int argc, ...) {
    va_list argv;
    size_t i;
    for (i = 0; i < vec->size; ++i) {
        va_start(argv, argc);
        gsl_vector_set(vec, i, fun(gsl_vector_get(vec, i), argv));
        va_end(argv);
    }
    return;
}

/**
 * @brief cumulative summation of a gsl_vector.
 *
 * @param vec
 * @param fun
 * @param argc
 * @param ...
 *
 * @return
 */
double cuscumsum(gsl_vector *vec, double (*fun) (double, va_list), int argc, ...) {
    va_list argv, tmp;
    double sum = 0;
    size_t i;
    va_start(argv, argc);
    for (i = 0; i < vec->size; ++i) {
        va_copy(tmp, argv);
        sum += fun(gsl_vector_get(vec, i), tmp);
    }
    va_end(tmp);
    va_end(argv);
    return sum;
}

int parse_title(char* str, char** output) {
    int idx = 0;
    char *token = strtok(str, " \t\r\n\v\f");
    while(token) {
        output[idx] = (char*)malloc(sizeof(char)*TITLE_ELEMENT_LEN);
        strcpy(output[idx++], token);
        token = strtok(NULL, " \t\r\n\v\f");
    }

    // int base = 0;
    // do {
    //     base += idx? strlen(output[idx-1]) + LEN_OF_NT : 0;
    //     output[idx] = (char*)malloc(sizeof(char)*TITLE_ELEMENT_LEN);
    // } while(sscanf(str + base, "%s \t\n", output[idx++]) != EOF);

    return idx-1;
}

/**
 * @brief
 *
 * @param input a string in format '1,2,3,4,5,6,...'.
 * @param vec
 *
 * @return
 */
int str_to_vector(char* input, gsl_vector** vec) {
    int idx = 0, count = 1;
    size_t i = 0;
    // gsl_vector_free(*vec);

    for (i = 0; i < strlen(input); ++i) {
        if (input[i] == ',') {
            ++count;
        }
    }

    *vec = gsl_vector_alloc(count);
    char *str = strtok(input, ",");
    while(str) {
        if (strcmp(str, "NA") == 0 || strcmp(str, "NAN") == 0 ||
            strcmp(str, "na") == 0 || strcmp(str, "nan") == 0) {
            return -1;
        }
        gsl_vector_set(*vec, idx++, atof(str));
        str = strtok(NULL, ",");
    }

    return count;
}

int parse_line(char* str, char* id, gsl_vector** inc1, gsl_vector** skp1,
               gsl_vector** inc2, gsl_vector** skp2,
               int* inclu_len, int* skip_len) {
    int idx = 0;
    char output[7][COLUMN_LEN];

    char *token = strtok(str, " \t\r\n\v\f");
    while(token) {
        strcpy(output[idx++], token);
        token = strtok(NULL, " \t\r\n\v\f");
    }

    // int base = 0;
    // do {
    //     base += idx? strlen(output[idx-1]) + LEN_OF_NT : 0;
    // } while(sscanf(str + base, "%s \t\n", output[idx++]) != EOF);

    strcpy(id, output[0]);
    if (str_to_vector(output[1], inc1) == -1 ||
        str_to_vector(output[2], skp1) == -1 ||
        str_to_vector(output[3], inc2) == -1 ||
        str_to_vector(output[4], skp2) == -1) {
        printf("An error occured. rMATS cannot handle missing value. Sample Id: %s\n", id);
        printf("Exiting.\n");
        exit(0);

    }
    *inclu_len = atoi(output[5]);
    *skip_len = atoi(output[6]);

    return idx-1;
}

int parse_file(const char* filename, diff_list_node* list, char** title_element_list) {
    FILE *ifp;
    size_t olen = MAX_LINE;
    char *str_line = (char*)malloc(sizeof(char)*olen);
    char title[TITLE_LEN], id[MAX_LINE];
    int row_num=0, inclu_len, skip_len;
    gsl_vector *inc1 = NULL, *skp1 = NULL, *inc2 = NULL, *skp2 = NULL;

    if ((ifp = fopen(filename, "r")) == NULL) {
        printf("Fail to open!");
        return 0;
    }
    fgets(title, TITLE_LEN, ifp);
    parse_title(title, title_element_list);

    while (getline(&str_line, &olen, ifp) != -1) {
        ++row_num;
        parse_line(str_line, id, &inc1, &skp1, &inc2, &skp2, &inclu_len, &skip_len);
        if (inc1->size != skp1->size || inc2->size != skp2->size) {
            printf("An error occured. The length of the pair of vector should be equal. Sample Id: %s\n", id);
            printf("Size of vector: %ld, %ld, %ld, %ld\n", inc1->size, skp1->size, inc2->size, skp2->size);
            printf("Exiting\n");
            exit(0);
        }
        gsl_vector_add(inc1, skp1), gsl_vector_add(inc2, skp2);

        // TODO | or || ?
        // According to original python code 'if (vecprod(vecadd(inc1,skp1))==0) | (vecprod(vecadd(inc2,skp2))==0):',
        // it's a bit arithmetic '|'. However,  it should be a logical 'or' in such senario.
        if (cumprod(inc1) == 0 || cumprod(inc2) == 0) {
            gsl_vector_sub(inc1, skp1), gsl_vector_sub(inc2, skp2);
            diff_append(list, diff_alloc(inc1, inc2, skp1, skp2, inclu_len, skip_len, 0, id));
        } else {
            gsl_vector_sub(inc1, skp1), gsl_vector_sub(inc2, skp2);

            // TODO According to original python code, this function will not change anything.
            // original comment: add 1 in both inclusion and skipping counts for robustness in small sample size
            gsl_vector_add_constant(inc1, 0.0), gsl_vector_add_constant(skp1, 0.0);
            gsl_vector_add_constant(inc2, 0.0), gsl_vector_add_constant(skp2, 0.0);
            diff_append(list, diff_alloc(inc1, inc2, skp1, skp2, inclu_len, skip_len, 1, id));
        }
    }
    free(str_line);
    fclose(ifp);

    return row_num;
}

void display_dvector(const gsl_vector* vec) {
    size_t i;
    for (i = 0; i < vec->size; ++i) {
        printf("%.10f ", gsl_vector_get(vec, i));
    }
}

/**
 * @brief handy function.
 *
 * @param i
 * @param argv
 *
 * @return
 */
double logit(double i) {
    if (i < 0.01) {
        i = 0.01;
    } else if (i > 0.99) {
        i = 0.99;
    }
    return log(i/(1-i));
}

/*
double log_and_minus_x(const double i, va_list argv) {
    double pa = va_arg(argv, double);
    return logit(i) - logit(pa);
}

double rev_log_and_minus_x(const double i, va_list argv) {
    double pa = va_arg(argv, double);
    return pow(M_E, i + log(pa));
}
*/

/**
 * @brief cumulative production of a gsl_vector.
 *
 * @param vec
 *
 * @return
 */
double cumprod(const gsl_vector* vec) {
    double res = 1;
    size_t i;
    for (i = 0; i < vec->size; ++i) {
        res *= gsl_vector_get(vec, i);
    }
    return res;
}

// for multivar, 1, 2
double sum_for_multivar(const double i, va_list argv) {
    double pa = va_arg(argv, double);
    return pow(logit(i)-logit(pa), 2);
}

// for multivar_der, 1_der, 2_der
double sum_for_multivar_der(const double i, va_list argv) {
    double pa = va_arg(argv, double);
    // return -2 * (logit(i) - logit(pa))/(pa*(1-pa));
    return 2 * (logit(i) - logit(pa))/(pa*pa-pa);
}

double myfunc_marginal_2der(const double x, const double I, const double S,
                            const double beta, const double var,
                            const int inclu_len, const int skip_len) {
    double tmp1, tmp2, tmp3, one_x = 1-x;
    double powx = pow(x,2), pow1_x = pow(one_x,2), pow_len = pow(inclu_len*x + skip_len*one_x,2);

    tmp1 = (((2 * x - 1) * (logit(beta) - logit(x)) - 1)/var - 1)/(powx*pow1_x);
    tmp2 = I * skip_len * (2*inclu_len*x+skip_len)/(powx * pow_len);
    tmp3 = S * inclu_len * (inclu_len+2*skip_len*one_x)/(pow1_x * pow_len);

    return tmp1 - tmp2 - tmp3;
}

// for marginal
double sum_for_marginal(const double i, va_list argv) {
    int *idx = va_arg(argv, int*);
    double beta = va_arg(argv, double);
    double I_ = gsl_vector_get(va_arg(argv, gsl_vector*), *idx);
    double S_ = gsl_vector_get(va_arg(argv, gsl_vector*), *idx);
    double var = va_arg(argv, double), new_psi, f1, f1_2der;
    int inclu_len = va_arg(argv, int), skip_len = va_arg(argv, int);
    *idx += 1;

    new_psi = inclu_len * i/(inclu_len * i + skip_len * (1 - i));
    f1 = I_ * log(new_psi) + S_ * log(1 - new_psi) -
         pow(logit(i) - logit(beta),2)/(2*var) - log(i) - log(1-i);
    f1_2der = fabs(myfunc_marginal_2der(i, I_, S_, beta,
                                        var, inclu_len, skip_len));

    return 0.5 * log(fabs(f1_2der) + 0.00001) - f1;
}

// for marginal_der
double sum_for_marginal_der(const double i, va_list argv) {
    int *idx = va_arg(argv, int*);
    double beta = va_arg(argv, double);
    double I_ = gsl_vector_get(va_arg(argv, gsl_vector*), *idx);
    double S_ = gsl_vector_get(va_arg(argv, gsl_vector*), *idx);
    double var = va_arg(argv, double);
    int inclu_len = va_arg(argv, int), skip_len = va_arg(argv, int);
    double f1_1der, f1_2der, f1_3der, one_i = 1-i;
    *idx += 1;
    double powi = pow(i,2), pow1_i = pow(one_i,2);

    f1_3der = (2 * i - 1)/(powi * pow1_i*beta*(1-beta)*var);
    f1_2der = myfunc_marginal_2der(i, I_, S_, beta, var, inclu_len, skip_len);
    f1_1der = (logit(i) - logit(beta))/(beta * (1-beta) * var);

    return 0.5 * f1_3der/f1_2der - f1_1der;
}

// function used to manipulate our linked list.
int diff_append(diff_list_node* header, odiff* data) {
    diff_list_node *tmp = (diff_list_node*)malloc(sizeof(diff_list_node));
    tmp->data = data;
    tmp->end = tmp;
    tmp->next = NULL;
    header->end->next = tmp;
    header->end->end = tmp;
    header->end = tmp;

    return 0;
}

int diff_insert(diff_list_node* header, odiff* data, int idx) {

    return 6;
}

int diff_get_next(diff_list_node* header, odiff* data) {

    return 6;
}

int diff_get_at(diff_list_node* header, odiff* data, int idx) {
    int i;
    diff_list_node* tmp = header;
    for (i = 0; i <= idx; ++i) {
        tmp = tmp->next;
    }
    *data = *(tmp->data);
    return 0;
}

void split_data_list(int lenofl, int batch_size, diff_list_node* list, batch_datum *ret) {
    int i = 0, group = lenofl/batch_size, carry = lenofl % batch_size;
    diff_list_node *node = list;
    odiff **datum = (odiff**)malloc(sizeof(odiff*)*(group+1));
    for (i = 0; i < group+1; ++i) {
        datum[i] = (odiff*)malloc(sizeof(odiff)*batch_size);
    }

    for (i = 0; i < lenofl; ++i) {
        node = node->next;
        datum[i/batch_size][i%batch_size] = *(node->data);
    }
    for (i = 0; i < group; ++i) {
        ret[i].batch_size = batch_size;
        ret[i].datum = (void**)&datum[i];
    }
    ret[group].batch_size = carry;
    ret[group].datum = (void**)&datum[group];

    return;
}

void mp_threadpool(int nthread, int ntask, void* (*func)(void *), void** datum, void **ret) {
    int i;
    omp_set_num_threads(nthread);

    #pragma omp parallel for private(i)
    for (i = 0; i < ntask; ++i) {
        ret[i] = (*func)(datum[i]);
    }

    return;
}

/**
 * @brief C wrapper of Fortran l_bfgs_b routine.
 *
 * @param n
 * @param m
 * @param x[]
 * @param l[]
 * @param u[]
 * @param nbd[]
 * @param fp
 * @param gp
 * @param factr
 * @param pgtol
 * @param iprint
 * @param maxiter
 * @param argc
 * @param ...
 *
 * @return
********************************************************************
c     --------------------------------------------------------------
c             DESCRIPTION OF THE VARIABLES IN L-BFGS-B
c     --------------------------------------------------------------
c
c     n is an INTEGER variable that must be set by the user to the
c       number of variables.  It is not altered by the routine.
c
c     m is an INTEGER variable that must be set by the user to the
c       number of corrections used in the limited memory matrix.
c       It is not altered by the routine.  Values of m < 3  are
c       not recommended, and large values of m can result in excessive
c       computing time. The range  3 <= m <= 20 is recommended.
c
c     x is a DOUBLE PRECISION array of length n.  On initial entry
c       it must be set by the user to the values of the initial
c       estimate of the solution vector.  Upon successful exit, it
c       contains the values of the variables at the best point
c       found (usually an approximate solution).
c
c     l is a DOUBLE PRECISION array of length n that must be set by
c       the user to the values of the lower bounds on the variables. If
c       the i-th variable has no lower bound, l(i) need not be defined.
c
c     u is a DOUBLE PRECISION array of length n that must be set by
c       the user to the values of the upper bounds on the variables. If
c       the i-th variable has no upper bound, u(i) need not be defined.
c
c     nbd is an INTEGER array of dimension n that must be set by the
c       user to the type of bounds imposed on the variables:
c       nbd(i)=0 if x(i) is unbounded,
c              1 if x(i) has only a lower bound,
c              2 if x(i) has both lower and upper bounds,
c              3 if x(i) has only an upper bound.
c
c     f is a DOUBLE PRECISION variable.  If the routine setulb returns
c       with task(1:2)= 'FG', then f must be set by the user to
c       contain the value of the function at the point x.
c
c     g is a DOUBLE PRECISION array of length n.  If the routine setulb
c       returns with taskb(1:2)= 'FG', then g must be set by the user to
c       contain the components of the gradient at the point x.
c
c     factr is a DOUBLE PRECISION variable that must be set by the user.
c       It is a tolerance in the termination test for the algorithm.
c       The iteration will stop when
c
c        (f^k - f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= factr*epsmch
c
c       where epsmch is the machine precision which is automatically
c       generated by the code. Typical values for factr on a computer
c       with 15 digits of accuracy in double precision are:
c       factr=1.d+12 for low accuracy;
c             1.d+7  for moderate accuracy;
c             1.d+1  for extremely high accuracy.
c       The user can suppress this termination test by setting factr=0.
c
c     pgtol is a double precision variable.
c       On entry pgtol >= 0 is specified by the user.  The iteration
c         will stop when
c
c                 max{|proj g_i | i = 1, ..., n} <= pgtol
c
c         where pg_i is the ith component of the projected gradient.
c       The user can suppress this termination test by setting pgtol=0.
c
c     wa is a DOUBLE PRECISION  array of length
c       (2mmax + 4)nmax + 12mmax^2 + 12mmax used as workspace.
c       This array must not be altered by the user.
c
c     iwa is an INTEGER  array of length 3nmax used as
c       workspace. This array must not be altered by the user.
c
c     task is a CHARACTER string of length 60.
c       On first entry, it must be set to 'START'.
c       On a return with task(1:2)='FG', the user must evaluate the
c         function f and gradient g at the returned value of x.
c       On a return with task(1:5)='NEW_X', an iteration of the
c         algorithm has concluded, and f and g contain f(x) and g(x)
c         respectively.  The user can decide whether to continue or stop
c         the iteration.
c       When
c         task(1:4)='CONV', the termination test in L-BFGS-B has been
c           satisfied;
c         task(1:4)='ABNO', the routine has terminated abnormally
c           without being able to satisfy the termination conditions,
c           x contains the best approximation found,
c           f and g contain f(x) and g(x) respectively;
c         task(1:5)='ERROR', the routine has detected an error in the
c           input parameters;
c       On exit with task = 'CONV', 'ABNO' or 'ERROR', the variable task
c         contains additional information that the user can print.
c       This array should not be altered unless the user wants to
c          stop the run for some reason.  See driver2 or driver3
c          for a detailed explanation on how to stop the run
c          by assigning task(1:4)='STOP' in the driver.
c
c     iprint is an INTEGER variable that must be set by the user.
c       It controls the frequency and type of output generated:
c        iprint<0    no output is generated;
c        iprint=0    print only one line at the last iteration;
c        0<iprint<99 print also f and |proj g| every iprint iterations;
c        iprint=99   print details of every iteration except n-vectors;
c        iprint=100  print also the changes of active set and final x;
c        iprint>100  print details of every iteration including x and g;
c       When iprint > 0, the file iterate.dat will be created to
c                        summarize the iteration.
c
c     csave  is a CHARACTER working array of length 60.
c
c     lsave is a LOGICAL working array of dimension 4.
c       On exit with task = 'NEW_X', the following information is
c         available:
c       lsave(1) = .true.  the initial x did not satisfy the bounds;
c       lsave(2) = .true.  the problem contains bounds;
c       lsave(3) = .true.  each variable has upper and lower bounds.
c
c     isave is an INTEGER working array of dimension 44.
c       On exit with task = 'NEW_X', it contains information that
c       the user may want to access:
c         isave(30) = the current iteration number;
c         isave(34) = the total number of function and gradient
c                         evaluations;
c         isave(36) = the number of function value or gradient
c                                  evaluations in the current iteration;
c         isave(38) = the number of free variables in the current
c                         iteration;
c         isave(39) = the number of active constraints at the current
c                         iteration;
c
c         see the subroutine setulb.f for a description of other
c         information contained in isave
c
c     dsave is a DOUBLE PRECISION working array of dimension 29.
c       On exit with task = 'NEW_X', it contains information that
c         the user may want to access:
c         dsave(2) = the value of f at the previous iteration;
c         dsave(5) = the machine precision epsmch generated by the code;
c         dsave(13) = the infinity norm of the projected gradient;
c
c         see the subroutine setulb.f for a description of other
c         information contained in dsave
c
c     --------------------------------------------------------------
c           END OF THE DESCRIPTION OF THE VARIABLES IN L-BFGS-B
c     --------------------------------------------------------------
c
c     << An example of subroutine 'timer' for AIX Version 3.2 >>
c
c     subroutine timer(ttime)
c     double precision ttime
c     integer itemp, integer mclock
c     itemp = mclock()
c     ttime = dble(itemp)*1.0d-2
c     return
c     end

*********************************************************************
 */
double l_bfgs_b_wrapper(integer n, integer m, doublereal x[], doublereal l[],
                        doublereal u[], integer nbd[],
                        double (*fp) (const double x[], va_list argv),
                        void (*gp) (const double x[], double res[], va_list argv),
                        doublereal factr, doublereal pgtol, integer iprint,
                        int maxfun, int maxiter, int argc, ...) {
    va_list argv, tmp;
    char task[SIXTY], csave[SIXTY];
    logical lsave[4];
    integer iwa[3*n], isave[44];
    doublereal dsave[29], wa[2*m*n+5*n+11*m*m+8*m];
    va_start(argv, argc);

    strcpy(task,"START");
    int i, count = 0, funcount = 0, maxls = 20;
    for(i=5;i<SIXTY;i++) {
        task[i]=' ';
    }

    doublereal f = 0, g[n];
    for (i = 0; i < n; ++i) {
        g[i] = 0;
    }

    // TODO # The minimization routine wants f and g at the current x.
    //      # Note that interruptions due to maxfun are postponed
    //      # until the completion of the current minimization iteration.
    //      # Overwrite f and g:
    do {
        setulb_(&n,&m,x,l,u,nbd,&f,g,&factr,&pgtol,wa,iwa,task,&iprint,
                csave,lsave,isave,dsave, &maxls);

        if (strncmp(task,"FG",2)==0) {

            va_copy(tmp, argv);
            f = fp(x, tmp);
            va_copy(tmp, argv);
            gp(x, g, tmp);
            ++funcount;

        } else if(strncmp(task,"NEW_X",5)==0){
            ++count;
        } else {
            break;
        }
    } while(count < maxiter && funcount < maxfun);

    f = fp(x, argv);
    va_end(argv);
    va_end(tmp);

    return f;
}

max_stars_repo_name: chunjie-sam-liu/rmats-turbo

max_issues_repo_name: chunjie-sam-liu/rmats-turbo

max_forks_repo_name: chunjie-sam-liu/rmats-turbo


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

#include <config.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_histogram.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_ieee_utils.h>

#include "urand.c"

void
test1d_resample (void)
{
  size_t i;
  int status = 0;

  gsl_histogram *h;

  gsl_ieee_env_setup ();

  h = gsl_histogram_calloc_uniform (10, 0.0, 1.0);

  gsl_histogram_increment (h, 0.1);
  gsl_histogram_increment (h, 0.2);
  gsl_histogram_increment (h, 0.2);
  gsl_histogram_increment (h, 0.3);

  {
    gsl_histogram_pdf *p = gsl_histogram_pdf_alloc (10);

    gsl_histogram *hh = gsl_histogram_calloc_uniform (100, 0.0, 1.0);

    gsl_histogram_pdf_init (p, h);

    for (i = 0; i < 100000; i++)
      {
        double u = urand();
        double x = gsl_histogram_pdf_sample (p, u);
        gsl_histogram_increment (hh, x);
      }

    for (i = 0; i < 100; i++)
      {
        double y = gsl_histogram_get (hh, i) / 2500;
        double x, xmax;
        size_t k;
        double ya;

        gsl_histogram_get_range (hh, i, &x, &xmax);

        gsl_histogram_find (h, x, &k);
        ya = gsl_histogram_get (h, k);

        if (ya == 0)
          {
            if (y != 0)
              {
                printf ("%d: %g vs %g\n", (int) i, y, ya);
                status = 1;
              }
          }
        else
          {
            double err = 1 / sqrt (gsl_histogram_get (hh, i));
            double sigma = fabs ((y - ya) / (ya * err));
            if (sigma > 3)
              {
                status = 1;
                printf ("%g vs %g err=%g sigma=%g\n", y, ya, err, sigma);
              }
          }
      }

    gsl_histogram_pdf_free (p) ;
    gsl_histogram_free (hh);

    gsl_test (status, "gsl_histogram_pdf_sample within statistical errors");
  }

  gsl_histogram_free (h);
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
#ifndef AMICI_SUNDIALS_MATRIX_WRAPPER_H
#define AMICI_SUNDIALS_MATRIX_WRAPPER_H

#include <sunmatrix/sunmatrix_band.h>   // SUNMatrix_Band
#include <sunmatrix/sunmatrix_dense.h>  // SUNMatrix_Dense
#include <sunmatrix/sunmatrix_sparse.h> // SUNMatrix_Sparse

#include <gsl/gsl-lite.hpp>

#include <vector>

#include "amici/vector.h"

namespace amici {

/**
 * @brief A RAII wrapper for SUNMatrix structs.
 *
 * This can create dense, sparse, or banded matrices using the respective
 * constructor.
 */
class SUNMatrixWrapper {
  public:
    SUNMatrixWrapper() = default;

    /**
     * @brief Create sparse matrix. See SUNSparseMatrix in sunmatrix_sparse.h
     * @param M Number of rows
     * @param N Number of columns
     * @param NNZ Number of nonzeros
     * @param sparsetype Sparse type
     */
    SUNMatrixWrapper(int M, int N, int NNZ, int sparsetype);

    /**
     * @brief Create dense matrix. See SUNDenseMatrix in sunmatrix_dense.h
     * @param M Number of rows
     * @param N Number of columns
     */
    SUNMatrixWrapper(int M, int N);

    /**
     * @brief Create banded matrix. See SUNBandMatrix in sunmatrix_band.h
     * @param M Number of rows and columns
     * @param ubw Upper bandwidth
     * @param lbw Lower bandwidth
     */
    SUNMatrixWrapper(int M, int ubw, int lbw);

    /**
     * @brief Create sparse matrix from dense or banded matrix. See
     * SUNSparseFromDenseMatrix and SUNSparseFromBandMatrix in
     * sunmatrix_sparse.h
     * @param A Wrapper for dense matrix
     * @param droptol tolerance for dropping entries
     * @param sparsetype Sparse type
     */
    SUNMatrixWrapper(const SUNMatrixWrapper &A, realtype droptol,
                     int sparsetype);

    /**
     * @brief Wrap existing SUNMatrix
     * @param mat
     */
    explicit SUNMatrixWrapper(SUNMatrix mat);

    ~SUNMatrixWrapper();

    /**
     * @brief Copy constructor
     * @param other
     */
    SUNMatrixWrapper(const SUNMatrixWrapper &other);

    /**
     * @brief Move constructor
     * @param other
     */
    SUNMatrixWrapper(SUNMatrixWrapper &&other);

    /**
     * @brief Copy assignment
     * @param other
     * @return
     */
    SUNMatrixWrapper &operator=(const SUNMatrixWrapper &other);

    /**
     * @brief Move assignment
     * @param other
     * @return
     */
    SUNMatrixWrapper &operator=(SUNMatrixWrapper &&other);

    /**
     * @brief Access raw data
     * @return raw data pointer
     */
    realtype *data() const;

    /**
     * @brief Get the wrapped SUNMatrix
     * @return SlsMat
     */
    SUNMatrix get() const;

    /**
     * @brief Get the number of rows
     * @return number
     */
    sunindextype rows() const;

    /**
     * @brief Get the number of columns
     * @return number
     */
    sunindextype columns() const;

    /**
     * @brief Get the number of non-zero elements (sparse matrices only)
     * @return number
     */
    sunindextype nonzeros() const;

    /**
     * @brief Get the index values of a sparse matrix
     * @return index array
     */
    sunindextype *indexvals() const;

    /**
     * @brief Get the index pointers of a sparse matrix
     * @return index array
     */
    sunindextype *indexptrs() const;

    /**
     * @brief Get the type of sparse matrix
     * @return index array
     */
    int sparsetype() const;

    /**
     * @brief reset data to zeroes
     */
    void reset();

    /**
     * @brief multiply with a scalar (in-place)
     * @param a scalar value to multiply matrix
     */
    void scale(realtype a);

    /**
     * @brief N_Vector interface for multiply
     * @param c output vector, may already contain values
     * @param b multiplication vector
     */
    void multiply(N_Vector c, const_N_Vector b) const;

    /**
     * @brief Perform matrix vector multiplication c += A*b
     * @param c output vector, may already contain values
     * @param b multiplication vector
     */
    void multiply(gsl::span<realtype> c, gsl::span<const realtype> b) const;

    /**
     * @brief Perform reordered matrix vector multiplication c += A[:,cols]*b
     * @param c output vector, may already contain values
     * @param b multiplication vector
     * @param cols int vector for column reordering
     * @param transpose bool transpose A before multiplication
     */
    void multiply(N_Vector c,
                  const N_Vector b,
                  gsl::span <const int> cols,
                  bool transpose) const;

    /**
     * @brief Perform reordered matrix vector multiplication c += A[:,cols]*b
     * @param c output vector, may already contain values
     * @param b multiplication vector
     * @param cols int vector for column reordering
     * @param transpose bool transpose A before multiplication
     */
    void multiply(gsl::span<realtype> c,
                  gsl::span<const realtype> b,
                  gsl::span <const int> cols,
                  bool transpose) const;

    /**
     * @brief Perform matrix matrix multiplication
              C[:, :] += A * B
              for sparse A, B, C
     * @param C output matrix, may already contain values
     * @param B multiplication matrix
     */
    void sparse_multiply(SUNMatrixWrapper *C,
                         SUNMatrixWrapper *B) const;

    /**
     * @brief Set to 0.0
     */
    void zero();

  private:
    void update_ptrs();

    /**
     * @brief CSC matrix to which all methods are applied
     */
    SUNMatrix matrix_ {nullptr};
    realtype *data_ptr_ {nullptr};
    sunindextype *indexptrs_ptr_ {nullptr};
    sunindextype *indexvals_ptr_ {nullptr};
};

} // namespace amici

namespace gsl {
/**
 * @brief Create span from SUNMatrix
 * @param nv
 * @return
 */
inline span<realtype> make_span(SUNMatrix m)
{
    switch (SUNMatGetID(m)) {
    case SUNMATRIX_DENSE:
        return span<realtype>(SM_DATA_D(m), SM_LDATA_D(m));
    case SUNMATRIX_SPARSE:
        return span<realtype>(SM_DATA_S(m), SM_NNZ_S(m));
    default:
        throw amici::AmiException("Unimplemented SUNMatrix type for make_span");
    }
}
} // namespace gsl

#endif // AMICI_SUNDIALS_MATRIX_WRAPPER_H

max_stars_repo_name: lcontento/AMICI

max_issues_repo_name: lcontento/AMICI

max_forks_repo_name: lcontento/AMICI


text:
#include <viaio/Vlib.h>
#include <viaio/VImage.h>

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

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include "gsl_utils.h"

#define ROUND(x) (int)(((x) >= 0) ? ((x) + 0.5) : ((x) - 0.5))

void dfs(VFloat** Dfs,gsl_matrix_float* pinvX,gsl_matrix_float* X,gsl_matrix_float* con,int numlags)
{
  int numcon = 1;
  int n = X->size2;
  /* these vals represent the matlab interval 2:n */
  int k1 = 1;
  int k2 = n-1;
  int i;

  /* X^T */
  gsl_matrix_float* transX = gsl_matrix_float_alloc(X->size2, X->size1);
  gsl_matrix_float_transpose_memcpy(transX,X);

  gsl_vector_float* CorX2 = gsl_vector_float_alloc(numcon);
  gsl_vector_float_set_zero(CorX2);

  /* cpinvX=contrast^T*pinvX */
  gsl_matrix_float* cpinvX = fmat_x_mat(con, pinvX, NULL);

  /* CovX0=cpinvX*cpinvX' */
  gsl_matrix_float* CovX0 = fmat_x_matT(cpinvX,cpinvX,NULL);

  if(numlags == 1) {
    /* CovX1=cpinvX(:,k1)*cpinvX(:,k1-1)' */
    gsl_matrix_float_view sub1 =
      gsl_matrix_float_submatrix(cpinvX,0,k1,cpinvX->size1,k2);
    gsl_matrix_float_view sub2 =
      gsl_matrix_float_submatrix(cpinvX,0,k1-1,cpinvX->size1,k2);
    gsl_matrix_float* CovX1 = fmat_x_matT(&sub1.matrix, &sub2.matrix, NULL);

    gsl_vector_float_view CovX1diag = gsl_matrix_float_diagonal(CovX1);
    gsl_vector_float_view CovX0diag = gsl_matrix_float_diagonal(CovX0);

    gsl_vector_float_div (&CovX1diag.vector, &CovX0diag.vector);
    gsl_vector_float_mul (&CovX1diag.vector, &CovX1diag.vector);

    gsl_vector_float_memcpy (CorX2, &CovX1diag.vector);
  }
  else {
    int lag;
    gsl_matrix_float_view sub1, sub2;
    gsl_matrix_float* CovX1 = gsl_matrix_float_alloc(cpinvX->size1,cpinvX->size1);
    for(lag=0;lag<numlags;lag++){
      sub1 = gsl_matrix_float_submatrix(cpinvX,0,0,cpinvX->size1,k2-lag);
      sub2 = gsl_matrix_float_submatrix(cpinvX,0,lag+1,cpinvX->size1,k2-lag);
      fmat_x_matT(&sub1.matrix, &sub2.matrix,CovX1);

      gsl_vector_float_view CovX1diag = gsl_matrix_float_diagonal(CovX1);
      gsl_vector_float_view CovX0diag = gsl_matrix_float_diagonal(CovX0);

      gsl_vector_float_div (&CovX1diag.vector, &CovX0diag.vector);
      gsl_vector_float_mul (&CovX1diag.vector, &CovX1diag.vector);

      gsl_vector_float_add(CorX2,&CovX1diag.vector);
    }
    gsl_matrix_float_free(CovX1);
  }

  float dfresid  = n-rank(transX);
  float dfmin = dfresid;

  for (i=0;i<numcon;i++) {
    (*Dfs)[i] = (VFloat)ROUND(dfresid / (1 + 2*gsl_vector_float_get(CorX2,i)));
    if ((*Dfs)[i]<dfmin) dfmin = (*Dfs)[i];
  }
  (*Dfs)[numcon] = dfresid;

  /*
    if(dfmin < 100)
    fprintf(stderr, "  warning, dfmin= %f\n",dfmin);
  */

  gsl_matrix_float_free(cpinvX);
  gsl_matrix_float_free(transX);
  gsl_matrix_float_free(CovX0);
  gsl_vector_float_free(CorX2);
}

max_stars_repo_name: zrajna/lipsia

max_issues_repo_name: zrajna/lipsia

max_forks_repo_name: zrajna/lipsia


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

#ifndef __GSL_MATRIX_SHORT_H__
#define __GSL_MATRIX_SHORT_H__

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_vector_short.h>

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

__BEGIN_DECLS

typedef struct
{
  size_t size1;
  size_t size2;
  size_t tda;
  short * data;
  gsl_block_short * block;
  int owner;
} gsl_matrix_short;

typedef struct
{
  gsl_matrix_short matrix;
} _gsl_matrix_short_view;

typedef _gsl_matrix_short_view gsl_matrix_short_view;

typedef struct
{
  gsl_matrix_short matrix;
} _gsl_matrix_short_const_view;

typedef const _gsl_matrix_short_const_view gsl_matrix_short_const_view;

/* Allocation */

GSL_EXPORT
gsl_matrix_short *
gsl_matrix_short_alloc (const size_t n1, const size_t n2);

GSL_EXPORT
gsl_matrix_short *
gsl_matrix_short_calloc (const size_t n1, const size_t n2);

GSL_EXPORT
gsl_matrix_short *
gsl_matrix_short_alloc_from_block (gsl_block_short * b,
                                   const size_t offset,
                                   const size_t n1,
                                   const size_t n2,
                                   const size_t d2);

GSL_EXPORT
gsl_matrix_short *
gsl_matrix_short_alloc_from_matrix (gsl_matrix_short * m,
                                    const size_t k1,
                                    const size_t k2,
                                    const size_t n1,
                                    const size_t n2);

GSL_EXPORT
gsl_vector_short *
gsl_vector_short_alloc_row_from_matrix (gsl_matrix_short * m,
                                        const size_t i);

GSL_EXPORT
gsl_vector_short *
gsl_vector_short_alloc_col_from_matrix (gsl_matrix_short * m,
                                        const size_t j);

GSL_EXPORT void gsl_matrix_short_free (gsl_matrix_short * m);

/* Views */

GSL_EXPORT
_gsl_matrix_short_view
gsl_matrix_short_submatrix (gsl_matrix_short * m,
                            const size_t i, const size_t j,
                            const size_t n1, const size_t n2);

GSL_EXPORT
_gsl_vector_short_view
gsl_matrix_short_row (gsl_matrix_short * m, const size_t i);

GSL_EXPORT
_gsl_vector_short_view
gsl_matrix_short_column (gsl_matrix_short * m, const size_t j);

GSL_EXPORT
_gsl_vector_short_view
gsl_matrix_short_diagonal (gsl_matrix_short * m);

GSL_EXPORT
_gsl_vector_short_view
gsl_matrix_short_subdiagonal (gsl_matrix_short * m, const size_t k);

GSL_EXPORT
_gsl_vector_short_view
gsl_matrix_short_superdiagonal (gsl_matrix_short * m, const size_t k);

GSL_EXPORT
_gsl_matrix_short_view
gsl_matrix_short_view_array (short * base,
                             const size_t n1,
                             const size_t n2);

GSL_EXPORT
_gsl_matrix_short_view
gsl_matrix_short_view_array_with_tda (short * base,
                                      const size_t n1,
                                      const size_t n2,
                                      const size_t tda);

GSL_EXPORT
_gsl_matrix_short_view
gsl_matrix_short_view_vector (gsl_vector_short * v,
                              const size_t n1,
                              const size_t n2);

GSL_EXPORT
_gsl_matrix_short_view
gsl_matrix_short_view_vector_with_tda (gsl_vector_short * v,
                                       const size_t n1,
                                       const size_t n2,
                                       const size_t tda);

GSL_EXPORT
_gsl_matrix_short_const_view
gsl_matrix_short_const_submatrix (const gsl_matrix_short * m,
                                  const size_t i, const size_t j,
                                  const size_t n1, const size_t n2);

GSL_EXPORT
_gsl_vector_short_const_view
gsl_matrix_short_const_row (const gsl_matrix_short * m,
                            const size_t i);

GSL_EXPORT
_gsl_vector_short_const_view
gsl_matrix_short_const_column (const gsl_matrix_short * m,
                               const size_t j);

GSL_EXPORT
_gsl_vector_short_const_view
gsl_matrix_short_const_diagonal (const gsl_matrix_short * m);

GSL_EXPORT
_gsl_vector_short_const_view
gsl_matrix_short_const_subdiagonal (const gsl_matrix_short * m,
                                    const size_t k);

GSL_EXPORT
_gsl_vector_short_const_view
gsl_matrix_short_const_superdiagonal (const gsl_matrix_short * m,
                                      const size_t k);

GSL_EXPORT
_gsl_matrix_short_const_view
gsl_matrix_short_const_view_array (const short * base,
                                   const size_t n1,
                                   const size_t n2);

GSL_EXPORT
_gsl_matrix_short_const_view
gsl_matrix_short_const_view_array_with_tda (const short * base,
                                            const size_t n1,
                                            const size_t n2,
                                            const size_t tda);

GSL_EXPORT
_gsl_matrix_short_const_view
gsl_matrix_short_const_view_vector (const gsl_vector_short * v,
                                    const size_t n1,
                                    const size_t n2);

GSL_EXPORT
_gsl_matrix_short_const_view
gsl_matrix_short_const_view_vector_with_tda (const gsl_vector_short * v,
                                             const size_t n1,
                                             const size_t n2,
                                             const size_t tda);

/* Operations */

GSL_EXPORT short   gsl_matrix_short_get(const gsl_matrix_short * m, const size_t i, const size_t j);
GSL_EXPORT void    gsl_matrix_short_set(gsl_matrix_short * m, const size_t i, const size_t j, const short x);

GSL_EXPORT short * gsl_matrix_short_ptr(gsl_matrix_short * m, const size_t i, const size_t j);
GSL_EXPORT const short * gsl_matrix_short_const_ptr(const gsl_matrix_short * m, const size_t i, const size_t j);

GSL_EXPORT void gsl_matrix_short_set_zero (gsl_matrix_short * m);
GSL_EXPORT void gsl_matrix_short_set_identity (gsl_matrix_short * m);
GSL_EXPORT void gsl_matrix_short_set_all (gsl_matrix_short * m, short x);

GSL_EXPORT int gsl_matrix_short_fread (FILE * stream, gsl_matrix_short * m) ;
GSL_EXPORT int gsl_matrix_short_fwrite (FILE * stream, const gsl_matrix_short * m) ;
GSL_EXPORT int gsl_matrix_short_fscanf (FILE * stream, gsl_matrix_short * m);
GSL_EXPORT int gsl_matrix_short_fprintf (FILE * stream, const gsl_matrix_short * m, const char * format);

GSL_EXPORT int gsl_matrix_short_memcpy(gsl_matrix_short * dest, const gsl_matrix_short * src);
GSL_EXPORT int gsl_matrix_short_swap(gsl_matrix_short * m1, gsl_matrix_short * m2);

GSL_EXPORT int gsl_matrix_short_swap_rows(gsl_matrix_short * m, const size_t i, const size_t j);
GSL_EXPORT int gsl_matrix_short_swap_columns(gsl_matrix_short * m, const size_t i, const size_t j);
GSL_EXPORT int gsl_matrix_short_swap_rowcol(gsl_matrix_short * m, const size_t i, const size_t j);
GSL_EXPORT int gsl_matrix_short_transpose (gsl_matrix_short * m);
GSL_EXPORT int gsl_matrix_short_transpose_memcpy (gsl_matrix_short * dest, const gsl_matrix_short * src);

GSL_EXPORT short gsl_matrix_short_max (const gsl_matrix_short * m);
GSL_EXPORT short gsl_matrix_short_min (const gsl_matrix_short * m);
GSL_EXPORT void gsl_matrix_short_minmax (const gsl_matrix_short * m, short * min_out, short * max_out);

GSL_EXPORT void gsl_matrix_short_max_index (const gsl_matrix_short * m, size_t * imax, size_t *jmax);
GSL_EXPORT void gsl_matrix_short_min_index (const gsl_matrix_short * m, size_t * imin, size_t *jmin);
GSL_EXPORT void gsl_matrix_short_minmax_index (const gsl_matrix_short * m, size_t * imin, size_t * jmin, size_t * imax, size_t * jmax);

GSL_EXPORT int gsl_matrix_short_isnull (const gsl_matrix_short * m);

GSL_EXPORT int gsl_matrix_short_add (gsl_matrix_short * a, const gsl_matrix_short * b);
GSL_EXPORT int gsl_matrix_short_sub (gsl_matrix_short * a, const gsl_matrix_short * b);
GSL_EXPORT int gsl_matrix_short_mul_elements (gsl_matrix_short * a, const gsl_matrix_short * b);
GSL_EXPORT int gsl_matrix_short_div_elements (gsl_matrix_short * a, const gsl_matrix_short * b);
GSL_EXPORT int gsl_matrix_short_scale (gsl_matrix_short * a, const double x);
GSL_EXPORT int gsl_matrix_short_add_constant (gsl_matrix_short * a, const double x);
GSL_EXPORT int gsl_matrix_short_add_diagonal (gsl_matrix_short * a, const double x);

/***********************************************************************/
/* The functions below are obsolete                                    */
/***********************************************************************/
GSL_EXPORT int gsl_matrix_short_get_row(gsl_vector_short * v, const gsl_matrix_short * m, const size_t i);
GSL_EXPORT int gsl_matrix_short_get_col(gsl_vector_short * v, const gsl_matrix_short * m, const size_t j);
GSL_EXPORT int gsl_matrix_short_set_row(gsl_matrix_short * m, const size_t i, const gsl_vector_short * v);
GSL_EXPORT int gsl_matrix_short_set_col(gsl_matrix_short * m, const size_t j, const gsl_vector_short * v);

/* inline functions if you are using GCC */

#ifdef HAVE_INLINE
extern inline
short
gsl_matrix_short_get(const gsl_matrix_short * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (i >= m->size1)
    {
      GSL_ERROR_VAL("first index out of range", GSL_EINVAL, 0) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_VAL("second index out of range", GSL_EINVAL, 0) ;
    }
#endif
  return m->data[i * m->tda + j] ;
}

extern inline
void
gsl_matrix_short_set(gsl_matrix_short * m, const size_t i, const size_t j, const short x)
{
#if GSL_RANGE_CHECK
  if (i >= m->size1)
    {
      GSL_ERROR_VOID("first index out of range", GSL_EINVAL) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_VOID("second index out of range", GSL_EINVAL) ;
    }
#endif
  m->data[i * m->tda + j] = x ;
}

extern inline
short *
gsl_matrix_short_ptr(gsl_matrix_short * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (i >= m->size1)
    {
      GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
    }
#endif
  return (short *) (m->data + (i * m->tda + j)) ;
}

extern inline
const short *
gsl_matrix_short_const_ptr(const gsl_matrix_short * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (i >= m->size1)
    {
      GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
    }
#endif
  return (const short *) (m->data + (i * m->tda + j)) ;
}

#endif

__END_DECLS

#endif /* __GSL_MATRIX_SHORT_H__ */

max_stars_repo_name: dynaryu/vaws

max_issues_repo_name: dynaryu/vaws

max_forks_repo_name: dynaryu/vaws