text:
//
//  emu.c
//  Created by Earl Lawrence on 11/10/16.
//  Modified E Chisari 05/10/17 for CCL
//  For details on the license, see ../LICENSE_COSMICEMU
//  in this repository.
//

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

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_errno.h>

#include "ccl.h"

#include "ccl_emu17_params.h"
#include "ccl_emu17.h"

// Sizes of stuff
static int m[2] = {111, 36}, neta=2808, peta[2]={7, 28}, rs=8, p=8;

// Kriging basis computed by emuInit
// Sizes of each basis will be peta[ee] and m[ee]
static double KrigBasis[2][28][111];

// Need to combine some things into a big thing.
static double beta[2][28][8];
static double w[2][28][111];
static double lamws[2][28];
static double lamz[2][28];

// Initialization to compute the kriging basis for the two parts
static void emuInit() {

    int ee, i, j, k, l;
    double cov;
    gsl_matrix *SigmaSim;
    gsl_vector *b;

    // Because of the structure of this emu, I need to do this horrible part first.
    // Fill in the correlation lenghths
    for(i=0; i<7; i++) {
        for(j=0; j<8; j++) {
            beta[0][i][j] = beta1[i][j];
        }
    }

    for(i=0; i<28; i++) {
        for(j=0; j<8; j++) {
            beta[1][i][j] = beta2[i][j];
        }
    }

    // Fill in the PC weights
    for(i=0; i<7; i++) {
        for(j=0; j<111; j++) {
            w[0][i][j] = w1[i][j];
        }
    }

    for(i=0; i<28; i++) {
        for(j=0; j<36; j++) {
            w[1][i][j] = w2[i][j];
        }
    }

    // Fill in the precisions
    for(i=0; i<7; i++) {
        lamws[0][i] = lamws1[i];
        lamz[0][i] = lamz1[i];
    }

    for(i=0; i<28; i++) {
        lamws[1][i] = lamws2[i];
        lamz[1][i] = lamz2[i];
    }

    // This emu has two parts: one that uses all m[0] of the data for the the first peta[0] components
    // and another that uses only the m[1] complete data for the next peta[1] components.
    for(ee=0; ee<2; ee++) {

        // Allocate some stuff
        SigmaSim = gsl_matrix_alloc(m[ee], m[ee]);
        b = gsl_vector_alloc(m[ee]);

        // Loop over the basis
        for(i=0; i<peta[ee]; i++) {

            // Loop over the number of simulations
            for(j=0; j<m[ee]; j++) {

                // Diagonal term
                gsl_matrix_set(SigmaSim, j, j, (1.0/lamz[ee][i]) + (1.0/lamws[ee][i]));

                // Off-diagonals
                for(k=0; k<j; k++) {

                    // compute the covariance
                    cov = 0.0;
                    for(l=0; l<p; l++) {
                        cov -= beta[ee][i][l]*pow(x[j][l] - x[k][l], 2.0);
                    } // for(l=0; l<p; l++)
                    cov = exp(cov) / lamz[ee][i];

                    // put the covariance where it belongs
                    gsl_matrix_set(SigmaSim, j, k, cov);
                    gsl_matrix_set(SigmaSim, k, j, cov);

                } // for(k=0; k<j; k++)

                // Vector for the PC weights
                gsl_vector_set(b, j, w[ee][i][j]);

            } // for(j=0; j<m[ee]; j++)

            // Cholesky and solve
            gsl_linalg_cholesky_decomp(SigmaSim);
            gsl_linalg_cholesky_svx(SigmaSim, b);

            // Put b where it belongs in the Kriging basis
            for(j=0; j<m[ee]; j++) {
                KrigBasis[ee][i][j] = gsl_vector_get(b, j);
            }

        } // for(i=0; i<peta[ee]; i++)

        // Clean this up
        gsl_matrix_free(SigmaSim);
        gsl_vector_free(b);

    } // for(ee=0; ee<2; ee+)

} // emuInit()

// Actual emulation
void ccl_pkemu(double *xstar, double **ystar, int* status, ccl_cosmology* cosmo) {

    static double inited=0;
    int ee, i, j, k;
    double wstar[peta[0]+peta[1]];
    double Sigmastar[2][peta[1]][m[0]];
    double ystaremu[neta];
    *ystar=(double *)malloc(sizeof(double)*NK_EMU);
    double ybyz[rs];
    double logc;
    double xstarstd[p];
    int zmatch;
    gsl_spline *zinterp = gsl_spline_alloc(gsl_interp_cspline, rs);
    gsl_interp_accel *accel = gsl_interp_accel_alloc();

    // Initialize if necessary
    if(inited==0) {
        emuInit();
        inited = 1;
    }

    // Transform w_a into (-w_0-w_a)^(1/4)
    xstar[6] = pow(-xstar[5]-xstar[6], 0.25);
    // Check the inputs to make sure we're interpolating.
    for(i=0; i<p; i++) {
        if((xstar[i] < xmin[i]) || (xstar[i] > xmax[i])) {
            switch(i) {
                case 0:
                    ccl_cosmology_set_status_message(cosmo,
                            "ccl_pkemu(): omega_m must be between %f and %f.\n",
                            xmin[i], xmax[i]);
                    break;
                case 1:
                    ccl_cosmology_set_status_message(cosmo,
                            "ccl_pkemu(): omega_b must be between %f and %f.\n",
                            xmin[i], xmax[i]);
                    break;
                case 2:
                    ccl_cosmology_set_status_message(cosmo,
                            "ccl_pkemu(): sigma8 must be between %f and %f.\n",
                            xmin[i], xmax[i]);
                    break;
                case 3:
                    ccl_cosmology_set_status_message(cosmo,
                            "ccl_pkemu(): h must be between %f and %f.\n",
                            xmin[i], xmax[i]);
                    break;
                case 4:
                    ccl_cosmology_set_status_message(cosmo,
                            "ccl_pkemu(): n_s must be between %f and %f.\n",
                            xmin[i], xmax[i]);
                    break;
                case 5:
                    ccl_cosmology_set_status_message(cosmo,
                            "ccl_pkemu(): w_0 must be between %f and %f.\n",
                            xmin[i], xmax[i]);
                    break;
                case 6:
                    ccl_cosmology_set_status_message(cosmo,
                            "ccl_pkemu(): (-w_0-w_a)^(1/4) must be between %f and %f.\n",
                            xmin[i], xmax[i]);
                    break;
                case 7:
                    ccl_cosmology_set_status_message(cosmo,
                            "ccl_pkemu(): omega_nu must be between %f and %f.\n",
                            xmin[i], xmax[i]);
                    break;
            }
            *status = CCL_ERROR_EMULATOR_BOUND;
            ccl_raise_exception(*status, cosmo->status_message);
	    gsl_spline_free(zinterp);
	    gsl_interp_accel_free(accel);
            return;
        }
    } // for(i=0; i<p; i++)
    if((xstar[p] < z[0]) || (xstar[p] > z[rs-1])) {
        ccl_cosmology_set_status_message(cosmo,
                "ccl_pkemu(): z must be between %f and %f.\n",
                z[0], z[rs-1]);
        *status = CCL_ERROR_EMULATOR_BOUND;
        ccl_raise_exception(*status, cosmo->status_message);
	gsl_spline_free(zinterp);
	gsl_interp_accel_free(accel);
        return;
    }

    // Standardize the inputs
    for(i=0; i<p; i++) {
        xstarstd[i] = (xstar[i] - xmin[i]) / xrange[i];
    }

    // compute the covariances between the new input and sims for all the PCs.
    for(ee=0; ee<2; ee++) {
        for(i=0; i<peta[ee]; i++) {
            for(j=0; j<m[ee]; j++) {
                logc = 0.0;
                for(k=0; k<p; k++) {
                    logc -= beta[ee][i][k]*pow(x[j][k] - xstarstd[k], 2.0);
                }
                Sigmastar[ee][i][j] = exp(logc) / lamz[ee][i];
            }
        }
    }

    // Compute wstar for the first chunk.
    for(i=0; i<peta[0]; i++) {
        wstar[i]=0.0;
        for(j=0; j<m[0]; j++) {
            wstar[i] += Sigmastar[0][i][j] * KrigBasis[0][i][j];
        }
    }
    // Compute wstar for the second chunk.
    for(i=0; i<peta[1]; i++) {
        wstar[peta[0]+i]=0.0;
        for(j=0; j<m[1]; j++) {
            wstar[peta[0]+i] += Sigmastar[1][i][j] * KrigBasis[1][i][j];
        }
    }

    /*
    for(i=0; i<peta[0]+peta[1]; i++) {
        printf("%f\n", wstar[i]);
    }
     */

    // Compute ystar, the new output
    for(i=0; i<neta; i++) {
        ystaremu[i] = 0.0;
        for(j=0; j<peta[0]+peta[1]; j++) {
            ystaremu[i] += K[i][j]*wstar[j];
        }
        ystaremu[i] = ystaremu[i]*sd + mean[i];

        // Convert to P(k)
        //ystaremu[i] = ystaremu[i] - 1.5*log10(mode[i % NK_EMU]);
        //ystaremu[i] = 2*M_PI*M_PI*pow(10, ystaremu[i]);
    }

    // Interpolate to the desired redshift
    // Natural cubic spline interpolation over z.

    // First check to see if the requested z is one of the training z.
    zmatch = -1;
    for(i=0; i<rs; i++) {
        if(xstar[p] == z[i]) {
            zmatch = rs-i-1;
        }
    }

    // z doesn't match a training z, interpolate
    if(zmatch == -1) {
        for(i=0; i<NK_EMU; i++) {
            for(j=0; j<rs; j++) {
                ybyz[rs-j-1] = ystaremu[j*NK_EMU+i];
            }
            gsl_spline_init(zinterp, z, ybyz, rs);
            (*ystar)[i] = gsl_spline_eval(zinterp, xstar[p], accel);
            gsl_interp_accel_reset(accel);
        }

    } else { //otherwise, copy in the emulated z without interpolating
        for(i=0; i<NK_EMU; i++) {
	  (*ystar)[i] = ystaremu[zmatch*NK_EMU + i];
        }
    }

    gsl_spline_free(zinterp);
    gsl_interp_accel_free(accel);

    // Convert to P(k)
    for(i=0; i<NK_EMU; i++) {
      (*ystar)[i] = (*ystar)[i] - 1.5*log10(mode[i]) + log10(2) + 2*log10(M_PI);
      (*ystar)[i] = pow(10, (*ystar)[i]);
    }
}

max_stars_repo_name: Russell-Jones-OxPhys/CCL

max_issues_repo_name: Russell-Jones-OxPhys/CCL

max_forks_repo_name: Russell-Jones-OxPhys/CCL


text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <stdarg.h>
#include <gbpLib.h>
#include <gbpMath.h>
#include <gbpCosmo_mass_functions.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>

// Scaled mass functions
double scaled_mass_function(double sigma, int mode, double *P) {
    double rval;
    if(SID_CHECK_BITFIELD_SWITCH(mode, MF_WATSON)) { // Watson et al. (2013) universal FoF mass function
        double A     = 0.282;
        double alpha = 2.163;
        double beta  = 1.406;
        double gamma = 1.210;
        if(SID_CHECK_BITFIELD_SWITCH(mode, MF_PASS_PARAMS)) {
            A     = P[0];
            alpha = P[1];
            beta  = P[2];
            gamma = P[3];
        }
        rval = A * (pow(beta / sigma, alpha) + 1.) * exp(-gamma / (sigma * sigma));
    } else if(SID_CHECK_BITFIELD_SWITCH(mode, MF_TIAMAT)) { // Poole et al. (2013) universal FoF mass function for Tiamat
        double A     = 0.03331;
        double alpha = 1.153;
        double beta  = 12.33;
        double gamma = 1.009;
        if(SID_CHECK_BITFIELD_SWITCH(mode, MF_PASS_PARAMS)) {
            A     = P[0];
            alpha = P[1];
            beta  = P[2];
            gamma = P[3];
        }
        rval = A * (pow(beta / sigma, alpha) + 1.) * exp(-gamma / (sigma * sigma));
    } else if(SID_CHECK_BITFIELD_SWITCH(mode, MF_JENKINS)) { // Jenkins et al. (2001)
        double A = 0.315;
        double B = 0.61;
        double C = 3.8;
        if(SID_CHECK_BITFIELD_SWITCH(mode, MF_PASS_PARAMS)) {
            A = P[0];
            B = P[1];
            C = P[2];
        }
        rval = A * exp(-pow((double)fabs((float)(take_ln(1. / sigma) + B)), C));
    } else if(SID_CHECK_BITFIELD_SWITCH(mode, MF_ST)) { // Sheth-Torman (1999)
        double delta_k = 1.686;
        double A       = 0.3222;
        double B       = 0.707;
        double C       = 0.3;
        if(SID_CHECK_BITFIELD_SWITCH(mode, MF_PASS_PARAMS)) {
            A = P[0];
            B = P[1];
            C = P[2];
        }
        rval = A * sqrt(2. * B / PI) * (delta_k / sigma) * exp(-B * delta_k * delta_k / (2. * sigma * sigma)) *
               (1. + pow(sigma * sigma / (B * delta_k * delta_k), C));
    } else if(SID_CHECK_BITFIELD_SWITCH(mode, MF_PS)) { // Press-Schechter (1974)
        double delta_k = 1.686;
        if(SID_CHECK_BITFIELD_SWITCH(mode, MF_PASS_PARAMS)) {
            delta_k = P[0];
        }
        rval = sqrt(2. / PI) * (delta_k / sigma) * exp(-delta_k * delta_k / (2. * sigma * sigma));
    } else
        SID_exit_error("A valid mass function was not specified with mode (%d) in scaled_mass_function().\n",
                       SID_ERROR_LOGIC, mode);
    return (rval);
}

ext: c

max_stars_repo_name: gbpoole/gbpCode

max_issues_repo_name: gbpoole/gbpCode

max_forks_repo_name: gbpoole/gbpCode


text:
#ifndef UTIL_H
#define UTIL_H

#include <vector>

#include <gsl/gsl_matrix.h>

class Matrices_keeper {

    private:
        std::vector<gsl_matrix*> matrices;
    public:
        Matrices_keeper();
        Matrices_keeper(const Matrices_keeper& keeper);
        Matrices_keeper& operator=(const Matrices_keeper& keeper);

        gsl_matrix* alloc_next(int s1, int s2);
        void dalloc_all();
        ~Matrices_keeper();

};

bool validate_inputs(char **args);
void usage();

void exit_on_status(int status, char *msg);
void exit_on_null(void *p, char *msg);

#endif // UTIL_H

ext: h

max_stars_repo_name: bainit/hextexsim

max_issues_repo_name: bainit/hextexsim

max_forks_repo_name: bainit/hextexsim


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

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_rng.h>

/* This is the TT800 twisted GSFR generator for 32 bit integers. It
   has been superceded by MT19937 (mt.c). The period is 2^800.

   This implementation is based on tt800.c, July 8th 1996 version by
   M. Matsumoto, email: matumoto@math.keio.ac.jp

   From: Makoto Matsumoto and Yoshiharu Kurita, "Twisted GFSR
   Generators II", ACM Transactions on Modelling and Computer
   Simulation, Vol. 4, No. 3, 1994, pages 254-266. */

static inline unsigned long int tt_get (void *vstate);
static double tt_get_double (void *vstate);
static void tt_set (void *state, unsigned long int s);

#define N 25
#define M 7

typedef struct
  {
    int n;
    unsigned long int x[N];
  }
tt_state_t;

static inline unsigned long int
tt_get (void *vstate)
{
  tt_state_t *state = (tt_state_t *) vstate;

  /* this is the magic vector, a */

  const unsigned long mag01[2] =
  {0x00000000, 0x8ebfd028UL};
  unsigned long int y;
  unsigned long int *const x = state->x;
  int n = state->n;

  if (n >= N)
    {
      int i;
      for (i = 0; i < N - M; i++)
	{
	  x[i] = x[i + M] ^ (x[i] >> 1) ^ mag01[x[i] % 2];
	}
      for (; i < N; i++)
	{
	  x[i] = x[i + (M - N)] ^ (x[i] >> 1) ^ mag01[x[i] % 2];
	};
      n = 0;
    }

  y = x[n];
  y ^= (y << 7) & 0x2b5b2500UL;		/* s and b, magic vectors */
  y ^= (y << 15) & 0xdb8b0000UL;	/* t and c, magic vectors */
  y &= 0xffffffffUL;	/* you may delete this line if word size = 32 */

  /* The following line was added by Makoto Matsumoto in the 1996
     version to improve lower bit's correlation.  Delete this line
     to use the code published in 1994.  */

  y ^= (y >> 16);	/* added to the 1994 version */

  state->n = n + 1;

  return y;
}

static double
tt_get_double (void * vstate)
{
  return tt_get (vstate) / 4294967296.0 ;
}

static void
tt_set (void *vstate, unsigned long int s)
{
  tt_state_t *state = (tt_state_t *) vstate;

  const tt_state_t init_state =
  {0,
   {0x95f24dabUL, 0x0b685215UL, 0xe76ccae7UL,
    0xaf3ec239UL, 0x715fad23UL, 0x24a590adUL,
    0x69e4b5efUL, 0xbf456141UL, 0x96bc1b7bUL,
    0xa7bdf825UL, 0xc1de75b7UL, 0x8858a9c9UL,
    0x2da87693UL, 0xb657f9ddUL, 0xffdc8a9fUL,
    0x8121da71UL, 0x8b823ecbUL, 0x885d05f5UL,
    0x4e20cd47UL, 0x5a9ad5d9UL, 0x512c0c03UL,
    0xea857ccdUL, 0x4cc1d30fUL, 0x8891a8a1UL,
    0xa6b7aadbUL}};

  if (s == 0)	/* default seed is given explicitly in the original code */
    {
      *state = init_state;
    }
  else
    {
      int i;

      state->n = 0;

      state->x[0] = s & 0xffffffffUL;

      for (i = 1; i < N; i++)
	state->x[i] = (69069 * state->x[i - 1]) & 0xffffffffUL;
    }

  return;
}

static const gsl_rng_type tt_type =
{"tt800",			/* name */
 0xffffffffUL,			/* RAND_MAX */
 0,			        /* RAND_MIN */
 sizeof (tt_state_t),
 &tt_set,
 &tt_get,
 &tt_get_double};

const gsl_rng_type *gsl_rng_tt800 = &tt_type;


text:
#include <stdio.h>
#include <gsl/gsl_blas.h>
#include "../../include/timing.h"

double* default_matrix;
double** results;

int max(int a, int b){
  if(a > b) return a;
  return b;
}

void save_res_to_csv(double** results, int max_size, int stride, int trials, char* file_name){
  FILE* csv = fopen(file_name, "w");
  if(csv == NULL) perror("Error creating file!");

  // header of the .csv
  fprintf(csv, "size,result\n" );
  // measurement results
  for(int i = 0; i < max_size; i++){
    for (int j = 0; j < trials; j++) {
      fprintf(csv, "%d,%lf\n", i * stride, results[i][j]);
    }
  }
  fclose(csv);

}

void init_matrix(int max_size, int stride){
  int actual_size = max_size * stride;
  default_matrix = (double*) malloc(actual_size* actual_size * sizeof(double));
  for (int i = 0; i < actual_size; i++) {
    for (int j = 0; j < actual_size; j++) {
      default_matrix[actual_size * i + j] = actual_size * i + j;
    }
  }
}

void init_results(int max_size, int trials){
  results = (double**) malloc(max_size * sizeof(double*));
  for (int i = 0; i < max_size; i++) {
    results[i] = (double*) malloc(trials * sizeof(double));
  }
}

void cleanup_matrix(int max_size, int stride){
  free(default_matrix);
}

void cleanup_results(int max_size){
  for (int i = 0; i < max_size; i++) {
    free(results[i]);
  }
  free(results);
}

double matrix_naive_test(int size){
  gsl_matrix* new_mat = gsl_matrix_alloc(size, size);
  gsl_matrix_view tmp_m = gsl_matrix_view_array(default_matrix, size, size);
  double res;
  for (int i = 0; i < size; i++) {
    for (int j = 0; j < size; j++) {
      gsl_matrix_set(new_mat, i, j, 0);
    }
  }
  double start = get_real_time();
  for(int j = 0; j< size; j++){
    for (int k = 0; k < size; k++) {
      for (int i = 0; i < size; i++) {
        double tmp_res = gsl_matrix_get(&tmp_m.matrix, i, k) * gsl_matrix_get(&tmp_m.matrix, k, j);
        gsl_matrix_set(new_mat, i, j, gsl_matrix_get(new_mat, i, j) + tmp_res);

      }
    }
  }
  double end = get_real_time();
  gsl_matrix_free(new_mat);
  return end - start;

}

double matrix_better_test(int size){
  gsl_matrix* new_mat = gsl_matrix_alloc(size, size);
  gsl_matrix_view tmp_m = gsl_matrix_view_array(default_matrix, size, size);
  double res;
  double start = get_real_time();
  for(int i = 0; i< size; i++){
    for (int j = 0; j < size; j++) {
      double tmp_res = 0;
      for (int k = 0; k < size; k++) {
        tmp_res += gsl_matrix_get(&tmp_m.matrix, i, k) * gsl_matrix_get(&tmp_m.matrix, k, j);
      }
      gsl_matrix_set(new_mat, i, j, tmp_res);
    }
  }
  double end = get_real_time();
  gsl_matrix_free(new_mat);
  return end - start;

}

double matrix_blas_test(int size){
  gsl_matrix_view tmp_m = gsl_matrix_view_array(default_matrix, size, size);
  double res;
  double start = get_real_time();
  gsl_blas_dsymm(CblasLeft, CblasUpper, 1, &tmp_m.matrix, &tmp_m.matrix, 0,&tmp_m.matrix);

  double end = get_real_time();
  return end - start;

}

void matrix_naive_analysis(int max_size, int stride, int trials){
  init_matrix(max_size, stride);
  for(int i = 1; i < max_size; i++){
    for (int j = 0; j < trials; j++) {
      results[i][j] = matrix_naive_test(i * stride);
    }
  }
  save_res_to_csv(results, max_size, stride, trials, "naive.csv");
  cleanup_matrix(max_size, stride);
}

void matrix_better_analysis(int max_size, int stride, int trials){
  init_matrix(max_size, stride);
  for(int i = 1; i < max_size; i++){
    for (int j = 0; j < trials; j++) {
      results[i][j] = matrix_better_test(i * stride);
    }
  }
  save_res_to_csv(results, max_size, stride, trials, "better.csv");
  cleanup_matrix(max_size, stride);
}

void matrix_blas_analysis(int max_size, int stride, int trials){
  init_matrix(max_size, stride);
  for(int i = 1; i < max_size; i++){
    for (int j = 0; j < trials; j++) {
      results[i][j] = matrix_blas_test(i * stride);
    }
  }
  save_res_to_csv(results, max_size, stride, trials, "blas.csv");
  cleanup_matrix(max_size, stride);
}

int
main (void)
{
  int max_size = 25;
  int stride = 20;
  int trials = 10;
  int results_size = max_size;
  init_results(results_size, trials);
  printf("Analyzing naive multiplication.\n" );
  matrix_naive_analysis(max_size, stride, trials);
  printf("Analyzing better multiplication.\n" );
  matrix_better_analysis(max_size, stride, trials);
  printf("Analyzing blas multiplication.\n" );
  matrix_blas_analysis(max_size, stride, trials);
  cleanup_results(results_size);

  return 0;
}


text:
/* histogram/ntuple.h
 *
 * Copyright (C) 2000 Simone Piccardi
 *
 * 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.
 *
 */

/* Jan/2001 Modified by Brian Gough. Minor changes for GSL */

#ifndef __GSL_NTUPLE_H__
#define __GSL_NTUPLE_H__

#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_histogram.h>

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

__BEGIN_DECLS

typedef struct {
    FILE * file;
    void * ntuple_data;
    size_t size;
} gsl_ntuple;

typedef struct {
  int (* function) (void * ntuple_data, void * params);
  void * params;
} gsl_ntuple_select_fn;

typedef struct {
  double (* function) (void * ntuple_data, void * params);
  void * params;
} gsl_ntuple_value_fn;

gsl_ntuple *
gsl_ntuple_open (char * filename, void * ntuple_data, size_t size);

gsl_ntuple *
gsl_ntuple_create (char * filename, void * ntuple_data, size_t size);

int gsl_ntuple_write (gsl_ntuple * ntuple);
int gsl_ntuple_read (gsl_ntuple * ntuple);

int gsl_ntuple_bookdata (gsl_ntuple * ntuple);  /* synonym for write */

int gsl_ntuple_project (gsl_histogram * h, gsl_ntuple * ntuple,
                        gsl_ntuple_value_fn *value_func,
                        gsl_ntuple_select_fn *select_func);

int gsl_ntuple_close (gsl_ntuple * ntuple);

__END_DECLS

#endif /* __GSL_NTUPLE_H__ */

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


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

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

/* The binomial distribution has the form,

   prob(k) =  n!/(k!(n-k)!) *  p^k (1-p)^(n-k) for k = 0, 1, ..., n

   This is the algorithm from Knuth */

unsigned int
gsl_ran_binomial (const gsl_rng * r, double p, unsigned int n)
{
  unsigned int i, a, b, k = 0;

  while (n > 10)	/* This parameter is tunable */
    {
      double X;
      a = 1 + (n / 2);
      b = 1 + n - a;

      X = gsl_ran_beta (r, (double) a, (double) b);

      if (X >= p)
	{
	  n = a - 1;
	  p /= X;
	}
      else
	{
	  k += a;
	  n = b - 1;
	  p = (p - X) / (1 - X);
	}
    }

  for (i = 0; i < n; i++)
    {
      double u = gsl_rng_uniform (r);
      if (u < p)
	k++;
    }

  return k;
}

double
gsl_ran_binomial_pdf (const unsigned int k, const double p,
		      const unsigned int n)
{
  if (k > n)
    {
      return 0 ;
    }
  else
    {
      double a = k;
      double b = n - k;
      double P;
      double Cnk = gsl_sf_choose (n, k) ;

      P = Cnk * pow (p, a) * pow (1 - p, b);

      return P;
    }
}


text:
/**
 *
 * @file qwrapper_dlaset2.c
 *
 *  PLASMA core_blas quark wrapper
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Hatem Ltaief
 * @date 2010-11-15
 * @generated d Tue Jan  7 11:44:59 2014
 *
 **/
#include <lapacke.h>
#include "common.h"

/***************************************************************************//**
 *
 **/
void QUARK_CORE_dlaset2(Quark *quark, Quark_Task_Flags *task_flags,
                       PLASMA_enum uplo, int M, int N,
                       double alpha, double *A, int LDA)
{
    DAG_CORE_LASET;
    QUARK_Insert_Task(quark, CORE_dlaset2_quark, task_flags,
        sizeof(PLASMA_enum),                &uplo,  VALUE,
        sizeof(int),                        &M,     VALUE,
        sizeof(int),                        &N,     VALUE,
        sizeof(double),         &alpha, VALUE,
        sizeof(double)*M*N,     A,      OUTPUT,
        sizeof(int),                        &LDA,   VALUE,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_dlaset2_quark = PCORE_dlaset2_quark
#define CORE_dlaset2_quark PCORE_dlaset2_quark
#endif
void CORE_dlaset2_quark(Quark *quark)
{
    PLASMA_enum uplo;
    int M;
    int N;
    double alpha;
    double *A;
    int LDA;

    quark_unpack_args_6(quark, uplo, M, N, alpha, A, LDA);
    CORE_dlaset2(uplo, M, N, alpha, A, LDA);
}


text:
//
//  syr2k.h
//  Linear Algebra Template Library
//
//  Created by Rodney James on 1/4/12.
//  Copyright (c) 2012 University of Colorado Denver. All rights reserved.
//

#ifndef _syr2k_h
#define _syr2k_h

/// @file syr2k.h Performs multiplication of two matrices resulting in a symmetric matrix.

#include <cctype>
#include "latl.h"

namespace LATL
{
   /// @brief Performs multiplcation of real matrices resulting in a symmetric matrix.
   ///
   /// For real matrices A and B, real symmetric matrix C, and real scalars alpha and beta
   ///
   ///        C := alpha*A*B'+alpha*B*A'+beta*C
   /// or
   ///
   ///        C := alpha*A'*B+alpha*B'*A+beta*C
   /// is computed.
   /// @return 0 if success.
   /// @return -i if the ith argument is invalid.
   /// @tparam real_t Floating point type.
   /// @param uplo Specifies whether the upper or lower triangular part of the symmetric matrix C
   /// is to be referenced:
   ///
   ///        if uplo = 'U' or 'u' then C is upper triangular,
   ///        if uplo = 'L' or 'l' then C is lower triangular.
   /// @param trans Specifies the operation to be perfomed as follows:
   ///
   ///        if trans = 'N' or 'n' then C := alpha*A*B'+alpha*B*A'+beta*C
   ///        if trans = 'T' or 't' then C := alpha*A'*B+alpha*B'*A+beta*C
   ///        if trans = 'C' or 'c' then C := alpha*A'*B+alpha*B'*A+beta*C
   /// @param n Specifies the order of the complex Hermitian matrix C.  n>=0
   /// @param k Specifies the other dimension of the complex matrices A and B.
   ///
   ///        if trans = 'N' or 'n' then A and B are n-by-k
   ///        if trans = 'T' or 't' then A and B are k-by-n
   ///        if trans = 'C' or 'c' then A and B are k-by-n
   /// @param alpha Real scalar.
   /// @param A Pointer to real matrix.
   /// @param ldA Column length of the matrix A. If trans = 'N' or 'n' ldA>=n, otherwise ldA>=k.
   /// @param B Pointer to real matrix.
   /// @param ldB Column length of the matrix B. If trans = 'N' or 'n' ldB>=n, otherwise ldB>=k.
   /// @param beta Real scalar.
   /// @param C Pointer to real symmetric n-by-n matrix C.
   /// Only the upper or lower triangular part of C is referenced, depending on the value of uplo above.
   /// @param ldC Column length of the matrix C.  ldC>=n
   /// @ingroup BLAS

   template <typename real_t>
   int SYR2K(char uplo, char trans, int_t n, int_t k, real_t alpha, real_t *A, int_t ldA, real_t *B, int_t ldB, real_t beta, real_t *C, int_t ldC)
   {
      using std::toupper;

      const real_t zero(0.0);
      const real_t one(1.0);
      int_t i,j,l;
      real_t *a,*b,*c,*at,*bt;
      real_t s,t;

      uplo=toupper(uplo);
      trans=toupper(trans);

      if((uplo!='U')&&(uplo!='L'))
         return -1;
      else if((trans!='N')&&(trans!='T')&&(trans!='C'))
         return -2;
      else if(n<0)
         return -3;
      else if(k<0)
         return -4;
      else if(ldA<((trans=='N')?n:k))
         return -7;
      else if(ldB<((trans=='N')?n:k))
         return -9;
      else if(ldC<n)
         return -12;
      else if((n==0)||(((alpha==zero)||(k==0))&&(beta==one)))
         return 0;

      if(alpha==zero)
      {
         if(uplo=='U')
         {
            c=C;
            for(j=0;j<n;j++)
            {
               for(i=0;i<=j;i++)
                  c[i]*=beta;
               c+=ldC;
            }

         }
         else
         {
            c=C;
            for(j=0;j<n;j++)
            {
               for(i=j;i<n;i++)
                  c[i]*=beta;
               c+=ldC;
            }
         }
      }
      else if(trans=='N')
      {
         if(uplo=='U')
         {
            c=C;
            for(j=0;j<n;j++)
            {
               for(i=0;i<=j;i++)
                  c[i]*=beta;
               a=A;
               b=B;
               for(l=0;l<k;l++)
               {
                  s=alpha*b[j];
                  t=alpha*a[j];
                  for(i=0;i<=j;i++)
                     c[i]+=s*a[i]+t*b[i];
                  a+=ldA;
                  b+=ldB;
               }
               c+=ldC;
            }
         }
         else
         {
            c=C;
            for(j=0;j<n;j++)
            {
               for(i=j;i<n;i++)
                  c[i]*=beta;
               a=A;
               b=B;
               for(l=0;l<k;l++)
               {
                  s=alpha*b[j];
                  t=alpha*a[j];
                  for(i=j;i<n;i++)
                     c[i]+=s*a[i]+t*b[i];
                  a+=ldA;
                  b+=ldB;
               }
               c+=ldC;
            }
         }
      }
      else
      {
         if(uplo=='U')
         {
            c=C;
            at=A;
            bt=B;
            for(j=0;j<n;j++)
            {
               a=A;
               b=B;
               for(i=0;i<=j;i++)
               {
                  s=zero;
                  t=zero;
                  for(l=0;l<k;l++)
                  {
                     s+=b[l]*at[l];
                     t+=a[l]*bt[l];
                  }
                  c[i]=alpha*t+alpha*s+beta*c[i];
                  a+=ldA;
                  b+=ldB;
               }
               at+=ldA;
               bt+=ldB;
               c+=ldC;
            }
         }
         else
         {
            c=C;
            at=A;
            bt=B;
            for(j=0;j<n;j++)
            {
               a=A+j*ldA;
               b=B+j*ldB;
               for(i=j;i<n;i++)
               {
                  s=zero;
                  t=zero;
                  for(l=0;l<k;l++)
                  {
                     s+=b[l]*at[l];
                     t+=a[l]*bt[l];
                  }
                  c[i]=alpha*t+alpha*s+beta*c[i];
                  a+=ldA;
                  b+=ldB;
               }
               at+=ldA;
               bt+=ldB;
               c+=ldC;
            }
         }
      }
      return 0;
   }

   /// @brief Performs multiplcation of complex matrices resulting in a symmetric matrix.
   ///
   /// For complex matrices A and B, complex Hermitian matrix C, and complex scalars alpha and beta
   ///
   ///        C := alpha*A*B.'+alpha*B*A.'+beta*C
   /// or
   ///
   ///        C := alpha*A.'*B+alpha*B.'*A+beta*C
   /// is computed.
   /// @return 0 if success.
   /// @return -i if the ith argument is invalid.
   /// @tparam real_t Floating point type.
   /// @param uplo Specifies whether the upper or lower triangular part of the symmetric matrix C
   /// is to be referenced:
   ///
   ///        if uplo = 'U' or 'u' then C is upper triangular,
   ///        if uplo = 'L' or 'l' then C is lower triangular.
   /// @param trans Specifies the operation to be perfomed as follows:
   ///
   ///        if trans = 'N' or 'n' then C := alpha*A*B.'+alpha*B*A.'+beta*C
   ///        if trans = 'T' or 'T' then C := alpha*A.'*B+alpha*B.'*A+beta*C
   /// @param n Specifies the order of the complex symmetric matrix C.  n>=0
   /// @param k Specifies the other dimension of the complex matrices A and B.
   ///
   ///        if trans = 'N' or 'n' then A and B are n-by-k
   ///        if trans = 'T' or 't' then A and B are k-by-n
   /// @param alpha Complex scalar.
   /// @param A Pointer to complex matrix.
   /// @param ldA Column length of the matrix A. If trans = 'N' or 'n' ldA>=n, otherwise ldA>=k.
   /// @param B Pointer to complex matrix.
   /// @param ldB Column length of the matrix B. If trans = 'N' or 'n' ldB>=n, otherwise ldB>=k.
   /// @param beta Complex scalar.
   /// @param C Pointer to complex symmetric n-by-n matrix C.
   /// Only the upper or lower triangular part of C is referenced, depending on the value of uplo above.
   /// @param ldC Column length of the matrix C.  ldC>=n
   /// @ingroup BLAS

   template <typename real_t>
   int SYR2K(char uplo, char trans, int_t n, int_t k, complex<real_t> alpha, complex<real_t> *A, int_t ldA, complex<real_t> *B, int_t ldB, complex<real_t> beta, complex<real_t> *C, int_t ldC)
   {
      using std::toupper;

      const complex<real_t> zero(0.0,0.0);
      const complex<real_t> one(1.0,0.0);
      int_t i,j,l;
      complex<real_t> *a,*b,*c,*at,*bt;
      complex<real_t> s,t;

      uplo=toupper(uplo);
      trans=toupper(trans);

      if((uplo!='U')&&(uplo!='L'))
         return -1;
      else if((trans!='N')&&(trans!='T'))
         return -2;
      else if(n<0)
         return -3;
      else if(k<0)
         return -4;
      else if(ldA<((trans=='N')?n:k))
         return -7;
      else if(ldB<((trans=='N')?n:k))
         return -9;
      else if(ldC<n)
         return -12;
      else if((n==0)||(((alpha==zero)||(k==0))&&(beta==one)))
         return 0;

      if(alpha==zero)
      {
         if(uplo=='U')
         {
            c=C;
            for(j=0;j<n;j++)
            {
               for(i=0;i<=j;i++)
                  c[i]*=beta;
               c+=ldC;
            }

         }
         else
         {
            c=C;
            for(j=0;j<n;j++)
            {
               for(i=j;i<n;i++)
                  c[i]*=beta;
               c+=ldC;
            }
         }
      }
      else if(trans=='N')
      {
         if(uplo=='U')
         {
            c=C;
            for(j=0;j<n;j++)
            {
               for(i=0;i<=j;i++)
                  c[i]*=beta;
               a=A;
               b=B;
               for(l=0;l<k;l++)
               {
                  s=alpha*b[j];
                  t=alpha*a[j];
                  for(i=0;i<=j;i++)
                     c[i]+=s*a[i]+t*b[i];
                  a+=ldA;
                  b+=ldB;
               }
               c+=ldC;
            }
         }
         else
         {
            c=C;
            for(j=0;j<n;j++)
            {
               for(i=j;i<n;i++)
                  c[i]*=beta;
               a=A;
               b=B;
               for(l=0;l<k;l++)
               {
                  s=alpha*b[j];
                  t=alpha*a[j];
                  for(i=j;i<n;i++)
                     c[i]+=s*a[i]+t*b[i];
                  a+=ldA;
                  b+=ldB;
               }
               c+=ldC;
            }
         }
      }
      else
      {
         if(uplo=='U')
         {
            c=C;
            at=A;
            bt=B;
            for(j=0;j<n;j++)
            {
               a=A;
               b=B;
               for(i=0;i<=j;i++)
               {
                  s=zero;
                  t=zero;
                  for(l=0;l<k;l++)
                  {
                     s+=b[l]*at[l];
                     t+=a[l]*bt[l];
                  }
                  c[i]=alpha*t+alpha*s+beta*c[i];
                  a+=ldA;
                  b+=ldB;
               }
               at+=ldA;
               bt+=ldB;
               c+=ldC;
            }
         }
         else
         {
            c=C;
            at=A;
            bt=B;
            for(j=0;j<n;j++)
            {
               a=A+j*ldA;
               b=B+j*ldB;
               for(i=j;i<n;i++)
               {
                  s=zero;
                  t=zero;
                  for(l=0;l<k;l++)
                  {
                     s+=b[l]*at[l];
                     t+=a[l]*bt[l];
                  }
                  c[i]=alpha*t+alpha*s+beta*c[i];
                  a+=ldA;
                  b+=ldB;
               }
               at+=ldA;
               bt+=ldB;
               c+=ldC;
            }
         }
      }
      return 0;
   }

#ifdef __latl_cblas
#include <cblas.h>

   template <> int SYR2K<float>(char uplo, char trans, int_t n, int_t k, float alpha, float *A, int_t ldA, float *B, int_t ldB, float beta, float *C, int_t ldC)
   {
      using std::toupper;
      uplo=toupper(uplo);
      trans=toupper(trans);
      if((uplo!='U')&&(uplo!='L'))
         return -1;
      else if((trans!='N')&&(trans!='T')&&(trans!='C'))
         return -2;
      else if(n<0)
         return -3;
      else if(k<0)
         return -4;
      else if(ldA<((trans=='N')?n:k))
         return -7;
      else if(ldB<((trans=='N')?n:k))
         return -9;
      else if(ldC<n)
         return -12;

      const CBLAS_UPLO Uplo=(uplo=='U')?CblasUpper:CblasLower;
      const CBLAS_TRANSPOSE Trans=(trans=='N')?CblasNoTrans:((trans=='T')?CblasTrans:CblasConjTrans);

      cblas_ssyr2k(CblasColMajor,Uplo,Trans,n,k,alpha,A,ldA,B,ldB,beta,C,ldC);

      return 0;
   }

   template <> int SYR2K<double>(char uplo, char trans, int_t n, int_t k, double alpha, double *A, int_t ldA, double *B, int_t ldB, double beta, double *C, int_t ldC)
   {
      using std::toupper;
      uplo=toupper(uplo);
      trans=toupper(trans);
      if((uplo!='U')&&(uplo!='L'))
         return -1;
      else if((trans!='N')&&(trans!='T')&&(trans!='C'))
         return -2;
      else if(n<0)
         return -3;
      else if(k<0)
         return -4;
      else if(ldA<((trans=='N')?n:k))
         return -7;
      else if(ldB<((trans=='N')?n:k))
         return -9;
      else if(ldC<n)
         return -12;

      const CBLAS_UPLO Uplo=(uplo=='U')?CblasUpper:CblasLower;
      const CBLAS_TRANSPOSE Trans=(trans=='N')?CblasNoTrans:((trans=='T')?CblasTrans:CblasConjTrans);

      cblas_dsyr2k(CblasColMajor,Uplo,Trans,n,k,alpha,A,ldA,B,ldB,beta,C,ldC);

      return 0;
   }

   template <> int SYR2K<float>(char uplo, char trans, int_t n, int_t k, complex<float> alpha, complex<float> *A, int_t ldA, complex<float> *B, int_t ldB, complex<float> beta, complex<float> *C, int_t ldC)
   {
      using std::toupper;
      uplo=toupper(uplo);
      trans=toupper(trans);
      if((uplo!='U')&&(uplo!='L'))
         return -1;
      else if((trans!='N')&&(trans!='T'))
         return -2;
      else if(n<0)
         return -3;
      else if(k<0)
         return -4;
      else if(ldA<((trans=='N')?n:k))
         return -7;
      else if(ldB<((trans=='N')?n:k))
         return -9;
      else if(ldC<n)
         return -12;

      const CBLAS_UPLO Uplo=(uplo=='U')?CblasUpper:CblasLower;
      const CBLAS_TRANSPOSE Trans=(trans=='N')?CblasNoTrans:((trans=='T')?CblasTrans:CblasConjTrans);

      cblas_csyr2k(CblasColMajor,Uplo,Trans,n,k,&alpha,A,ldA,B,ldB,&beta,C,ldC);

      return 0;
   }

   template <> int SYR2K<double>(char uplo, char trans, int_t n, int_t k, complex<double> alpha, complex<double> *A, int_t ldA, complex<double> *B, int_t ldB, complex<double> beta, complex<double> *C, int_t ldC)
   {
      using std::toupper;
      uplo=toupper(uplo);
      trans=toupper(trans);
      if((uplo!='U')&&(uplo!='L'))
         return -1;
      else if((trans!='N')&&(trans!='T'))
         return -2;
      else if(n<0)
         return -3;
      else if(k<0)
         return -4;
      else if(ldA<((trans=='N')?n:k))
         return -7;
      else if(ldB<((trans=='N')?n:k))
         return -9;
      else if(ldC<n)
         return -12;

      const CBLAS_UPLO Uplo=(uplo=='U')?CblasUpper:CblasLower;
      const CBLAS_TRANSPOSE Trans=(trans=='N')?CblasNoTrans:((trans=='T')?CblasTrans:CblasConjTrans);

      cblas_zsyr2k(CblasColMajor,Uplo,Trans,n,k,&alpha,A,ldA,B,ldB,&beta,C,ldC);

      return 0;
   }

#endif
}

#endif

max_stars_repo_name: langou/latl

max_issues_repo_name: langou/latl

max_forks_repo_name: langou/latl


text:
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

/*
 Contains the PERDIF program that implements personalized diffusions for
 recommendations.

 Code by: Dimitris Berberidis and Athanasios N. Nikolakopoulos
 University of Minnesota 2018-2019
*/

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

#include <inttypes.h>
#include <math.h>
#include <stdbool.h>
#include <stddef.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
//#include <cblas.h>
#include <mkl.h>

#include "perdif_fit.h"
#include "rec_defs.h"

#if DEBUG
static double cost_func(double *, double *, double *, sz_med_t);
#endif

static void matrix_matrix_product(double *C, double *A, double *B, int, int,  int);
static double *unfold_matrix(d_mat_t);
static void constr_QP_with_PG(double *, double *, double *, sz_med_t);
static void grammian_matrix(double *, double *, int, int);
static double max_abs_dif(double *, double *, sz_long_t);
// static void project_to_simplex(double *, sz_med_t);
// static void project_to_simplex_thanos(double *, double *, sz_med_t length);
static void simplexproj_Duchi(double *y, double *x, const unsigned int length);
static void matvec(double *, double *, double *, int, int);
static double frob_norm(double *, sz_med_t);

///////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////

void k_simplex(double *theta, d_mat_t P, int max_walk, double *dummy, bool rmse) {

  if (!rmse) {
    // adjust for BPR  metric
    for (int i = 1; i < P.num_rows; i++) {
      for (int j = 0; j < P.num_cols; j++)
        P.val[i][j] = P.val[0][j] - P.val[i][j];
    }
  }

  // Unfold P for BLAS compatibility
  double *P_unfold = unfold_matrix(P);

  // Prepare parameters for constraned QP

  // first the Hessian matrix
  double A[max_walk * max_walk];
  grammian_matrix(A, P_unfold, (int)P.num_rows, (int)P.num_cols);

  // then the linear part
  double b[max_walk];

  if (rmse) {
    for (int i = 0; i < max_walk; i++)
      b[i] = -2.0 * (P_unfold[i]);
  } else {
    for (int i = 0; i < max_walk; i++) {
      b[i] = 0.0;
      for (int j = 0; j < P.num_rows; j++) {
        b[i] += P.val[i][j];
      }
      b[i] *= -2.0;
    }
  }

  constr_QP_with_PG(theta, A, b, (sz_med_t)max_walk);

  // free
  free(P_unfold);
}

// Simple search returns the single best step ( i.e. theta = [0 0 0 1 0 0 .... 0
// 0 ] )

void single_best(double *theta, d_mat_t P, int max_walk, double *dummy, bool rmse) {

  if (!rmse) {
    // adjust for BPR  metric
    for (int i = 1; i < P.num_rows; i++) {
      for (int j = 0; j < P.num_cols; j++)
        P.val[i][j] = P.val[0][j] - P.val[i][j];
    }
  }

  for (int k = 0; k < max_walk; k++)
    theta[k] = 0.0;

  double loss[max_walk];

  // Compute loss of each k-step and find smallest
  int best_k = 0;
  double min_loss = 1.0e5;
  for (int k = 0; k < max_walk; k++) {
    loss[k] += pow(1.0 - P.val[0][k], 2);
    if (rmse) {
      for (int j = 1; j < P.num_rows; j++)
        loss[k] += P.val[j][k] * P.val[j][k];
    } else {
      for (int j = 1; j < P.num_rows; j++)
        loss[k] += pow(1.0 - P.val[j][k], 2);
    }
    if (loss[k] < min_loss) {
      min_loss = loss[k];
      best_k = k;
    }
  }

  theta[best_k] = 1.0;
}

void ppr_grid(double *theta, d_mat_t P, int max_walk, double *ppr_c, bool rmse) {

  if (!rmse) {
    // adjust for BPR  metric
    for (int i = 1; i < P.num_rows; i++) {
      for (int j = 0; j < P.num_cols; j++)
        P.val[i][j] = P.val[0][j] - P.val[i][j];
    }
  }

  // Unfold P for BLAS compatibility
  double *P_unfold = unfold_matrix(P);

  // Candidate diffusions
  double *C = (double *)malloc(P.num_rows * GRID_POINTS * sizeof(double));
  matrix_matrix_product(C, P_unfold, ppr_c, P.num_rows, max_walk, GRID_POINTS);

  // find the best one
  double loss[GRID_POINTS];
  int best_point = 0;
  double min_loss = 1.0e5;
  for (int i = 0; i < GRID_POINTS; i++) {
    // first the regularizer
    loss[i] = 0.0;
    // then the ratings
    if (rmse) {
      loss[i] += pow(1.0 - C[i], 2);
      for (int j = 1; j < P.num_rows; j++)
        loss[i] += pow(C[j * GRID_POINTS + i], 2);
    } else {
      for (int j = 0; j < P.num_rows; j++)
        loss[i] += pow(1.0 - C[j * GRID_POINTS + i], 2);
    }

    if (loss[i] < min_loss) {
      min_loss = loss[i];
      best_point = i;
    }
  }

  // return the best one
  for (int k = 0; k < max_walk; k++)
    theta[k] = ppr_c[k * GRID_POINTS + best_point];

  // free
  free(P_unfold);
  free(C);
}

void hk_grid(double *theta, d_mat_t P, int max_walk, double *hk_c, bool rmse) {

  if (!rmse) {
    // adjust for BPR  metric
    for (int i = 1; i < P.num_rows; i++) {
      for (int j = 0; j < P.num_cols; j++)
        P.val[i][j] = P.val[0][j] - P.val[i][j];
    }
  }

  // Unfold P for BLAS compatibility
  double *P_unfold = unfold_matrix(P);

  double *C = (double *)malloc(P.num_rows * GRID_POINTS * sizeof(double));
  matrix_matrix_product(C, P_unfold, hk_c, P.num_rows, max_walk, GRID_POINTS);

  // find the best one
  double loss[GRID_POINTS];
  int best_point = 0;
  double min_loss = 1.0e5;
  for (int i = 0; i < GRID_POINTS; i++) {
    // first the regularizer
    loss[i] = 0.0;

    // then the ratings
    if (rmse) {
      loss[i] += pow(1.0 - C[i], 2);
      for (int j = 1; j < P.num_rows; j++)
        loss[i] += pow(C[j * GRID_POINTS + i], 2);
    } else {
      for (int j = 0; j < P.num_rows; j++)
        loss[i] += pow(1.0 - C[j * GRID_POINTS + i], 2);
    }

    if (loss[i] < min_loss) {
      min_loss = loss[i];
      best_point = i;
    }
  }

  // return the best one
  for (int k = 0; k < max_walk; k++)
    theta[k] = hk_c[k * GRID_POINTS + best_point];

  // free
  free(P_unfold);
  free(C);
}

void dictionary_single(double *theta, d_mat_t P, int max_walk, double *dict, bool rmse) {

  if (!rmse) {
    // adjust for BPR  metric
    for (int i = 1; i < P.num_rows; i++) {
      for (int j = 0; j < P.num_cols; j++)
        P.val[i][j] = P.val[0][j] - P.val[i][j];
    }
  }

  // Unfold P for BLAS compatibility
  double *P_unfold = unfold_matrix(P);

  // Prepare matrices
  double *PD = (double *)malloc(P.num_rows * GRID_POINTS * sizeof(double));
  matrix_matrix_product(PD, P_unfold, dict, (int)P.num_rows, (int)P.num_cols,
                        GRID_POINTS);
  double DD[GRID_POINTS * GRID_POINTS];
  grammian_matrix(DD, dict, max_walk, GRID_POINTS);
  double A[GRID_POINTS * GRID_POINTS];
  grammian_matrix(A, PD, (int)P.num_rows, GRID_POINTS);

  double g_D[GRID_POINTS];
  for (int i = 0; i < GRID_POINTS; i++) {
    g_D[i] = 0.0;
  }

  if (rmse) {
    for (int i = 0; i < GRID_POINTS; i++)
      g_D[i] += -2.0 * PD[i];
  } else {
    for (int i = 0; i < GRID_POINTS; i++) {
      for (int j = 0; j < P.num_rows; j++)
        g_D[i] += -2.0 * PD[j * GRID_POINTS + i];
    }
  }

  double theta_dict[GRID_POINTS];

  // call QP solver for dictionary coefficients
  constr_QP_with_PG(theta_dict, A, g_D, (sz_med_t)GRID_POINTS);

  // construct diffusion from dictionary coefficients
  matvec(theta, dict, theta_dict, max_walk, GRID_POINTS);

  // free
  free(PD);
  free(P_unfold);
}

void dictionary(double *theta, d_mat_t P, int max_walk, double *dict, bool rmse) {

  if (!rmse) {
    // adjust for BPR  metric
    for (int i = 1; i < P.num_rows; i++) {
      for (int j = 0; j < P.num_cols; j++)
        P.val[i][j] = P.val[0][j] - P.val[i][j];
    }
  }

  // Unfold P for BLAS compatibility
  double *P_unfold = unfold_matrix(P);

  // Prepare matrices
  double *PD = (double *)malloc(P.num_rows * 2 * GRID_POINTS * sizeof(double));
  matrix_matrix_product(PD, P_unfold, dict, (int)P.num_rows, (int)P.num_cols,
                        2 * GRID_POINTS);
  double DD[4 * GRID_POINTS * GRID_POINTS];
  grammian_matrix(DD, dict, max_walk, 2 * GRID_POINTS);
  double A[4 * GRID_POINTS * GRID_POINTS];
  grammian_matrix(A, PD, (int)P.num_rows, 2 * GRID_POINTS);

  double g_D[2 * GRID_POINTS];
  for (int i = 0; i < 2 * GRID_POINTS; i++) {
    g_D[i] = 0.0;
  }

  if (rmse) {
    for (int i = 0; i < 2 * GRID_POINTS; i++)
      g_D[i] += -2.0 * PD[i];
  } else {
    for (int i = 0; i < 2 * GRID_POINTS; i++) {
      for (int j = 0; j < P.num_rows; j++)
        g_D[i] += -2.0 * PD[j * 2 * GRID_POINTS + i];
    }
  }

  double theta_dict[2 * GRID_POINTS];

  // call QP solver for dictionary coefficients
  constr_QP_with_PG(theta_dict, A, g_D, (sz_med_t)2 * GRID_POINTS);

  // construct diffusion from dictionary coefficients
  matvec(theta, dict, theta_dict, max_walk, 2 * GRID_POINTS);

  // free
  free(PD);
  free(P_unfold);
}

// Returns PPR coefficients for given points on a grid of parameter a
double *ppr_coefficients(int max_walk) {

  int grid_points = GRID_POINTS;

  double *coef = (double *)malloc(max_walk * grid_points * sizeof(double));
  double width = 1.0 / (double)(grid_points + 2);

  double col_sum[grid_points];
  for (int j = 0; j < grid_points; j++)
    col_sum[j] = 0.0;

  for (int i = 0; i < max_walk; i++) {
    for (int j = 0; j < grid_points; j++) {
      // pagerank weights per step
      coef[i * grid_points + j] = (i == 0)
                                      ? (1 - (0.01 + width * ((double)j)))
                                      : (1 - (0.01 + width * ((double)j))) *
                                            pow(0.01 + width * ((double)j), i);
      col_sum[j] += coef[i * grid_points + j];
    }
  }

  // Normalize coefficients so tha they sum to 1
  for (int i = 0; i < max_walk; i++) {
    for (int j = 0; j < grid_points; j++)
      coef[i * grid_points + j] /= col_sum[j];
  }

  return coef;
}

// Returns HK coefficients for given points on a grid of parameter \tau
double *hk_coefficients(int max_walk) {

  int grid_points = GRID_POINTS;

  double *coef = (double *)malloc(max_walk * grid_points * sizeof(double));
  double width = (double)MAX_HK_COEF / (double)(grid_points + 2);

  // compute factorials from 1 to max_walk
  sz_long_t factorial[max_walk];
  factorial[0] = 1;
  for (int j = 1; j < max_walk; j++)
    factorial[j] = factorial[j - 1] * j;

  double col_sum[grid_points];
  for (int j = 0; j < grid_points; j++)
    col_sum[j] = 0.0;

  for (int i = 0; i < max_walk; i++) {
    for (int j = 0; j < grid_points; j++) {
      coef[i * grid_points + j] =
          pow(1.0 + width * ((double)j), i) / (double)factorial[i];
      col_sum[j] += coef[i * grid_points + j];
    }
  }
  // Normalize coefficients so that they sum to 1
  for (int i = 0; i < max_walk; i++) {
    for (int j = 0; j < grid_points; j++)
      coef[i * grid_points + j] /= col_sum[j];
  }

  return coef;
}

// Returns dictionary coefficients by concatenating PPR and HK
double *dict_coefficients(int max_walk) {
  int grid_points = GRID_POINTS;

  double *coef = (double *)malloc(2 * max_walk * grid_points * sizeof(double));

  double *hk = hk_coefficients(max_walk);
  double *ppr = ppr_coefficients(max_walk);

  for (int i = 0; i < max_walk; i++) {
    for (int j = 0; j < grid_points; j++) {
      coef[2 * i * grid_points + j] = ppr[i * grid_points + j];
    }
    for (int j = grid_points; j < 2 * grid_points; j++) {
      coef[2 * i * grid_points + j] = hk[i * grid_points + j - grid_points];
    }
  }

  // free temporary coefs
  free(ppr);
  free(hk);
  return coef;
}

///////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////

// AUXILIARY FUNCTIONS

// Interface for CBLAS mutrix matrix product
// C =A*B
// A : m x k
// B : k x n
static void matrix_matrix_product(double *C, double *A, double *B, int m, int k,
                                  int n) {

  for (int i = 0; i < m * n; i++)
    C[i] = 0.0f;

  cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, n, k, 1.0f, A, k, B,
              n, 0.0f, C, n);
}

// Unfold matrix from A[i[j] to A[i*M + j]
static double *unfold_matrix(d_mat_t A) {

  double *A_u = (double *)malloc(A.num_cols * A.num_rows * sizeof(double));

  for (sz_long_t i = 0; i < A.num_rows; i++)
    memcpy(A_u + i * A.num_cols, A.val[i], A.num_cols * sizeof(double));

  return A_u;
}

// Solving constrained quadratic minimization via projected gradient descent
// The following function returns x =arg min {x^T*A*x +x^T*B} s.t. x in Prob.
// Simplex
static void constr_QP_with_PG(double *x, double *A, double *b, sz_med_t K) {
  double inf_norm, step_size;
  double x_temp[K];
  double x_prev[K];

  step_size = STEPSIZE_2 / frob_norm(A, K);

  // Initialize to uniform
  for (sz_med_t i = 0; i < K; i++)
    x[i] = 1.0f / (double)K;

  sz_med_t iter = 0;
  memcpy(x_prev, x, K * sizeof(double));
  // int count=0;
  do {
    iter++;
    // step_size = 1 / ((double)iter) * step_size;
    // step_size = 1/(sqrt(iter));//*step_size;
    // step_size = 1/(double)iter;
    // Take gradient step
    matvec(x_temp, A, x, K, K);

    for (sz_med_t j = 0; j < K; j++)
      x[j] -= step_size * (2.0f * x_temp[j] + b[j]);

    // project to feasible set
    simplexproj_Duchi(x, x, K);

#if DEBUG
    printf("\n COST: ");
    printf(" %lf ", cost_func(A, b, x, K));
#endif

    inf_norm = max_abs_dif(x_prev, x, (sz_long_t)K);

    memcpy(x_prev, x, K * sizeof(double));

  } while (iter < MAXIT_GD && inf_norm > GD_TOL_2);

  // printf("\n Counter = %d", count);
  // printf("\n Optimization finished after: %"PRIu32" iterations\n", (uint32_t)
  // iter);
}

// Grammian matrix G =A'*A using CBLAS
// A : m x n
static void grammian_matrix(double *G, double *A, int m, int n) {

  for (int i = 0; i < n * n; i++)
    G[i] = 0.0f;

  double *A_copy = (double *)malloc(m * n * sizeof(double));

  memcpy(A_copy, A, m * n * sizeof(double));

  cblas_dgemm(CblasRowMajor, CblasTrans, CblasNoTrans, n, n, m, 1.0f, A_copy, n,
              A, n, 0.0f, G, n);

  free(A_copy);
}

// Interface for CBLAS matrix vector product
// A : M x N
static void matvec(double *y, double *A, double *x, int M, int N) {

  // for (int i = 0; i < M; i++)
  // {
  // 	y[i] = 0.0f;
  // }

  cblas_dgemv(CblasRowMajor, CblasNoTrans, M, N, 1.0f, A, N, x, 1, 0.0f, y, 1);
}

//

/* Algorithm using partitioning with respect to a pivot, chosen randomly,
as given by Duchi et al. in "Efficient Projections onto the l1-Ball for
Learning in High Dimensions" */
static void simplexproj_Duchi(double *y, double *x, const unsigned int length) {
  double *auxlower = (x == y ? (double *)malloc(length * sizeof(double)) : x);
  double *auxupper = (double *)malloc(length * sizeof(double));
  double *aux = auxlower;
  double pivot;
  int auxlowerlength = 0;
  int auxupperlength = 1;
  int upperlength;
  int auxlength;
  int i = 0;
  int pospivot = (int)(rand() / (((double)RAND_MAX + 1.0) / length));
  double tauupper;
  double tau = (pivot = y[pospivot]) - 1.0;
  while (i < pospivot)
    if (y[i] < pivot)
      auxlower[auxlowerlength++] = y[i++];
    else {
      auxupper[auxupperlength++] = y[i];
      tau += (y[i++] - tau) / auxupperlength;
    }
  i++;
  while (i < length)
    if (y[i] < pivot)
      auxlower[auxlowerlength++] = y[i++];
    else {
      auxupper[auxupperlength++] = y[i];
      tau += (y[i++] - tau) / auxupperlength;
    }
  if (tau < pivot) {
    upperlength = auxupperlength;
    tauupper = tau;
    auxlength = auxlowerlength;
  } else {
    tauupper = 0.0;
    upperlength = 0;
    aux = auxupper + 1;
    auxlength = auxupperlength - 1;
  }
  while (auxlength > 0) {
    pospivot = (int)(rand() / (((double)RAND_MAX + 1.0) / auxlength));
    if (upperlength == 0)
      tau = (pivot = aux[pospivot]) - 1.0;
    else
      tau = tauupper + ((pivot = aux[pospivot]) - tauupper) / (upperlength + 1);
    i = 0;
    auxlowerlength = 0;
    auxupperlength = 1;
    while (i < pospivot)
      if (aux[i] < pivot)
        auxlower[auxlowerlength++] = aux[i++];
      else {
        auxupper[auxupperlength++] = aux[i];
        tau += (aux[i++] - tau) / (upperlength + auxupperlength);
      }
    i++;
    while (i < auxlength)
      if (aux[i] < pivot)
        auxlower[auxlowerlength++] = aux[i++];
      else {
        auxupper[auxupperlength++] = aux[i];
        tau += (aux[i++] - tau) / (upperlength + auxupperlength);
      }
    if (tau < pivot) {
      upperlength += auxupperlength;
      tauupper = tau;
      auxlength = auxlowerlength;
      aux = auxlower;
    } else {
      aux = auxupper + 1;
      auxlength = auxupperlength - 1;
    }
  }
  for (i = 0; i < length; i++)
    x[i] = (y[i] > tau ? y[i] - tauupper : 0.0000000001);
  double sum = 0.0;
  for (i = 0; i < length; i++)
    sum += x[i];
  for (i = 0; i < length; i++)
    x[i] /= sum;
  if (x == y)
    free(auxlower);
  free(auxupper);
}

// Infinity norm
static double max_abs_dif(double *a, double *b, sz_long_t len) {
  double dif, max = 0.0;

  for (sz_long_t i = 0; i < len; i++) {
    dif = fabs(a[i] - b[i]);
    max = (dif > max) ? dif : max;
  }

  return max;
}

#if DEBUG
// Evaluates quadratic with Hessian A and linear part b at x
static double cost_func(double *A, double *b, double *x, sz_med_t len) {

  double quad = 0.0f, lin = 0.0f;

  for (sz_med_t i = 0; i < len; i++) {
    for (sz_med_t j = 0; j < len; j++) {
      quad += A[i * len + j] * x[i] * x[j];
    }
    lin += b[i] * x[i];
  }
  return quad + lin;
}
#endif

// frobenious norm of double-valued square matrix
static double frob_norm(double *A, sz_med_t dim) {
  double norm = 0.0f;

  for (sz_med_t i = 0; i < dim * dim; i++) {
    norm += pow(A[i], 2.0f);
  }

  return sqrt(norm);
}

max_stars_repo_name: nikolakopoulos/Personalized-Diffusions

max_issues_repo_name: nikolakopoulos/Personalized-Diffusions

max_forks_repo_name: nikolakopoulos/Personalized-Diffusions


text:
#ifndef SMSVM_E_H
#define SMSVM_E_H

#include "Matrix.h"
#include "LADMM.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 <math.h>

//This class solves problem: min_{x,omega} \sum_i max(1- b_i^T(A_ix+ omega), 0)+ P(x) by LADMM
// where P(x)=\frac{lambda2}{2}\|x\|_2^2 +lambda1\|x\|_1.

// Let X= (x; omega), M= (b.*A b), b_i^T(A_ix+ omega)= M_iX
// phi_i(x)= 0, h_i(x)= max(1- x,0), g(X)= P(x).

template<typename L, typename D>
class SMSVM_e: public LADMM<L, D>
{
private:

D lambda1;

D lambda2;

Matrix<L,D> my_M;

D val_lambda_f;

protected:

public:

  SMSVM_e(const char* Matrix_file,D val_lambda1, D val_lambda2)
  :LADMM<L,D>(),my_M(Matrix_file)
  {
    lambda1=val_lambda1;
    lambda2=val_lambda2;
    this->lambda_f=0;
    my_M.matrix_modified_to_svm(my_M);

  }

  L get_n(){return my_M.nfeatures;}
  L get_m(){return my_M.nsamples;}

  inline D gradient_of_f_j(D x1, L i){
      return 0;
  }

  inline D value_of_g_j(D x, L i){
  	    if(i== this->m_1-1){
  	    	return 0;
		  }
		else{
           return lambda2*x*x/2+lambda1*fabs(x);
       }
   }

  inline D value_of_f_j(D x1, L i){return 0;}

  inline D value_of_h_j(D x, L i){
  	if(x<= 1.0){
  		return 1.0- x;
	  }
	  else{
	  	return 0;
	  }
  }

  inline D prox_of_g_j(D x1,D x2, L i){
  	    if(i== this->m_1 - 1){
  	    	return x1;
		  }
		else{
            D new_x;
        if(x1*x2> lambda1)
         new_x=(x1*x2- lambda1)/(x2+lambda2);
        else if(x1*x2< -lambda1)
         new_x=(x1*x2+ lambda1)/(x2+ lambda2);
        else
         new_x=0;
        return new_x;
        }
      }

// compute argmin_x h(x)+ x2/2*(x- x1)*(x- x1)
  inline D prox_of_h_j(D x1, D x2, L j){
  	if (x1+ 1.0/x2< 1.0){
  		return x1+ 1.0/x2;
   }
	  else if (x1> 1.0){
	  	return x1;
	  }
	  else{
	  	return 1.0;;
	  }
	}

  inline void set_matrix_M(){
    this->data_M=my_M;
  }

  inline void set_matrix_A(){
    this->data_A.nsamples=0;
    this->data_A.nfeatures=this->data_M.nfeatures;
    this->data_A.nnz=0;
    this->data_A.ptr.resize(1,0);
    this->data_A.ptr_t.resize(this->data_M.nfeatures+1,0);
  }

  void LADMM_solver(D beta_0,  D rho,vector<D> & x0,vector<D> & y0, vector<D> & lambda0, L max_nb_outer, L p_N_1, string filename1, D time){
    this->ADMM_solve_with_Linear(beta_0, rho,x0,y0, lambda0, max_nb_outer,  p_N_1, filename1, time);
}

};

#endif


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

#ifndef __GSL_MATRIX_INT_H__
#define __GSL_MATRIX_INT_H__

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

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

__BEGIN_DECLS

typedef struct
{
  size_t size1;
  size_t size2;
  size_t tda;
  int * data;
  gsl_block_int * block;
  int owner;
} gsl_matrix_int;

typedef struct
{
  gsl_matrix_int matrix;
} _gsl_matrix_int_view;

typedef _gsl_matrix_int_view gsl_matrix_int_view;

typedef struct
{
  gsl_matrix_int matrix;
} _gsl_matrix_int_const_view;

typedef const _gsl_matrix_int_const_view gsl_matrix_int_const_view;

/* Allocation */

gsl_matrix_int *
gsl_matrix_int_alloc (const size_t n1, const size_t n2);

gsl_matrix_int *
gsl_matrix_int_calloc (const size_t n1, const size_t n2);

gsl_matrix_int *
gsl_matrix_int_alloc_from_block (gsl_block_int * b,
                                   const size_t offset,
                                   const size_t n1,
                                   const size_t n2,
                                   const size_t d2);

gsl_matrix_int *
gsl_matrix_int_alloc_from_matrix (gsl_matrix_int * m,
                                    const size_t k1,
                                    const size_t k2,
                                    const size_t n1,
                                    const size_t n2);

gsl_vector_int *
gsl_vector_int_alloc_row_from_matrix (gsl_matrix_int * m,
                                        const size_t i);

gsl_vector_int *
gsl_vector_int_alloc_col_from_matrix (gsl_matrix_int * m,
                                        const size_t j);

void gsl_matrix_int_free (gsl_matrix_int * m);

/* Views */

_gsl_matrix_int_view
gsl_matrix_int_submatrix (gsl_matrix_int * m,
                            const size_t i, const size_t j,
                            const size_t n1, const size_t n2);

_gsl_vector_int_view
gsl_matrix_int_row (gsl_matrix_int * m, const size_t i);

_gsl_vector_int_view
gsl_matrix_int_column (gsl_matrix_int * m, const size_t j);

_gsl_vector_int_view
gsl_matrix_int_diagonal (gsl_matrix_int * m);

_gsl_vector_int_view
gsl_matrix_int_subdiagonal (gsl_matrix_int * m, const size_t k);

_gsl_vector_int_view
gsl_matrix_int_superdiagonal (gsl_matrix_int * m, const size_t k);

_gsl_vector_int_view
gsl_matrix_int_subrow (gsl_matrix_int * m, const size_t i,
                         const size_t offset, const size_t n);

_gsl_vector_int_view
gsl_matrix_int_subcolumn (gsl_matrix_int * m, const size_t j,
                            const size_t offset, const size_t n);

_gsl_matrix_int_view
gsl_matrix_int_view_array (int * base,
                             const size_t n1,
                             const size_t n2);

_gsl_matrix_int_view
gsl_matrix_int_view_array_with_tda (int * base,
                                      const size_t n1,
                                      const size_t n2,
                                      const size_t tda);

_gsl_matrix_int_view
gsl_matrix_int_view_vector (gsl_vector_int * v,
                              const size_t n1,
                              const size_t n2);

_gsl_matrix_int_view
gsl_matrix_int_view_vector_with_tda (gsl_vector_int * v,
                                       const size_t n1,
                                       const size_t n2,
                                       const size_t tda);

_gsl_matrix_int_const_view
gsl_matrix_int_const_submatrix (const gsl_matrix_int * m,
                                  const size_t i, const size_t j,
                                  const size_t n1, const size_t n2);

_gsl_vector_int_const_view
gsl_matrix_int_const_row (const gsl_matrix_int * m,
                            const size_t i);

_gsl_vector_int_const_view
gsl_matrix_int_const_column (const gsl_matrix_int * m,
                               const size_t j);

_gsl_vector_int_const_view
gsl_matrix_int_const_diagonal (const gsl_matrix_int * m);

_gsl_vector_int_const_view
gsl_matrix_int_const_subdiagonal (const gsl_matrix_int * m,
                                    const size_t k);

_gsl_vector_int_const_view
gsl_matrix_int_const_superdiagonal (const gsl_matrix_int * m,
                                      const size_t k);

_gsl_vector_int_const_view
gsl_matrix_int_const_subrow (const gsl_matrix_int * m, const size_t i,
                               const size_t offset, const size_t n);

_gsl_vector_int_const_view
gsl_matrix_int_const_subcolumn (const gsl_matrix_int * m, const size_t j,
                                  const size_t offset, const size_t n);

_gsl_matrix_int_const_view
gsl_matrix_int_const_view_array (const int * base,
                                   const size_t n1,
                                   const size_t n2);

_gsl_matrix_int_const_view
gsl_matrix_int_const_view_array_with_tda (const int * base,
                                            const size_t n1,
                                            const size_t n2,
                                            const size_t tda);

_gsl_matrix_int_const_view
gsl_matrix_int_const_view_vector (const gsl_vector_int * v,
                                    const size_t n1,
                                    const size_t n2);

_gsl_matrix_int_const_view
gsl_matrix_int_const_view_vector_with_tda (const gsl_vector_int * v,
                                             const size_t n1,
                                             const size_t n2,
                                             const size_t tda);

/* Operations */

void gsl_matrix_int_set_zero (gsl_matrix_int * m);
void gsl_matrix_int_set_identity (gsl_matrix_int * m);
void gsl_matrix_int_set_all (gsl_matrix_int * m, int x);

int gsl_matrix_int_fread (FILE * stream, gsl_matrix_int * m) ;
int gsl_matrix_int_fwrite (FILE * stream, const gsl_matrix_int * m) ;
int gsl_matrix_int_fscanf (FILE * stream, gsl_matrix_int * m);
int gsl_matrix_int_fprintf (FILE * stream, const gsl_matrix_int * m, const char * format);

int gsl_matrix_int_memcpy(gsl_matrix_int * dest, const gsl_matrix_int * src);
int gsl_matrix_int_swap(gsl_matrix_int * m1, gsl_matrix_int * m2);

int gsl_matrix_int_swap_rows(gsl_matrix_int * m, const size_t i, const size_t j);
int gsl_matrix_int_swap_columns(gsl_matrix_int * m, const size_t i, const size_t j);
int gsl_matrix_int_swap_rowcol(gsl_matrix_int * m, const size_t i, const size_t j);
int gsl_matrix_int_transpose (gsl_matrix_int * m);
int gsl_matrix_int_transpose_memcpy (gsl_matrix_int * dest, const gsl_matrix_int * src);

int gsl_matrix_int_max (const gsl_matrix_int * m);
int gsl_matrix_int_min (const gsl_matrix_int * m);
void gsl_matrix_int_minmax (const gsl_matrix_int * m, int * min_out, int * max_out);

void gsl_matrix_int_max_index (const gsl_matrix_int * m, size_t * imax, size_t *jmax);
void gsl_matrix_int_min_index (const gsl_matrix_int * m, size_t * imin, size_t *jmin);
void gsl_matrix_int_minmax_index (const gsl_matrix_int * m, size_t * imin, size_t * jmin, size_t * imax, size_t * jmax);

int gsl_matrix_int_equal (const gsl_matrix_int * a, const gsl_matrix_int * b);

int gsl_matrix_int_isnull (const gsl_matrix_int * m);
int gsl_matrix_int_ispos (const gsl_matrix_int * m);
int gsl_matrix_int_isneg (const gsl_matrix_int * m);
int gsl_matrix_int_isnonneg (const gsl_matrix_int * m);

int gsl_matrix_int_add (gsl_matrix_int * a, const gsl_matrix_int * b);
int gsl_matrix_int_sub (gsl_matrix_int * a, const gsl_matrix_int * b);
int gsl_matrix_int_mul_elements (gsl_matrix_int * a, const gsl_matrix_int * b);
int gsl_matrix_int_div_elements (gsl_matrix_int * a, const gsl_matrix_int * b);
int gsl_matrix_int_scale (gsl_matrix_int * a, const double x);
int gsl_matrix_int_add_constant (gsl_matrix_int * a, const double x);
int gsl_matrix_int_add_diagonal (gsl_matrix_int * a, const double x);

/***********************************************************************/
/* The functions below are obsolete                                    */
/***********************************************************************/
int gsl_matrix_int_get_row(gsl_vector_int * v, const gsl_matrix_int * m, const size_t i);
int gsl_matrix_int_get_col(gsl_vector_int * v, const gsl_matrix_int * m, const size_t j);
int gsl_matrix_int_set_row(gsl_matrix_int * m, const size_t i, const gsl_vector_int * v);
int gsl_matrix_int_set_col(gsl_matrix_int * m, const size_t j, const gsl_vector_int * v);
/***********************************************************************/

/* inline functions if you are using GCC */

INLINE_DECL int   gsl_matrix_int_get(const gsl_matrix_int * m, const size_t i, const size_t j);
INLINE_DECL void    gsl_matrix_int_set(gsl_matrix_int * m, const size_t i, const size_t j, const int x);
INLINE_DECL int * gsl_matrix_int_ptr(gsl_matrix_int * m, const size_t i, const size_t j);
INLINE_DECL const int * gsl_matrix_int_const_ptr(const gsl_matrix_int * m, const size_t i, const size_t j);

#ifdef HAVE_INLINE
INLINE_FUN
int
gsl_matrix_int_get(const gsl_matrix_int * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_VAL("first index out of range", GSL_EINVAL, 0) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_VAL("second index out of range", GSL_EINVAL, 0) ;
        }
    }
#endif
  return m->data[i * m->tda + j] ;
}

INLINE_FUN
void
gsl_matrix_int_set(gsl_matrix_int * m, const size_t i, const size_t j, const int x)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_VOID("first index out of range", GSL_EINVAL) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_VOID("second index out of range", GSL_EINVAL) ;
        }
    }
#endif
  m->data[i * m->tda + j] = x ;
}

INLINE_FUN
int *
gsl_matrix_int_ptr(gsl_matrix_int * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
        }
    }
#endif
  return (int *) (m->data + (i * m->tda + j)) ;
}

INLINE_FUN
const int *
gsl_matrix_int_const_ptr(const gsl_matrix_int * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
        }
    }
#endif
  return (const int *) (m->data + (i * m->tda + j)) ;
}

#endif

__END_DECLS

#endif /* __GSL_MATRIX_INT_H__ */

max_stars_repo_licenses: Net-SNMP Xnet

max_issues_repo_licenses: Net-SNMP Xnet

max_forks_repo_licenses: Net-SNMP Xnet


text:
/* Copyright (c) 2014, Giuseppe Argentieri <giuseppe.argentieri@ts.infn.it>

 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright notice,
 * this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright notice,
 * this list of conditions and the following disclaimer in the documentation
 * and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 * POSSIBILITY OF SUCH DAMAGE.
 *
 */

/*
 * =====================================================================================
 *
 *       Filename:  Img_ss.c
 *
 *    Description:
 *
 *        Version:  1.0
 *        Created:  15/04/2013 16:44:51
 *       Revision:  none
 *       Compiler:  gcc
 *
 *         Author:  Giuseppe Argentieri (), argenti@ts.infn.it
 *   Organization:
 *
 * =====================================================================================
 */

#include "funcs.h"
/* #include <gsl/gsl_sf_expint.h> */

int im_gss ( void* params, double* val, double* error )
{
	struct f_params* pars = (struct f_params*) params ;
	double o_c, b, O, o_1, alpha ;
	assign_p ( pars, &o_c, &b, &O, &o_1 ) ;
	alpha = pars->alpha ;

	double e1, e2, E1, E2 ;
	double err1, err2, ERR1, ERR2 ;

	double beta1 = O + o_1 ;
	double beta2 = o_1 - O ;
	double mu = 1/o_c ;

	expi ( beta1*mu, &e1, &err1 ) ;
	expi_plus ( beta1*mu, &E1, &ERR1 ) ;
	expi_plus ( beta2*mu, &E2, &ERR2 ) ;
	expi ( beta2*mu, &e2, &err2 ) ;

/* 	e1 = - gsl_sf_expint_E1 ( beta1*mu ) ;
 * 	e2 = - gsl_sf_expint_E1 ( beta2*mu ) ;
 * 	E1 = gsl_sf_expint_Ei ( beta1*mu ) ;
 * 	E2 = gsl_sf_expint_Ei ( beta2*mu ) ;
 */

	double imgss = (alpha/4)*(beta1*(e1*exp(beta1*mu)-E1*exp(-beta1*mu)) +
				  beta2*(E2*exp(-beta2*mu)-e2*exp(beta2*mu))) ;
	*val = imgss ;

	double err = (alpha/4)*(beta1*(err1*exp(beta1*mu))+ERR1*exp(-beta1*mu) +
				  beta2*(ERR2*exp(-beta2*mu)+err2*exp(beta2*mu))) ;
	*error = err ;

	return 0;
}