text:
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <gsl/gsl_linalg.h>

int DEBUG = 1;
int DEBUG_SOLUTION = 1;
int STORE_TIMING = 0;

FILE* decompFile;
FILE* solveFile;

double* generateVector(int size)
{
    double* vec = calloc(sizeof(double), size);
    for (int i = 0; i < size; ++i){
        double val = (double)rand() / RAND_MAX; // double [0, 1]
        vec[i] = val;
    }
    return vec;
}

void printMatrix(double* m, int size)
{
    for (int i = 0; i < size; ++i){
        for (int j = 0; j < size; ++j){
            printf("%g    ", m[j + i * size]);
        }
        printf("\n");
    }
    printf("\n");
}

void printVector(double* v, int size)
{
    for (int i = 0; i < size; ++i){
        printf("%g", v[i]);
        printf("\n");
    }
    printf("\n");
}

double* copyVector(double* src, int size)
{
    double* vec = calloc(sizeof(double), size);
    for (int i = 0; i < size; ++i){
        vec[i] = src[i];
    }
    return vec;
}

void checkSolution(double* m, double* b, gsl_vector x, int size)
{
    const double epsilon = 0.0001;
    int failures = 0;

    for (int i = 0; i < size; ++i){
        double sum = 0;
        for (int j = 0; j < size; ++j){
            sum += m[j + i * size] * x.data[j];
            // printf("%g * %g = %g\n", m[j + i * size], x.data[j], m[j + i * size] * x.data[j]);
        }
        if (abs(sum - b[i]) > epsilon){
            failures++;
            if (DEBUG_SOLUTION){
                printf("Wrong solution was calculated: %g --- %g\n", sum, b[i]);
                printf("\n");
            }
        }
    }
    if (DEBUG_SOLUTION){
        if (failures == 0){
            printf("All calculated values were correct");
        } else {
            printf("Total %n errors in calculated value", failures);
        }
        printf("\n");
    }
}

void decomp(gsl_matrix *A, gsl_permutation* p, int *signum, int size)
{
    clock_t start = clock();
    gsl_linalg_LU_decomp(A, p, signum);
    clock_t end = clock();

    double t = ((double) (end - start)) / CLOCKS_PER_SEC;
    if (STORE_TIMING) fprintf (decompFile,"%d %g\n", size, t);
}

void solve(const gsl_matrix* LU, const gsl_permutation* p, const gsl_vector* b, gsl_vector* x, int size)
{
    clock_t start = clock();
    gsl_linalg_LU_solve(LU, p, b, x);
    clock_t end = clock();

    double t = ((double) (end - start)) / CLOCKS_PER_SEC;
    if (STORE_TIMING) fprintf (solveFile,"%d %g\n", size, t);
}

void calculate(int n)
{
    double* a_data = generateVector(n * n);
    double* b_data = generateVector(n);

    double* originalMatrix = copyVector(a_data, n * n);

    gsl_matrix_view m = gsl_matrix_view_array(a_data, n, n);
    gsl_vector_view b = gsl_vector_view_array(b_data, n);

    if (DEBUG){
        printf("Matrix M:\n");
        printMatrix(a_data, n);

        printf("Vector B:\n");
        printVector(b_data, n);
    }
    gsl_vector *x = gsl_vector_alloc(n);
    int s;

    gsl_permutation *p = gsl_permutation_alloc(n);

    decomp(&m.matrix, p, &s, n);
    solve(&m.matrix, p, &b.vector, x, n);

    if (DEBUG){
        printf("Solution X:\n");
        gsl_vector_fprintf(stdout, x, "%g");
        printf("\n");
    }
    checkSolution(originalMatrix, b_data, *x, n);

    gsl_permutation_free(p);
    gsl_vector_free(x);
}

int main(int argc, char* argv[])
{
    if (argc != 2){
        return -1;
    }
    char *pEnd;
    int n = strtol(argv[1], &pEnd, 0);

    srand(time(NULL));
    decompFile = fopen("decomp.txt", "w");
    solveFile = fopen("solve.txt", "w");

    if (n == 0){
        return -1;
    } else if (n < 0){
        DEBUG = 0;
        DEBUG_SOLUTION = 0;
        STORE_TIMING = 1;
        for (int i = 10; i <= 1000; i+=10){
            calculate(i);
        }
    } else {
        calculate(n);
    }

    return 0;
}


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

/* L D L^T decomposition of a symmetric banded positive semi-definite matrix */

#include <config.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>

static double symband_norm1(const gsl_matrix * A);
static int ldlt_band_Ainv(CBLAS_TRANSPOSE_t TransA, gsl_vector * x, void * params);

/*
gsl_linalg_ldlt_band_decomp()
  L D L^T decomposition of a square symmetric positive semi-definite banded
matrix

Inputs: A - matrix in symmetric banded format, N-by-ndiag where N is the size of
            the matrix and ndiag is the number of nonzero diagonals.

Notes:
1) The matrix D is stored in the first column of A;
the first subdiagonal of L in the second column and so on.

2) If ndiag > 1, the 1-norm of A is stored in A(N,ndiag) on output

3) At each diagonal element, the matrix is factored as

A(j:end,j:end) = [ A11 A21^T ] = [ 1 0 ] [ alpha 0 ] [ 1 v^T ]
                 [ A21 A22   ]   [ v L ] [   0   D ] [ 0 L^T ]

where:

alpha = A(j,j)
v = A(j+1:end, j) / alpha
A22 = L D L^T + alpha v v^T

So we start at A(1,1) and work right. Pseudo-code is:

loop j = 1, ..., N
  alpha = A(j,j)
  A(j+1:end, j) := A(j+1:end, j) / alpha     (DSCAL)
  A(j+1:end, j+1:end) -= alpha v v^T         (DSYR)

Due to the banded structure, v has at most p non-zero elements, where
p is the lower bandwidth
*/

int
gsl_linalg_ldlt_band_decomp(gsl_matrix * A)
{
  const size_t N = A->size1;     /* size of matrix */
  const size_t ndiag = A->size2; /* number of diagonals in band, including main diagonal */

  if (ndiag > N)
    {
      GSL_ERROR ("invalid matrix dimensions", GSL_EBADLEN);
    }
  else
    {
      const size_t p = ndiag - 1; /* lower bandwidth */
      const int kld = (int) GSL_MAX(1, p);
      double Anorm;
      size_t j;

      /* check for quick return */
      if (ndiag == 1)
        return GSL_SUCCESS;

      /*
       * calculate 1-norm of A and store in lower right of matrix, which is not accessed
       * by rest of routine. gsl_linalg_ldlt_band_rcond() will use this later. If
       * A is diagonal, there is no empty slot to store the 1-norm, so the rcond routine
       * will have to compute it.
       */
      Anorm = symband_norm1(A);
      gsl_matrix_set(A, N - 1, p, Anorm);

      for (j = 0; j < N - 1; ++j)
        {
          double ajj = gsl_matrix_get(A, j, 0);
          size_t lenv;

          if (ajj == 0.0)
            {
              GSL_ERROR("matrix is singular", GSL_EDOM);
            }

          /* number of elements in v, which will normally be p, unless we
           * are in lower right corner of matrix */
          lenv = GSL_MIN(p, N - j - 1);

          if (lenv > 0)
            {
              gsl_vector_view v = gsl_matrix_subrow(A, j, 1, lenv);
              gsl_matrix_view m = gsl_matrix_submatrix(A, j + 1, 0, lenv, lenv);

              gsl_blas_dscal(1.0 / ajj, &v.vector);

              m.matrix.tda = kld;
              gsl_blas_dsyr(CblasUpper, -ajj, &v.vector, &m.matrix);
            }
        }

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_ldlt_band_solve (const gsl_matrix * LDLT,
                            const gsl_vector * b,
                            gsl_vector * x)
{
  if (LDLT->size1 != b->size)
    {
      GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
    }
  else if (LDLT->size1 != x->size)
    {
      GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
    }
  else
    {
      int status;

      /* copy x <- b */
      gsl_vector_memcpy (x, b);

      status = gsl_linalg_ldlt_band_svx(LDLT, x);

      return status;
    }
}

int
gsl_linalg_ldlt_band_svx (const gsl_matrix * LDLT, gsl_vector * x)
{
  if (LDLT->size1 != x->size)
    {
      GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
    }
  else
    {
      gsl_vector_const_view diag = gsl_matrix_const_column(LDLT, 0);

      /* solve for z using forward-substitution, L z = b */
      cblas_dtbsv(CblasColMajor, CblasLower, CblasNoTrans, CblasUnit,
                  (int) LDLT->size1, (int) (LDLT->size2 - 1), LDLT->data, LDLT->tda,
                  x->data, x->stride);

      /* solve for y, D y = z */
      gsl_vector_div(x, &diag.vector);

      /* perform back-substitution, L^T x = y */
      cblas_dtbsv(CblasColMajor, CblasLower, CblasTrans, CblasUnit,
                  (int) LDLT->size1, (int) (LDLT->size2 - 1), LDLT->data, LDLT->tda,
                  x->data, x->stride);

      return GSL_SUCCESS;
    }
}

/*
gsl_linalg_ldlt_band_unpack()
  Unpack symmetric banded format matrix LDLT into
larger matrix L and diagonal vector D
*/

int
gsl_linalg_ldlt_band_unpack (const gsl_matrix * LDLT, gsl_matrix * L, gsl_vector * D)
{
  const size_t N = LDLT->size1;

  if (N != L->size1)
    {
      GSL_ERROR("L matrix does not match LDLT dimensions", GSL_EBADLEN);
    }
  else if (L->size1 != L->size2)
    {
      GSL_ERROR("L matrix is not square", GSL_ENOTSQR);
    }
  else if (N != D->size)
    {
      GSL_ERROR("D vector does not match LDLT dimensions", GSL_EBADLEN);
    }
  else
    {
      const size_t p = LDLT->size2 - 1; /* lower bandwidth */
      gsl_vector_const_view diag = gsl_matrix_const_column(LDLT, 0);
      gsl_vector_view diagL = gsl_matrix_diagonal(L);
      size_t i;

      /* copy diagonal entries */
      gsl_vector_memcpy(D, &diag.vector);

      /* copy subdiagonals into L */
      for (i = 1; i <= p; ++i)
        {
          gsl_vector_const_view v = gsl_matrix_const_subcolumn(LDLT, i, 0, N - i);
          gsl_vector_view w = gsl_matrix_subdiagonal(L, i);
          gsl_vector_memcpy(&w.vector, &v.vector);
        }

      /* set main diagonal of L */
      gsl_vector_set_all(&diagL.vector, 1.0);

      /* zero out remaining subdiagonals */
      for (i = p + 1; i < N; ++i)
        {
          gsl_vector_view w = gsl_matrix_subdiagonal(L, i);
          gsl_vector_set_zero(&w.vector);
        }

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_ldlt_band_rcond (const gsl_matrix * LDLT, double * rcond, gsl_vector * work)
{
  const size_t N = LDLT->size1;

  if (work->size != 3 * N)
    {
      GSL_ERROR ("work vector must have length 3*N", GSL_EBADLEN);
    }
  else
    {
      int status;
      const size_t ndiag = LDLT->size2;
      double Anorm;    /* ||A||_1 */
      double Ainvnorm; /* ||A^{-1}||_1 */

      if (ndiag == 1)
        {
          /* diagonal matrix, compute 1-norm since it has not been stored */
          Anorm = symband_norm1(LDLT);
        }
      else
        {
          /* 1-norm is stored in A(N, ndiag) by gsl_linalg_ldlt_band_decomp() */
          Anorm = gsl_matrix_get(LDLT, N - 1, ndiag - 1);
        }

      *rcond = 0.0;

      /* return if matrix is singular */
      if (Anorm == 0.0)
        return GSL_SUCCESS;

      status = gsl_linalg_invnorm1(N, ldlt_band_Ainv, (void *) LDLT, &Ainvnorm, work);
      if (status)
        return status;

      if (Ainvnorm != 0.0)
        *rcond = (1.0 / Anorm) / Ainvnorm;

      return GSL_SUCCESS;
    }
}

/* compute 1-norm of symmetric banded matrix */
static double
symband_norm1(const gsl_matrix * A)
{
  const size_t N = A->size1;
  const size_t ndiag = A->size2; /* number of diagonals in band, including main diagonal */
  double value;

  if (ndiag == 1)
    {
      /* diagonal matrix */
      gsl_vector_const_view v = gsl_matrix_const_column(A, 0);
      CBLAS_INDEX_t idx = gsl_blas_idamax(&v.vector);
      value = gsl_vector_get(&v.vector, idx);
    }
  else
    {
      size_t j;

      value = 0.0;
      for (j = 0; j < N; ++j)
        {
          size_t ncol = GSL_MIN(ndiag, N - j); /* number of elements in column j below and including main diagonal */
          gsl_vector_const_view v = gsl_matrix_const_subrow(A, j, 0, ncol);
          double sum = gsl_blas_dasum(&v.vector);
          size_t k, l;

          /* sum now contains the absolute sum of elements below and including main diagonal for column j; we
           * have to add the symmetric elements above the diagonal */
          k = j;
          l = 1;
          while (k > 0 && l < ndiag)
            {
              double Akl = gsl_matrix_get(A, --k, l++);
              sum += fabs(Akl);
            }

          value = GSL_MAX(value, sum);
        }
    }

  return value;
}

/* x := A^{-1} x = A^{-t} x, A = L D L^T */
static int
ldlt_band_Ainv(CBLAS_TRANSPOSE_t TransA, gsl_vector * x, void * params)
{
  gsl_matrix * LDLT = (gsl_matrix * ) params;
  gsl_vector_const_view diag = gsl_matrix_const_column(LDLT, 0);

  (void) TransA; /* unused parameter warning */

  /* compute x := L^{-1} x */
  cblas_dtbsv(CblasColMajor, CblasLower, CblasNoTrans, CblasUnit,
              (int) LDLT->size1, (int) (LDLT->size2 - 1), LDLT->data, LDLT->tda,
              x->data, x->stride);

  /* compute x := D^{-1} x */
  gsl_vector_div(x, &diag.vector);

  /* compute x := L^{-T} x */
  cblas_dtbsv(CblasColMajor, CblasLower, CblasTrans, CblasUnit,
              (int) LDLT->size1, (int) (LDLT->size2 - 1), LDLT->data, LDLT->tda,
              x->data, x->stride);

  return GSL_SUCCESS;
}

ext: c

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
#ifndef ALM_L_KATYUSHA_H
#define ALM_L_KATYUSHA_H

#include "L_Katyusha.h"

#include <string>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <stdio.h>      /* printf */
#include <time.h>
#include <fstream>
#include <algorithm>
#include <iomanip>
#include <ctime>
#include <sstream>
//#include "cmd_line.h"

//This class implements the method IPALM_L_Katyusha

/*
The optimization problem to solve is:

min f_0(x)+\sum_{i=1}^m h_i(f_i(x))+P(x)

x is a d-dimensional vector.

// Each subproblem solves problem of the form
f_0(x)+ \sum_{i=1}^m h_{\beta_s}^i(f_i(x);\lambda_s^i) +g(x)+ \beta_s/2\|x-x_s\|^2
 by L_Katyuhsa
*/

template<typename L, typename D>
class ALM_L_Katyusha: public L_Katyusha<L, D>
{
private:

protected:

  std::vector<D> x_s;

  std::vector<D> old_x_s;

  std::vector<D> lambda_s;

  std::vector<D> old_lambda_s;

  std::vector<D> x_tmp;

  std::vector<D> g_tmp;

  std::vector<D> g_tmp_x;
  std::vector<D> g_tmp_w;
  std::vector<D> valuefi;

  std::vector<D> is_feasiblity_constraint;

  D beta_s;

  D epsilon_s;

  D m_s;

  D m_0;

  L m;

  L d;

  D eta;

  D rho;

  D tau_s;

  std::vector<D> L_phi;

  D function_value;

  D infeasibility;

  L print_every_N_ALM;

  D running_time_ALM;

  L nb_outer_iters;

  ofstream samp_ALM;

public:

  D mu_g;

  D lambda1;

  D lambda2;

  virtual inline void set_Li_Lf(){}

  virtual inline D value_of_P(vector<D> &){return D(NULL);}
  virtual inline void prox_of_P(D, vector<D> &, vector<D> &){} //prox_of_P(L,x,y) computes y=argmin{P(u)+L/2||u-x||^2}
  virtual inline D prox_of_h_star_j(D,D, L){return D(NULL);} //prox_of_h_star_j(x,beta,j) computes argmin{h*(y)+beta/2(y-x)^2}
  virtual inline D value_of_h_star_j(D, L){return D(NULL);}
  virtual inline void compute_gradient_f0(vector<D> &,vector<D> &){} //compute_gradient_f0(x,g) computes g=nabla f0(x)
  virtual inline D compute_gradient_and_value_f_i(vector<D> &,vector<D> &, L){return D(NULL);}//compute_gradient_and_value_f_i(i,x,g) computes g=nabla f_i(x) and return f_i(x)

  virtual inline D value_of_h_j(D, L){return D(NULL);}
  virtual inline D distance_to_domain_of_h_j(D,L){return D(NULL);}

  virtual inline D value_of_f0(vector<D> &){return D(NULL);} //f0(x)
  virtual inline D value_of_f_i(vector<D> &,L){return D(NULL);} // f_i(x)
  virtual inline void set_dimension(){} // set the value of m and d
  virtual inline void set_is_feasibility_constraint(){}

  ALM_L_Katyusha()
  : L_Katyusha<L,D>()
  {

  }

  D value_of_phi_i(L i){
    D res= 0;
    if (i==0)
    res=value_of_f0(this->x);
    else
    {
      D fi=value_of_f_i(this->x,i);
      i=i-1;
      D lambdasi=lambda_s[i];
      D tmp= prox_of_h_star_j(fi/beta_s+ lambdasi,beta_s,i);
      res=fi*tmp- value_of_h_star_j(tmp,i)- beta_s/2*(tmp- lambdasi)*(tmp- lambdasi);
    }
    return this->nsamples*res;
  }

  void compute_full_gradient(vector<D> & x, vector<D> & gx){
    compute_gradient_f0(x,gx);
    for(L i=0;i<m;i++){
          D fi=compute_gradient_and_value_f_i(x,g_tmp,i+1);
          D tmp= prox_of_h_star_j(fi/beta_s+ lambda_s[i],beta_s,i);
          for(L j=0;j<d;j++)
            gx[j]+=g_tmp[j]*tmp;
    }
  }

  void compute_batch_delta_gradient(){
    for (L j=0; j< this->batch_size; j++){
      L s= this->batch_i[j];
      if(s==0){
        compute_gradient_f0(this->x,g_tmp_x);
        compute_gradient_f0(this->w,g_tmp_w);
        add_up(1.,1.,this->theta_S[s]);
      }
      else{
        L i=s-1;
        D fi_x=compute_gradient_and_value_f_i(this->x,g_tmp_x,s);
        D tmp_x= prox_of_h_star_j(fi_x/beta_s+ lambda_s[i],beta_s,i);
        D fi_w=compute_gradient_and_value_f_i(this->w,g_tmp_w,s);
        D tmp_w= prox_of_h_star_j(fi_w/beta_s+ lambda_s[i],beta_s,i);
        add_up(tmp_x,tmp_w,this->theta_S[s]);
      }
    }
  }

  void  add_up(D fx, D fw, D c){
      for(L i=0;i<d;i++)
        this->batch_delta_gradient[i]+=(fx*g_tmp_x[i]-fw*g_tmp_w[i])*c;
    }

  inline D value_of_g(){
    D res= value_of_P(this->x);
    for(L i=0;i<d;i++)
      res+= 0.5*beta_s*(this->x[i]- x_s[i])*(this->x[i]- x_s[i]);
    return res;
  }

  void prox_of_g(D stepsz, vector<D> & x, vector<D> & g, vector<D> & nextx){
    D tmp0=beta_s+stepsz;
    for(L i=0;i<d;i++)
     x_tmp[i]=(beta_s*x_s[i]+stepsz*x[i]-g[i])/tmp0;
    prox_of_P(tmp0, x_tmp, nextx);
  }

  void compute_x(){
    for (L i=0; i< d; i++){
      x_s[i]= this->x[i];
    }
    for (L i= 0; i< m; i++){
      valuefi[i]=value_of_f_i(x_s,i+1);
    }
  }

  void update_lambda(){
    for(L j=0;j<m;j++)
    {
      D fi= valuefi[j];
      lambda_s[j]= prox_of_h_star_j(fi/beta_s+ lambda_s[j],beta_s,j);
    }
  }

  void compute_function_value(){
    function_value=value_of_f0(x_s);
    D res= 0;
    for (L i= 0; i< m; i++){
      D fi=value_of_f_i(x_s,i+1);
      if(is_feasiblity_constraint[i]==1){
        D tmp=distance_to_domain_of_h_j(fi,i);
        res+=tmp*tmp;
      }
      else
       function_value+=value_of_h_j(fi,i);
    }
    function_value+=value_of_P(x_s);
    infeasibility=sqrt(res);
  }

inline void compute_m0(D beta0, D val_eta, D val_rho, L val_tau){
  cout<<"beta0="<<beta0<<"val_eta="<<val_eta<<"; val_rho="<<val_rho<<"; val_tau="<<val_tau<<endl;
  m_s= this->nsamples/val_tau*1e+6;

  //D tmp7=1-val_tau/this->n*sqrt(beta0/(max_Lf_s+max_M_s/beta0+beta0));
  //cout<<"tmp7: "<<tmp7<<"; "<<val_tau/(this->n+0.)*sqrt(beta0/(max_Lf_s+max_M_s/beta0+beta0))<<endl;
  //m_0=(2*log(val_rho)+log(val_eta)+log(2))/log(tmp7);

  cout<<" m_s="<<m_s<<endl;
}

void Initialize(D beta_0, D epsilon_0,  D val_eta, D val_rho,L val_tau, vector<D> & x0,vector<D> & y0){
  cout<<"start initializing ALM"<<endl;
  set_dimension();
  this->nsamples=m+1;
  this->nfeatures=d;
  this->tau=val_tau;

  beta_s=beta_0;
  tau_s= 1;
  epsilon_s=epsilon_0;
  compute_m0(beta_s,val_eta,val_rho,val_tau);

  eta=val_eta;
  rho=val_rho;

  x_s.resize(d,0);
  old_x_s.resize(d,0);
  lambda_s.resize(m,0);
  old_lambda_s.resize(m,0);
  for(L i=0;i<d;i++){
    x_s[i]=x0[i];
    old_x_s[i]= x0[i];
  }
  for(L j=0;j<m;j++){
    lambda_s[j]=y0[j];
    old_lambda_s[j]= y0[j];
  }
  g_tmp.clear();
  g_tmp.resize(d,0);
  g_tmp_x.clear();
  g_tmp_x.resize(d,0);
  g_tmp_w.clear();
  g_tmp_w.resize(d,0);
  valuefi.clear();
  valuefi.resize(m,0);
  x_tmp.clear();
  x_tmp.resize(d,0);
  is_feasiblity_constraint.resize(m,0);

  set_is_feasibility_constraint();
  cout<<"Initialization ALM finished!"<<endl;

}

void reset_everything(){
  epsilon_s*=rho;
  beta_s*=eta;
}

void update_m_s(L val_tau){
  D tmp= 0;
  D tmpx= 0;
  D tmpy= 0;
  D tmpy2= 0;
  for(L j=0;j<m;j++)
  {
    D tmp_p= prox_of_h_star_j(valuefi[j]/beta_s/eta+ lambda_s[j], beta_s/eta, j);
    tmpy+=(tmp_p- lambda_s[j])*(tmp_p- lambda_s[j]);
    tmp+= (lambda_s[j]- old_lambda_s[j])*(lambda_s[j]- old_lambda_s[j]);
    tmpy2+= (beta_s*eta*lambda_s[j]- beta_s*old_lambda_s[j])*(beta_s*eta*lambda_s[j]- beta_s*old_lambda_s[j]);
  }
  for (L i= 0; i< d;i++){
    tmpx= (x_s[i]- old_x_s[i])*(x_s[i]- old_x_s[i]);
  }
  D tmp4= 2*epsilon_s+ beta_s*tmp;
  D tmp5= (1- eta)*beta_s/2*tmpy;
  D tmp6= 0;
  D tmp7= sqrt(tmp)*sqrt(tmpy2);
  D tmp1= log(tmp4+ tmp5+ tmp6+ tmp7+ 2/(2*eta- 1)*tau_s*beta_s/2*tmpx);
  D tmp2= log(epsilon_s)+ log(rho)- log(2);
  D tmp3= sqrt(beta_s/this->sumLi);
  m_s= ceil((tmp1- tmp2)/log(2)*4*max(1.0+this->nsamples,1/tmp3)/val_tau);
  cout<<"here here="<<tmp4<<" tmp5="<<tmp5<<" tmp6="<<tmp6<<" tmp7="<<tmp7<<" tmp1="<<tmp1<<" tmp2="<<tmp2<<" tmp3="<<tmp3<<endl;
}

inline void compute_and_record_res(){
  if(nb_outer_iters%print_every_N_ALM==0){
    compute_function_value();
    cout<<setprecision(9)<<"Iteration: "<<nb_outer_iters<<"; time="<<running_time_ALM<< "; function value="<<function_value<<"; infeasibility= "<< infeasibility<<endl;
    samp_ALM<<setprecision(9)<<nb_outer_iters<<" "<<running_time_ALM<<" "<< function_value<<" "<< infeasibility<<endl;
  }
}

void ALM_solve_with_L_Katyusha(D beta_0, D epsilon_0,  D eta, D rho,vector<D> & x0,vector<D> & y0, L val_tau, L max_nb_outer, L p_N_1, L p_N_2, string filename1, string filename2, D time){
  Initialize(beta_0, epsilon_0, eta, rho,val_tau, x0, y0);
  nb_outer_iters=0;
  filename1= "results/ALM_"+filename1;
  samp_ALM.open(filename1.c_str());
  running_time_ALM=0;
  print_every_N_ALM=p_N_1;
  this->set_print_every_N(p_N_2);
  compute_and_record_res();

 while(nb_outer_iters<max_nb_outer){
  D start= std::clock();
  cout<<"m_s= "<< ceil(m_s/this->nsamples*val_tau)<<"; beta_s="<<beta_s<<"; epsilon_s="<<epsilon_s<<endl;
  this->loopless2(x_s, filename2, mu_g+ tau_s*beta_s,  ceil(m_s/this->nsamples*val_tau), epsilon_s, val_tau, 1, 1, 0, 1);
  for(L i=0;i<d;i++){
    old_x_s[i]=x_s[i];
  }
  for(L i=0;i<m;i++){
    old_lambda_s[i]=lambda_s[i];
  }
  compute_x();
  update_lambda();
  update_m_s(val_tau);
  nb_outer_iters++;
  running_time_ALM+=( std::clock() - start ) / (double) CLOCKS_PER_SEC;
  compute_and_record_res();
  start = std::clock();
  reset_everything();
  running_time_ALM+=( std::clock() - start ) / (double) CLOCKS_PER_SEC;
    if (running_time_ALM> time){
      break;
    }
  }
}

};

#endif /* MIN_SMOOTH_CONVEX_H */

ext: h


text:
/* -*- c -*- */
/*
 * Copyright 2002-2018 Free Software Foundation, Inc.
 *
 * This file is part of GNU Radio
 *
 * GNU Radio 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, or (at your option)
 * any later version.
 *
 * GNU Radio 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 GNU Radio; see the file COPYING.  If not, write to
 * the Free Software Foundation, Inc., 51 Franklin Street,
 * Boston, MA 02110-1301, USA.
 */

/*
 * generate MMSE FIR interpolation table values
 */

#include <math.h>
#include <assert.h>
#include <gsl/gsl_integration.h>

#define	MU	0.5			/* the MU for which we're computing coeffs */

#define	Ts	(1.0)			/* sampling period 		*/
#define	B	(1.0/(4*Ts))		/* one-sided signal bandwidth 	*/
//#define	B	(1.0/(8./3*Ts))		/* one-sided signal bandwidth 	*/

static unsigned global_n;
static double  *global_h;
double		global_mu = MU;
double		global_B  = B;

gsl_integration_workspace *global_gsl_int_workspace = NULL;

/*
 * This function computes the difference squared between the ideal
 * interpolator frequency response at frequency OMEGA and the
 * approximation defined by the FIR coefficients in global_h[]
 *
 * See eqn (9-7), "Digital Communication Receivers", Meyr, Moeneclaey
 * and Fechtel, Wiley, 1998.
 */

static double
integrand (double omega, void *params)
{
  double real_ideal;
  double real_approx;
  double real_diff;
  double imag_ideal;
  double imag_approx;
  double imag_diff;

  int	 i, n;
  int	 I1;
  double *h;

  real_ideal = cos (omega * Ts * global_mu);
  imag_ideal = sin (omega * Ts * global_mu);

  n = global_n;
  h = global_h;
  I1 = -(n / 2);

  real_approx = 0;
  imag_approx = 0;

  for (i = 0; i < n; i++){
    real_approx += h[i] * cos (-omega * Ts * (i + I1));
    imag_approx += h[i] * sin (-omega * Ts * (i + I1));
  }

  real_diff = real_ideal - real_approx;
  imag_diff = imag_ideal - imag_approx;

  return real_diff * real_diff + imag_diff * imag_diff;
}

/*
 * Integrate the difference squared over all frequencies of interest.
 */
double
c_fcn (double *x, int n)
{
  gsl_function F;
  double result, error;

  F.function = integrand;
  F.params = NULL;

  assert ((n & 1) == 0);	/* assert n is even */
  global_n = n;
  global_h = x;
  gsl_integration_qag(&F, -2 * M_PI * global_B, 2 * M_PI * global_B,
                      0.0, 1e-12, 1000, GSL_INTEG_GAUSS61,
                      global_gsl_int_workspace, &result, &error);
  return result;
}

/* this is the interface expected by the calling fortran code */

double
objective (double x[], int *ndim)
{
  return c_fcn (x, *ndim);
}

static double
si (double x)
{
  if (fabs (x) < 1e-9)
    return 1.0;

  return sin(x) / x;
}

/*
 * starting guess for optimization
 */
void initpt (double x[], int ndim)
{
  int i;
  for (i = 0; i < ndim; i++){
    x[i] = si (M_PI * ((double) (i - ndim/2) + global_mu));
  }
}


text:
/* $Id$         */
/*--------------------------------------------------------------------*/
/*;  Copyright (C) 2003-2013                                          */
/*;  Associated Universities, Inc. Washington DC, USA.                */
/*;                                                                   */
/*;  This program is free software; you can redistribute it and/or    */
/*;  modify it under the terms of the GNU General Public License as   */
/*;  published by the Free Software Foundation; either version 2 of   */
/*;  the License, or (at your option) any later version.              */
/*;                                                                   */
/*;  This program is distributed in the hope that it will be useful,  */
/*;  but WITHOUT ANY WARRANTY; without even the implied warranty of   */
/*;  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the    */
/*;  GNU General Public License for more details.                     */
/*;                                                                   */
/*;  You should have received a copy of the GNU General Public        */
/*;  License along with this program; if not, write to the Free       */
/*;  Software Foundation, Inc., 675 Massachusetts Ave, Cambridge,     */
/*;  MA 02139, USA.                                                   */
/*;                                                                   */
/*;Correspondence about this software should be addressed as follows: */
/*;         Internet email: bcotton@nrao.edu.                         */
/*;         Postal address: William Cotton                            */
/*;                         National Radio Astronomy Observatory      */
/*;                         520 Edgemont Road                         */
/*;                         Charlottesville, VA 22903-2475 USA        */
/*--------------------------------------------------------------------*/
#ifndef OBITFFT_H
#define OBITFFT_H

/* FFTW3 ?  */
#if HAVE_FFTW3==1
#include <fftw3.h>
#elif HAVE_FFTW==1  /* else Try FFTW 2 */
#include <fftw.h>
#include <rfftw.h>
#elif HAVE_GSL==1  /* Else try GSL version */
#include <gsl/gsl_fft_complex_float.h>
#include <gsl/gsl_fft_halfcomplex_float.h>
#include <gsl/gsl_fft_real_float.h>
#include <gsl/gsl_fft_complex.h>
#include <gsl/gsl_fft_halfcomplex.h>
#include <gsl/gsl_fft_real.h>
#endif /* HAVE_GSL */
#include "Obit.h"
#include "ObitThread.h"
#include "ObitFArray.h"
#include "ObitCArray.h"

/*-------- Obit: Merx mollis mortibus nuper ------------------*/
/**
 * \file ObitFFT.h
 * ObitFFT Fast Fourier Transform class definition.
 *
 * This class is derived from the #Obit class.
 *
 * This class is for performing FFT on memory resident data.
 * This implementation uses the FFTW package if available, else gsl
 *
 * \section FFTOrder  Data order
 * Data passed to and from ObitFFT routines are as an ObitFArray or ObitCArray
 * which are stored in column major (Fortran) order.
 * Data are passed and returned in "center-at-the-edge" (i.e. unnatural) order
 * and there is NO transpose of the array axes.
 * In the half complex form, the first axis is nx/2+1 where nx is the number of
 * real elements.
 * Only even numbers of elements on each axis will work well.
 *
 * \section ObitFFTaccess Creators and Destructors
 * An ObitFFT can be created using newObitFFT which allows specifying
 * a name for the object, and the type, size and direction of the transform.
 *
 * A copy of a pointer to an ObitFFT should always be made using the
 * #ObitFFTRef function which updates the reference count in the object.
 * Then whenever freeing an ObitFFT or changing a pointer, the function
 * #ObitFFTUnref will decrement the reference count and destroy the object
 * when the reference count hits 0.
 * There is no explicit destructor.
 *
 * \section Multi-threaded operation
 * If this module is compiled with the -DHAVE_FFTW3THREADS==1 option
 * and linked with the -lfftw3f_threads option, multithreaded operation
 * is enabled.
 * NOTE: this should NOT be used inside ObitThreads as the threading will
 * collide.
 * The maximum number of threads to be used in subsequent ObitFFT objects is
 * set using #ObitFFTNThreads; setting nThreads to 1 turns off threading.
 * After all ObitFFT objects are destroyed, calling #ObitFFTClearThreads
 * will clean up after the threading.
 */

/*-------------- enumerations -------------------------------------*/
/**
 * \enum ObitFFTdir
 * enum for FFT direction
 * This specifies the status of the connection a disk resident data.
 */
enum obitFFTdir {
  /** Sign of exponent in transform = -1, real to complex */
  OBIT_FFT_Forward,
  /** Sign of exponent in transform = +1, complex to real  */
  OBIT_FFT_Reverse
}; /* end enum obitFFTdir */
/** typedef for enum for ObitIO object status. */
typedef enum obitFFTdir ObitFFTdir;

/**
 * \enum ObitFFTtype
 * enum for FFT type, complex to complex or real to complex/
 * This specifies the status of the connection a disk resident data.
 */
enum obitFFTtype {
  /** Full complex to complex transforms */
  OBIT_FFT_FullComplex,
  /** Real to half complex or reverse.  */
  OBIT_FFT_HalfComplex
}; /* end enum obitFFTtype */
/** typedef for enum for ObitIO object status. */
typedef enum obitFFTtype ObitFFTtype;

/*--------------Class definitions-------------------------------------*/
/** ObitFFT Class structure. */
typedef struct {
#include "ObitFFTDef.h"   /* this class definition */
} ObitFFT;

/*----------------- Macroes ---------------------------*/
/**
 * Macro to unreference (and possibly destroy) an ObitFFT
 * returns a ObitFFT*.
 * in = object to unreference
 */
#define ObitFFTUnref(in) ObitUnref (in)

/**
 * Macro to reference (update reference count) an ObitFFT.
 * returns a ObitFFT*.
 * in = object to reference
 */
#define ObitFFTRef(in) ObitRef (in)

/**
 * Macro to determine if an object is the member of this or a
 * derived class.
 * Returns TRUE if a member, else FALSE
 * in = object to reference
 */
#define ObitFFTIsA(in) ObitIsA (in, ObitFFTGetClass())

/*---------------Public functions---------------------------*/
/** Public: Class initializer. */
void ObitFFTClassInit (void);

/** Public: Constructor. */
ObitFFT* newObitFFT (gchar* name, ObitFFTdir dir, ObitFFTtype type,
		     olong rank, olong *dim);
/** Typedef for definition of class pointer structure */
typedef ObitFFT* (*newObitFFTFP) (gchar *name, ObitFFTdir dir,
				  ObitFFTtype type,
				  olong rank, olong *dim);

/** Public: ClassInfo pointer */
gconstpointer ObitFFTGetClass (void);

/** Public: Suggest efficient size for a transform */
olong ObitFFTSuggestSize (olong length);

/** Public: Enable threading */
void ObitFFTNThreads (olong nThreads);

/** Public: Cleanup threading */
void ObitFFTClearThreads (void);

/** Public: Real to half Complex. */
void ObitFFTR2C (ObitFFT *in, ObitFArray *inArray, ObitCArray *outArray);
typedef void (*ObitFFTR2CFP) (ObitFFT *in, ObitFArray *inArray,
			    ObitCArray *outArray);

/** Public: Half Complex to Real. */
void ObitFFTC2R (ObitFFT *in, ObitCArray *inArray, ObitFArray *outArray);
typedef void (*ObitFFTC2RFP) (ObitFFT *in, ObitCArray *inArray,
			    ObitFArray *outArray);

/** Public: Full Complex to Complex. */
void ObitFFTC2C (ObitFFT *in, ObitCArray *inArray, ObitCArray *outArray);
typedef void (*ObitFFTC2CFP) (ObitFFT *in, ObitCArray *inArray,
			    ObitCArray *outArray);

/*----------- ClassInfo Structure -----------------------------------*/
/**
 * ClassInfo Structure.
 * Contains class name, a pointer to any parent class
 * (NULL if none) and function pointers.
 */
typedef struct  {
#include "ObitFFTClassDef.h"
} ObitFFTClassInfo;

#endif /* OBITFFT_H */

max_stars_repo_name: sarrvesh/Obit

max_issues_repo_name: sarrvesh/Obit

max_forks_repo_name: sarrvesh/Obit


text:
//--------------------------------------------------------------------------------------------------
//
//	WIN-BLUETOOTH
//
//--------------------------------------------------------------------------------------------------
//
// The MIT License (MIT)
//
// Permission is hereby granted, free of charge, to any person obtaining a copy of this software
// and associated documentation files (the "Software"), to deal in the Software without
// restriction, including without limitation the rights to use, copy, modify, merge, publish,
// distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in all copies or
// substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING
// BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
// DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
//
//--------------------------------------------------------------------------------------------------
//
// Copyright (c) 2018 Nic Holthaus
//
//--------------------------------------------------------------------------------------------------
//
// ATTRIBUTION:
// https://github.com/martinmoene/gsl-lite/releases
//
//--------------------------------------------------------------------------------------------------
//
/// @file	obexResponse.h
/// @brief	OBEX server response classes
//
//--------------------------------------------------------------------------------------------------

#pragma once
#ifndef obexResponse_h__
#define obexResponse_h__

//-------------------------
//	INCLUDES
//-------------------------

#include <QObject>
#include <gsl-lite.h>
#include <obexHeader.h>

//-------------------------
//	FORWARD DECLARATIONS
//-------------------------

class QDataStream;

//--------------------------------------------------------------------------------------------------
//	OBEXResponse
//--------------------------------------------------------------------------------------------------
class OBEXResponse
{
	Q_GADGET

public:

	friend QDataStream& operator>>(QDataStream& in, OBEXResponse& response);

	// OBEX operation response codes. See OBEX standard chapter 3.2.1
	class Code
	{
	public:

		enum class Enum : quint8
		{
			INVALID					= 0x00,
			CONTINUE				= 0x90,
			SUCCESS					= 0xA0,
			SERVICEUNAVAILABLE		= 0xD3,
		};

		Code(quint8 value);
		operator quint8() const;
		bool isFinal() const;
		bool operator==(const Code& other) const;
		bool operator!=(const Code& other) const;

	private:

		Code::Enum	m_code;
		bool		m_isFinal = false;
	};

public:

	OBEXResponse() = default;
	virtual ~OBEXResponse() = default;
	virtual gsl::string_span data() = 0;	// returns a span representing the response-specific data storage
	virtual quint16 packetLength() = 0;		// packet length value from the response data
	virtual bool validateAndFixup() = 0;	// called at the end of operator>>. Used to validate code fields, and convert
											// large types from big to little endian. Returns true on success.
	bool isValid() const;

protected:

	bool m_valid = false;
	std::vector<OBEXHeader>	m_optionalHeaders;
};

QDataStream& operator>>(QDataStream& in, OBEXResponse& response);

//--------------------------------------------------------------------------------------------------
//	OBEXConnectResponse
//--------------------------------------------------------------------------------------------------

class OBEXConnectResponse : public OBEXResponse
{
public:

	OBEXConnectResponse() = default;
	virtual ~OBEXConnectResponse() = default;
	virtual gsl::string_span data() override;
	virtual quint16 packetLength() override;	// the return value will generally not be valid until after data has been streamed into the class.
	virtual bool validateAndFixup() override;
	quint16 maxPacketLength() const;

private:

#pragma pack(push, 1)
	struct Data
	{
		OBEXResponse::Code::Enum    code			= Code::Enum::INVALID;
		quint16						length			= 0;
		quint8						version			= 0;
		quint8						flags			= 0;
		quint16						maxPacketLength = 0;
	} m_data;
#pragma pack(pop)
};

//--------------------------------------------------------------------------------------------------
//	OBEXDisconnectResponse
//--------------------------------------------------------------------------------------------------

class OBEXDisconnectResponse : public OBEXResponse
{
public:

	OBEXDisconnectResponse() = default;
	virtual ~OBEXDisconnectResponse() = default;

	virtual gsl::string_span data() override;
	virtual quint16 packetLength() override;
	virtual bool validateAndFixup() override;

#pragma pack(push, 1)
	struct Data
	{
		OBEXResponse::Code::Enum    code = Code::Enum::SUCCESS;
		quint16						length = 0;
	} m_data;
#pragma pack(pop)
};

//--------------------------------------------------------------------------------------------------
//	OBEXPutResponse
//--------------------------------------------------------------------------------------------------

class OBEXPutResponse : public OBEXResponse
{
public:

	OBEXPutResponse() = default;
	virtual ~OBEXPutResponse() = default;

	virtual gsl::string_span data() override;
	virtual quint16 packetLength() override;
	virtual bool validateAndFixup() override;
	bool continueSending() const;

protected:

#pragma pack(push, 1)
	struct Data
	{
		OBEXResponse::Code::Enum    code = Code::Enum::CONTINUE;	// assume success until proven otherwise
		quint16						length = 0;
	} m_data;
#pragma pack(pop)

};
#endif // obexResponse_h__

max_stars_repo_name: nholthaus/win-bluetooth

max_issues_repo_name: nholthaus/win-bluetooth

max_forks_repo_name: nholthaus/win-bluetooth


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

#include "demo_fn.h"
#include "demo_fn.c"

int
main (void)
{
  int status;
  int iter = 0, max_iter = 100;
  const gsl_root_fdfsolver_type *T;
  gsl_root_fdfsolver *s;
  double x0, x = 5.0, r_expected = sqrt (5.0);
  gsl_function_fdf FDF;
  struct quadratic_params params = {1.0, 0.0, -5.0};

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

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

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

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

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

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

  gsl_root_fdfsolver_free (s);
  return status;
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


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

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

/* The Type I Gumbel distribution has the form,

   p(x) dx = a b exp(-(b exp(-ax) + ax)) dx

   and the Type II Gumbel distribution has the form,

   p(x) dx = b a x^-(a+1) exp(-b x^-a)) dx

 */

double
gsl_ran_gumbel1 (const gsl_rng * r, const double a, const double b)
{
  double x = gsl_rng_uniform_pos (r);

  double z = (log(b) - log(-log(x))) / a;

  return z;
}

double
gsl_ran_gumbel1_pdf (const double x, const double a, const double b)
{
  double p = a * b *  exp (-(b * exp(-a * x) + a * x));
  return p;
}

double
gsl_ran_gumbel2 (const gsl_rng * r, const double a, const double b)
{
  double x = gsl_rng_uniform_pos (r);

  double z = pow(-b / log(x), 1/a);

  return z;
}

double
gsl_ran_gumbel2_pdf (const double x, const double a, const double b)
{
  if (x <= 0)
    {
      return 0 ;
    }
  else
    {
      double p = b * a *  pow(x,-(a+1)) * exp (-b * pow(x, -a));
      return p;
    }
}

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


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

#include "ccl.h"

//Spline creator
//n     -> number of points
//x     -> x-axis
//y     -> f(x)-axis
//y0,yf -> values of f(x) to use beyond the interpolation range
ccl_f1d_t *ccl_f1d_t_new(int n,double *x,double *y,double y0,double yf)
{
  ccl_f1d_t *spl=malloc(sizeof(ccl_f1d_t));
  if(spl==NULL)
    return NULL;

  spl->spline=gsl_spline_alloc(gsl_interp_cspline,n);
  if (spl->spline == NULL) {
    free(spl);
    return NULL;
  }
  int parstatus=gsl_spline_init(spl->spline,x,y,n);
  if(parstatus) {
    gsl_spline_free(spl->spline);
    free(spl);
    return NULL;
  }

  spl->x0=x[0];
  spl->xf=x[n-1];
  spl->y0=y0;
  spl->yf=yf;

  return spl;
}

//Evaluates spline at x checking for bound errors
double ccl_f1d_t_eval(ccl_f1d_t *spl,double x)
{
  if(x<=spl->x0)
    return spl->y0;
  else if(x>=spl->xf)
    return spl->yf;
  else {
    double y;
    int stat=gsl_spline_eval_e(spl->spline,x,NULL,&y);
    if (stat!=GSL_SUCCESS) {
      ccl_raise_gsl_warning(stat, "ccl_utils.c: ccl_splin_eval():");
      return NAN;
    }
    return y;
  }
}

//Spline destructor
void ccl_f1d_t_free(ccl_f1d_t *spl)
{
  if (spl != NULL) {
    gsl_spline_free(spl->spline);
  }
  free(spl);
}

max_stars_repo_name: benediktdiemer/CCL

max_issues_repo_name: benediktdiemer/CCL

max_forks_repo_name: benediktdiemer/CCL


text:
#include "../include/paralleltt.h"
#include <math.h>

#include <stdlib.h>
#include <stdio.h>
#include <mpi.h>
#include <cblas.h>
#include <string.h>
#include <lapacke.h>

#define min(a,b) ((a)>(b)?(b):(a))
#define abs(a) ((a)>(0)?(a):(-a))

#ifndef HEAD
#define HEAD (int) 0
#endif

matrix_tt* matrix_tt_init(const int m, const int n)
{
    matrix_tt* A = (matrix_tt*) malloc(sizeof(matrix_tt));
    A->m = m;
    A->n = n;
    A->transpose = 0;
    A->offset = 0;
    A->lda = m;
    A->X_size = (long) m*n;
    A->X = (double*) malloc(A->X_size*sizeof(double));

    return A;
}

void matrix_tt_wrap_update(matrix_tt* A, int m, int n, double* X)
{
    A->m = m;
    A->n = n;
    A->transpose = 0;
    A->offset = 0;
    A->lda = m;
    A->X_size = (long) m*n;
    A->X = X;
}

matrix_tt* matrix_tt_wrap(int m, int n, double* X)
{
    matrix_tt* A = (matrix_tt*) malloc(sizeof(matrix_tt));
    A->m = m;
    A->n = n;
    A->transpose = 0;
    A->offset = 0;
    A->lda = m;
    A->X_size = (long) m*n;
    A->X = X;
    return A;
}

matrix_tt* matrix_tt_copy(const matrix_tt* A)
{
    matrix_tt* B = (matrix_tt*) malloc(sizeof(matrix_tt));
    int m = A->m; int n = A->n; int X_size = A->X_size;
    B->m = m;
    B->n = n;
    B->transpose = A->transpose;
    B->offset = A->offset;
    B->lda = A->lda;
    B->X_size = X_size;
    B->X = (double*) malloc(X_size * sizeof(double));
    memcpy(B->X, A->X, X_size * sizeof(double));

    return B;
}

// Copy from matrix A to matrix B
void matrix_tt_copy_data(matrix_tt* B, const matrix_tt* A)
{
    int BT = (B->transpose == 0) ? 0 : 1;
    int AT = (A->transpose == 0) ? 0 : 1;
    if (BT == AT){
        if ((B->m != A->m) || (B->n != A->n)){
            printf("matrix_tt_copy_data: B (%d x %d) is not the same size as A (%d x %d)\n", B->n, B->m, A->n, A->m);
            return;
        }
        for (int jj = 0; jj < B->n; ++jj){
            memcpy((B->X) + (B->offset) + (B->lda)*jj, (A->X) + (A->offset) + jj*(A->lda), (B->m)*sizeof(double));
        }
    }
    else{
        if ((B->m != A->n) || (B->n != A->m)){
            printf("matrix_tt_copy_data: B (%d x %d, transpose = %d) is not the same size as A (%d x %d, transpose = %d)\n",
                    B->n, B->m, B->transpose, A->n, A->m, A->transpose);
            return;
        }
        for (int jj = 0; jj < B->n; ++jj){
            double* BX = (B->X) + (B->offset) + (B->lda)*jj;
            double* AX = (A->X) + (A->offset) + jj;
            for (int ii = 0; ii < B->m; ++ii){
                BX[ii] = AX[ii*(A->lda)];
            }
        }
    }
}

matrix_tt* submatrix_copy(const matrix_tt* A)
{
    matrix_tt* B = (matrix_tt*) malloc(sizeof(matrix_tt));
    int m = A->m; int n = A->n; int X_size = m*n;
    B->m = m;
    B->n = n;
    B->transpose = A->transpose;
    B->offset = 0;
    B->lda = m;
    B->X_size = X_size;
    B->X = (double*) malloc(X_size * sizeof(double));
    matrix_tt_copy_data(B, A);

    return(B);
}

// NOTE: reshape and submatrix will not work well together.
void matrix_tt_reshape(int m, int n, matrix_tt* A)
{
    long new_size = (long) m*n;
    if (new_size > A->X_size){
        printf("Reshape failed: m*n=%ld is too large for X_size=%ld\n", (long) new_size, A->X_size);
    }
    else{
        A->m = m;
        A->n = n;
        A->offset = 0;
        A->lda = m;
    }
}

// In matlab notation, returns A[ii1:ii2-1, jj1:jj2-1]
matrix_tt* submatrix(const matrix_tt* A, int ii1, int ii2, int jj1, int jj2)
{
    int mA = A->m; int nA = A->n; int lda = A->lda;
    if ((ii1 < 0) || (ii2 < ii1) || (mA < ii2) || (jj1 < 0) || (jj2 < jj1) || (nA < jj2)){
        printf("Cannot take the A[%d:%d,%d:%d] of a %d by %d matrix (Matlab notation)\n", ii1, ii2-1, jj1, jj2-1, mA, nA);
    }
    matrix_tt* B = (matrix_tt*) malloc(sizeof(matrix_tt));

    B->m = ii2 - ii1;
    B->n = jj2 - jj1;
    B->transpose = A->transpose;
    B->offset = ii1 + jj1*lda;
    B->lda = lda;
    B->X_size = A->X_size;
    B->X = A->X + A->offset;

    return B;
}

// Updates the indices of the submatrix A
void submatrix_update(matrix_tt* A, int ii1, int ii2, int jj1, int jj2)
{
    A->m = ii2 - ii1;
    A->n = jj2 - jj1;
    A->offset = ii1 + jj1*(A->lda);
}

// Get an element of the matrix
double matrix_tt_element(const matrix_tt* A, int ii, int jj)
{
    return A->X[A->offset + ii + (A->lda)*jj];
}

void matrix_tt_free(matrix_tt* A)
{
    free(A->X); A->X = NULL;
    free(A);
    return;
}

void matrix_tt_print(const matrix_tt* A, int just_the_matrix)
{
    int m = A->m; int n = A->n; int T = A->transpose;
    int offset = A->offset; int lda = A->lda;

    if (!just_the_matrix){
        printf("Size of matrix: m=%d, n=%d\n", m, n);
        printf("Transpose=%d, offset=%d, lda=%d\n", T, offset, lda);
        printf("X_size=%ld\n", A->X_size);
        printf("X");
        if (T){ printf("^T"); }
        printf(" is \n");
    }

    int mloop = T ? m : n;
    int nloop = T ? n : m;
    printf("[");
    for (int jj = 0; jj < nloop; ++jj){
        for (int ii = 0; ii < mloop; ++ii){
            if (!T){ printf("%f",A->X[lda*ii + jj + offset]); }
            else{    printf("%f",A->X[lda*jj + ii + offset]); }
            if (ii < mloop-1){ printf(", "); }
        }

        if (jj == nloop - 1){
            printf("]\n");
        }
        else{
            printf(",\n ");
        }
    }
}

double frobenius_norm(matrix_tt* A){
    double norm = 0;
    if (A->m == A->lda){
        norm = cblas_dnrm2(A->m * A->n, A->offset + A->X, 1);
    }
    else{
        for (int ii = 0; ii < A->n; ++ii){
            double tmp = cblas_dnrm2(A->m, A->offset + (A->lda)*ii + A->X, 1);
            norm = sqrt(norm*norm + tmp*tmp);
        }
    }
    return norm;
}

// Performs the operation C = alpha*A*B + beta*C, where A and B are appropriately transposed
int matrix_tt_dgemm(const matrix_tt* A, const matrix_tt* B, matrix_tt* C, const double alpha, const double beta)
{
    int mA = A->m; int nA = A->n;
    int mB = B->m; int nB = B->n;
    int m_check = C->m; int n_check = C->n; int k_check;

    int m; int n; int k;
    CBLAS_TRANSPOSE TransA;
    CBLAS_TRANSPOSE TransB;
    if (A->transpose == 0){
        TransA = CblasNoTrans;
        m = mA;
        k = nA;
    }
    else{
        TransA = CblasTrans;
        m = nA;
        k = mA;
    }

    if (B->transpose == 0){
        TransB = CblasNoTrans;
        k_check = mB;
        n = nB;
    }
    else{
        TransB = CblasTrans;
        k_check = nB;
        n = mB;
    }

    if ((m != m_check) || (n != n_check) || (k != k_check)){
        printf("Dimensions of input arrays do not match:\n");
        printf("A: m=%d, n=%d, tranpose=%d\n",mA,nA,A->transpose);
        printf("B: m=%d, n=%d, tranpose=%d\n",mB,nB,B->transpose);
        printf("C: m=%d, n=%d, tranpose=%d\n",m_check,n_check,C->transpose);
        return 1;
    }

    C->transpose = 0;

    cblas_dgemm(CblasColMajor, TransA, TransB,
                m, n, k,
                alpha,
                A->X + A->offset, A->lda,
                B->X + B->offset, B->lda,
                beta,
                C->X + C->offset, C->lda);

    return 0;
}

// Solves the least squares problem Ax = b
// At the end, it just copies the result into x. This probably isn't the most efficient way to do this
void matrix_tt_dgels(matrix_tt* x, matrix_tt* A, const matrix_tt* b)
{
    if (x->transpose != 0){
        printf("matrix_tt_least_squares: x cannot be transposed for dgels\n");
        return;
    }

    if (b->transpose != 0){
        printf("matrix_tt_least_squares: b cannot be transposed for dgels\n");
        return;
    }

    int nrhs = x->n;
    char TransA;
    int m_check;
    int n_check;

    if (A->transpose == 0){
        TransA = 'N';
        m_check = A->m;
        n_check = A->n;
    }
    else{
        TransA = 'T';
        m_check = A->n;
        n_check = A->m;
    }

    if (n_check != x->m){
        printf("right index of A (%d) does not match left index of x (%d)\n", n_check, x->m);
    }
    if (m_check != b->m){
        printf("left index of A (%d) does not match left index of b (%d)\n", m_check, b->m);
    }
    if (x->n != b->n){
        printf("right index of x (%d) does not match right index of b (%d)\n", x->n, b->n);
    }

    int m = A->m;
    int n = A->n;

    int info = LAPACKE_dgels(LAPACK_COL_MAJOR, TransA, m, n, nrhs, A->X + A->offset, A->lda, b->X + b->offset, b->lda);

    matrix_tt* result = submatrix(b, 0, x->m, 0, x->n);
    matrix_tt_copy_data(x, result);
    free(result);
}

// Input:
//   Q - the matrix you want to take the QR
// Output:
//   Q - the Q matrix of QR
//   R - (input is NULL) ? nothing : the R matrix of QR

// NOTE: Currently assumes R and Q are not transposed
int matrix_tt_truncated_qr(matrix_tt* Q, matrix_tt* R, int r)
{
    int Qm = Q->m; int Qn = Q->n; double* QX = Q->X;
    int Qtranspose = Q->transpose; long Qoffset = Q->offset; int Qlda = Q->lda;

    lapack_int info;

    if (Qtranspose == 0){
        if (r > Qn)
        {
            printf("ERROR: (matrix_tt_truncated_qr) r = %d is too large for Q->n = %d \n", r, Qn);
            return 1;
        }
        // The QR step
        double* tau = (double*) malloc(sizeof(double)*min(Qm,Qn));
        info = LAPACKE_dgeqrf(LAPACK_COL_MAJOR, Qm, Qn, QX + Qoffset, Qlda, tau);

        if (R != NULL){
            int Rm = R->m; int Rn = R->n; double* RX = R->X;
            int Rtranspose = R->transpose; long Roffset = R->offset; int Rlda = R->lda;

            if ((Rm != r) || (Rn != r)){
                printf("ERROR: (matrix_tt_truncated_qr) R->m = %d or R->n = %d is not equal to r = %d\n", Rm, Rn, r);
                return 1;
            }
            if (Rtranspose != 0){
                printf("ERROR: (matrix_tt_truncated_qr) Only defined for R->transpose = 0\n");
                return 1;
            }

            for (int jj = 0; jj < r; ++jj){
                int r_col_length = (Qm < jj+1) ? (Qm) : jj+1;
                memcpy(RX + Roffset + jj*Rlda, QX + Qoffset + jj*Qlda, r_col_length*sizeof(double));
                memset(RX + Roffset + jj*Rlda + r_col_length, 0, (r-r_col_length)*sizeof(double));
            }
        }

        if (r > Qm){
            LAPACKE_dorgqr(LAPACK_COL_MAJOR, Qm, Qm, Qm, QX + Qoffset, Qlda, tau);
            matrix_tt* Q_sub = submatrix(Q, 0, Qm, Qm, r);
            matrix_tt_fill_zeros(Q_sub);
            free(Q_sub);
        }
        else{
            LAPACKE_dorgqr(LAPACK_COL_MAJOR, Qm, r, r, QX + Qoffset, Qlda, tau);
        }

        Q->n = r;

        free(tau);       tau = NULL;
    }
    else{
        if (r > Qm)
        {
            printf("ERROR: (matrix_tt_truncated_qr) r = %d is too large for Q->n = %d \n", r, Qn);
            return 1;
        }
        // The LQ step
        printf("Performing LQ\n");
        double* tau = (double*) malloc(sizeof(double)*min(Qm,Qn));
        info = LAPACKE_dgelqf(LAPACK_COL_MAJOR, Qm, Qn, QX + Qoffset, Qlda, tau);
        printf("\nRight after LQ, Q = \n");
        matrix_tt_print(Q, 1);

        if (R != NULL){
            int Rm = R->m; int Rn = R->n; double* RX = R->X;
            int Rtranspose = R->transpose; long Roffset = R->offset; int Rlda = R->lda;

            if ((Rm != r) || (Rn != r)){
                printf("ERROR: (matrix_tt_truncated_qr) R->m = %d or R->n = %d is not equal to r = %d\n", Rm, Rn, r);
                return 1;
            }
            if (Rtranspose == 0){
                printf("ERROR: (matrix_tt_truncated_qr) Only defined for R->transpose = Q->transpose\n");
                return 1;
            }

            int ncols = (r > Qm) ? Qm : r;
            for (int jj = 0; jj < r; ++jj){
                if (jj < Qm){
                    memcpy(RX + Roffset + jj*Rlda + jj, QX + Qoffset + jj*Qlda + jj, (Qm - jj)*sizeof(double));
                    memset(RX + Roffset + jj*Rlda, 0, jj*sizeof(double));
                }
            }

        }

        if (r > Qn){
            printf("Doing first dorglq\n");
            LAPACKE_dorglq(LAPACK_COL_MAJOR, Qn, Qn, Qn, QX + Qoffset, Qlda, tau);
            matrix_tt* Q_sub = submatrix(Q, Qn, r, 0, Qn);
            matrix_tt_fill_zeros(Q_sub);
            free(Q_sub);
        }
        else{
            printf("Doing second dorglq\n");
            LAPACKE_dorglq(LAPACK_COL_MAJOR, Qm, Qn, Qm, QX + Qoffset, Qlda, tau);
        }

        Q->m = r;

        free(tau);       tau = NULL;
    }

    return 0;
}/**/

void matrix_tt_group_reduce(MPI_Comm comm, int rank, matrix_tt* A, matrix_tt* buf, int head, int* group_ranks, int nranks)
{
    int in_the_group = 0;
    for (int ii = 0; ii < nranks; ++ii){
        if (rank == group_ranks[ii]){
            in_the_group = 1;
        }
    }

    if (in_the_group){
        int buf_assigned = 0;
        if (!buf){
            buf_assigned = 1;
            buf = matrix_tt_init(A->m, A->n);
        }

        if ((buf->m < A->m) || (buf->n < A->n)){
            printf("buf (%d x %d) must be larger than A (%d x %d)\n", buf->m, buf->n, A->m, A->n);
        }

        int send = -1;
        int recv = -1;
        while (nranks > 1){
            int half_nranks = (nranks + 1) / 2;
            for (int ii = 0; ii < half_nranks; ++ii){
                if (ii + half_nranks >= nranks){
                    recv = group_ranks[ii];
                    send = -1;
                }
                else if (group_ranks[ii + half_nranks] == head){
                    recv = head;
                    send = group_ranks[ii];
                }
                else{
                    recv = group_ranks[ii];
                    send = group_ranks[ii + half_nranks];
                }

                if (rank == send){
                    matrix_tt_send(comm, A, recv);
                }
                else if ((rank == recv) && (send != -1)){
                    matrix_tt_recv(comm, buf, send);
                    for (int jj = 0; jj < A->n; ++jj){
//                        MPI_Recv(buf->X + buf->offset, A->m, MPI_DOUBLE, send, 0, comm, MPI_STATUS_IGNORE);
                        long A_offset = A->offset + jj * (A->lda);
                        long buf_offset = buf->offset +  jj * (buf->lda);
                        for (int kk = 0; kk < A->m; ++kk){
                            A->X[A_offset + kk] = A->X[A_offset + kk] + buf->X[buf_offset + kk];
                        }
                    }
                }
                group_ranks[ii] = recv;
            }
            nranks = half_nranks;
        }

        if (buf_assigned){
            matrix_tt_free(buf);
        }
    }

}

// buf just has to be size A->lda x 1 or more
void matrix_tt_reduce(MPI_Comm comm, int rank, matrix_tt* A, matrix_tt* buf, int head)
{
    if (buf == NULL){
        matrix_tt* buf = matrix_tt_init(A->m, 1);
        matrix_tt_reduce(comm, rank, A, buf, head);
        matrix_tt_free(buf);
    }
    else{
        int size;
        MPI_Comm_size(comm, &size);

        if ((head < 0) || (head >= size)){
            printf("matrix_tt_reduce: head = %d is not a valid rank for size = %d\n", head, size);
        }

        if (buf->lda < A->m){
            printf("buf->lda = %d must be larger than A->m = %d\n", buf->lda, A->m);
        }

        int* activated_ranks = (int*) calloc(size, sizeof(int));
        for (int ii = 0; ii < size; ++ii){
            activated_ranks[ii] = ii;
        }
        int send = -1;
        int recv = -1;
        while (size > 1){
            int half_size = (size + 1) / 2;
            for (int ii = 0; ii < half_size; ++ii){
                if (ii + half_size >= size){
                    recv = activated_ranks[ii];
                    send = -1;
                }
                else if (activated_ranks[ii + half_size] == head){
                    recv = head;
                    send = activated_ranks[ii];
                }
                else{
                    recv = activated_ranks[ii];
                    send = activated_ranks[ii + half_size];
                }

                if (rank == send){
                    for (int jj = 0; jj < A->n; ++jj){
                        MPI_Send(A->X + A->offset + A->lda * jj, A->m, MPI_DOUBLE, recv, 0, comm);
                    }
                }
                else if ((rank == recv) && (send != -1)){
                    for (int jj = 0; jj < A->n; ++jj){
                        MPI_Recv(buf->X + buf->offset, A->m, MPI_DOUBLE, send, 0, comm, MPI_STATUS_IGNORE);
                        long A_offset = A->offset + jj * (A->lda);
                        for (int kk = 0; kk < A->m; ++kk){
                            A->X[A_offset + kk] = A->X[A_offset + kk] + buf->X[buf->offset + kk];
                        }
                    }
                }
                activated_ranks[ii] = recv;
            }
            size = half_size;
        }
        free(activated_ranks);
    }
}

void matrix_tt_allreduce(MPI_Comm comm, matrix_tt* A)
{
    for (int jj = 0; jj < A->n; ++jj){
        MPI_Allreduce(MPI_IN_PLACE, A->X + A->offset + jj*(A->lda), A->m, MPI_DOUBLE, MPI_SUM, comm);
    }
}

void matrix_tt_broadcast(MPI_Comm comm, matrix_tt* buf)
{
    int count = buf->m * buf->n;
    MPI_Bcast(buf->X, count, MPI_DOUBLE, HEAD, comm);
}

void matrix_tt_send(MPI_Comm comm, matrix_tt* buf, int dest){
    MPI_Datatype buf_type;
    MPI_Type_vector(buf->n, buf->m, buf->lda, MPI_DOUBLE, &buf_type);
    MPI_Type_commit(&buf_type);
    MPI_Send(buf->X + buf->offset, 1, buf_type, dest, 0, comm);
    MPI_Type_free(&buf_type);
}

void matrix_tt_recv(MPI_Comm comm, matrix_tt* buf, int source){
    MPI_Datatype buf_type;
    MPI_Type_vector(buf->n, buf->m, buf->lda, MPI_DOUBLE, &buf_type);
    MPI_Type_commit(&buf_type);
    MPI_Recv(buf->X + buf->offset, 1, buf_type, source, 0, comm, MPI_STATUS_IGNORE);
    MPI_Type_free(&buf_type);
}

// Column khatri-rao product C = A x B
void khatri_rao(const matrix_tt* restrict A, const matrix_tt* restrict B, matrix_tt* restrict C){
    int n = A->n;
    int m_A = A->m;

    int n_B = B->n;
    int m_B = B->m;

    int n_C = C->n;
    int m_C = C->m;

    if ((n != n_B) || (n != n_C) || (m_A*m_B != m_C)){
        printf("khatri_rao: the sizes dont work!\n");
        printf("A is %d x %d\n", m_A, n);
        printf("B is %d x %d\n", m_B, n_B);
        printf("C is %d x %d\n", m_C, n_C);
    }

    for (int kk = 0; kk < n; ++kk){
        for (int ii = 0; ii < m_B; ++ii){
            double Bii = B->X[kk*m_B + ii];
            for (int jj = 0; jj < m_A; ++jj){
                C->X[kk*m_C + ii*m_B + jj] = Bii * A->X[kk*m_A + jj];
            }
        }
    }
}

// Recursively KR multiply the list of matrices A
void list_khatri_rao(int d_kr, matrix_tt** A, matrix_tt* Omega){
    if (d_kr == 1){
        for (int ii = 0; ii < A[0]->X_size; ++ii){
            Omega->X[ii] = A[0]->X[ii];
        }
        return;
    }

    if (d_kr == 2){
        khatri_rao(A[0], A[1], Omega);
    }
    else{
        int m_B = 1;
        int n_B = Omega->n;
        for (int ii = 0; ii < d_kr-1; ++ii){
            m_B = m_B * A[ii]->m;
        }

        matrix_tt* B = matrix_tt_init(m_B, n_B);
        list_khatri_rao(d_kr - 1, A, B);
        khatri_rao(B, A[d_kr - 1], Omega);
        matrix_tt_free(B); B = NULL;
    }

}

void matrix_tt_dlarnv(matrix_tt* A){
    int r1 = rand()%4096, r2 = rand()%4096, r3 = rand()%4096, r4 = rand()%4096;
    int iseed[4] = {r1, r2, r3, r4+(r4%2 == 0?1:0)};
    LAPACKE_dlarnv(3, iseed, (A->n) * (A->m), A->X);
}

void matrix_tt_fill_zeros(matrix_tt* A)
{
    for (int jj = 0; jj < A->n; ++jj){
        memset(A->X + A->offset + jj*A->lda, 0, A->m * sizeof(double));
    }
}

// Takes a matrix in row major and transforms it to column major
void row_to_col_major(const matrix_tt* mat_row, matrix_tt* mat_col)
{
    int m = mat_col->m;
    int n = mat_col->n;
    if ((m != mat_row->n) || (n != mat_row->m)){
        printf("The dimensions don't work out! mat_row should have the dimensions of mat_col transposed\n");
        printf("mat_row->m = %d, mat_row->n = %d\n", mat_row->m, mat_row->n);
        printf("mat_col->m = %d, mat_col->n = %d\n", mat_col->m, mat_col->n);
    }

    for (int ii = 0; ii < m; ++ii){
        for (int jj = 0; jj < n; ++jj){
            mat_col->X[jj * (mat_col->lda) + ii + mat_col->offset] = mat_row->X[ii * (mat_row->lda) + jj + mat_row->offset];
        }
    }
}


text:
#pragma once
#include "utility/types.h"
#include "utility/c++std/optional.h"
#include "utility/c++std/string_view.h"
#include <fmt/format.h>
#include <gsl/gsl>

namespace cws80 {

class NkScreen;
class NativeUI;
struct Bank;
struct Program;

///
class UIController {
public:
    virtual ~UIController() {}

    virtual NkScreen &screen() const = 0;
    virtual NativeUI &native_ui() const = 0;
    virtual const Program &program() const = 0;
    virtual uint bank_number() const = 0;
    virtual uint program_number() const = 0;
    virtual gsl::span<const std::string> program_names() const = 0;
    virtual i32 get_parameter(uint idx) const = 0;
    virtual bool set_parameter(uint idx, i32 val) = 0;
    virtual f32 get_f32_parameter(uint idx) const = 0;
    virtual bool set_f32_parameter(uint idx, f32 val) = 0;
    virtual void currently_editing_parameter(cxx::optional<uint> idx) = 0;
    virtual void request_next_bank() = 0;
    virtual void request_bank(uint num) = 0;
    virtual void request_program(uint num) = 0;
    virtual void request_rename_program(cxx::string_view name) = 0;
    virtual void request_init_program() = 0;
    virtual void request_write_program() = 0;
    virtual void dlg_rename_program() = 0;
    virtual void dlg_load_bank() = 0;
    virtual void dlg_save_bank() = 0;
    virtual void send_piano_events(const i8 events[128]) = 0;
    virtual const std::string &led_message() = 0;
    virtual void led_message(std::string msg) = 0;
    virtual void led_priority_message(f64 timeout, std::string msg) = 0;
    virtual const std::string &status_message() = 0;
    virtual void status_message(std::string msg) = 0;
    //
    template <class A, class... As>
    void led_fmt(const char *fmt, const A &x, const As &... xs)
    {
        led_message(fmt::format(fmt, x, xs...));
    }
    template <class A, class... As>
    void led_priority_fmt(f64 timeout, const char *fmt, const A &x, const As &... xs)
    {
        led_priority_message(timeout, fmt::format(fmt, x, xs...));
    }
    //
    template <class A, class... As>
    void status_fmt(const char *fmt, const A &x, const As &... xs)
    {
        status_message(fmt::format(fmt, x, xs...));
    }
    //
    enum {
        debugging_input = true,
        debugging_logic = true,
    };
    inline void debug_input(std::string str)
    {
        if (debugging_input)
            status_message(std::move(str));
    }
    inline void debug_logic(std::string str)
    {
        if (debugging_logic)
            status_message(std::move(str));
    }
    template <class A, class... As>
    void debug_input_fmt(const char *fmt, const A &x, const As &... xs)
    {
        if (debugging_input)
            status_fmt(fmt, x, xs...);
    }
    template <class A, class... As>
    void debug_logic_fmt(const char *fmt, const A &x, const As &... xs)
    {
        if (debugging_logic)
            status_fmt(fmt, x, xs...);
    }
};

}  // namespace cws80


text:
/* gsl_sf_hermite.h
 *
 * Copyright (C) 2011-2014 Konrad Griessinger
 *
 * 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.
 */

/*----------------------------------------------------------------------*
 * (konradg(at)gmx.net)                                                 *
 *----------------------------------------------------------------------*/

#ifndef __GSL_SF_HERMITE_H__
#define __GSL_SF_HERMITE_H__

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

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

GSL_FUN int gsl_sf_hermite_prob_e(const int n, const double x, gsl_sf_result * result);
GSL_FUN double gsl_sf_hermite_prob(const int n, const double x);
GSL_FUN int gsl_sf_hermite_prob_deriv_e(const int m, const int n, const double x, gsl_sf_result * result);
GSL_FUN double gsl_sf_hermite_prob_deriv(const int m, const int n, const double x);
GSL_FUN int gsl_sf_hermite_e(const int n, const double x, gsl_sf_result * result);
GSL_FUN double gsl_sf_hermite(const int n, const double x);
GSL_FUN int gsl_sf_hermite_deriv_e(const int m, const int n, const double x, gsl_sf_result * result);
GSL_FUN double gsl_sf_hermite_deriv(const int m, const int n, const double x);
GSL_FUN int gsl_sf_hermite_func_e(const int n, const double x, gsl_sf_result * result);
GSL_FUN double gsl_sf_hermite_func(const int n, const double x);
GSL_FUN int gsl_sf_hermite_func_fast_e(const int n, const double x, gsl_sf_result * result);
GSL_FUN double gsl_sf_hermite_func_fast(const int n, const double x);
GSL_FUN int gsl_sf_hermite_prob_array(const int nmax, const double x, double * result_array);
GSL_FUN int gsl_sf_hermite_prob_array_deriv(const int m, const int nmax, const double x, double * result_array);
GSL_FUN int gsl_sf_hermite_prob_deriv_array(const int mmax, const int n, const double x, double * result_array);
GSL_FUN int gsl_sf_hermite_prob_series_e(const int n, const double x, const double * a, gsl_sf_result * result);
GSL_FUN double gsl_sf_hermite_prob_series(const int n, const double x, const double * a);
GSL_FUN int gsl_sf_hermite_array(const int nmax, const double x, double * result_array);
GSL_FUN int gsl_sf_hermite_array_deriv(const int m, const int nmax, const double x, double * result_array);
GSL_FUN int gsl_sf_hermite_deriv_array(const int mmax, const int n, const double x, double * result_array);
GSL_FUN int gsl_sf_hermite_series_e(const int n, const double x, const double * a, gsl_sf_result * result);
GSL_FUN double gsl_sf_hermite_series(const int n, const double x, const double * a);
GSL_FUN int gsl_sf_hermite_func_array(const int nmax, const double x, double * result_array);
GSL_FUN int gsl_sf_hermite_func_series_e(const int n, const double x, const double * a, gsl_sf_result * result);
GSL_FUN double gsl_sf_hermite_func_series(const int n, const double x, const double * a);
GSL_FUN int gsl_sf_hermite_func_der_e(const int m, const int n, const double x, gsl_sf_result * result);
GSL_FUN double gsl_sf_hermite_func_der(const int m, const int n, const double x);
GSL_FUN int gsl_sf_hermite_prob_zero_e(const int n, const int s, gsl_sf_result * result);
GSL_FUN double gsl_sf_hermite_prob_zero(const int n, const int s);
GSL_FUN int gsl_sf_hermite_zero_e(const int n, const int s, gsl_sf_result * result);
GSL_FUN double gsl_sf_hermite_zero(const int n, const int s);
GSL_FUN int gsl_sf_hermite_func_zero_e(const int n, const int s, gsl_sf_result * result);
GSL_FUN double gsl_sf_hermite_func_zero(const int n, const int s);

#ifndef GSL_DISABLE_DEPRECATED

GSL_FUN int gsl_sf_hermite_phys_e(const int n, const double x, gsl_sf_result * result);
GSL_FUN double gsl_sf_hermite_phys(const int n, const double x);
GSL_FUN int gsl_sf_hermite_phys_der_e(const int m, const int n, const double x, gsl_sf_result * result);
GSL_FUN double gsl_sf_hermite_phys_der(const int m, const int n, const double x);
GSL_FUN int gsl_sf_hermite_phys_array(const int nmax, const double x, double * result_array);
GSL_FUN int gsl_sf_hermite_phys_series_e(const int n, const double x, const double * a, gsl_sf_result * result);
GSL_FUN double gsl_sf_hermite_phys_series(const int n, const double x, const double * a);
GSL_FUN int gsl_sf_hermite_phys_array_der(const int m, const int nmax, const double x, double * result_array);
GSL_FUN int gsl_sf_hermite_phys_der_array(const int mmax, const int n, const double x, double * result_array);
GSL_FUN int gsl_sf_hermite_phys_zero_e(const int n, const int s, gsl_sf_result * result);
GSL_FUN double gsl_sf_hermite_phys_zero(const int n, const int s);

GSL_FUN int gsl_sf_hermite_prob_array_der(const int m, const int nmax, const double x, double * result_array);
GSL_FUN int gsl_sf_hermite_prob_der_array(const int mmax, const int n, const double x, double * result_array);
GSL_FUN int gsl_sf_hermite_prob_der_e(const int m, const int n, const double x, gsl_sf_result * result);
GSL_FUN double gsl_sf_hermite_prob_der(const int m, const int n, const double x);

#endif /* !GSL_DISABLE_DEPRECATED */

__END_DECLS

#endif /* __GSL_SF_HERMITE_H__ */

max_stars_repo_name: mgreter/astrometrylib

max_issues_repo_name: zzpwahaha/VimbaCamJILA

max_forks_repo_name: zzpwahaha/VimbaCamJILA


text:
/*
 * Copyright (c) 1997-1999 Massachusetts Institute of Technology
 *
 * 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., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 */

/* fftw.h -- system-wide definitions */
/* $Id: fftw-int.h,v 1.39 1999/02/19 17:22:00 athena Exp $ */

#ifndef FFTW_INT_H
#define FFTW_INT_H
#include <config.h>
#include <fftw.h>

#ifdef __cplusplus
extern "C" {
#else
#endif				/* __cplusplus */

/****************************************************************************/
/*                            Private Functions                             */
/****************************************************************************/

extern fftw_twiddle *fftw_create_twiddle(int n, const fftw_codelet_desc *d);
extern void fftw_destroy_twiddle(fftw_twiddle *tw);

extern void fftw_strided_copy(int, fftw_complex *, int, fftw_complex *);
extern void fftw_executor_simple(int, const fftw_complex *, fftw_complex *,
				 fftw_plan_node *, int, int);

extern fftwnd_plan fftwnd_create_plan_aux(int rank, const int *n,
					  fftw_direction dir, int flags);
extern fftw_plan *fftwnd_new_plan_array(int rank);
extern fftw_plan *fftwnd_create_plans_generic(fftw_plan *plans,
					      int rank, const int *n,
					      fftw_direction dir, int flags);
extern fftw_plan *fftwnd_create_plans_specific(fftw_plan *plans,
					       int rank, const int *n,
					       const int *n_after,
					       fftw_direction dir, int flags,
					       fftw_complex *in, int istride,
					       fftw_complex *out, int ostride);
extern int fftwnd_work_size(int rank, const int *n, int flags, int ncopies);

extern void fftwnd_aux(fftwnd_plan p, int cur_dim,
		       fftw_complex *in, int istride,
		       fftw_complex *out, int ostride,
		       fftw_complex *work);
extern void fftwnd_aux_howmany(fftwnd_plan p, int cur_dim,
			       int howmany,
			       fftw_complex *in, int istride, int idist,
			       fftw_complex *out, int ostride, int odist,
			       fftw_complex *work);

/* wisdom prototypes */
enum fftw_wisdom_category {
     FFTW_WISDOM, RFFTW_WISDOM
};

extern int fftw_wisdom_lookup(int n, int flags, fftw_direction dir,
			      enum fftw_wisdom_category category,
			      int istride, int ostride,
			      enum fftw_node_type *type,
			      int *signature, int replace_p);
extern void fftw_wisdom_add(int n, int flags, fftw_direction dir,
			    enum fftw_wisdom_category cat,
			    int istride, int ostride,
			    enum fftw_node_type type,
			    int signature);

/* Private planner functions: */
extern double fftw_estimate_node(fftw_plan_node *p);
extern fftw_plan_node *fftw_make_node_notw(int size,
					const fftw_codelet_desc *config);
extern fftw_plan_node *fftw_make_node_real2hc(int size,
					const fftw_codelet_desc *config);
extern fftw_plan_node *fftw_make_node_hc2real(int size,
					const fftw_codelet_desc *config);
extern fftw_plan_node *fftw_make_node_twiddle(int n,
					 const fftw_codelet_desc *config,
					      fftw_plan_node *recurse,
					      int flags);
extern fftw_plan_node *fftw_make_node_hc2hc(int n,
					    fftw_direction dir,
					 const fftw_codelet_desc *config,
					    fftw_plan_node *recurse,
					    int flags);
extern fftw_plan_node *fftw_make_node_generic(int n, int size,
					      fftw_generic_codelet *codelet,
					      fftw_plan_node *recurse,
					      int flags);
extern fftw_plan_node *fftw_make_node_rgeneric(int n, int size,
					       fftw_direction dir,
					       fftw_rgeneric_codelet * codelet,
					       fftw_plan_node *recurse,
					       int flags);
extern int fftw_factor(int n);
extern fftw_plan_node *fftw_make_node(void);
extern fftw_plan fftw_make_plan(int n, fftw_direction dir,
				fftw_plan_node *root, int flags,
				enum fftw_node_type wisdom_type,
				int wisdom_signature);
extern void fftw_use_plan(fftw_plan p);
extern void fftw_use_node(fftw_plan_node *p);
extern void fftw_destroy_plan_internal(fftw_plan p);
extern fftw_plan fftw_pick_better(fftw_plan p1, fftw_plan p2);
extern fftw_plan fftw_lookup(fftw_plan *table, int n, int flags);
extern void fftw_insert(fftw_plan *table, fftw_plan this_plan, int n);
extern void fftw_make_empty_table(fftw_plan *table);
extern void fftw_destroy_table(fftw_plan *table);
extern void fftw_complete_twiddle(fftw_plan_node *p, int n);

extern fftw_plan_node *fftw_make_node_rader(int n, int size,
					    fftw_direction dir,
					    fftw_plan_node *recurse,
					    int flags);
extern fftw_rader_data *fftw_rader_top;

/****************************************************************************/
/*                           Floating Point Types                           */
/****************************************************************************/

/*
 * We use these definitions to make it easier for people to change
 * FFTW to use long double and similar types. You shouldn't have to
 * change this just to use float or double.
 */

/*
 * Change this if your floating-point constants need to be expressed
 * in a special way.  For example, if fftw_real is long double, you
 * will need to append L to your fp constants to make them of the
 * same precision.  Do this by changing "x" below to "x##L".
 */
#define FFTW_KONST(x) ((fftw_real) x)

#define FFTW_TRIG_SIN sin
#define FFTW_TRIG_COS cos
typedef double FFTW_TRIG_REAL;	/* the argument type for sin and cos */

#define FFTW_K2PI FFTW_KONST(6.2831853071795864769252867665590057683943388)

/****************************************************************************/
/*                               gcc/x86 hacks                              */
/****************************************************************************/

/*
 * gcc 2.[78].x and x86 specific hacks.  These macros align the stack
 * pointer so that the double precision temporary variables in the
 * codelets will be aligned to a multiple of 8 bytes (*way* faster on
 * pentium and pentiumpro)
 */
#ifdef __GNUC__
#ifdef __i386__
#ifdef FFTW_ENABLE_I386_HACKS
#ifndef FFTW_ENABLE_FLOAT
#define FFTW_USING_I386_HACKS
#define HACK_ALIGN_STACK_EVEN() {                        \
     if ((((long) (__builtin_alloca(0))) & 0x7)) __builtin_alloca(4);    \
}

#define HACK_ALIGN_STACK_ODD() {                         \
     if (!(((long) (__builtin_alloca(0))) & 0x7)) __builtin_alloca(4);   \
}

#ifdef FFTW_DEBUG_ALIGNMENT
#define ASSERT_ALIGNED_DOUBLE() {                        \
     double __foo;                                       \
     if ((((long) &__foo) & 0x7)) abort();                \
}
#endif

#endif
#endif
#endif
#endif

#ifndef HACK_ALIGN_STACK_EVEN
#define HACK_ALIGN_STACK_EVEN()
#endif
#ifndef HACK_ALIGN_STACK_ODD
#define HACK_ALIGN_STACK_ODD()
#endif
#ifndef ASSERT_ALIGNED_DOUBLE
#define ASSERT_ALIGNED_DOUBLE()
#endif

/****************************************************************************/
/*                                  Timers                                  */
/****************************************************************************/

/*
 * Here, you can use all the nice timers available in your machine.
 */

/*
 *
 Things you should define to include your own clock:

 fftw_time -- the data type used to store a time

 extern fftw_time fftw_get_time(void);
 -- a function returning the current time.  (We have
 implemented this as a macro in most cases.)

 extern fftw_time fftw_time_diff(fftw_time t1, fftw_time t2);
 -- returns the time difference (t1 - t2).
 If t1 < t2, it may simply return zero (although this
 is not required).  (We have implemented this as a macro
 in most cases.)

 extern double fftw_time_to_sec(fftw_time t);
 -- returns the time t expressed in seconds, as a double.
 (Implemented as a macro in most cases.)

 FFTW_TIME_MIN -- a double-precision macro holding the minimum
 time interval (in seconds) for accurate time measurements.
 This should probably be at least 100 times the precision of
 your clock (we use even longer intervals, to be conservative).
 This will determine how long the planner takes to measure
 the speeds of different possible plans.

 Bracket all of your definitions with an appropriate #ifdef so that
 they will be enabled on your machine.  If you do add your own
 high-precision timer code, let us know (at fftw@theory.lcs.mit.edu).

 Only declarations should go in this file.  Any function definitions
 that you need should go into timer.c.
 */

/*
 * define a symbol so that we know that we have the fftw_time_diff
 * function/macro (it did not exist prior to FFTW 1.2)
 */
#define FFTW_HAS_TIME_DIFF

/**********************************************
 *              SOLARIS
 **********************************************/
#if defined(HAVE_GETHRTIME)

/* we use the nanosecond virtual timer */
#ifdef HAVE_SYS_TIME_H
#include <sys/time.h>
#endif

typedef hrtime_t fftw_time;

#define fftw_get_time() gethrtime()
#define fftw_time_diff(t1,t2) ((t1) - (t2))
#define fftw_time_to_sec(t) ((double) t / 1.0e9)

/*
 * a measurement is valid if it runs for at least
 * FFTW_TIME_MIN seconds.
 */
#define FFTW_TIME_MIN (1.0e-4)	/* for Solaris nanosecond timer */
#define FFTW_TIME_REPEAT 8

/**********************************************
 *        Pentium time stamp counter
 **********************************************/
#elif defined(__GNUC__) && defined(__i386__) && defined(FFTW_ENABLE_PENTIUM_TIMER)

/*
 * Use internal Pentium register (time stamp counter). Resolution
 * is 1/FFTW_CYCLES_PER_SEC seconds (e.g. 5 ns for Pentium 200 MHz).
 * (This code was contributed by Wolfgang Reimer)
 */

#ifndef FFTW_CYCLES_PER_SEC
#error "Must define FFTW_CYCLES_PER_SEC in fftw/config.h to use the Pentium cycle counter"
#endif

typedef unsigned long long fftw_time;

static __inline__ fftw_time read_tsc()
{
     struct {
	  long unsigned lo, hi;
     } counter;
     long unsigned sav_eax, sav_edx;
   __asm__("movl %%eax,%0":"=m"(sav_eax));
   __asm__("movl %%edx,%0":"=m"(sav_edx));
     __asm__("rdtsc");
   __asm__("movl %%eax,%0":"=m"(counter.lo));
   __asm__("movl %%edx,%0":"=m"(counter.hi));
   __asm__("movl %0,%%eax": : "m"(sav_eax):"eax");
   __asm__("movl %0,%%edx": : "m"(sav_edx):"edx");
     return *(fftw_time *) & counter;
}

#define fftw_get_time()  read_tsc()
#define fftw_time_diff(t1,t2) ((t1) - (t2))
#define fftw_time_to_sec(t) (((double) (t)) / FFTW_CYCLES_PER_SEC)
#define FFTW_TIME_MIN (1.0e-4)	/* for Pentium TSC register */

/************* generic systems having gettimeofday ************/
#elif defined(HAVE_GETTIMEOFDAY) || defined(HAVE_BSDGETTIMEOFDAY)
#ifdef HAVE_SYS_TIME_H
#include <sys/time.h>
#endif
#ifdef HAVE_UNISTD_H
#include <unistd.h>
#endif
#define FFTW_USE_GETTIMEOFDAY

typedef struct timeval fftw_time;

extern fftw_time fftw_gettimeofday_get_time(void);
extern fftw_time fftw_gettimeofday_time_diff(fftw_time t1, fftw_time t2);
#define fftw_get_time() fftw_gettimeofday_get_time()
#define fftw_time_diff(t1, t2) fftw_gettimeofday_time_diff(t1, t2)
#define fftw_time_to_sec(t) ((double)(t).tv_sec + (double)(t).tv_usec * 1.0E-6)

#ifndef FFTW_TIME_MIN
/* this should be fine on any system claiming a microsecond timer */
#define FFTW_TIME_MIN (1.0e-2)
#endif

/**********************************************
 *              MACINTOSH
 **********************************************/
#elif defined(HAVE_MAC_TIMER)

/*
 * By default, use the microsecond-timer in the Mac Time Manager.
 * Alternatively, by changing the following #if 1 to #if 0, you
 * can use the nanosecond timer available *only* on PCI PowerMacs.
 */
#ifndef HAVE_MAC_PCI_TIMER	/* use time manager */

/*
 * Use Macintosh Time Manager routines (maximum resolution is about 20
 * microseconds).
 */
typedef struct fftw_time_struct {
     unsigned long hi, lo;
} fftw_time;

extern fftw_time get_Mac_microseconds(void);

#define fftw_get_time() get_Mac_microseconds()

/* define as a function instead of a macro: */
extern fftw_time fftw_time_diff(fftw_time t1, fftw_time t2);

#define fftw_time_to_sec(t) ((t).lo * 1.0e-6 + 4294967295.0e-6 * (t).hi)

/* very conservative, since timer should be accurate to 20e-6: */
/* (although this seems not to be the case in practice) */
#define FFTW_TIME_MIN (5.0e-2)	/* for MacOS Time Manager timer */

#else				/* use nanosecond timer */

/* Use the nanosecond timer available on PCI PowerMacs. */

#include <DriverServices.h>

typedef AbsoluteTime fftw_time;
#define fftw_get_time() UpTime()
#define fftw_time_diff(t1,t2) SubAbsoluteFromAbsolute(t1,t2)
#define fftw_time_to_sec(t) (AbsoluteToNanoseconds(t).lo * 1.0e-9)

/* Extremely conservative minimum time: */
/* for MacOS PCI PowerMac nanosecond timer */
#define FFTW_TIME_MIN (5.0e-3)

#endif				/* use nanosecond timer */

/**********************************************
 *              WINDOWS
 **********************************************/
#elif defined(HAVE_WIN32_TIMER)

#include <time.h>

typedef unsigned long fftw_time;
extern unsigned long GetPerfTime(void);
extern double GetPerfSec(double ticks);

#define fftw_get_time() GetPerfTime()
#define fftw_time_diff(t1,t2) ((t1) - (t2))
#define fftw_time_to_sec(t) GetPerfSec(t)

#define FFTW_TIME_MIN (5.0e-2)	/* for Win32 timer */

/**********************************************
 *              CRAY
 **********************************************/
#elif defined(_CRAYMPP)		/* Cray MPP system */

double SECONDR(void);		/*
				 * I think you have to link with -lsci to
				 * get this
				 */

typedef double fftw_time;
#define fftw_get_time() SECONDR()
#define fftw_time_diff(t1,t2) ((t1) - (t2))
#define fftw_time_to_sec(t) (t)

#define FFTW_TIME_MIN (1.0e-1)	/* for Cray MPP SECONDR timer */

/**********************************************
 *          VANILLA UNIX/ISO C SYSTEMS
 **********************************************/
/* last resort: use good old Unix clock() */
#else

#include <time.h>

typedef clock_t fftw_time;

#ifndef CLOCKS_PER_SEC
#ifdef sun
/* stupid sunos4 prototypes */
#define CLOCKS_PER_SEC 1000000
extern long clock(void);
#else				/* not sun, we don't know CLOCKS_PER_SEC */
#error Please define CLOCKS_PER_SEC
#endif
#endif

#define fftw_get_time() clock()
#define fftw_time_diff(t1,t2) ((t1) - (t2))
#define fftw_time_to_sec(t) (((double) (t)) / CLOCKS_PER_SEC)

/*
 * ***VERY*** conservative constant: this says that a
 * measurement must run for 200ms in order to be valid.
 * You had better check the manual of your machine
 * to discover if it can do better than this
 */
#define FFTW_TIME_MIN (2.0e-1)	/* for default clock() timer */

#endif				/* UNIX clock() */

/* take FFTW_TIME_REPEAT measurements... */
#ifndef FFTW_TIME_REPEAT
#define FFTW_TIME_REPEAT 4
#endif

/* but do not run for more than TIME_LIMIT seconds while measuring one FFT */
#ifndef FFTW_TIME_LIMIT
#define FFTW_TIME_LIMIT 2.0
#endif

#ifdef __cplusplus
}				/* extern "C" */

#endif				/* __cplusplus */

#endif				/* FFTW_INT_H */

max_stars_repo_name: jlprieur/jlplib

max_issues_repo_name: jlprieur/jlplib

max_forks_repo_name: jlprieur/jlplib


text:
/***************************************************************************
    File                 : muParserScripting.h
    Project              : QtiPlot
    --------------------------------------------------------------------

    Copyright            : (C) 2006 by Ion Vasilief, Knut Franke
    Email (use @ for *)  : ion_vasilief*yahoo.fr, knut.franke*gmx.de
    Description          : Evaluate mathematical expressions using muParser

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

/***************************************************************************
 *                                                                         *
 *  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 MUPARSER_SCRIPTING_H
#define MUPARSER_SCRIPTING_H

#include "ScriptingEnv.h"
#include "Script.h"
#include "muParserScript.h"

#include <muParser.h>
#include "math.h"
#include <gsl/gsl_sf.h>
#include <gsl/gsl_cdf.h>
#include <q3asciidict.h>

//! TODO
class muParserScripting: public ScriptingEnv
{
  Q_OBJECT

  public:
    static const char *langName;
    muParserScripting(ApplicationWindow *parent) : ScriptingEnv(parent, langName) { d_initialized=true; }
    static ScriptingEnv *constructor(ApplicationWindow *parent) { return new muParserScripting(parent); }

    bool isRunning() const { return true; }
    Script *newScript(const QString &code, QObject *context, const QString &name="<input>")
    {
      return new muParserScript(this, code, context, name);
    }

    // we do not support global variables
    bool setQObject(QObject*, const char*) { return false; }
    bool setInt(int, const char*) { return false; }
    bool setDouble(double, const char*) { return false; }

    const QStringList mathFunctions() const;
    const QString mathFunctionDoc (const QString &name) const;

    struct mathFunction
    {
      char *name;
      int numargs;
      double (*fun1)(double);
      double (*fun2)(double,double);
      double (*fun3)(double,double,double);
      char *description;
    };
    static const mathFunction math_functions[];

  private:
    static double mod(double x, double y)
      { return fmod(x,y); }
    static double bessel_J0(double x)
      { return gsl_sf_bessel_J0 (x); }
    static double bessel_J1(double x)
      { return gsl_sf_bessel_J1 (x); }
    static double bessel_Jn(double x, double n)
      { return gsl_sf_bessel_Jn ((int)n, x); }
    static double bessel_Yn(double x, double n)
      { return gsl_sf_bessel_Yn ((int)n, x); }
    static double bessel_Jn_zero(double n, double s)
      { return gsl_sf_bessel_zero_Jnu(n, (unsigned int) s); }
    static double bessel_Y0(double x)
      { return gsl_sf_bessel_Y0 (x); }
    static double bessel_Y1(double x)
      { return gsl_sf_bessel_Y1 (x); }
    static double beta(double a, double b)
      { return gsl_sf_beta (a,b); }
    static double erf(double x)
      { return gsl_sf_erf (x); }
    static double erfc(double x)
      { return gsl_sf_erfc (x); }
    static double erf_Z(double x)
      { return gsl_sf_erf_Z (x); }
    static double erf_Q(double x)
      { return gsl_sf_erf_Q (x); }
    static double gamma(double x)
      { return gsl_sf_gamma (x); }
    static double lngamma(double x)
      { return gsl_sf_lngamma (x); }
    static double hazard(double x)
      { return gsl_sf_hazard (x); }
	static double lambert_W0(double x)
	   { return gsl_sf_lambert_W0(x); }
	static double lambert_Wm1(double x)
	   { return gsl_sf_lambert_Wm1(x); }
	static double ttable(double x, double n)
	   { return gsl_cdf_tdist_Pinv(x, n); }
};

class EmptySourceError : public mu::ParserError
{
	public:
		EmptySourceError() {}
};

#endif

max_stars_repo_name: hoehnp/SpaceDesignTool

max_issues_repo_name: hoehnp/SpaceDesignTool

max_forks_repo_name: hoehnp/SpaceDesignTool


text:
//
//  parameter.h
//  EpiGenMCMC
//
//  Created by Lucy Li on 08/04/2016.
//  Copyright (c) 2016 Lucy Li, Imperial College London. All rights reserved.
//

#ifndef EpiGenMCMC_parameter_h
#define EpiGenMCMC_parameter_h

#include <fstream>
#include <sstream>
#include <vector>
#include <string>
#include <algorithm>
#include <cmath>
#include <limits>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_cdf.h>

class Parameter {
    std::vector <double> parameter_values;
    std::vector <std::string> parameter_names;
    std::vector <bool> estimate;
    std::vector <std::string> transform;
    std::vector <std::string> prior;
    std::vector <double> prior_par_1;
    std::vector <double> prior_par_2;
    std::vector <double> prior_par_3;
    std::vector <std::string> proposal;
    std::vector <double> proposal_sd;
    std::vector <double> proposal_lower;
    std::vector <double> proposal_upper;
    std::vector <double> accepted;
    std::vector <double> rejected;
    std::vector <double> acceptance_rate;
    std::vector <int> params_to_estim;
    int total_params;
    int curr_param_to_estimate;
    double old_param_value;
    double optimal_acceptance;
    double lower_acceptance;
    double upper_acceptance;
    int adapt_every;
    int max_adapt_times;
    bool stop_adapting;
    double get_prior(double, int) const;
    double get_transform(double, std::string, bool) const;
public:
    Parameter();
    Parameter(const Parameter &);
    Parameter (std::string);
    bool is_estim(int) const;
    int get_curr_estim() const;
    void set_next_param();
    int get_total_params() const;
    int get_total_params_to_estim() const;
    int get_estim_index(int) const;
    double get(std::string) const;
    double get(int) const;
    double get_lower(int) const;
    double get_upper(int) const;
    void set(int, double);
    std::vector<double> get_values_vector() const;
    std::vector <std::string> get_names_vector() const;
    std::string getname(int) const;
    void reset();
    void accept();
    void reject();
    double get_acceptance() const;
    void stop_adapt();
    void start_adapt();
    void adapt();
    double propose(gsl_rng *);
    double get_prior_ratio() const;
    double get_prior_all() const;
    bool param_exists(std::string) const;
    void transform_param(int, bool);
    double get_lognormal_sd (double, double) const;
    double get_lognormal_mean (double, double) const;
};

struct MCMCoptions {
    int particles;
    int iterations;
    int log_every;
    int pfilter_every;
    double pfilter_threshold;
    int which_likelihood;
    int num_trees;
    std::string log_filename;
    std::string traj_filename;
    std::string model;
    int total_dt;
    double sim_dt;
    int num_groups;
    int seed;
    bool verbose;
    bool save_traj;
    bool use_lhs;
    int lhs_divides;
    int lhs_iterations;
    int num_threads;
    double heat_factor;
    int heat_length;
    double cool_rate;
    gsl_rng ** rng;
    MCMCoptions();
    MCMCoptions(const MCMCoptions&);
    MCMCoptions(std::string);

};

#endif

max_stars_repo_name: lucymli/EpiGenMCMC

max_issues_repo_name: lucymli/EpiGenMCMC

max_forks_repo_name: lucymli/EpiGenMCMC