text:
/*
 *  PIXON
 *  A Pixon-based method for reconstructing velocity-delay map in reverberation mapping.
 *
 *  Yan-Rong Li, liyanrong@mail.ihep.ac.cn
 *
 */

/*!
 *  \file mathfun.c
 *  \brief mathematic functions.
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <float.h>
#include <lapacke.h>
#include <cblas.h>

#include "mathfun.h"

#define EPS (DBL_MIN)

/*!
 * This function calculates matrix multiply C(nxn) = A(nxn) * B(nxn).
 */
void multiply_mat(double * a, double *b, double *c, int n)
{
  cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, n, n, n, 1.0f
                             , a, n, b, n, 0.0f, c, n);
}

/*!
 * This function calculates matrix multiply C(nxn) = A^T(nxn) * B(nxn).
 */
void multiply_mat_transposeA(double * a, double *b, double *c, int n)
{
  cblas_dgemm(CblasRowMajor, CblasTrans, CblasNoTrans, n, n, n, 1.0f
                             , a, n, b, n, 0.0f, c, n);
}

/*!
 * This function calculates matrix multiply C(nxn) = A(nxn) * B^T(nxn).
 */
void multiply_mat_transposeB(double * a, double *b, double *c, int n)
{
  cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasTrans, n, n, n, 1.0f
                             , a, n, b, n, 0.0f, c, n);
}

/*!
 * This function calculates matrix multiply Y(n) = A^T(nxn) * X(n).
 */
void multiply_matvec_transposeA(double *a, double *x, int n, double *y)
{
  cblas_dgemv(CblasRowMajor, CblasTrans, n, n, 1.0f, a, n, x, 1, 0.0f, y, 1);
}

/*!
 * This function calculates matrix multiply Y(n) = A(nxn) * X(n).
 */
void multiply_matvec(double *a, double *x, int n, double *y)
{
  cblas_dgemv(CblasRowMajor, CblasNoTrans, n, n, 1.0f, a, n, x, 1, 0.0f, y, 1);
}

/*!
 * This function calculates  Y(m) = A(m, n) * X(n).
 */
void multiply_matvec_MN(double * a, int m, int n, double *x, double *y)
{
  cblas_dgemv(CblasRowMajor, CblasNoTrans, m, n, 1.0f, a, n, x, 1, 0.0f, y, 1);
}

/* C(m*n) = A(m*k) * B(k*n) */
void multiply_mat_MN(double * a, double *b, double *c, int m, int n, int k)
{
  cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, n, k, 1.0f
                             , a, k, b, n, 0.0f, c, n);
}

/* C(m*n) = A^T(m*k) * B(k*n) */
void multiply_mat_MN_transposeA(double * a, double *b, double *c, int m, int n, int k)
{
  cblas_dgemm(CblasRowMajor, CblasTrans, CblasNoTrans, m, n, k, 1.0f
                             , a, m, b, n, 0.0f, c, n);
}
/* C(m*n) = A(m*k) * B^T(k*n) */
void multiply_mat_MN_transposeB(double * a, double *b, double *c, int m, int n, int k)
{
  cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasTrans, m, n, k, 1.0f
                             , a, k, b, k, 0.0f, c, n);
}

/* A(mxm)^-1 * B(mxn), store the output in B
 * note that A will be changed on exit. */
int multiply_mat_MN_inverseA(double * a, double *b, int m, int n)
{
  int * ipiv, info;
  ipiv=malloc(m*sizeof(int));

  info=LAPACKE_dgetrf(LAPACK_ROW_MAJOR, m, m, a, m, ipiv);
  if(info!=0)
  {
    printf("multiply_mat_MN_inverseA 1.\n this usually caused by improper nc.\n increase the low limit of nc");
    exit(0);
    return info;
  }
  info = LAPACKE_dgetrs(LAPACK_ROW_MAJOR, 'N', m, n, a, m, ipiv, b, n);
  if(info!=0)
  {
    printf("multiply_mat_MN_inverseA 2\n this usually caused by improper nc.\n increase the low limit of nc");
    exit(0);
    return info;
  }

  free(ipiv);
  return info;
}

/*!
 * This functions calculate A^-1(nxn).
 * A is a generic matrix.
 */
void inverse_mat(double * a, int n, int *info)
{
  int * ipiv;
  ipiv=malloc(n*sizeof(int));

//  dgetrf_(&n, &n, a, &n, ipiv, info);
//  dgetri_(&n, a, &n, ipiv, work, &lwork, info);

  LAPACKE_dgetrf(LAPACK_ROW_MAJOR, n, n, a, n, ipiv);
  LAPACKE_dgetri(LAPACK_ROW_MAJOR, n, a, n, ipiv);

  free(ipiv);
  return;
}

/*!
 * This functions calculate A^-1(nxn).
 * A is a symmetric matrix.
 */
void inverse_symat(double * a, int n, int *info)
{
  int * ipiv, i, j;
  ipiv=malloc(n*sizeof(int));

  LAPACKE_dsytrf(LAPACK_ROW_MAJOR, 'U', n, a, n, ipiv);
  LAPACKE_dsytri(LAPACK_ROW_MAJOR, 'U', n, a, n, ipiv);

  /* fill up the lower triangle */
  for(i=0; i<n; i++)
    for(j=0; j<i; j++)
      a[i*n+j] = a[j*n+i];

  free(ipiv);
  return;
}

/*!
 * This functions calculate A^-1(nxn).
 * A is a postive-definite symmetrix matrix.
 */
void inverse_pomat(double * a, int n, int *info)
{
  int i, j;

  LAPACKE_dpotrf(LAPACK_ROW_MAJOR, 'U', n, a, n);
  LAPACKE_dpotri(LAPACK_ROW_MAJOR, 'U', n, a, n);

  /* fill up the lower triangle */
  for(i=0; i<n; i++)
    for(j=0; j<i; j++)
      a[i*n+j] = a[j*n+i];
  return;
}

/*!
 * This function calculate eigenvectors and eigenvalues of matrix.
 */
void eigen_sym_mat(double *a, int n, double *val, int *info)
{
    char jobz='V', uplo='U';

/* store the eigenvectors  in a by rows.
 * store the eigenvalues in val in ascending order.
 */
//    dsyev_(&jobz, &uplo, &n, a, &n, val, work, &lwork, info);

/* store the eigenvectors  in a by columns.
 * store the eigenvalues in val in ascending order.
 */
    LAPACKE_dsyev(LAPACK_ROW_MAJOR, jobz, uplo, n, a, n, val);
    return;
}

/*!
 * This function calculate A(nxn) = X^T(1xn)*X(1xn)
 */
void multiply_vec2mat(double * x, double * a, int n)
{
//  cblas_dsyr(CblasRowMajor, CblasUpper, n, 1.0f, x, 1, a, n);
  int i, j;
  for(i=0; i<n; i++)
    for(j=0; j<=i; j++)
    {
      a[i*n+j] = a[j*n+i] = x[i]*x[j];
    }
}

/*!
 * This function calculates determinant of matrix A. \n
 * There are two versions in the internet, the main difference lies at dealing with the sign of det. \n
 * The version II is verifed to be \b INCORRECT. \n
 * Note that LAPACK is written in Fortran, the indix diffes by 1 with that in C.
 * LAPACK version 3.5.0
 * ********************************
 * Version I:
 * \code{.sh}
 *  det = 1.0;
 *  for(i=0; i<n; i++)
 *  {
 *   det *= a[i*n+i];
 *   if (ipiv[i] != i+1)
 *   {
 *     det = -det;
 *   }
 *  }
 * \endcode
 * ********************************
 * Version II:
 * \code{.sh}
 *  det = 1.0;
 *  for(i=0; i<n; i++)
 *  {
 *   det *= a[i*n+i];
 *   if (ipiv[i] != i+1)
 *   {
 *     ipiv[ipiv[i]-1] = ipiv[i];
 *     det = -det;
 *   }
 *  }
 * \endcode
 * ********************************
 */
double det_mat(double *a, int n, int *info)
{
  int * ipiv;
  int i;
  double det;
  ipiv=malloc(n*sizeof(int));

  /* LU decomposition */
//  dgetrf_(&n, &n, a, &n, ipiv, info);
  *info=LAPACKE_dgetrf(LAPACK_ROW_MAJOR, n, n, a, n, ipiv);
  if(*info!=0)
  {
    printf("# Error, Wrong in det_mat!\n");
    exit(-1);
  }

  det = 1.0;
  for(i=0; i<n; i++)
  {
    printf("%d\n", ipiv[i]);
    det *= a[i*n+i];
    if (ipiv[i] != i+1) // note that LAPACK is written in C, the indix diffes by 1 with C.
    {
      det = -det;
    }
  }
  //det=fabs(det);
  free(ipiv);
  return det;
}

/*!
 * This function calculates logarithm determinant of matrix.
 */
double lndet_mat(double *a, int n, int *info)
{
  int * ipiv;
  int i;
  double lndet;
  ipiv=malloc(n*sizeof(int));

  /* LU factorization */
//  dgetrf_(&n, &n, a, &n, ipiv, info);
  *info=LAPACKE_dgetrf(LAPACK_ROW_MAJOR, n, n, a, n, ipiv);
  if(*info!=0)
  {
    printf("Wrong!\n");
    exit(-1);
  }

  lndet = 0.0;
  for(i=0; i<n; i++)
  {
    lndet += log(a[i*n+i]);
  }
  free(ipiv);
  return lndet;
}

/*!
 * This function performs Cholesky decomposition of matrix into upper triangle matrixs
 */
void Chol_decomp_U(double *a, int n, int *info)
{
  int i,j;
  char uplo = 'U';   // decomposite as A = U^T*U, i.e., upper triangle matrix.
//  dpotrf_(&uplo, &n, a, &n, info);
  *info=LAPACKE_dpotrf(LAPACK_ROW_MAJOR, uplo, n, a, n);
  if(*info<0)
  {
    fprintf(stderr, "The %d-th argument had an illegal value!\n", *info);
//    exit(-1);
    return;
  }
  else if (*info>0)
  {
    fprintf(stderr, "The leading minor of order %d is not positive definite, and the factorization could not be completed.\n", *info);
//    exit(-1);
    return;
  }

  // only the upper triangle is referenced by dpotrf output,
  // so the strictly lower triangle must be set to zero
  for(i=0;i<n;i++)
    for(j=0;j<i;j++)
      a[i*n+j] = 0.0;
  return;
}

/*!
 * This function performs Cholesky decomposition of matrix into lower triangle matrixs
 */
void Chol_decomp_L(double *a, int n, int *info)
{
  int i,j;
  char uplo = 'L';   // decomposite as A = L*L^T, i.e., lower triangle matrix.
//  dpotrf_(&uplo, &n, a, &n, info);
  *info=LAPACKE_dpotrf(LAPACK_ROW_MAJOR, uplo, n, a, n);
  if(*info<0)
  {
    fprintf(stderr, "The %d-th argument had an illegal value!\n", *info);
//    exit(-1);
    return;
  }
  else if (*info>0)
  {
    fprintf(stderr, "The leading minor of order %d is not positive definite, and the factorization could not be completed.\n", *info);
//    exit(-1);
    return;
  }

  // only the lower triangle is referenced by dpotrf output,
  // so the strictly upper triangle must be set to zero
  for(i=0;i<n;i++)
    for(j=i+1;j<n;j++)
      a[i*n+j] = 0.0;
  return;
}

/*
 * semiseparable matrix
 */
void compute_semiseparable_drw(double *t, int n, double a1, double c1, double *sigma, double syserr, double *W, double *D, double *phi)
{
  int i;
  double S, A;
  phi[0] = 0.0;
  for(i=1; i<n; i++)
  {
    phi[i] = exp(-c1 * (t[i] - t[i-1]));
  }

  S = 0.0;
  A = sigma[0]*sigma[0] + syserr*syserr + a1;
  D[0] = A;
  W[0] = 1.0/D[0];
  for(i=1; i<n; i++)
  {
    S = phi[i]*phi[i] * (S + D[i-1]*W[i-1]*W[i-1]);
    A = sigma[i]*sigma[i] + syserr*syserr + a1;
    D[i] = A - a1 * a1 * S;
    W[i] = 1.0/D[i] * (1.0 - a1*S);
  }
}
/*
 * z = C^-1 x y
 *
 * y is a vector
 */
void multiply_matvec_semiseparable_drw(double *y, double  *W, double *D, double *phi, int n, double a1, double *z)
{
  int i;
  double f, g;

  // forward substitution
  f = 0.0;
  z[0] = y[0];
  for(i=1; i<n;i++)
  {
    f = phi[i] * (f + W[i-1] * z[i-1]);
    z[i] = y[i] - a1*f;
  }

  //backward substitution
  g = 0.0;
  z[n-1] = z[n-1]/D[n-1];
  for(i=n-2; i>=0; i--)
  {
    g = phi[i+1] *(g + a1*z[i+1]);
    z[i] = z[i]/D[i] - W[i]*g;
  }
}
/*
 * Z = C^-1 x Y
 *
 * Y is an (nxm) matrix.
 * Note that Y is row-major
 */
void multiply_mat_semiseparable_drw(double *Y, double  *W, double *D, double *phi, int n, int m, double a1, double *Z)
{
  int i, j;
  double f, g;

  // forward substitution
  for(j=0; j<m; j++)
  {
    f = 0.0;
    Z[0*m+j] = Y[0*m+j];
    for(i=1; i<n;i++)
    {
      f = phi[i] * (f + W[i-1] * Z[(i-1)*m + j]);
      Z[i*m+j] = Y[i*m+j] - a1*f;
    }
  }

  //backward substitution
  for(j=0; j<m; j++)
  {
    g = 0.0;
    Z[(n-1)*m+j] = Z[(n-1)*m+j]/D[n-1];
    for(i=n-2; i>=0; i--)
    {
      g = phi[i+1] *(g + a1*Z[(i+1)*m+j]);
      Z[i*m+j] = Z[i*m+j]/D[i] - W[i]*g;
    }
  }
}

/*
 * Z = C^-1 x Y^T
 *
 * Y is an (mxn) matrix.
 * Note that Y is row-major
 */
void multiply_mat_transposeB_semiseparable_drw(double *Y, double  *W, double *D, double *phi, int n, int m, double a1, double *Z)
{
  int i, j;
  double f, g;

  // forward substitution
  for(j=0; j<m; j++)
  {
    f = 0.0;
    Z[0*m+j] = Y[0+j*n];
    for(i=1; i<n;i++)
    {
      f = phi[i] * (f + W[i-1] * Z[(i-1)*m + j]);
      Z[i*m+j] = Y[i+j*n] - a1*f;
    }
  }

  //backward substitution
  for(j=0; j<m; j++)
  {
    g = 0.0;
    Z[(n-1)*m+j] = Z[(n-1)*m+j]/D[n-1];
    for(i=n-2; i>=0; i--)
    {
      g = phi[i+1] *(g + a1*Z[(i+1)*m+j]);
      Z[i*m+j] = Z[i*m+j]/D[i] - W[i]*g;
    }
  }
}

/**
 *  calculate A^-1.
 *
 *  A = LxDxL^T, L = I + tril(UxW^T), D is a diagonal matrix.
 *
 *  M = LxD^1/2,  A = MxM^T,  A^-1 = (M^T)^-1xM^-1.
 */
void inverse_semiseparable(double *t, int n, double a1, double c1, double *sigma,
                           double syserr, double *W, double *D, double *phi,
                           double *A, double *lndet)
{
  int i, j;

  compute_semiseparable_drw(t, n, a1, c1, sigma, syserr, W, D, phi);

  *lndet = 0.0;
  for(i=0; i<n; i++)
  {
    A[i*n + i] = 1.0 * sqrt(D[i]);
    for(j=0; j<i; j++)
    {
      A[i*n + j] = a1 * (exp(-c1*(t[i]-t[j]))*W[j]) * sqrt(D[j]);
    }

    *lndet += log(D[i]);
  }
  LAPACKE_dtrtri(LAPACK_ROW_MAJOR, 'L', 'N', n, A, n);

  LAPACKE_dlauum(LAPACK_ROW_MAJOR, 'L', n, A, n);

  /* fill up upper triangle */
  for(i=0; i<n; i++)
    for(j=i+1; j<n; j++)
      A[i*n+j] = A[j*n+i];
  return;
}

/**
 *  caclulate inverse of a DRW semiseparable matrix, which
 *  is a tridiagonal matrix.
 */
void inverse_semiseparable_uv(double *t, int n, double a1, double c1, double *A)
{
  int i, j;
  double b1, b2, dt;

  dt = t[1] - t[0];
  b1 = 1.0/( a1 * (exp(-c1*dt) - exp(c1*dt)) );
  A[0] = - b1 * exp(c1*dt);
  A[1] = A[1*n+0] = b1;
  for(i=2; i<n; i++)
    A[i] = 0.0;

  for(i=1; i<n-1; i++)
  {
    for(j=0; j<i-1; j++)
      A[i*n+j] = A[j*n+i] = 0.0;

    for(j=i+1; j<n; j++)
      A[i*n+j] = A[j*n+i] = 0.0;

    dt = t[i+1] - t[i];
    b2 = 1.0/( a1 * (exp(-c1*dt) - exp(c1*dt)) );
    dt = t[i+1] - t[i-1];
    A[i*n+i] = -b1 * b2 *  a1 * (exp(-c1*dt) - exp(c1*dt));

    A[i*n+(i+1)] = A[(i+1)*n + i] = b2;

    b1 = b2;
  }

  i = n-1;
  dt = t[i] - t[i-1];
  A[i*n+i] = - b1 * exp(c1 * dt);

  for(j=0; j<n-2; j++)
    A[i*n+j] = 0.0;
  return;
}

/**
 *  calculate Q = [S^-1+N^-1]^-1, where S is a DRW symmetric semiseparable matrix
 *  and N is a diagonal matrix.
 *
 *  S^-1 + N^-1 is a tridiagonal symmetric matrix.
 *
 *  [S^-1 + N^-1]^-1 is a semiseparable symmetric matrix.
 *
 */
void compute_inverse_semiseparable_plus_diag(double *t, int n, double a1, double c1,
                double *sigma, double syserr, double *u, double *v,
                double *W, double *D, double *phi, double *work)
{
  int i;
  double *a, *b, f1, f2;

  a = work;
  b = work + n;

  phi[0] = 0.0;
  for(i=1; i<n; i++)
  {
    phi[i] = exp(-c1 * (t[i] - t[i-1]));
  }

  /* first S^-1 + N^-1 */
  f1 = (phi[1]*phi[1] - 1.0) + EPS;
  b[0] = phi[1]/( a1 * f1);
  a[0] = - 1.0/( a1 * f1) + 1.0/(sigma[0]*sigma[0] + syserr*syserr);

  for(i=1; i<n-1; i++)
  {
    f2 = (phi[i+1]*phi[i+1] -1.0) + EPS;
    b[i] = phi[i+1]/( a1 *  f2);
    a[i] = -1.0/a1 * (1.0 + 1.0/f1 + 1.0/f2)  + 1.0/(sigma[i]*sigma[i] + syserr*syserr);

    f1 = f2;
  }
  i = n-1;
  a[i] = - 1.0 / (a1 * f1) + 1.0/(sigma[i]*sigma[i] + syserr*syserr);;

  /* now inverse [S^-1+N^-1]^-1*/
  v[0] = 1.0 ;
  v[1] = -a[0]/b[0] * phi[1];
  for(i=2; i<n; i++)
  {
    v[i] = -(a[i-1]*v[i-1] * phi[i]  + b[i-2]*v[i-2] * phi[i]*phi[i-1] )/b[i-1];
  }

  u[n-1] = 1.0/(b[n-2]*v[n-2] * phi[n-1] + a[n-1]*v[n-1]);
  for(i=n-2; i>0; i--)
  {
    u[i] = (1.0 - b[i]*v[i]*u[i+1] * phi[i+1] )/(a[i]*v[i]+b[i-1]*v[i-1]*phi[i]);
  }
  u[0] = (1.0-b[0]*v[0]*u[1] * phi[1] )/(a[0]*v[0]);

  /* calculate W, D */
  D[0] = u[0]*v[0];
  W[0] = 1.0/D[0];
  for(i=1; i<n; i++)
  {
    /*S = phi[i]*phi[i]*(S + D[i-1]*W[i-1]*W[i-1]);
    A = u[i]*v[i];
    D[i] = A - u[i]*u[i] * S;
    W[i] = 1.0/D[i] * (v[i] - u[i]*S);*/

    D[i] = -u[i]/u[i-1]/b[i-1] * phi[i];
    W[i] = 1.0/u[i];
    //printf("%e %e %e %e %e %e\n", u[i], v[i], D[i], W[i], a[i], b[i]);
  }
  return;
}

/*
 * z = C^1/2 x y
 *
 * C is a semiseparable matrix with (u,v) representation, y is a vector
 *
 * C = LxDxL^T,   L=I + tril(UxW^T)
 * C^1/2 = LxD^1/2
 */
void multiply_matvec_semiseparable_uv(double *y, double *u, double  *W, double *D, double *phi,
                                      int n, double *z)
{
  int i;
  double f;

  f = 0.0;
  z[0] = sqrt(D[0]) * y[0];
  for(i=1; i<n; i++)
  {
    f = phi[i] * (f + W[i-1] * sqrt(D[i-1]) * y[i-1]);
    z[i] = sqrt(D[i]) * y[i] + u[i] * f;
  }
  return;
}

/*!
 * This function display matrix on the screen.
 */
void display_mat(double *a, int m, int n)
{
    int i, j;
    for(i=0;i<m;i++)
    {
        for(j=0;j<n;j++)
        {
            printf("%e\t", a[i*n+j]);
        }
        printf("\n");
    }
}

/*!
 * This function allocates memory for matrix.
 */
double ** matrix_malloc(int n1, int n2)
{
  double ** mat;
  int i;

  if(!(mat = malloc(n1*sizeof(double*))))
  {
    fprintf(stderr, "Unable to allocate the matrix!\n");
    exit(-1);
  }

  for(i=0; i<n1; i++)
  {
    if(!(mat[i] = malloc(n2*sizeof(double))))
    {
      fprintf(stderr, "Unable to allocate the matrix!\n");
      exit(-1);
    }
  }
  return mat;
}

/*!
 * This function allocates memory for array.
 */
double * array_malloc(int n)
{
  double *array;

  if(!(array = malloc(n*sizeof(double))))
  {
    fprintf(stderr, "Unable to allocate the matrix!\n");
    exit(-1);
  }

  return array;
}

/*!
 * This function is to test the functions defined in mathfun.c.
 */
void test_mathfun()
{
  double *A, det;
  int n = 3, info;
  A = malloc(n*n*sizeof(double));

  A[0*n+0] = 1.0;
  A[0*n+1] = 10.3;
  A[0*n+2] = 6.3;

  A[1*n+0] = -8.3;
  A[1*n+1] = -8.0;
  A[1*n+2] = 5.3;

  A[2*n+0] = 0.3;
  A[2*n+1] = -9.3;
  A[2*n+2] = -3.0;

  det = det_mat(A, n, &info);
  printf("%f\n", det);
}

/* used for qsort sorting, ascending order */
int compare(const void* a, const void* b)
{
	if (*(double *)a > *(double *)b)
		return 1;
	else if (*(double *)a < *(double *)b)
		return -1;
	return 0;
}

max_stars_repo_name: LiyrAstroph/PyCALI

max_issues_repo_name: LiyrAstroph/pyCALI

max_forks_repo_name: LiyrAstroph/pyCALI


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

#pragma once

#include <Eigen/Geometry>
#include <gsl/gsl>

namespace mage
{
    inline void ToEuler(const Eigen::Quaternionf& quat, float& xRad, float& yRad, float& zRad)
    {
        Eigen::Matrix3f eigenMat = quat.toRotationMatrix();
        Eigen::Vector3f eulerAngles = eigenMat.eulerAngles(0, 1, 2);
        xRad = eulerAngles.x();
        yRad = eulerAngles.y();
        zRad = eulerAngles.z();
    }

    inline Eigen::Quaternionf FromEuler(float xRad, float yRad, float zRad)
    {
        return Eigen::AngleAxisf(xRad, Eigen::Vector3f::UnitX())
            * Eigen::AngleAxisf(yRad, Eigen::Vector3f::UnitY())
            * Eigen::AngleAxisf(zRad, Eigen::Vector3f::UnitZ());
    }

    template<int N>
    inline Eigen::Map<const Eigen::Matrix<float, N, 1>> ToCMap(gsl::span<const float, N> values)
    {
        return Eigen::Map<const Eigen::Matrix<float, N, 1>>{ values.data() };
    }

    template<int N>
    inline Eigen::Map<Eigen::Matrix<float, N, 1>> ToMap(gsl::span<float, N> values)
    {
        return Eigen::Map<const Eigen::Matrix<float, N, 1>>{ values.data() };
    }

    template<typename T>
    inline typename T::ConstMapType ToCMap(const T& element)
    {
        return T::ConstMapType{ element.data() };
    }

    template<typename T>
    inline typename T::MapType ToMap(T& element)
    {
        return T::MapType{ element.data() };
    }

}

ext: h

max_stars_repo_name: syntheticmagus/mageslam

max_issues_repo_name: syntheticmagus/mageslam

max_forks_repo_name: syntheticmagus/mageslam


text:
//This takes the output of STFT, X, and a transform matrix, T,
//and outputs a T*X or X*T', depending on orientation of X.
//X is either WxF or FxW.
//T must be BxF, as from get_spectrogram_T_mat.
//Y is either WxB or BxW.

#include <stdio.h>
#include <math.h>
#include <cblas.h>
//#include <time.h>

#ifdef __cplusplus
namespace codee {
extern "C" {
#endif

int apply_spectrogram_T_mat_s (float *Y, const float *X, const float *T, const size_t W, const size_t F, const size_t B);
int apply_spectrogram_T_mat_d (double *Y, const double *X, const double *T, const size_t W, const size_t F, const size_t B);

int apply_spectrogram_T_mat_s (float *Y, const float *X, const float *T, const size_t W, const size_t F, const size_t B)
{
    if (B<1u) { fprintf(stderr,"error in apply_spectrogram_T_mat_s: B (num output freq bands) must be positive\n"); return 1; }
    if (F<1u) { fprintf(stderr,"error in apply_spectrogram_T_mat_s: F (num STFT freqs) must be positive\n"); return 1; }

    if (B*W<1200000u)
    {
        //struct timespec tic, toc; clock_gettime(CLOCK_REALTIME,&tic);
        for (size_t w=W; w>0u; --w, X+=F, T-=B*F)
        {
            for (size_t b=B; b>0u; --b, X-=F, ++Y)
            {
                float sm = 0.0f;
                //size_t f = 0u;
                //while (f<F && *T==0.0f) { ++f; }
                //X += f; T += f;
                //while (f<F) { sm += *X * *T; ++X; ++T; ++f; }

                for (size_t f=F; f>0u; --f, ++X, ++T) { sm += *X * *T; }
                *Y = sm;
            }
        }
        //clock_gettime(CLOCK_REALTIME,&toc); fprintf(stderr,"fmaf: elapsed time = %.6f ms\n",(toc.tv_sec-tic.tv_sec)*1e3+(toc.tv_nsec-tic.tv_nsec)/1e6);
    }
    else
    {
        //struct timespec tic, toc; clock_gettime(CLOCK_REALTIME,&tic);
        cblas_sgemm(CblasRowMajor,CblasNoTrans,CblasTrans,(int)W,(int)B,(int)F,1.0f,X,(int)F,T,(int)F,0.0f,Y,(int)B);
        //clock_gettime(CLOCK_REALTIME,&toc); fprintf(stderr,"sgemm: elapsed time = %.6f ms\n",(toc.tv_sec-tic.tv_sec)*1e3+(toc.tv_nsec-tic.tv_nsec)/1e6);
    }

    return 0;
}

int apply_spectrogram_T_mat_d (double *Y, const double *X, const double *T, const size_t W, const size_t F, const size_t B)
{
    if (B<1u) { fprintf(stderr,"error in apply_spectrogram_T_mat_d: B (num output freq bands) must be positive\n"); return 1; }
    if (F<1u) { fprintf(stderr,"error in apply_spectrogram_T_mat_d: F (num STFT freqs) must be positive\n"); return 1; }

    if (B*W<12000u)
    {
        for (size_t w=W; w>0u; --w, X+=F, T-=B*F)
        {
            for (size_t b=B; b>0u; --b, X-=F, ++Y)
            {
                double sm = 0.0;
                //size_t f = F;
                //while (*T==0.0f && f>0u) { ++X; ++T; --f; }
                //while (*T!=0.0f && f>0u) { sm += *X**T; ++X; ++T; --f; }
                //X -= f; T -= f;
                //for (size_t f=F; f>0u; --f, ++X, ++T) { sm = fma(*X,*T,sm); }
                for (size_t f=F; f>0u; --f, ++X, ++T) { sm += *X * *T; }
                *Y = sm;
            }
        }
    }
    else
    {
        cblas_dgemm(CblasRowMajor,CblasNoTrans,CblasTrans,(int)W,(int)B,(int)F,1.0,X,(int)F,T,(int)F,0.0,Y,(int)B);
    }

    return 0;
}

#ifdef __cplusplus
}
}
#endif

ext: c


text:
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "sphere.h"
#include <petsc.h>
#undef __FUNCT__
#define __FUNCT__ "CartesianToSpherical"
PetscErrorCode CartesianToSpherical(Vec xyz, Vec sph)
{
  PetscErrorCode ierr;
  PetscScalar *xyzPt, *sphPt;
  PetscInt size;
  PetscReal x, y, z, r, theta, phi;
  PetscFunctionBegin;
  //get size of xyz vec
  ierr = VecGetSize(xyz, &size);CHKERRQ(ierr);

  ierr = VecGetArrayRead(xyz, &xyzPt);CHKERRQ(ierr);
  ierr = VecGetArray(sph, &sphPt);CHKERRQ(ierr);
  for(PetscInt k=0; k<size/3; ++k) {
    //get xyz points
    x = xyzPt[3*k+0];
    y = xyzPt[3*k+1];
    z = xyzPt[3*k+2];
    //compute transformation
    r = sqrt(x*x + y*y + z*z);
    theta = acos(z/r);
    phi = atan2(y,x);
    //set values in sph vec
    sphPt[3*k+0] = r;
    sphPt[3*k+1] = 0;
    sphPt[3*k+2] = phi;
    if(r>0) sphPt[3*k+1] = theta;

  }
  ierr = VecRestoreArray(sph, &sphPt);CHKERRQ(ierr);
  ierr = VecRestoreArrayRead(xyz, &xyzPt);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}
void cartesianToSpherical(SPoint *point) {

  double r, theta, phi;

  r = sqrt((point->x1)*(point->x1) + (point->x2)*(point->x2) + (point->x3)*(point->x3));
  theta = acos(point->x3/r);
  phi = atan2((point->x2),(point->x1));

  point->x1 = r;
  point->x2 = 0;
  point->x3 = phi;

  if(r > 0) point->x2 = theta;

  point->type = 's';
}

void sphericalToCartesian(SPoint *point) {
  double r = point->x1;
  double theta = point->x2;
  double phi = point->x3;

  double x = r*sin(theta)*cos(phi);
  double y = r*sin(theta)*sin(phi);
  double z = r*cos(theta);

  point->x1 = x;
  point->x2 = y;
  point->x3 = z;
}

max_stars_repo_name: tom-klotz/ellipsoid-solvation

max_issues_repo_name: tom-klotz/ellipsoid-solvation

max_forks_repo_name: tom-klotz/ellipsoid-solvation


text:
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

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

#ifdef VRNA_WITH_GSL
#include <gsl/gsl_multimin.h>
#endif

#include "ViennaRNA/eval.h"
#include "ViennaRNA/fold_vars.h"
#include "ViennaRNA/constraints/hard.h"
#include "ViennaRNA/constraints/soft.h"
#include "ViennaRNA/fold.h"
#include "ViennaRNA/part_func.h"
#include "ViennaRNA/utils/basic.h"
#include "ViennaRNA/params/basic.h"
#include "ViennaRNA/perturbation_fold.h"

static void
calculate_probability_unpaired(vrna_fold_compound_t *vc,
                               double               *probability)
{
  int         length = vc->length;
  FLT_OR_DBL  *probs = vc->exp_matrices->probs;
  int         *iidx = vc->iindx;
  int         i, j;

  for (i = 0; i <= length; ++i)
    probability[i] = 1;

  for (i = 1; i <= length; ++i)
    for (j = i + 1; j <= length; ++j) {
      probability[i]  -= probs[iidx[i] - j];
      probability[j]  -= probs[iidx[i] - j];
    }
}

#if 0
static double
calculate_norm(double *vector,
               int    length)
{
  double  sum = 0;
  int     i;

  for (i = 1; i <= length; ++i)
    sum += vector[i] * vector[i];

  return sqrt(sum);
}

#endif

static void
addSoftConstraint(vrna_fold_compound_t  *vc,
                  const double          *epsilon,
                  int                   length)
{
  /* remove previous soft constraints */
  vrna_sc_init(vc);

  /* prepare vector of unpaired constraints in kcal/mol */
  FLT_OR_DBL *constraints = (FLT_OR_DBL *)vrna_alloc(sizeof(FLT_OR_DBL) * (length + 1));

  memcpy(constraints + 1, epsilon + 1, sizeof(FLT_OR_DBL) * length);

  /* add new soft constraints */
  vrna_sc_set_up(vc, (const FLT_OR_DBL *)constraints, VRNA_OPTION_DEFAULT);

  free(constraints);
}

static double
evaluate_objective_function_contribution(double value,
                                         int    objective_function)
{
  if (objective_function == VRNA_OBJECTIVE_FUNCTION_QUADRATIC)
    return value * value;

  if (objective_function == VRNA_OBJECTIVE_FUNCTION_ABSOLUTE)
    return fabs(value);

  assert(0);
  return 0;
}

static double
evaluate_perturbation_vector_score(vrna_fold_compound_t *vc,
                                   const double         *epsilon,
                                   const double         *q_prob_unpaired,
                                   double               sigma_squared,
                                   double               tau_squared,
                                   int                  objective_function)
{
  double  ret   = 0;
  double  ret2  = 0.;
  double  *p_prob_unpaired;
  int     i;
  int     length = vc->length;

  /* calculate pairing probabilty in the pertubated energy model */
  p_prob_unpaired = vrna_alloc(sizeof(double) * (length + 1));

  addSoftConstraint(vc, epsilon, length);

  vc->params->model_details.compute_bpp     = 1;
  vc->exp_params->model_details.compute_bpp = 1;

  /* get new (constrained) MFE to scale pf computations properly */
  double mfe = (double)vrna_mfe(vc, NULL);

  vrna_exp_params_rescale(vc, &mfe);

  vrna_pf(vc, NULL);

  calculate_probability_unpaired(vc, p_prob_unpaired);

  vrna_sc_remove(vc);

  for (i = 1; i <= length; ++i) {
    /* add penalty for pertubation energies */
    ret += evaluate_objective_function_contribution(epsilon[i], objective_function) / tau_squared;

    /* add penalty for mismatches between observed and predicted probabilities */
    if (q_prob_unpaired[i] >= 0) /* ignore positions with missing data */
      ret2 += evaluate_objective_function_contribution(p_prob_unpaired[i] - q_prob_unpaired[i],
                                                       objective_function) / sigma_squared;
  }

  vrna_message_info(stderr, "Score: pertubation: %g\tdiscrepancy: %g", ret, ret2);
  free(p_prob_unpaired);

  return ret + ret2;
}

static void
pairing_probabilities_from_restricted_pf(vrna_fold_compound_t *vc,
                                         const double         *epsilon,
                                         double               *prob_unpaired,
                                         double               **conditional_prob_unpaired)
{
  int length = vc->length;
  int i;

  addSoftConstraint(vc, epsilon, length);
  vc->params->model_details.compute_bpp     = 1;
  vc->exp_params->model_details.compute_bpp = 1;

  /* get new (constrained) MFE to scale pf computations properly */
  double mfe = (double)vrna_mfe(vc, NULL);

  vrna_exp_params_rescale(vc, &mfe);

  vrna_pf(vc, NULL);

  calculate_probability_unpaired(vc, prob_unpaired);

#ifdef _OPENMP
#pragma omp parallel for private(i)
#endif
  for (i = 1; i <= length; ++i) {
    vrna_fold_compound_t  *restricted_vc;
    char                  *hc_string;
    unsigned int          constraint_options = VRNA_CONSTRAINT_DB
                                               | VRNA_CONSTRAINT_DB_PIPE
                                               | VRNA_CONSTRAINT_DB_DOT
                                               | VRNA_CONSTRAINT_DB_X
                                               | VRNA_CONSTRAINT_DB_ANG_BRACK
                                               | VRNA_CONSTRAINT_DB_RND_BRACK;

    hc_string = vrna_alloc(sizeof(char) * (length + 1));
    memset(hc_string, '.', length);
    hc_string[i - 1] = 'x';

    restricted_vc = vrna_fold_compound(vc->sequence,
                                       &(vc->exp_params->model_details),
                                       VRNA_OPTION_DEFAULT);
    vrna_constraints_add(restricted_vc, hc_string, constraint_options);
    free(hc_string);

    vrna_exp_params_subst(restricted_vc, vc->exp_params);

    vrna_pf(restricted_vc, NULL);
    calculate_probability_unpaired(restricted_vc, conditional_prob_unpaired[i]);

    restricted_vc->sc = NULL;
    vrna_fold_compound_free(restricted_vc);
  }

  vrna_sc_remove(vc);
}

static void
pairing_probabilities_from_sampling(vrna_fold_compound_t  *vc,
                                    const double          *epsilon,
                                    int                   sample_size,
                                    double                *prob_unpaired,
                                    double                **conditional_prob_unpaired,
                                    unsigned int          options)
{
  char    **samples, **ptr;
  int     length, i, j;
  double  mfe;

  length = vc->length;
  addSoftConstraint(vc, epsilon, length);

  vc->params->model_details.compute_bpp     = 0;
  vc->exp_params->model_details.compute_bpp = 0;

  /* get new (constrained) MFE to scale pf computations properly */
  mfe = (double)vrna_mfe(vc, NULL);
  vrna_exp_params_rescale(vc, &mfe);

  vrna_pf(vc, NULL);

  samples = vrna_pbacktrack_num(vc,
                                (unsigned int)sample_size,
                                options);

  for (ptr = samples; (*ptr); ptr++) {
    for (i = length; i > 0; i--) {
      if ((*ptr)[i - 1] == '.') {
        ++prob_unpaired[i];

        for (j = length; j > 0; j--)
          if ((*ptr)[j - 1] == '.')
            ++conditional_prob_unpaired[i][j];
      }
    }

    free(*ptr);
  }

  free(samples);

  for (i = 1; i <= length; ++i) {
    if (prob_unpaired[i])
      for (j = 1; j <= length; ++j)
        conditional_prob_unpaired[i][j] /= prob_unpaired[i];

    prob_unpaired[i] /= sample_size;

    assert(prob_unpaired[i] >= 0 && prob_unpaired[i] <= 1);
  }

  vrna_sc_remove(vc);
}

static void
allocateProbabilityArrays(double  **unpaired,
                          double  ***conditional_unpaired,
                          int     length)
{
  int i;

  *unpaired             = vrna_alloc(sizeof(double) * (length + 1));
  *conditional_unpaired = vrna_alloc(sizeof(double *) * (length + 1));

  for (i = 1; i <= length; ++i)
    (*conditional_unpaired)[i] = vrna_alloc(sizeof(double) * (length + 1));
}

static void
freeProbabilityArrays(double  *unpaired,
                      double  **conditional_unpaired,
                      int     length)
{
  int i;

  free(unpaired);
  for (i = 1; i <= length; ++i)
    free(conditional_unpaired[i]);
  free(conditional_unpaired);
}

static void
evaluate_perturbation_vector_gradient(vrna_fold_compound_t  *vc,
                                      const double          *epsilon,
                                      const double          *q_prob_unpaired,
                                      double                sigma_squared,
                                      double                tau_squared,
                                      int                   objective_function,
                                      int                   sample_size,
                                      double                *gradient)
{
  double  *p_prob_unpaired;
  double  **p_conditional_prob_unpaired;
  int     i, mu;
  int     length  = vc->length;
  double  kT      = vc->exp_params->kT / 1000;

  allocateProbabilityArrays(&p_prob_unpaired, &p_conditional_prob_unpaired, length);

  if (sample_size > 0) {
    pairing_probabilities_from_sampling(vc,
                                        epsilon,
                                        sample_size,
                                        p_prob_unpaired,
                                        p_conditional_prob_unpaired,
                                        VRNA_PBACKTRACK_DEFAULT);
  } else if (sample_size < 0) {
    pairing_probabilities_from_sampling(vc,
                                        epsilon,
                                        -sample_size,
                                        p_prob_unpaired,
                                        p_conditional_prob_unpaired,
                                        VRNA_PBACKTRACK_NON_REDUNDANT);
  } else {
    pairing_probabilities_from_restricted_pf(vc,
                                             epsilon,
                                             p_prob_unpaired,
                                             p_conditional_prob_unpaired);
  }

  for (mu = 1; mu <= length; ++mu) {
    double sum = 0;

    if (objective_function == VRNA_OBJECTIVE_FUNCTION_QUADRATIC) {
      for (i = 1; i <= length; ++i) {
        if (q_prob_unpaired[i] < 0) /* ignore positions with missing data */
          continue;

        sum += (p_prob_unpaired[i] - q_prob_unpaired[i])
               * p_prob_unpaired[i] * (p_prob_unpaired[mu] - p_conditional_prob_unpaired[i][mu])
               / sigma_squared;
      }

      gradient[mu] = 2 * (epsilon[mu] / tau_squared + sum / kT);
    } else if (objective_function == VRNA_OBJECTIVE_FUNCTION_ABSOLUTE) {
      for (i = 1; i <= length; ++i)
        if (q_prob_unpaired[i] >= 0 && p_prob_unpaired[i] != q_prob_unpaired[i]) {
          sum += (p_prob_unpaired[i] * (p_prob_unpaired[mu] - p_conditional_prob_unpaired[i][mu])) /
                 kT
                 / sigma_squared
                 * (p_prob_unpaired[i] > q_prob_unpaired[i] ? 1. : -1.);
        }

      if (epsilon[mu])
        sum += (epsilon[mu] > 0 ? 1. : -1.) / tau_squared;

      gradient[mu] = sum;
    }
  }

  freeProbabilityArrays(p_prob_unpaired, p_conditional_prob_unpaired, length);
}

#ifdef VRNA_WITH_GSL
typedef struct parameters_gsl {
  vrna_fold_compound_t  *vc;
  const double          *q_prob_unpaired;
  double                sigma_squared;
  double                tau_squared;
  int                   objective_function;
  int                   sample_size;
} parameters_gsl;

static double
f_gsl(const gsl_vector  *x,
      void              *params)
{
  parameters_gsl *p = params;

  return evaluate_perturbation_vector_score(p->vc,
                                            x->data,
                                            p->q_prob_unpaired,
                                            p->sigma_squared,
                                            p->tau_squared,
                                            p->objective_function);
}

static void
df_gsl(const gsl_vector *x,
       void             *params,
       gsl_vector       *df)
{
  parameters_gsl *p = params;

  gsl_vector_set(df, 0, 0);
  evaluate_perturbation_vector_gradient(p->vc,
                                        x->data,
                                        p->q_prob_unpaired,
                                        p->sigma_squared,
                                        p->tau_squared,
                                        p->objective_function,
                                        p->sample_size,
                                        df->data);
}

static void
fdf_gsl(const gsl_vector  *x,
        void              *params,
        double            *f,
        gsl_vector        *g)
{
  *f = f_gsl(x, params);
  df_gsl(x, params, g);
}

#endif /* VRNA_WITH_GSL */

PUBLIC void
vrna_sc_minimize_pertubation(vrna_fold_compound_t *vc,
                             const double         *q_prob_unpaired,
                             int                  objective_function,
                             double               sigma_squared,
                             double               tau_squared,
                             int                  algorithm,
                             int                  sample_size,
                             double               *epsilon,
                             double               initialStepSize,
                             double               minStepSize,
                             double               minImprovement,
                             double               minimizerTolerance,
                             progress_callback    callback)
{
  int                                   iteration_count = 0;
  const int                             max_iterations  = 100;
  int                                   length          = vc->length;

#ifdef VRNA_WITH_GSL
  const gsl_multimin_fdfminimizer_type  *minimizer_type = 0;

  struct {
    int                                   type;
    const gsl_multimin_fdfminimizer_type  *gsl_type;
  } algorithms[] =
  { { VRNA_MINIMIZER_CONJUGATE_FR,
      gsl_multimin_fdfminimizer_conjugate_fr },
    { VRNA_MINIMIZER_CONJUGATE_PR,
      gsl_multimin_fdfminimizer_conjugate_pr },
    { VRNA_MINIMIZER_VECTOR_BFGS,
      gsl_multimin_fdfminimizer_vector_bfgs },
    { VRNA_MINIMIZER_VECTOR_BFGS2,
      gsl_multimin_fdfminimizer_vector_bfgs2 },
    { VRNA_MINIMIZER_STEEPEST_DESCENT,
      gsl_multimin_fdfminimizer_steepest_descent },
    { 0,
      NULL } };
  int i;
  for (i = 0; algorithms[i].type; ++i)
    if (algorithms[i].type == algorithm) {
      minimizer_type = algorithms[i].gsl_type;
      break;
    }

  if (minimizer_type) {
    parameters_gsl            parameters;
    gsl_multimin_function_fdf fdf;
    gsl_multimin_fdfminimizer *minimizer;
    gsl_vector                *vector;

    int                       status;

    parameters.vc                 = vc;
    parameters.q_prob_unpaired    = q_prob_unpaired;
    parameters.sigma_squared      = sigma_squared;
    parameters.tau_squared        = tau_squared;
    parameters.objective_function = objective_function;
    parameters.sample_size        = sample_size;

    fdf.n       = length + 1;
    fdf.f       = &f_gsl;
    fdf.df      = &df_gsl;
    fdf.fdf     = &fdf_gsl;
    fdf.params  = (void *)&parameters;

    minimizer = gsl_multimin_fdfminimizer_alloc(minimizer_type, length + 1);
    vector    = gsl_vector_calloc(length + 1);

    /* gsl_multimin_fdfminimizer_set(minimizer, &fdf, vector, 0.01, 1e-4); */
    gsl_multimin_fdfminimizer_set(minimizer, &fdf, vector, initialStepSize, minimizerTolerance);

    if (callback)
      callback(0, minimizer->f, minimizer->x->data);

    do {
      ++iteration_count;
      status = gsl_multimin_fdfminimizer_iterate(minimizer);

      if (callback)
        callback(iteration_count, minimizer->f, minimizer->x->data);

      if (status)
        break;

      status = gsl_multimin_test_gradient(minimizer->gradient, minimizerTolerance);
    } while (status == GSL_CONTINUE && iteration_count < max_iterations);

    memcpy(epsilon, minimizer->x->data, sizeof(double) * (length + 1));

    gsl_multimin_fdfminimizer_free(minimizer);
    gsl_vector_free(vector);

    return;
  }

#endif /* VRNA_WITH_GSL */

  double        improvement;
  const double  min_improvement = minImprovement;

  double        *new_epsilon  = vrna_alloc(sizeof(double) * (length + 1));
  double        *gradient     = vrna_alloc(sizeof(double) * (length + 1));

  double        score = evaluate_perturbation_vector_score(vc,
                                                           epsilon,
                                                           q_prob_unpaired,
                                                           sigma_squared,
                                                           tau_squared,
                                                           objective_function);

  if (callback)
    callback(0, score, epsilon);

  do {
    double  new_score;
    double  step_size;

    ++iteration_count;

    evaluate_perturbation_vector_gradient(vc,
                                          epsilon,
                                          q_prob_unpaired,
                                          sigma_squared,
                                          tau_squared,
                                          objective_function,
                                          sample_size,
                                          gradient);

    /*    step_size = 0.5 / calculate_norm(gradient, length);*/
    step_size = initialStepSize;

    do {
      int i;
      for (i = 1; i <= length; ++i)
        new_epsilon[i] = epsilon[i] - step_size * gradient[i];

      new_score = evaluate_perturbation_vector_score(vc,
                                                     new_epsilon,
                                                     q_prob_unpaired,
                                                     sigma_squared,
                                                     tau_squared,
                                                     objective_function);
      improvement = 1 - new_score / score;
      step_size   /= 2;
    } while ((improvement < min_improvement) && (step_size >= minStepSize));

    if (new_score > score)
      break;

    if (callback)
      callback(iteration_count, new_score, new_epsilon);

    score = new_score;
    memcpy(epsilon, new_epsilon, sizeof(double) * (length + 1));
  } while (improvement >= min_improvement && iteration_count < max_iterations);

  free(gradient);
  free(new_epsilon);
}

max_stars_repo_name: tsjzz/ViennaRNA

max_issues_repo_name: tsjzz/ViennaRNA

max_forks_repo_name: tsjzz/ViennaRNA


text:
/* Copyright (c) 2016 PaddlePaddle Authors. All Rights Reserve.

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License. */

#pragma once
#ifdef PADDLE_USE_MKLML
#include <mkl_cblas.h>
#include <mkl_lapacke.h>
#include <mkl_vml_functions.h>
#endif

#ifdef PADDLE_USE_MKL
#include <mkl.h>
#include <mkl_lapacke.h>
#endif

#ifdef PADDLE_USE_ATLAS
extern "C" {
#include <cblas.h>
#include <clapack.h>
}
#endif

#ifdef PADDLE_USE_OPENBLAS
#include <cblas.h>
#include <lapacke.h>
#endif

#ifndef LAPACK_FOUND
extern "C" {
#include <cblas.h>
int LAPACKE_sgetrf(int matrix_layout, int m, int n, float* a, int lda,
                   int* ipiv);
int LAPACKE_dgetrf(int matrix_layout, int m, int n, double* a, int lda,
                   int* ipiv);
int LAPACKE_sgetri(int matrix_layout, int n, float* a, int lda,
                   const int* ipiv);
int LAPACKE_dgetri(int matrix_layout, int n, double* a, int lda,
                   const int* ipiv);
}
#endif

#include <cmath>

#include "paddle/framework/tensor.h"
#include "paddle/platform/device_context.h"
#include "paddle/platform/enforce.h"

namespace paddle {
namespace operators {
namespace math {

// Support continuous memory now
// If transA = N, and transB = N
// Then matrixA: M * K, matrixB: K * N matrixC : M * N
// For more detailed info, please refer to
// http://www.netlib.org/lapack/explore-html/d4/de2/sgemm_8f.html
template <typename Place, typename T>
void gemm(const platform::DeviceContext& context, const CBLAS_TRANSPOSE transA,
          const CBLAS_TRANSPOSE transB, const int M, const int N, const int K,
          const T alpha, const T* A, const T* B, const T beta, T* C);

// matrix multiply with continuous memory
template <typename Place, typename T>
void matmul(const platform::DeviceContext& context,
            const framework::Tensor& matrix_a, bool trans_a,
            const framework::Tensor& matrix_b, bool trans_b, T alpha,
            framework::Tensor* matrix_out, T beta);

}  // namespace math
}  // namespace operators
}  // namespace paddle

max_stars_repo_name: AI-books/Paddle

max_issues_repo_name: AI-books/Paddle

max_forks_repo_name: AI-books/Paddle


text:
#ifndef __INC_OPTSTRUCT__
#define __INC_OPTSTRUCT__

#include <gsl/gsl_matrix.h>

/**
 * @file optstruct.h
 * \brief defines the optstruct which holds all the mundane options and dimensions etc
 */

#define POWEREXPCOVFN 1
#define MATERN32 2
#define MATERN52 3

/**
 * \struct optstruct
 * \brief holds the main parameters etc needed by all the emulator fns
 *
 * this is the most important structure after the modelstruct, this contains
 * all the dimensions of xmodel, training_vec etc
 *
 */
typedef struct optstruct{
	/**
	 * the number of parameters in the covariance fn
	 */
	int nthetas;
	/**
	 * the number of parameters in the model to be emulated, this sets the dimensionality
	 * of the estimation space and the "hardness" of the estimation process
	 */
	int nparams;
	/**
	 * the number of locations in the model space that we have training values at
	 * another contributing factor to how complex the problem becomes, the covariance matrix
	 * scales as (nmodel_points**2) and we need to invert this every step of the estimation
	 * process
	 */
	int nmodel_points;
	/**
	 * how many points to evaluate the emulated model at
	 */
	int nemulate_points;

	/**
	 * what order (if any) should the regression model be
	 * 0 -> a single constant value (a0)
	 * 1 -> a constant plus a slope vector (a0 + a1*x)
	 * 2 -> (a0 + a1*x + a2*x^2)
	 * 3 -> (a0 + a1*x + a2*x^2 + a3*x^3)
	 * >3 (not supported)
	 */
	int regression_order;

	/**
	 * indirectly controls the shape of the regression model, be very careful with this one
	 * this is set correctly by calling setup_regression
	 */
	int nregression_fns;

	/**
	 * fixed nugget?
	 * allow the user to try and supply a nugget which will not be optimized over, but fixed
	 */
	int fixed_nugget_mode;
	double fixed_nugget;

	/**
	 * set which cov fn to use, can be one of
	 * POWEREXPCOVFN
	 * MATERN32
	 * MATERN52
	 *
	 * this is determined by setup_cov_fn
	 */
	int cov_fn_index;

	// set this to not zero if you want to use the length scales set by the data
	int use_data_scales;

	/** this holds the ranges for the optimisation routine*/
	gsl_matrix* grad_ranges;

} optstruct;

void free_optstruct(optstruct *opts);
void copy_optstruct(optstruct *dst, optstruct* src);
void dump_optstruct(FILE *fptr, optstruct* opts);
void load_optstruct(FILE *fptr, optstruct* opts);

void setup_cov_fn(optstruct *opts);
void setup_regression(optstruct *opts);

#include "modelstruct.h"

void setup_optimization_ranges(optstruct* options, modelstruct* the_model);

#endif

max_stars_repo_name: jackdawjackdaw/emulator

max_issues_repo_name: MADAI/MADAIEmulator

max_forks_repo_name: MADAI/MADAIEmulator


text:
#pragma once

#include <queue>
#include <vector>

#include <gsl/gsl>

#define EMPTY_NODE '\0'

class minqueue_node {
    public:
    minqueue_node(char value, unsigned int score);
    ~minqueue_node();

    unsigned int score;
    char value;

    minqueue_node* left;
    minqueue_node* right;
};

struct cmp_minqueue_nodes {
    bool operator()(
        const gsl::not_null<minqueue_node*> left, const gsl::not_null<minqueue_node*> right);
};

class minqueue {
    std::priority_queue<
        minqueue_node*, std::vector<minqueue_node*>, cmp_minqueue_nodes> nodes;

    public:
    minqueue();
    ~minqueue();

    void push(gsl::not_null<minqueue_node*> node);
    minqueue_node* pop();

    bool empty();
    int size();
};

max_stars_repo_name: joshleeb/Huffman

max_issues_repo_name: joshleeb/Huffman

max_forks_repo_name: joshleeb/Huffman


text:
#include "cosmocalc.h"
#include "weaklens.h"

cosmocalcData cosmoData;
weaklensData wlData;

#include <gsl/gsl_errno.h>
void turn_off_gsl_errs(void) {
  gsl_set_error_handler_off();
}

max_stars_repo_name: beckermr/cosmocalc

max_issues_repo_name: beckermr/cosmocalc

max_forks_repo_name: beckermr/cosmocalc


text:
//This is a simple example meant to test your capabilities with basic library linking functionality. All you have to do is return the output to modelica.

#include <stdio.h>
#include <gsl/gsl_sf_bessel.h>

double sf_bessel(double x)
{
	double y = gsl_sf_bessel_J0 (x);
	return y;
}

//The main function has been provided solely for the sake of testing the function. It will be removed when the problem set is formally passed on.
/*
int main(void)
{
	double x = 5;
	double a = sf_bessel(x);

	printf("a = %f \n",a);
	return 0;
}
*/

max_stars_repo_name: RITUait/OpenModelica

max_issues_repo_name: RITUait/OpenModelica

max_forks_repo_name: RITUait/OpenModelica


text:
/* bspline/bspline.c
 *
 * Copyright (C) 2006, 2007, 2008, 2009 Patrick Alken
 * Copyright (C) 2008 Rhys Ulerich
 *
 * 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_errno.h>
#include <gsl/gsl_bspline.h>

/*
 * This module contains routines related to calculating B-splines.
 * The algorithms used are described in
 *
 * [1] Carl de Boor, "A Practical Guide to Splines", Springer
 *     Verlag, 1978.
 *
 * The bspline_pppack_* internal routines contain code adapted from
 *
 * [2] "PPPACK - Piecewise Polynomial Package",
 *     http://www.netlib.org/pppack/
 *
 */

#include "bspline.h"

/*
gsl_bspline_alloc()
  Allocate space for a bspline workspace. The size of the
workspace is O(5k + nbreak)

Inputs: k      - spline order (cubic = 4)
        nbreak - number of breakpoints

Return: pointer to workspace
*/

gsl_bspline_workspace *
gsl_bspline_alloc (const size_t k, const size_t nbreak)
{
  if (k == 0)
    {
      GSL_ERROR_NULL ("k must be at least 1", GSL_EINVAL);
    }
  else if (nbreak < 2)
    {
      GSL_ERROR_NULL ("nbreak must be at least 2", GSL_EINVAL);
    }
  else
    {
      gsl_bspline_workspace *w;

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

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

      w->k = k;
      w->km1 = k - 1;
      w->nbreak = nbreak;
      w->l = nbreak - 1;
      w->n = w->l + k - 1;

      w->knots = gsl_vector_alloc (w->n + k);
      if (w->knots == 0)
	{
	  free (w);
	  GSL_ERROR_NULL ("failed to allocate space for knots vector",
			  GSL_ENOMEM);
	}

      w->deltal = gsl_vector_alloc (k);
      if (w->deltal == 0)
	{
	  gsl_vector_free (w->knots);
	  free (w);
	  GSL_ERROR_NULL ("failed to allocate space for deltal vector",
			  GSL_ENOMEM);
	}

      w->deltar = gsl_vector_alloc (k);
      if (w->deltar == 0)
	{
	  gsl_vector_free (w->deltal);
	  gsl_vector_free (w->knots);
	  free (w);
	  GSL_ERROR_NULL ("failed to allocate space for deltar vector",
			  GSL_ENOMEM);
	}

      w->B = gsl_vector_alloc (k);
      if (w->B == 0)
	{
	  gsl_vector_free (w->deltar);;
	  gsl_vector_free (w->deltal);
	  gsl_vector_free (w->knots);
	  free (w);
	  GSL_ERROR_NULL
	    ("failed to allocate space for temporary spline vector",
	     GSL_ENOMEM);
	}

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

/*
gsl_bspline_deriv_alloc()
  Allocate space for a bspline derivative workspace. The size of the
workspace is O(2k^2)

Inputs: k      - spline order (cubic = 4)

Return: pointer to workspace
*/

gsl_bspline_deriv_workspace *
gsl_bspline_deriv_alloc (const size_t k)
{
  if (k == 0)
    {
      GSL_ERROR_NULL ("k must be at least 1", GSL_EINVAL);
    }
  else
    {
      gsl_bspline_deriv_workspace *dw;

      dw =
	(gsl_bspline_deriv_workspace *)
	malloc (sizeof (gsl_bspline_deriv_workspace));

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

      dw->A = gsl_matrix_alloc (k, k);
      if (dw->A == 0)
	{
	  free (dw);
	  GSL_ERROR_NULL
	    ("failed to allocate space for derivative work matrix",
	     GSL_ENOMEM);
	}

      dw->dB = gsl_matrix_alloc (k, k + 1);
      if (dw->dB == 0)
	{
	  gsl_matrix_free (dw->A);
	  free (dw);
	  GSL_ERROR_NULL
	    ("failed to allocate space for temporary derivative matrix",
	     GSL_ENOMEM);
	}

      dw->k = k;

      return dw;
    }
}				/* gsl_bspline_deriv_alloc() */

/* Return number of coefficients */
size_t
gsl_bspline_ncoeffs (gsl_bspline_workspace * w)
{
  return w->n;
}

/* Return order */
size_t
gsl_bspline_order (gsl_bspline_workspace * w)
{
  return w->k;
}

/* Return number of breakpoints */
size_t
gsl_bspline_nbreak (gsl_bspline_workspace * w)
{
  return w->nbreak;
}

/* Return the location of the i-th breakpoint*/
double
gsl_bspline_breakpoint (size_t i, gsl_bspline_workspace * w)
{
  size_t j = i + w->k - 1;
  return gsl_vector_get (w->knots, j);
}

/*
gsl_bspline_free()
  Free a gsl_bspline_workspace.

Inputs: w - workspace to free

Return: none
*/

void
gsl_bspline_free (gsl_bspline_workspace * w)
{
  RETURN_IF_NULL (w);
  gsl_vector_free (w->knots);
  gsl_vector_free (w->deltal);
  gsl_vector_free (w->deltar);
  gsl_vector_free (w->B);
  free (w);
}				/* gsl_bspline_free() */

/*
gsl_bspline_deriv_free()
  Free a gsl_bspline_deriv_workspace.

Inputs: dw - workspace to free

Return: none
*/

void
gsl_bspline_deriv_free (gsl_bspline_deriv_workspace * dw)
{
  RETURN_IF_NULL (dw);
  gsl_matrix_free (dw->A);
  gsl_matrix_free (dw->dB);
  free (dw);
}				/* gsl_bspline_deriv_free() */

/*
gsl_bspline_knots()
  Compute the knots from the given breakpoints:

   knots(1:k) = breakpts(1)
   knots(k+1:k+l-1) = breakpts(i), i = 2 .. l
   knots(n+1:n+k) = breakpts(l + 1)

where l is the number of polynomial pieces (l = nbreak - 1) and
   n = k + l - 1
(using matlab syntax for the arrays)

The repeated knots at the beginning and end of the interval
correspond to the continuity condition there. See pg. 119
of [1].

Inputs: breakpts - breakpoints
        w        - bspline workspace

Return: success or error
*/

int
gsl_bspline_knots (const gsl_vector * breakpts, gsl_bspline_workspace * w)
{
  if (breakpts->size != w->nbreak)
    {
      GSL_ERROR ("breakpts vector has wrong size", GSL_EBADLEN);
    }
  else
    {
      size_t i;			/* looping */

      for (i = 0; i < w->k; i++)
	gsl_vector_set (w->knots, i, gsl_vector_get (breakpts, 0));

      for (i = 1; i < w->l; i++)
	{
	  gsl_vector_set (w->knots, w->k - 1 + i,
			  gsl_vector_get (breakpts, i));
	}

      for (i = w->n; i < w->n + w->k; i++)
	gsl_vector_set (w->knots, i, gsl_vector_get (breakpts, w->l));

      return GSL_SUCCESS;
    }
}				/* gsl_bspline_knots() */

/*
gsl_bspline_knots_uniform()
  Construct uniformly spaced knots on the interval [a,b] using
the previously specified number of breakpoints. 'a' is the position
of the first breakpoint and 'b' is the position of the last
breakpoint.

Inputs: a - left side of interval
        b - right side of interval
        w - bspline workspace

Return: success or error

Notes: 1) w->knots is modified to contain the uniformly spaced
          knots

       2) The knots vector is set up as follows (using octave syntax):

          knots(1:k) = a
          knots(k+1:k+l-1) = a + i*delta, i = 1 .. l - 1
          knots(n+1:n+k) = b
*/

int
gsl_bspline_knots_uniform (const double a, const double b,
			   gsl_bspline_workspace * w)
{
  size_t i;			/* looping */
  double delta;			/* interval spacing */
  double x;

  delta = (b - a) / (double) w->l;

  for (i = 0; i < w->k; i++)
    gsl_vector_set (w->knots, i, a);

  x = a + delta;
  for (i = 0; i < w->l - 1; i++)
    {
      gsl_vector_set (w->knots, w->k + i, x);
      x += delta;
    }

  for (i = w->n; i < w->n + w->k; i++)
    gsl_vector_set (w->knots, i, b);

  return GSL_SUCCESS;
}				/* gsl_bspline_knots_uniform() */

/*
gsl_bspline_eval()
  Evaluate the basis functions B_i(x) for all i. This is
a wrapper function for gsl_bspline_eval_nonzero() which
formats the output in a nice way.

Inputs: x - point for evaluation
        B - (output) where to store B_i(x) values
            the length of this vector is
            n = nbreak + k - 2 = l + k - 1 = w->n
        w - bspline workspace

Return: success or error

Notes: The w->knots vector must be initialized prior to calling
       this function (see gsl_bspline_knots())
*/

int
gsl_bspline_eval (const double x, gsl_vector * B, gsl_bspline_workspace * w)
{
  if (B->size != w->n)
    {
      GSL_ERROR ("vector B not of length n", GSL_EBADLEN);
    }
  else
    {
      size_t i;			/* looping */
      size_t istart;		/* first non-zero spline for x */
      size_t iend;		/* last non-zero spline for x, knot for x */
      int error;		/* error handling */

      /* find all non-zero B_i(x) values */
      error = gsl_bspline_eval_nonzero (x, w->B, &istart, &iend, w);
      if (error)
	{
	  return error;
	}

      /* store values in appropriate part of given vector */
      for (i = 0; i < istart; i++)
	gsl_vector_set (B, i, 0.0);

      for (i = istart; i <= iend; i++)
	gsl_vector_set (B, i, gsl_vector_get (w->B, i - istart));

      for (i = iend + 1; i < w->n; i++)
	gsl_vector_set (B, i, 0.0);

      return GSL_SUCCESS;
    }
}				/* gsl_bspline_eval() */

/*
gsl_bspline_eval_nonzero()
  Evaluate all non-zero B-spline functions at point x.
These are the B_i(x) for i in [istart, iend].
Always B_i(x) = 0 for i < istart and for i > iend.

Inputs: x      - point at which to evaluate splines
        Bk     - (output) where to store B-spline values (length k)
        istart - (output) B-spline function index of
                 first non-zero basis for given x
        iend   - (output) B-spline function index of
                 last non-zero basis for given x.
                 This is also the knot index corresponding to x.
        w      - bspline workspace

Return: success or error

Notes: 1) the w->knots vector must be initialized before calling
          this function

       2) On output, B contains

             [B_{istart,k}, B_{istart+1,k},
             ..., B_{iend-1,k}, B_{iend,k}]

          evaluated at the given x.
*/

int
gsl_bspline_eval_nonzero (const double x, gsl_vector * Bk, size_t * istart,
			  size_t * iend, gsl_bspline_workspace * w)
{
  if (Bk->size != w->k)
    {
      GSL_ERROR ("Bk vector length does not match order k", GSL_EBADLEN);
    }
  else
    {
      size_t i;			/* spline index */
      size_t j;			/* looping */
      int flag = 0;		/* interval search flag */
      int error = 0;		/* error flag */

      i = bspline_find_interval (x, &flag, w);
      error = bspline_process_interval_for_eval (x, &i, flag, w);
      if (error)
	{
	  return error;
	}

      *istart = i - w->k + 1;
      *iend = i;

      bspline_pppack_bsplvb (w->knots, w->k, 1, x, *iend, &j, w->deltal,
			     w->deltar, Bk);

      return GSL_SUCCESS;
    }
}				/* gsl_bspline_eval_nonzero() */

/*
gsl_bspline_deriv_eval()
  Evaluate d^j/dx^j B_i(x) for all i, 0 <= j <= nderiv.
This is a wrapper function for gsl_bspline_deriv_eval_nonzero()
which formats the output in a nice way.

Inputs: x      - point for evaluation
        nderiv - number of derivatives to compute, inclusive.
        dB     - (output) where to store d^j/dx^j B_i(x)
                 values. the size of this matrix is
                 (n = nbreak + k - 2 = l + k - 1 = w->n)
                 by (nderiv + 1)
        w      - bspline derivative workspace

Return: success or error

Notes: 1) The w->knots vector must be initialized prior to calling
          this function (see gsl_bspline_knots())

       2) based on PPPACK's bsplvd
*/

int
gsl_bspline_deriv_eval (const double x, const size_t nderiv, gsl_matrix * dB,
			gsl_bspline_workspace * w,
			gsl_bspline_deriv_workspace * dw)
{
  if (dB->size1 != w->n)
    {
      GSL_ERROR ("dB matrix first dimension not of length n", GSL_EBADLEN);
    }
  else if (dB->size2 < nderiv + 1)
    {
      GSL_ERROR
	("dB matrix second dimension must be at least length nderiv+1",
	 GSL_EBADLEN);
    }
  else if (dw->k < w->k)
    {
      GSL_ERROR ("derivative workspace is too small", GSL_EBADLEN);
    }
  else
    {
      size_t i;			/* looping */
      size_t j;			/* looping */
      size_t istart;		/* first non-zero spline for x */
      size_t iend;		/* last non-zero spline for x, knot for x */
      int error;		/* error handling */

      /* find all non-zero d^j/dx^j B_i(x) values */
      error =
	gsl_bspline_deriv_eval_nonzero (x, nderiv, dw->dB, &istart, &iend, w,
					dw);
      if (error)
	{
	  return error;
	}

      /* store values in appropriate part of given matrix */
      for (j = 0; j <= nderiv; j++)
	{
	  for (i = 0; i < istart; i++)
	    gsl_matrix_set (dB, i, j, 0.0);

	  for (i = istart; i <= iend; i++)
	    gsl_matrix_set (dB, i, j, gsl_matrix_get (dw->dB, i - istart, j));

	  for (i = iend + 1; i < w->n; i++)
	    gsl_matrix_set (dB, i, j, 0.0);
	}

      return GSL_SUCCESS;
    }
}				/* gsl_bspline_deriv_eval() */

/*
gsl_bspline_deriv_eval_nonzero()
  At point x evaluate all requested, non-zero B-spline function
derivatives and store them in dB.  These are the
d^j/dx^j B_i(x) with i in [istart, iend] and j in [0, nderiv].
Always d^j/dx^j B_i(x) = 0 for i < istart and for i > iend.

Inputs: x      - point at which to evaluate splines
        nderiv - number of derivatives to request, inclusive
        dB     - (output) where to store dB-spline derivatives
                 (size k by nderiv + 1)
        istart - (output) B-spline function index of
                 first non-zero basis for given x
        iend   - (output) B-spline function index of
                 last non-zero basis for given x.
                 This is also the knot index corresponding to x.
        w      - bspline derivative workspace

Return: success or error

Notes: 1) the w->knots vector must be initialized before calling
          this function

       2) On output, dB contains

            [[B_{istart,  k}, ..., d^nderiv/dx^nderiv B_{istart  ,k}],
             [B_{istart+1,k}, ..., d^nderiv/dx^nderiv B_{istart+1,k}],
             ...
             [B_{iend-1,  k}, ..., d^nderiv/dx^nderiv B_{iend-1,  k}],
             [B_{iend,    k}, ..., d^nderiv/dx^nderiv B_{iend,    k}]]

          evaluated at x.  B_{istart, k} is stored in dB(0,0).
          Each additional column contains an additional derivative.

       3) Note that the zero-th column of the result contains the
          0th derivative, which is simply a function evaluation.

       4) based on PPPACK's bsplvd
*/

int
gsl_bspline_deriv_eval_nonzero (const double x, const size_t nderiv,
				gsl_matrix * dB, size_t * istart,
				size_t * iend, gsl_bspline_workspace * w,
				gsl_bspline_deriv_workspace * dw)
{
  if (dB->size1 != w->k)
    {
      GSL_ERROR ("dB matrix first dimension not of length k", GSL_EBADLEN);
    }
  else if (dB->size2 < nderiv + 1)
    {
      GSL_ERROR
	("dB matrix second dimension must be at least length nderiv+1",
	 GSL_EBADLEN);
    }
  else if (dw->k < w->k)
    {
      GSL_ERROR ("derivative workspace is too small", GSL_EBADLEN);
    }
  else
    {
      size_t i;			/* spline index */
      size_t j;			/* looping */
      int flag = 0;		/* interval search flag */
      int error = 0;		/* error flag */
      size_t min_nderivk;

      i = bspline_find_interval (x, &flag, w);
      error = bspline_process_interval_for_eval (x, &i, flag, w);
      if (error)
	{
	  return error;
	}

      *istart = i - w->k + 1;
      *iend = i;

      bspline_pppack_bsplvd (w->knots, w->k, x, *iend,
			     w->deltal, w->deltar, dw->A, dB, nderiv);

      /* An order k b-spline has at most k-1 nonzero derivatives
         so we need to zero all requested higher order derivatives */
      min_nderivk = GSL_MIN_INT (nderiv, w->k - 1);
      for (j = min_nderivk + 1; j <= nderiv; j++)
	{
	  for (i = 0; i < w->k; i++)
	    {
	      gsl_matrix_set (dB, i, j, 0.0);
	    }
	}

      return GSL_SUCCESS;
    }
}				/* gsl_bspline_deriv_eval_nonzero() */

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

/*
bspline_find_interval()
  Find knot interval such that t_i <= x < t_{i + 1}
where the t_i are knot values.

Inputs: x    - x value
        flag - (output) error flag
        w    - bspline workspace

Return: i (index in w->knots corresponding to left limit of interval)

Notes: The error conditions are reported as follows:

       Condition                        Return value        Flag
       ---------                        ------------        ----
       x < t_0                               0               -1
       t_i <= x < t_{i+1}                    i                0
       t_i < x = t_{i+1} = t_{n+k-1}         i                0
       t_{n+k-1} < x                       l+k-1             +1
*/

static inline size_t
bspline_find_interval (const double x, int *flag, gsl_bspline_workspace * w)
{
  size_t i;

  if (x < gsl_vector_get (w->knots, 0))
    {
      *flag = -1;
      return 0;
    }

  /* find i such that t_i <= x < t_{i+1} */
  for (i = w->k - 1; i < w->k + w->l - 1; i++)
    {
      const double ti = gsl_vector_get (w->knots, i);
      const double tip1 = gsl_vector_get (w->knots, i + 1);

      if (tip1 < ti)
	{
	  GSL_ERROR ("knots vector is not increasing", GSL_EINVAL);
	}

      if (ti <= x && x < tip1)
	break;

      if (ti < x && x == tip1 && tip1 == gsl_vector_get (w->knots, w->k + w->l
							 - 1))
	break;
    }

  if (i == w->k + w->l - 1)
    *flag = 1;
  else
    *flag = 0;

  return i;
}				/* bspline_find_interval() */

/*
bspline_process_interval_for_eval()
  Consumes an x location, left knot from bspline_find_interval, flag
from bspline_find_interval, and a workspace.  Checks that x lies within
the splines' knots, enforces some endpoint continuity requirements, and
avoids divide by zero errors in the underlying bspline_pppack_* functions.
*/
static inline int
bspline_process_interval_for_eval (const double x, size_t * i, const int flag,
				   gsl_bspline_workspace * w)
{
  if (flag == -1)
    {
      GSL_ERROR ("x outside of knot interval", GSL_EINVAL);
    }
  else if (flag == 1)
    {
      if (x <= gsl_vector_get (w->knots, *i) + GSL_DBL_EPSILON)
	{
	  *i -= 1;
	}
      else
	{
	  GSL_ERROR ("x outside of knot interval", GSL_EINVAL);
	}
    }

  if (gsl_vector_get (w->knots, *i) == gsl_vector_get (w->knots, *i + 1))
    {
      GSL_ERROR ("knot(i) = knot(i+1) will result in division by zero",
		 GSL_EINVAL);
    }

  return GSL_SUCCESS;
}				/* bspline_process_interval_for_eval */

/********************************************************************
 * PPPACK ROUTINES
 *
 * The routines herein deliberately avoid using the bspline workspace,
 * choosing instead to pass all work areas explicitly.  This allows
 * others to more easily adapt these routines to low memory or
 * parallel scenarios.
 ********************************************************************/

/*
bspline_pppack_bsplvb()
  calculates the value of all possibly nonzero b-splines at x of order
jout = max( jhigh , (j+1)*(index-1) ) with knot sequence t.

Parameters:
   t      - knot sequence, of length left + jout , assumed to be
            nondecreasing.  assumption t(left).lt.t(left + 1).
            division by zero  will result if t(left) = t(left+1)
   jhigh  -
   index  - integers which determine the order jout = max(jhigh,
            (j+1)*(index-1))  of the b-splines whose values at x
            are to be returned.  index  is used to avoid
            recalculations when several columns of the triangular
            array of b-spline values are needed (e.g., in  bsplpp
            or in  bsplvd ).  precisely,

            if  index = 1 ,
               the calculation starts from scratch and the entire
               triangular array of b-spline values of orders
               1,2,...,jhigh  is generated order by order , i.e.,
               column by column .

            if  index = 2 ,
               only the b-spline values of order j+1, j+2, ..., jout
               are generated, the assumption being that biatx, j,
               deltal, deltar are, on entry, as they were on exit
               at the previous call.

            in particular, if jhigh = 0, then jout = j+1, i.e.,
            just the next column of b-spline values is generated.
   x      - the point at which the b-splines are to be evaluated.
   left   - an integer chosen (usually) so that
            t(left) .le. x .le. t(left+1).
   j      - (output) a working scalar for indexing
   deltal - (output) a working area which must be of length at least jout
   deltar - (output) a working area which must be of length at least jout
   biatx  - (output) array of length jout, with  biatx(i)
            containing the value at  x  of the polynomial of order
            jout which agrees with the b-spline b(left-jout+i,jout,t)
            on the interval (t(left), t(left+1)) .

Method:
   the recurrence relation

                      x - t(i)              t(i+j+1) - x
      b(i,j+1)(x) = -----------b(i,j)(x) + ---------------b(i+1,j)(x)
                    t(i+j)-t(i)            t(i+j+1)-t(i+1)

   is used (repeatedly) to generate the (j+1)-vector  b(left-j,j+1)(x),
   ...,b(left,j+1)(x)  from the j-vector  b(left-j+1,j)(x),...,
   b(left,j)(x), storing the new values in  biatx  over the old. the
   facts that

      b(i,1) = 1  if  t(i) .le. x .lt. t(i+1)

   and that

      b(i,j)(x) = 0  unless  t(i) .le. x .lt. t(i+j)

   are used. the particular organization of the calculations follows
   algorithm (8) in chapter x of [1].

Notes:

   (1) This is a direct translation of PPPACK's bsplvb routine with
       j, deltal, deltar rewritten as input parameters and
       utilizing zero-based indexing.

   (2) This routine contains no error checking.  Please use routines
       like gsl_bspline_eval().
*/

static void
bspline_pppack_bsplvb (const gsl_vector * t,
		       const size_t jhigh,
		       const size_t index,
		       const double x,
		       const size_t left,
		       size_t * j,
		       gsl_vector * deltal,
		       gsl_vector * deltar, gsl_vector * biatx)
{
  size_t i;			/* looping */
  double saved;
  double term;

  if (index == 1)
    {
      *j = 0;
      gsl_vector_set (biatx, 0, 1.0);
    }

  for ( /* NOP */ ; *j < jhigh - 1; *j += 1)
    {
      gsl_vector_set (deltar, *j, gsl_vector_get (t, left + *j + 1) - x);
      gsl_vector_set (deltal, *j, x - gsl_vector_get (t, left - *j));

      saved = 0.0;

      for (i = 0; i <= *j; i++)
	{
	  term = gsl_vector_get (biatx, i) / (gsl_vector_get (deltar, i)
					      + gsl_vector_get (deltal,
								*j - i));

	  gsl_vector_set (biatx, i,
			  saved + gsl_vector_get (deltar, i) * term);

	  saved = gsl_vector_get (deltal, *j - i) * term;
	}

      gsl_vector_set (biatx, *j + 1, saved);
    }

  return;
}				/* gsl_bspline_pppack_bsplvb */

/*
bspline_pppack_bsplvd()
  calculates value and derivs of all b-splines which do not vanish at x

Parameters:
   t      - the knot array, of length left+k (at least)
   k      - the order of the b-splines to be evaluated
   x      - the point at which these values are sought
   left   - an integer indicating the left endpoint of the interval
            of interest. the k b-splines whose support contains the
            interval (t(left), t(left+1)) are to be considered.
            it is assumed that t(left) .lt. t(left+1)
            division by zero will result otherwise (in  bsplvb).
            also, the output is as advertised only if
            t(left) .le. x .le. t(left+1) .
   deltal - a working area which must be of length at least k
   deltar - a working area which must be of length at least k
   a      - an array of order (k,k), to contain b-coeffs of the
            derivatives of a certain order of the k b-splines
            of interest.
   dbiatx - an array of order (k,nderiv). its entry (i,m) contains
            value of (m)th derivative of (left-k+i)-th b-spline
            of order k for knot sequence  t, i=1,...,k, m=0,...,nderiv.
   nderiv - an integer indicating that values of b-splines and
            their derivatives up to AND INCLUDING the nderiv-th
            are asked for. (nderiv is replaced internally by the
            integer mhigh in (1,k) closest to it.)

Method:
   values at x of all the relevant b-splines of order k,k-1,..., k+1-nderiv
   are generated via bsplvb and stored temporarily in dbiatx.  then, the
   b-coeffs of the required derivatives of the b-splines of interest are
   generated by differencing, each from the preceeding one of lower order,
   and combined with the values of b-splines of corresponding order in
   dbiatx  to produce the desired values .

Notes:

   (1) This is a direct translation of PPPACK's bsplvd routine with
       deltal, deltar rewritten as input parameters (to later feed them
       to bspline_pppack_bsplvb) and utilizing zero-based indexing.

   (2) This routine contains no error checking.
*/

static void
bspline_pppack_bsplvd (const gsl_vector * t,
		       const size_t k,
		       const double x,
		       const size_t left,
		       gsl_vector * deltal,
		       gsl_vector * deltar,
		       gsl_matrix * a,
		       gsl_matrix * dbiatx, const size_t nderiv)
{
  int i, ideriv, il, j, jlow, jp1mid, kmm, ldummy, m, mhigh;
  double factor, fkmm, sum;

  size_t bsplvb_j;
  gsl_vector_view dbcol = gsl_matrix_column (dbiatx, 0);

  mhigh = GSL_MIN_INT (nderiv, k - 1);
  bspline_pppack_bsplvb (t, k - mhigh, 1, x, left, &bsplvb_j, deltal, deltar,
			 &dbcol.vector);
  if (mhigh > 0)
    {
      /* the first column of dbiatx always contains the b-spline
         values for the current order. these are stored in column
         k-current order before bsplvb is called to put values
         for the next higher order on top of it.  */
      ideriv = mhigh;
      for (m = 1; m <= mhigh; m++)
	{
	  for (j = ideriv, jp1mid = 0; j < (int) k; j++, jp1mid++)
	    {
	      gsl_matrix_set (dbiatx, j, ideriv,
			      gsl_matrix_get (dbiatx, jp1mid, 0));
	    }
	  ideriv--;
	  bspline_pppack_bsplvb (t, k - ideriv, 2, x, left, &bsplvb_j, deltal,
				 deltar, &dbcol.vector);
	}

      /* at this point,  b(left-k+i, k+1-j)(x) is in  dbiatx(i,j)
         for i=j,...,k-1 and j=0,...,mhigh. in particular, the
         first column of dbiatx is already in final form. to obtain
         corresponding derivatives of b-splines in subsequent columns,
         generate their b-repr. by differencing, then evaluate at x. */
      jlow = 0;
      for (i = 0; i < (int) k; i++)
	{
	  for (j = jlow; j < (int) k; j++)
	    {
	      gsl_matrix_set (a, j, i, 0.0);
	    }
	  jlow = i;
	  gsl_matrix_set (a, i, i, 1.0);
	}

      /* at this point, a(.,j) contains the b-coeffs for the j-th of the
         k b-splines of interest here. */
      for (m = 1; m <= mhigh; m++)
	{
	  kmm = k - m;
	  fkmm = (float) kmm;
	  il = left;
	  i = k - 1;

	  /* for j=1,...,k, construct b-coeffs of (m)th  derivative
	     of b-splines from those for preceding derivative by
	     differencing and store again in  a(.,j) . the fact that
	     a(i,j) = 0  for i .lt. j  is used. */
	  for (ldummy = 0; ldummy < kmm; ldummy++)
	    {
	      factor =
		fkmm / (gsl_vector_get (t, il + kmm) -
			gsl_vector_get (t, il));
	      /* the assumption that t(left).lt.t(left+1) makes
	         denominator in factor nonzero. */
	      for (j = 0; j <= i; j++)
		{
		  gsl_matrix_set (a, i, j, factor * (gsl_matrix_get (a, i, j)
						     - gsl_matrix_get (a,
								       i - 1,
								       j)));
		}
	      il--;
	      i--;
	    }

	  /* for i=1,...,k, combine b-coeffs a(.,i) with b-spline values
	     stored in dbiatx(.,m) to get value of (m)th  derivative
	     of i-th b-spline (of interest here) at x, and store in
	     dbiatx(i,m). storage of this value over the value of a
	     b-spline of order m there is safe since the remaining
	     b-spline derivatives of the same order do not use this
	     value due to the fact that a(j,i) = 0 for j .lt. i . */
	  for (i = 0; i < (int) k; i++)
	    {
	      sum = 0;
	      jlow = GSL_MAX_INT (i, m);
	      for (j = jlow; j < (int) k; j++)
		{
		  sum +=
		    gsl_matrix_get (a, j, i) * gsl_matrix_get (dbiatx, j, m);
		}
	      gsl_matrix_set (dbiatx, i, m, sum);
	    }
	}
    }

  return;

}				/* bspline_pppack_bsplvd */

max_stars_repo_name: ruslankuzmin/julia

max_issues_repo_name: ruslankuzmin/julia

max_forks_repo_name: ruslankuzmin/julia


text:
// The various memcached commands that we support.

#include "commands.h"
#include "connections.h"
#include "locking.h"
#include "protocol.h"
#include "server.h"
#include "threads.h"
#include "utils.h"

#include <assert.h>
#include <gsl/gsl_randist.h>
#include <stdio.h>
#include <stdlib.h>

static const char* default_key = "skeleton";

/* // process a memcached get(s) command. (we don't support CAS). This function */
/* // performs the request parsing and setup of backend RPC's. */
/* void process_get_command(conn *c, token_t *tokens, size_t ntokens, */
/*                                 bool return_cas) { */
/* 	char *key; */
/* 	size_t nkey; */
/* 	int i = 0; */
/* 	item *it; */
/* 	token_t *key_token = &tokens[KEY_TOKEN]; */
/* 	char *suffix; */
/* 	worker_thread_t *t = c->thread; */
/* 	memcached_t *mc; */
/*  */
/* 	assert(c != NULL); */
/*  */
/* 	key  = key_token->value; */
/* 	nkey = key_token->length; */
/*  */
/* 	if (config.use_dist) { */
/* 		long size = config.dist_arg1 + gsl_ran_gaussian(config.r, config.dist_arg2); */
/* 		if (config.verbose > 1) { */
/* 			fprintf(stderr, "allocated blob: %ld\n", size); */
/* 		} */
/* 		c->mem_blob = malloc(sizeof(char) * size); */
/* 	} */
/*  */
/* 	if(nkey > KEY_MAX_LENGTH) { */
/* 		error_response(c, "CLIENT_ERROR bad command line format"); */
/* 		return; */
/* 	} */
/*  */
/* 	// lookup key-value. */
/* 	it = item_get(key, nkey); */
/*  */
/* 	// hit. */
/* 	if (it) { */
/* 		if (i >= c->isize && !conn_expand_items(c)) { */
/* 			item_remove(it); */
/* 			error_response(c, "SERVER_ERROR out of memory writing get response"); */
/* 			return; */
/* 		} */
/* 		// add item to remembered list (i.e., we've taken ownership of them */
/* 		// through refcounting and later must release them once we've */
/* 		// written out the iov associated with them). */
/* 		item_update(it); */
/* 		*(c->ilist + i) = it; */
/* 		i++; */
/* 	} */
/*  */
/* 	// make sure it's a single get */
/* 	key_token++; */
/* 	if (key_token->length != 0 || key_token->value != NULL) { */
/* 		error_response(c, "SERVER_ERROR only support single `get`"); */
/* 		return; */
/* 	} */
/*  */
/* 	// update our rememberd reference set. */
/* 	c->icurr = c->ilist; */
/* 	c->ileft = i; */
/*  */
/* 	// setup RPC calls. */
/* 	for (i = 0; i < t->memcache_used; i++) { */
/* 		mc = t->memcache[i]; */
/* 		if (!conn_add_msghdr(mc) != 0) { */
/* 			error_response(mc, "SERVER_ERROR out of memory preparing response"); */
/* 			return; */
/* 		} */
/* 		memcache_get(mc, c, default_key); */
/* 	} */
/* 	conn_set_state(c, conn_rpc_wait); */
/* } */

// complete the response to a get request.
void finish_get_command(conn *c) {
	item *it;
	int i;

	// setup all items for writing out.
	for (i = 0; i < c->ileft; i++) {
		it = *(c->ilist + i);
		if (it) {
			// Construct the response. Each hit adds three elements to the
			// outgoing data list:
			//   "VALUE <key> <flags> <data_length>\r\n"
			//   "<data>\r\n"
			// The <data> element is stored on the connection item list, not on
			// the iov list.
			if (!conn_add_iov(c, "VALUE ", 6) ||
				 !conn_add_iov(c, ITEM_key(it), it->nkey) ||
				 !conn_add_iov(c, ITEM_suffix(it), it->nsuffix + it->nbytes)) {
				item_remove(it);
				error_response(c, "SERVER_ERROR out of memory writing get response");
				return;
			}

			if (config.verbose > 1) {
				fprintf(stderr, ">%d sending key %s\n", c->sfd, ITEM_key(it));
			}
		} else {
			fprintf(stderr, "ERROR corrupted ilist!\n");
			exit(1);
		}
	}

	if (config.verbose > 1) {
		fprintf(stderr, ">%d END\n", c->sfd);
	}

	if (!conn_add_iov(c, "END\r\n", 5) != 0) {
		error_response(c, "SERVER_ERROR out of memory writing get response");
	} else {
		conn_set_state(c, conn_mwrite);
	}
}

// process a memcached get(s) command. (we don't support CAS).
void process_get_command(conn *c, token_t *tokens, size_t ntokens,
                                bool return_cas) {
	char *key;
	size_t nkey;
	int i = 0;
	item *it;
	token_t *key_token = &tokens[KEY_TOKEN];
	char *suffix;

	assert(c != NULL);

	if (config.alloc && c->mem_blob == NULL) {
		long size = config.alloc_mean + gsl_ran_gaussian(c->thread->r, config.alloc_stddev);
		size = size <= 0 ? 10 : size;
		if (config.verbose > 0) {
			fprintf(stderr, "allocated blob: %ld\n", size);
		}

		c->mem_blob = malloc(sizeof(char) * size);
		c->mem_free_delay = 0;

		if (config.rtt_delay) {
			double r = config.rtt_mean + gsl_ran_gaussian(c->thread->r, config.rtt_stddev);
			if (r >= config.rtt_cutoff) {
				int wait = r / 100;
				if (config.verbose > 0) {
					fprintf(stderr, "delay: %d\n", wait);
				}
				c->mem_free_delay = wait;
				conn_set_state(c, conn_mwrite);
			}
		}
	}

	// process the whole command line, (only part of it may be tokenized right now)
	do {
		// process all tokenized keys at this stage.
		while(key_token->length != 0) {
			key = key_token->value;
			nkey = key_token->length;

			if(nkey > KEY_MAX_LENGTH) {
				error_response(c, "CLIENT_ERROR bad command line format");
				return;
			}

			// lookup key-value.
			it = item_get(key, nkey);

			// hit.
			if (it) {
				if (i >= c->isize && !conn_expand_items(c)) {
					item_remove(it);
					break;
				}

				// Construct the response. Each hit adds three elements to the
				// outgoing data list:
				//   "VALUE <key> <flags> <data_length>\r\n"
				//   "<data>\r\n"
				// The <data> element is stored on the connection item list, not on
				// the iov list.
				if (!conn_add_iov(c, "VALUE ", 6) != 0 ||
				    !conn_add_iov(c, ITEM_key(it), it->nkey) != 0 ||
				    !conn_add_iov(c, ITEM_suffix(it), it->nsuffix + it->nbytes) != 0) {
					item_remove(it);
					break;
				}

				if (config.verbose > 1) {
					fprintf(stderr, ">%d sending key %s\n", c->sfd, key);
				}

				// add item to remembered list (i.e., we've taken ownership of them
				// through refcounting and later must release them once we've
				// written out the iov associated with them).
				item_update(it);
				*(c->ilist + i) = it;
				i++;
			}

			key_token++;
		}

		/*
		 * If the command string hasn't been fully processed, get the next set
		 * of tokens.
		 */
		if(key_token->value != NULL) {
			ntokens = tokenize_command(key_token->value, tokens, MAX_TOKENS);
			key_token = tokens;
		}

	} while(key_token->value != NULL);

	c->icurr = c->ilist;
	c->ileft = i;

	if (config.verbose > 1) {
		fprintf(stderr, ">%d END\n", c->sfd);
	}

	// If the loop was terminated because of out-of-memory, it is not reliable
	// to add END\r\n to the buffer, because it might not end in \r\n. So we
	// send SERVER_ERROR instead.
	if (key_token->value != NULL || !conn_add_iov(c, "END\r\n", 5) != 0) {
		error_response(c, "SERVER_ERROR out of memory writing get response");
	} else {
		conn_set_state(c, conn_mwrite);
	}
}

// process a memcached set command.
void process_update_command(conn *c, token_t *tokens,
                            const size_t ntokens,
                            int comm, bool handle_cas) {
	int vlen;
	assert(c != NULL);

	if (tokens[KEY_TOKEN].length > KEY_MAX_LENGTH ||
	    !safe_strtol(tokens[4].value, (int32_t *)&vlen)) {
		error_response(c, "CLIENT_ERROR bad command line format");
		return;
	}

	if (vlen < 0) {
		error_response(c, "CLIENT_ERROR bad command line format");
		return;
	}

	// setup value to be read
	c->sbytes = vlen + 2; // for \r\n consumption.
	conn_set_state(c, conn_read_value);
}

// process a memcached stat command.
void process_stat_command(conn *c, token_t *tokens, const size_t ntokens) {
	mutex_lock(&c->stats->lock);

	// just for debugging right now
	fprintf(stderr, "STAT client_id %d\n", c->stats->client_id);
	fprintf(stderr, "STAT total_connections %d\n", c->stats->total_connections);
	fprintf(stderr, "STAT live_connections %d\n", c->stats->live_connections);
	fprintf(stderr, "STAT requests %d\n", c->stats->requests);

	mutex_unlock(&c->stats->lock);
	conn_set_state(c, conn_new_cmd);
}

max_stars_repo_name: dterei/synthetic-client

max_issues_repo_name: dterei/synthetic-client

max_forks_repo_name: dterei/synthetic-client


text:
/****************************************************************************
 *                                                                          *
 *  Author : lukasz.iwaszkiewicz@gmail.com                                  *
 *  ~~~~~~~~                                                                *
 *  License : see COPYING file for details.                                 *
 *  ~~~~~~~~~                                                               *
 ****************************************************************************/

#pragma once
#include <chrono>
#include <cstdint>
#include <gsl/gsl>
#include <iostream>
#include <vector>

namespace tp {
class CanFrame;

/**
 * @brief The TimeProvider struct
 * TODO if I'm not mistaken, those classes ought to have 100µs resolutiuon instead of 1ms (1000µs).
 */
struct ChronoTimeProvider {
        long operator() () const
        {
                using namespace std::chrono;
                system_clock::time_point now = system_clock::now ();
                auto duration = now.time_since_epoch ();
                return duration_cast<milliseconds> (duration).count ();
        }
};

struct CoutPrinter {
        template <typename T> void operator() (T &&a) { std::cout << std::forward<T> (a) << std::endl; }
};

/**
 * Buffer for all input and ouptut operations
 */
using IsoMessage = std::vector<uint8_t>;

} // namespace tp

max_stars_repo_name: semasquare/cpp-can-isotp

max_issues_repo_name: semasquare/cpp-can-isotp

max_forks_repo_name: semasquare/cpp-can-isotp


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

#ifndef __GSL_MATRIX_DOUBLE_H__
#define __GSL_MATRIX_DOUBLE_H__

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/block/gsl_check_range.h>
#include <gsl/vector/gsl_vector_double.h>

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

__BEGIN_DECLS

typedef struct
{
  size_t size1;		//nr rows
  size_t size2;		//nr cols
  size_t tda;		//space between two consecutive rows
  double * data;	//a pointer to the data
  gsl_block * block;//the block pointer.
  int owner;		//this is 1 if the current object owns the memory block, 0 otherwise
} gsl_matrix;

typedef struct
{
  gsl_matrix matrix;
} _gsl_matrix_view;

typedef _gsl_matrix_view gsl_matrix_view;

typedef struct
{
  gsl_matrix matrix;
} _gsl_matrix_const_view;

typedef const _gsl_matrix_const_view gsl_matrix_const_view;

/* Allocation */

GSL_FUNC gsl_matrix *
gsl_matrix_alloc (const size_t n1, const size_t n2);

GSL_FUNC gsl_matrix *
gsl_matrix_calloc (const size_t n1, const size_t n2);

GSL_FUNC gsl_matrix *
gsl_matrix_alloc_from_block (gsl_block * b,
                                   const size_t offset,
                                   const size_t n1,
                                   const size_t n2,
                                   const size_t d2);

GSL_FUNC gsl_matrix *
gsl_matrix_alloc_from_matrix (gsl_matrix * m,
                                    const size_t k1,
                                    const size_t k2,
                                    const size_t n1,
                                    const size_t n2);

GSL_FUNC gsl_vector *
gsl_vector_alloc_row_from_matrix (gsl_matrix * m,
                                        const size_t i);

GSL_FUNC gsl_vector *
gsl_vector_alloc_col_from_matrix (gsl_matrix * m,
                                        const size_t j);

GSL_FUNC void gsl_matrix_free (gsl_matrix * m);

/* Views */

GSL_FUNC _gsl_matrix_view
gsl_matrix_submatrix (gsl_matrix * m,
                            const size_t i, const size_t j,
                            const size_t n1, const size_t n2);

GSL_FUNC _gsl_vector_view
gsl_matrix_row (gsl_matrix * m, const size_t i);

GSL_FUNC _gsl_vector_view
gsl_matrix_column (gsl_matrix * m, const size_t j);

GSL_FUNC _gsl_vector_view
gsl_matrix_diagonal (gsl_matrix * m);

GSL_FUNC _gsl_vector_view
gsl_matrix_subdiagonal (gsl_matrix * m, const size_t k);

GSL_FUNC _gsl_vector_view
gsl_matrix_superdiagonal (gsl_matrix * m, const size_t k);

GSL_FUNC _gsl_matrix_view
gsl_matrix_view_array (double * base,
                             const size_t n1,
                             const size_t n2);

GSL_FUNC _gsl_matrix_view
gsl_matrix_view_array_with_tda (double * base,
                                      const size_t n1,
                                      const size_t n2,
                                      const size_t tda);

GSL_FUNC _gsl_matrix_view
gsl_matrix_view_vector (gsl_vector * v,
                              const size_t n1,
                              const size_t n2);

GSL_FUNC _gsl_matrix_view
gsl_matrix_view_vector_with_tda (gsl_vector * v,
                                       const size_t n1,
                                       const size_t n2,
                                       const size_t tda);

GSL_FUNC _gsl_matrix_const_view
gsl_matrix_const_submatrix (const gsl_matrix * m,
                                  const size_t i, const size_t j,
                                  const size_t n1, const size_t n2);

GSL_FUNC _gsl_vector_const_view
gsl_matrix_const_row (const gsl_matrix * m,
                            const size_t i);

GSL_FUNC _gsl_vector_const_view
gsl_matrix_const_column (const gsl_matrix * m,
                               const size_t j);

GSL_FUNC _gsl_vector_const_view
gsl_matrix_const_diagonal (const gsl_matrix * m);

GSL_FUNC _gsl_vector_const_view
gsl_matrix_const_subdiagonal (const gsl_matrix * m,
                                    const size_t k);

GSL_FUNC _gsl_vector_const_view
gsl_matrix_const_superdiagonal (const gsl_matrix * m,
                                      const size_t k);

GSL_FUNC _gsl_matrix_const_view
gsl_matrix_const_view_array (const double * base,
                                   const size_t n1,
                                   const size_t n2);

GSL_FUNC _gsl_matrix_const_view
gsl_matrix_const_view_array_with_tda (const double * base,
                                            const size_t n1,
                                            const size_t n2,
                                            const size_t tda);

GSL_FUNC _gsl_matrix_const_view
gsl_matrix_const_view_vector (const gsl_vector * v,
                                    const size_t n1,
                                    const size_t n2);

GSL_FUNC _gsl_matrix_const_view
gsl_matrix_const_view_vector_with_tda (const gsl_vector * v,
                                             const size_t n1,
                                             const size_t n2,
                                             const size_t tda);

/* Operations */

GSL_FUNC double   gsl_matrix_get(const gsl_matrix * m, const size_t i, const size_t j);
GSL_FUNC void    gsl_matrix_set(gsl_matrix * m, const size_t i, const size_t j, const double x);

GSL_FUNC double * gsl_matrix_ptr(gsl_matrix * m, const size_t i, const size_t j);
GSL_FUNC const double * gsl_matrix_const_ptr(const gsl_matrix * m, const size_t i, const size_t j);

GSL_FUNC void gsl_matrix_set_zero (gsl_matrix * m);
GSL_FUNC void gsl_matrix_set_identity (gsl_matrix * m);
GSL_FUNC void gsl_matrix_set_all (gsl_matrix * m, double x);

GSL_FUNC int gsl_matrix_fread (FILE * stream, gsl_matrix * m) ;
GSL_FUNC int gsl_matrix_fwrite (FILE * stream, const gsl_matrix * m) ;
GSL_FUNC int gsl_matrix_fscanf (FILE * stream, gsl_matrix * m);
GSL_FUNC int gsl_matrix_fprintf (FILE * stream, const gsl_matrix * m, const char * format);

GSL_FUNC int gsl_matrix_memcpy(gsl_matrix * dest, const gsl_matrix * src);
GSL_FUNC int gsl_matrix_swap(gsl_matrix * m1, gsl_matrix * m2);

GSL_FUNC int gsl_matrix_swap_rows(gsl_matrix * m, const size_t i, const size_t j);
GSL_FUNC int gsl_matrix_swap_columns(gsl_matrix * m, const size_t i, const size_t j);
GSL_FUNC int gsl_matrix_swap_rowcol(gsl_matrix * m, const size_t i, const size_t j);
GSL_FUNC int gsl_matrix_transpose (gsl_matrix * m);
GSL_FUNC int gsl_matrix_transpose_memcpy (gsl_matrix * dest, const gsl_matrix * src);

GSL_FUNC double gsl_matrix_max (const gsl_matrix * m);
GSL_FUNC double gsl_matrix_min (const gsl_matrix * m);
GSL_FUNC void gsl_matrix_minmax (const gsl_matrix * m, double * min_out, double * max_out);

GSL_FUNC void gsl_matrix_max_index (const gsl_matrix * m, size_t * imax, size_t *jmax);
GSL_FUNC void gsl_matrix_min_index (const gsl_matrix * m, size_t * imin, size_t *jmin);
GSL_FUNC void gsl_matrix_minmax_index (const gsl_matrix * m, size_t * imin, size_t * jmin, size_t * imax, size_t * jmax);

GSL_FUNC int gsl_matrix_isnull (const gsl_matrix * m);
GSL_FUNC int gsl_matrix_ispos (const gsl_matrix * m);
GSL_FUNC int gsl_matrix_isneg (const gsl_matrix * m);

GSL_FUNC int gsl_matrix_add (gsl_matrix * a, const gsl_matrix * b);
GSL_FUNC int gsl_matrix_add_scaled (gsl_matrix * a, const gsl_matrix * b, double scaleA, double scaleB);
GSL_FUNC int gsl_matrix_sub (gsl_matrix * a, const gsl_matrix * b);
GSL_FUNC int gsl_matrix_mul_elements (gsl_matrix * a, const gsl_matrix * b);
GSL_FUNC int gsl_matrix_div_elements (gsl_matrix * a, const gsl_matrix * b);
GSL_FUNC int gsl_matrix_scale (gsl_matrix * a, const double x);
GSL_FUNC int gsl_matrix_add_constant (gsl_matrix * a, const double x);
GSL_FUNC int gsl_matrix_add_diagonal (gsl_matrix * a, const double x);

/***********************************************************************/
/* The functions below are obsolete                                    */
/***********************************************************************/
GSL_FUNC int gsl_matrix_get_row(gsl_vector * v, const gsl_matrix * m, const size_t i);
GSL_FUNC int gsl_matrix_get_col(gsl_vector * v, const gsl_matrix * m, const size_t j);
GSL_FUNC int gsl_matrix_set_row(gsl_matrix * m, const size_t i, const gsl_vector * v);
GSL_FUNC int gsl_matrix_set_col(gsl_matrix * m, const size_t j, const gsl_vector * v);

/* inline functions if you are using GCC */

#ifdef HAVE_INLINE
extern inline
double
gsl_matrix_get(const gsl_matrix * 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_set(gsl_matrix * m, const size_t i, const size_t j, const double 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
double *
gsl_matrix_ptr(gsl_matrix * 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 (double *) (m->data + (i * m->tda + j)) ;
}

extern inline
const double *
gsl_matrix_const_ptr(const gsl_matrix * 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 double *) (m->data + (i * m->tda + j)) ;
}

#endif

__END_DECLS

#endif /* __GSL_MATRIX_DOUBLE_H__ */

max_stars_repo_name: MontyThibault/centre-of-mass-awareness

max_issues_repo_name: MontyThibault/centre-of-mass-awareness

max_forks_repo_name: MontyThibault/centre-of-mass-awareness