text:
/*
 * Copyright 2008-2016 Jan Gasthaus
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

/**
 * Functions for computing and tabulating the generalized stirling numbers
 * of type (-1,-d,0).
 *
 */

#ifndef STIRLING_H_
#define STIRLING_H_
#include <cmath>
#include <vector>
#include <string>
#include <iostream>
#include <cassert>
#include <gsl/gsl_sf_gamma.h>

namespace gatsby { namespace libplump {

typedef std::vector<double> d_vec;
typedef std::vector<d_vec> d_vec_vec;

/**
 * computes log(exp(a) + exp(b)) while avoiding numerical
 * instabilities.
 */
inline double fast_logsumexp(double a, double b) {
  if (a>b) {
    return log(1+exp(b-a)) + a;
  } else {
    return log(1+exp(a-b)) + b;
  }
}

/**
 * computes log(exp(a) - exp(b)) while avoiding numerical
 * instabilities.
 */
inline double fast_logminusexp(double a, double b) {
  if (a>b) {
    return log(1-exp(b-a)) + a;
  } else {
    return log(1-exp(a-b)) + b;
  }
}

/**
 * Compute S_d(c,t) recursively by directly applying the recursion
 * S_d(c,t) = S_d(c-1,t-1) + (c-1 - t*d)S_d(c-1,t)
 *
 * This does a _lot_ of duplicate work and should only be used for debugging!
 */
double gen_stirling_recursive(double d, int c, int t);

/**
 * Compute S_d(c,t) recursively in log space by directly applying the recursion
 * S_d(c,t) = S_d(c-1,t-1) + (c-1 - t*d)S_d(c-1,t)
 *
 * This does a _lot_ of duplicate work and should only be used for debugging!
 */
double log_gen_stirling_recursive(double d, int c, int t);

/**
 * Compute S_d(c,t) directly in log space using the recursion
 * S_d(c,t) = S_d(c-1,t-1) + (c-1 - t*d)S_d(c-1,t)
 */
double log_gen_stirling_direct(double d, int c, int t);

double log_gen_stirling_ratio(double d, int c, int t);

d_vec_vec log_gen_stirling_table(double d, int c);

void log_gen_stirling_table_extend(double d, int c, d_vec_vec& table);

double log_get_stirling_from_table(d_vec_vec& table, int c, int t);

class stirling_generator_recompute_log {
  public:
    stirling_generator_recompute_log(double d, int c, int t);

    double ratio(int c, int t);

    static std::string statsToString();

  private:
    double d;
};

class stirling_generator_fast_log {

  public:
    stirling_generator_fast_log(double d, int c, int t);

    double get(int c, int t);

    double ratio(int c, int t);

    static std::string statsToString();

  private:
    void incC();
    void decC();

    double d;
    int current_c;
    d_vec row, col, prev_col;
};

/**
 * Class to encapsulate the generation of ratios of stirling numbers
 * for removing customers from restaurants.
 *
 * One of these objects should be constructed for each restaurant,
 * providing the current number of tables, and the total number of customers.
 */
class stirling_generator_full_log {
  public:
    stirling_generator_full_log(double d, int c, int t);

    double ratio(int c, int t);

    double getLog(int c, int t);

    static std::string statsToString();

  private:

    static int global_c_max, num_ratio_calls, num_extends, num_construct;

    d_vec_vec table;
    int c_max;
    double d;
};

inline double log_get_stirling_from_table(d_vec_vec& table, int c, int t) {
  // c and t must be non-negative
  assert(c >= 0 && t>= 0);

  if ((c == 1  && t == 1) || (c == 0 && t == 0))
    return 0;
  if (c == 0 || t == 0)
    return -INFINITY;
  if (t > c)
    return -INFINITY;
  if(c == t)
    return 0;
  if (t > (int)table.size() || (c - t) > (int)table[t-1].size()) {
    std::cout<< c << ", " << t << std::endl;
  }
  return table[t-1][c-t-1];
}

inline double log_stirling_asymptotic(double d, int c, int t) {
  return   gsl_sf_lngamma(c)
         - gsl_sf_lngamma(1 - d)
         - gsl_sf_lngamma(t)
         - (t-1) * log(d)
         - d * log(c);
}

// class log_gen_stirling_table {
//     private:
//         std::vector<std::vector<double> > table;
//         size_t c_cur, t_cur;
//
//         /*
//          * Grow the arrays to the required size.
//          */
//         void _grow(size_t c_max, size_t t_max) {
//             for(int i=t_cur;i<t_max;i++) {
//                 std::vector<double> tmp(c_cur,0);
//                 table.push_back(std::vector<double>());
//             }
//         }
//
//         /*
//          * Tabulate up to the specified limits
//          */
//         void _tabulate(size_t c_max, size_t t_max) {
//             for t in xrange(1,max_t):
//         for c in xrange(t,max_c):
//             a[c,t] = a[c-1,t-1] + (c-1-d*t)*a[c-1,t]
//
//         }
//
//     public:
//         double d;
//         gen_stirling_table(double d) : d(d), c_cur(0), t_cur(0)  {}
//
// };

}} // namespace gatsby::libplump

#endif /* STIRLING_H_ */

max_stars_repo_name: jgasthaus/libPLUMP

max_issues_repo_name: jgasthaus/libPLUMP

max_forks_repo_name: jgasthaus/libPLUMP


text:
#ifndef __MMA__
#define __MMA__

#include <petsc.h>

/* -----------------------------------------------------------------------------
Authors: Niels Aage
 Copyright (C) 2013-2020,
This MMA implementation is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version.
This Module is distributed in the hope that it will be useful,implementation
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with this Module; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
-------------------------------------------------------------------------- */
class MMA {
  public:
    // Construct using defaults subproblem penalization
    MMA(PetscInt n, PetscInt m, Vec x);
    // User defined subproblem penalization
    MMA(PetscInt n, PetscInt m, Vec x, PetscScalar* a, PetscScalar* c, PetscScalar* d);
    // Initialize with restart from itr
    MMA(PetscInt n, PetscInt m, PetscInt itr, Vec xo1, Vec xo2, Vec U, Vec L);
    // Initialize with restart and specify subproblem parameters
    MMA(PetscInt n, PetscInt m, PetscInt itr, Vec xo1, Vec xo2, Vec U, Vec L, PetscScalar* a, PetscScalar* c,
        PetscScalar* d);
    // Destructor
    ~MMA();

    // Set and solve a subproblem: return new xval
    PetscErrorCode Update(Vec xval, Vec dfdx, PetscScalar* gx, Vec* dgdx, Vec xmin, Vec xmax);

    // Return necessary data for possible restart
    PetscErrorCode Restart(Vec xo1, Vec xo2, Vec U, Vec L);

    // Set the aggresivity of the moving asymptotes
    PetscErrorCode SetAsymptotes(PetscScalar init, PetscScalar decrease, PetscScalar increase);

    // do/don't add convexity approx to constraints: default=false
    PetscErrorCode ConstraintModification(PetscBool conMod) {
        constraintModification = conMod;
        return 0;
    };

    // val=0: default, val=1: increase robustness, i.e
    // control the spacing between L < alp < x < beta < U,
    PetscErrorCode SetRobustAsymptotesType(PetscInt val);

    // Sets outer movelimits on all primal design variables
    // This is often requires to prevent the solver from oscilating
    PetscErrorCode SetOuterMovelimit(PetscScalar Xmin, PetscScalar Xmax, PetscScalar movelim, Vec x, Vec xmin,
                                     Vec xmax);

    // Return KKT residual norms (norm2 and normInf)
    PetscErrorCode KKTresidual(Vec xval, Vec dfdx, PetscScalar* gx, Vec* dgdx, Vec xmin, Vec xmax, PetscScalar* norm2,
                               PetscScalar* normInf);

    // Inf norm on diff between two vectors: SHOULD NOT BE HERE - USE BASIC
    // PETSc!!!!!
    PetscScalar DesignChange(Vec x, Vec xold);

  private:
    // Set up the MMA subproblem based on old x's and xval
    PetscErrorCode GenSub(Vec xval, Vec dfdx, PetscScalar* gx, Vec* dgdx, Vec xmin, Vec xmax);

    // Interior point solver for the subproblem
    PetscErrorCode SolveDIP(Vec xval);

    // Compute primal vars based on dual solution
    PetscErrorCode XYZofLAMBDA(Vec x);

    // Dual gradient
    PetscErrorCode DualGrad(Vec x);

    // Dual Hessian
    PetscErrorCode DualHess(Vec x);

    // Dual line search
    PetscErrorCode DualLineSearch();

    // Dual residual
    PetscScalar DualResidual(Vec x, PetscScalar epsi);

    // Problem size and iteration counter
    PetscInt n, m, k;

    // "speed-control" for the asymptotes
    PetscScalar asyminit, asymdec, asyminc;

    // do/don't add convexity constraint approximation in subproblem
    PetscBool constraintModification; // default = FALSE

    // Bool specifying if non lin constraints are included or not
    PetscBool NonLinConstraints;

    // 0: (default) span between alp L x U beta,
    // 1: increase the span for further robustness
    PetscInt RobustAsymptotesType;

    // Local vectors: penalty numbers for subproblem
    PetscScalar *a, *c, *d;

    // Local vectors: elastic variables
    PetscScalar* y;
    PetscScalar  z;

    // Local vectors: Lagrange multipliers:
    PetscScalar *lam, *mu, *s;

    // Global: Asymptotes, bounds, objective approx., constraint approx.
    Vec L, U, alpha, beta, p0, q0, *pij, *qij;

    // Local: subproblem constant terms, dual gradient, dual hessian
    PetscScalar *b, *grad, *Hess;

    // Global: Old design variables
    Vec xo1, xo2;

    // Math helpers
    PetscErrorCode Factorize(PetscScalar* K, PetscInt nn);
    PetscErrorCode Solve(PetscScalar* K, PetscScalar* x, PetscInt nn);
    PetscScalar    Min(PetscScalar d1, PetscScalar d2);
    PetscScalar    Max(PetscScalar d1, PetscScalar d2);
    PetscInt       Min(PetscInt d1, PetscInt d2);
    PetscInt       Max(PetscInt d1, PetscInt d2);
    PetscScalar    Abs(PetscScalar d1);
};

#endif

ext: h

max_stars_repo_name: TTitscher/BitJug

max_issues_repo_name: TTitscher/BitJug

max_forks_repo_name: TTitscher/BitJug


text:
#pragma once
#include <vector>
#include <gsl/gsl>
#include "halley/utils/utils.h"

namespace Halley
{
	class NetworkPacketBase
	{
	public:
		size_t copyTo(gsl::span<gsl::byte> dst) const;
		size_t getSize() const;
		gsl::span<const gsl::byte> getBytes() const;

		NetworkPacketBase(NetworkPacketBase&& other) = delete;
		NetworkPacketBase& operator=(NetworkPacketBase&& other) = delete;

	protected:
		NetworkPacketBase();
		NetworkPacketBase(gsl::span<const gsl::byte> data, size_t prePadding);

		size_t dataStart;
		std::vector<gsl::byte> data;
	};

	class OutboundNetworkPacket : public NetworkPacketBase
	{
	public:
		OutboundNetworkPacket(const OutboundNetworkPacket& other);
		explicit OutboundNetworkPacket(OutboundNetworkPacket&& other) noexcept;
		explicit OutboundNetworkPacket(gsl::span<const gsl::byte> data);
		explicit OutboundNetworkPacket(const Bytes& data);

		void addHeader(gsl::span<const gsl::byte> src);

		template <typename T>
		void addHeader(const T& h)
		{
			addHeader(gsl::as_bytes(gsl::span<const T>(&h, 1)));
		}

		OutboundNetworkPacket& operator=(OutboundNetworkPacket&& other) noexcept;
	};

	class InboundNetworkPacket : public NetworkPacketBase
	{
	public:
		InboundNetworkPacket();
		explicit InboundNetworkPacket(InboundNetworkPacket&& other);
		explicit InboundNetworkPacket(gsl::span<const gsl::byte> data);
		void extractHeader(gsl::span<gsl::byte> dst);

		template <typename T>
		void extractHeader(T& h)
		{
			extractHeader(gsl::as_writeable_bytes(gsl::span<T>(&h, 1)));
		}

		InboundNetworkPacket& operator=(InboundNetworkPacket&& other);
	};
}

max_stars_repo_name: sunhay/halley

max_issues_repo_name: sunhay/halley

max_forks_repo_name: sunhay/halley


text:
#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sum.h>

gsl_sum_levin_u_workspace *
gsl_sum_levin_u_alloc (size_t n)
{
  gsl_sum_levin_u_workspace * w;

  if (n == 0)
    {
      GSL_ERROR_VAL ("length n must be positive integer", GSL_EDOM, 0);
    }

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

  if (w == NULL)
    {
      GSL_ERROR_VAL ("failed to allocate struct", GSL_ENOMEM, 0);
    }

  w->q_num = (double *) malloc (n * sizeof (double));

  if (w->q_num == NULL)
    {
      free(w) ; /* error in constructor, prevent memory leak */

      GSL_ERROR_VAL ("failed to allocate space for q_num", GSL_ENOMEM, 0);
    }

  w->q_den = (double *) malloc (n * sizeof (double));

  if (w->q_den == NULL)
    {
      free (w->q_num);
      free (w) ; /* error in constructor, prevent memory leak */

      GSL_ERROR_VAL ("failed to allocate space for q_den", GSL_ENOMEM, 0);
    }

  w->dq_num = (double *) malloc (n * n * sizeof (double));

  if (w->dq_num == NULL)
    {
      free (w->q_den);
      free (w->q_num);
      free(w) ; /* error in constructor, prevent memory leak */

      GSL_ERROR_VAL ("failed to allocate space for dq_num", GSL_ENOMEM, 0);
    }

  w->dq_den = (double *) malloc (n * n * sizeof (double));

  if (w->dq_den == NULL)
    {
      free (w->dq_num);
      free (w->q_den);
      free (w->q_num);
      free (w) ; /* error in constructor, prevent memory leak */

      GSL_ERROR_VAL ("failed to allocate space for dq_den", GSL_ENOMEM, 0);
    }

  w->dsum = (double *) malloc (n * sizeof (double));

  if (w->dsum == NULL)
    {
      free (w->dq_den);
      free (w->dq_num);
      free (w->q_den);
      free (w->q_num);
      free (w) ; /* error in constructor, prevent memory leak */

      GSL_ERROR_VAL ("failed to allocate space for dsum", GSL_ENOMEM, 0);
    }

  w->size = n;
  w->terms_used = 0;
  w->sum_plain = 0;

  return w;
}

void
gsl_sum_levin_u_free (gsl_sum_levin_u_workspace * w)
{
  free (w->dsum);
  free (w->dq_den);
  free (w->dq_num);
  free (w->q_den);
  free (w->q_num);
  free (w);
}

ext: c

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
#include <stdio.h>       /* standard I/O routines                 */
#include <stdlib.h>
#include <string.h>
#include <stdio.h>       /* standard I/O routines                 */
#include <omp.h>
#include <mpi.h>

#include <sys/time.h>
#include "utils.h"
//#include <cblas.h>
#ifdef PAPI
#include "cscs_papi.h"
#endif

#define NN 8000;
#define NTIMES 10

#if 0
extern "C"{
void dgemm_ (char *transa,
		char *transb,
		int *m, int *n, int *k,
		double *alpha, double *A, int *lda,
		double *B, int *ldb,
		double *beta, double *C, int *ldc) ;
}
#endif

double mysecond()
{
	struct timeval tp;
	struct timezone tzp;
	int i;

	i = gettimeofday(&tp,&tzp);
	return ( (double) tp.tv_sec + (double) tp.tv_usec * 1.e-6 );
}

int main (int argc, char *argv[])
{

	int N = NN;
	int M = NN;
	int K = NN;
	int ii;

	/* Find out my identity in the default communicator */

	int Ncpu = 0;

	if ( argc == 2 )
	{
		N    = atoi(argv[1]);
		M = N;
		K = N;
	}

	double  alpha =  1.;
	double  beta  = -1.;

	int lda    = M;
	int ldb    = K;
	int ldc    = M;

	int size_A = K*lda;
	int size_B = N*ldb;
	int size_C = N*ldc;

	double *A  = (double*) malloc(sizeof(double)*size_A);
	if (A == 0) printf("Could not allocate A.\n");

	double *B  = (double*) malloc(sizeof(double)*size_B);
	if (B == 0) printf("Could not allocate B.\n");

	double *C  = (double*) malloc(sizeof(double)*size_C);
	if (C == 0) printf("Could not allocate C.\n");

	double *Cg = (double*) malloc(sizeof(double)*size_C);
	if (Cg == 0) printf("Could not allocate Cg.\n");

	fill(A,  size_A,  31.);
	eye (B,     ldb,   N );
	fill(C,  size_C,  31.);
	fill(Cg, size_C,  31.);

	int t;

	int nthreads, tid, procs, maxthreads, inpar, dynamic, nested;

	char transa = 'n';
	char transb = 'n';

	/* Get environment information */
#ifdef _OPENMP
	procs      = omp_get_num_procs();
	nthreads   = omp_get_num_threads();
	maxthreads = omp_get_max_threads();
	inpar      = omp_in_parallel();
	dynamic    = omp_get_dynamic();
	nested     = omp_get_nested();
	/* Print environment information */
	printf("Number of processors          = %d\n", procs);
	printf("Number of threads             = %d\n", nthreads);
	printf("Max threads                   = %d\n", maxthreads);
	printf("In parallel?                  = %d\n", inpar);
	printf("Dynamic threads enabled?      = %d\n", dynamic);
	printf("Nested parallelism supported? = %d\n", nested);
#endif

	double otime = 0.;
	otime -= mysecond();

	for (ii = 0; ii < NTIMES; ++ii)
		dgemm_(&transa, &transb,
				&M, &N, &K,
				&alpha, A, &lda,
				B, &ldb,
				&beta, Cg, &ldc);

	otime += mysecond();
	otime /= NTIMES;

	printf("Size = %d, Gflops Max: %f  Min: %f  Avg:  %f\n", N, 2.*M*N*K/otime/1e9, 2.*M*N*K/otime/1e9, 2.*M*N*K/otime/1e9);

	free(A);
	free(B);
	free(C);
	free(Cg);

	exit(0);
}


text:
//g++ open_blas_measure.cpp -I /include/ -L/lib/ -lopenblas -lpthread
// #include <cblas.h>
#include <stdio.h>
#include <iostream>

#ifndef OPEN_BLAS
#define OPEN_BLAS

extern "C"{
  int dgemm_(char *, char *, int *, int *, int *, double *, double *, int *,
              double *, int *, double *, double *, int *);
}

/* Assuming row major order of matrices.

    This function is used to call open_blas multiplication. It takes the paramters given and
    puts the result of multiplication in *C.

    To be used in multiplication by open_blas
*/
void call_open_blas(double *A, int rowA, double *B, int colB, int commonSide, double *C) {
    char t = 'N';
    double alpha =1.0;
    double beta = 0.0;
    dgemm_(&t, &t, &rowA, &colB, &commonSide, &alpha, A, &rowA, B, &commonSide, &beta, C, &rowA);
}

#endif

max_stars_repo_path: Subtask 2/include/impl/OBM.h

max_issues_repo_path: Subtask 2/include/impl/OBM.h

max_forks_repo_path: Subtask 2/include/impl/OBM.h


text:
/* specfunc/gsl_sf_fermi_dirac.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 *
 * 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.
 */

/* Author:  G. Jungman */

#ifndef __GSL_SF_FERMI_DIRAC_H__
#define __GSL_SF_FERMI_DIRAC_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

/* Complete Fermi-Dirac Integrals:
 *
 *  F_j(x)   := 1/Gamma[j+1] Integral[ t^j /(Exp[t-x] + 1), {t,0,Infinity}]
 *
 *
 * Incomplete Fermi-Dirac Integrals:
 *
 *  F_j(x,b) := 1/Gamma[j+1] Integral[ t^j /(Exp[t-x] + 1), {t,b,Infinity}]
 */

/* Complete integral F_{-1}(x) = e^x / (1 + e^x)
 *
 * exceptions: GSL_EUNDRFLW
 */
GSL_FUN int     gsl_sf_fermi_dirac_m1_e(const double x, gsl_sf_result * result);
GSL_FUN double     gsl_sf_fermi_dirac_m1(const double x);

/* Complete integral F_0(x) = ln(1 + e^x)
 *
 * exceptions: GSL_EUNDRFLW
 */
GSL_FUN int     gsl_sf_fermi_dirac_0_e(const double x, gsl_sf_result * result);
GSL_FUN double     gsl_sf_fermi_dirac_0(const double x);

/* Complete integral F_1(x)
 *
 * exceptions: GSL_EUNDRFLW, GSL_EOVRFLW
 */
GSL_FUN int     gsl_sf_fermi_dirac_1_e(const double x, gsl_sf_result * result);
GSL_FUN double     gsl_sf_fermi_dirac_1(const double x);

/* Complete integral F_2(x)
 *
 * exceptions: GSL_EUNDRFLW, GSL_EOVRFLW
 */
GSL_FUN int     gsl_sf_fermi_dirac_2_e(const double x, gsl_sf_result * result);
GSL_FUN double     gsl_sf_fermi_dirac_2(const double x);

/* Complete integral F_j(x)
 * for integer j
 *
 * exceptions: GSL_EUNDRFLW, GSL_EOVRFLW
 */
GSL_FUN int     gsl_sf_fermi_dirac_int_e(const int j, const double x, gsl_sf_result * result);
GSL_FUN double     gsl_sf_fermi_dirac_int(const int j, const double x);

/* Complete integral F_{-1/2}(x)
 *
 * exceptions: GSL_EUNDRFLW, GSL_EOVRFLW
 */
GSL_FUN int     gsl_sf_fermi_dirac_mhalf_e(const double x, gsl_sf_result * result);
GSL_FUN double     gsl_sf_fermi_dirac_mhalf(const double x);

/* Complete integral F_{1/2}(x)
 *
 * exceptions: GSL_EUNDRFLW, GSL_EOVRFLW
 */
GSL_FUN int     gsl_sf_fermi_dirac_half_e(const double x, gsl_sf_result * result);
GSL_FUN double     gsl_sf_fermi_dirac_half(const double x);

/* Complete integral F_{3/2}(x)
 *
 * exceptions: GSL_EUNDRFLW, GSL_EOVRFLW
 */
GSL_FUN int     gsl_sf_fermi_dirac_3half_e(const double x, gsl_sf_result * result);
GSL_FUN double     gsl_sf_fermi_dirac_3half(const double x);

/* Incomplete integral F_0(x,b) = ln(1 + e^(b-x)) - (b-x)
 *
 * exceptions: GSL_EUNDRFLW, GSL_EDOM
 */
GSL_FUN int     gsl_sf_fermi_dirac_inc_0_e(const double x, const double b, gsl_sf_result * result);
GSL_FUN double     gsl_sf_fermi_dirac_inc_0(const double x, const double b);

__END_DECLS

#endif /* __GSL_SF_FERMI_DIRAC_H__ */

max_stars_repo_name: berkus/music-cs

max_issues_repo_name: berkus/music-cs

max_forks_repo_name: berkus/music-cs


text:
/*
* -----------------------------------------------------------------
*  thrm_lib.h
*  Thermochemistry Library
*  Version: 2.0
*  Last Update: Nov 3, 2019
*
*  Programmer: Americo Barbosa da Cunha Junior
*              americo.cunhajr@gmail.com
* -----------------------------------------------------------------
*  Copyright (c) 2010-2019, Americo Barbosa da Cunha Junior
*  All rights reserved.
* -----------------------------------------------------------------
*  This is the header file for THRM_LIB module, a computational
*  library with thermochemistry routines.
* -----------------------------------------------------------------
*/

#ifndef __THRM_LIB_H__
#define __THRM_LIB_H__

#include <gsl/gsl_vector.h>
#include <nvector/nvector_serial.h>

/*
* ------------------------------------------------------------
*  Types: struct thrm_struct, thrm_wrk
* ------------------------------------------------------------
* This structure contains fields of the Chemkin workspace.
*
* last update: Mar 5, 2009
* ------------------------------------------------------------
*/

typedef struct thrm_struct
{
    int    *iwk;   /* integer work vector */
    double *rwk;   /* real work vector */
    char   *cwk;   /* character work vector */
    int leniwk;    /* iwk dimension */
    int lenrwk;    /* rwk dimension */
    int lencwk;    /* cwk dimension */
    int n_e;       /* # of chemical elements */
    int n_s;       /* # of elementary species */
    int n_r;       /* # of elementary reactions */
} thrm_wrk;
/*------------------------------------------------------------*/

/*
*------------------------------------------------------------
* function prototypes
*------------------------------------------------------------
*/

void thrm_title();

void thrm_mech(int n_e,
               int n_s,
               int n_r,
               char *mech);

int thrm_input(int *steps,
               double *t0,
               double *delta_t,
               double *Tmax,
               double *Tmin,
               double *tol);

thrm_wrk* thrm_alloc();

void thrm_free(void **thrm_bl);

int thrm_init(thrm_wrk *thrm);

void thrm_composition(thrm_wrk *thrm,
                      gsl_vector *phi);

void thrm_page_heading0(FILE *file);

double thrm_h2T(thrm_wrk *thrm,
                double *phi,
                double Tmax,
                double Tmin,
                double tol);

void thrm_temp_meanvar(int Np,
                       gsl_vector **ps,
                       thrm_wrk *thrm,
                       double Tmax,
                       double Tmin,
                       double tol,
                       double *mean,
                       double *var);

double thrm_rrsum(gsl_vector *phi);

double thrm_mfsum(gsl_vector *phi);

int thrm_mfsign(gsl_vector *phi);

int thrm_eqs(double t,
             N_Vector phi,
             N_Vector Rphi,
             void *thrm_data);

int thrm_eqs2(double t,
              N_Vector phi,
              N_Vector Sphi,
              void *thrm_data);

#endif /* __THRM_LIB_H__ */

max_stars_repo_name: americocunhajr/CRFlowLib

max_issues_repo_name: americocunhajr/CRFlowLib

max_forks_repo_name: americocunhajr/CRFlowLib


text:
#ifndef VolViz_AtomicCache_h
#define VolViz_AtomicCache_h

#include <gsl/gsl>

#include <atomic>

#pragma clang diagnostic ignored "-Wpadded"

namespace VolViz {

/// Class represening a cached value, that has a very basic thread safety.
/// The thread saftey is as follows: the cache'd value must not be accessed
/// from different threads (at least without external synchronization), but
/// the cache can be marked dirty from any thread.
/// A typical use case is to cache a derived property that depends on some
/// shared (beween threads) resources (e.g. AtomicWrapper<> objects) and the
/// cached value is frequently accessed but the dependent resource does not
/// change so often.
///
/// @tparam T the cached type
template <class T> class AtomicCache {
public:
  /// Creates a dirty cache with the given fetch operation.
  /// @tparam FetchOp the fetch operation. It takes no parameter and returns an
  /// object convertible to T.
  /// @note At construction no fetch operation is issued.
  template <class FetchOp>
  AtomicCache(FetchOp fetchOperation)
      : fetchOperation_(fetchOperation) {
    Expects(fetchOperation_);
    dirtyFlag_.clear();
  }

  /// Marks the cache as dirty, i.e. a fetch operations must be perfomed on the
  /// next read.
  /// @note this method is thread safe.
  inline void markAsDirty() noexcept { dirtyFlag_.clear(); }

  /// Returns the cached object. If the cache is dirty, a fetch operations is
  /// performed first.
  /// @note This method must not be called from more than one thread.
  inline operator T const &() const noexcept {
    if (dirtyFlag_.test_and_set()) return value_;

    Expects(fetchOperation_);
    value_ = std::move(fetchOperation_());
    return value_;
  }

private:
  /// The caced value
  mutable T value_;

  /// The atomic dirty flag
  mutable std::atomic_flag dirtyFlag_;

  /// The fetch operation
  std::function<T()> const fetchOperation_;
};

} // namespace VolViz

#endif // VolViz_AtomicCache_h

max_stars_repo_name: ithron/VolViz

max_issues_repo_name: ithron/VolViz

max_forks_repo_name: ithron/VolViz


text:
/* linalg/householdercomplex.c
 *
 * Copyright (C) 2001 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_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_complex_math.h>

#include "gsl_linalg.h"

gsl_complex
gsl_linalg_complex_householder_transform (gsl_vector_complex * v)
{
  /* replace v[0:n-1] with a householder vector (v[0:n-1]) and
     coefficient tau that annihilate v[1:n-1] */

  const size_t n = v->size ;

  if (n == 1)
    {
      gsl_complex alpha = gsl_vector_complex_get (v, 0) ;
      double absa = gsl_complex_abs (alpha);
      double beta_r = - (GSL_REAL(alpha) >= 0 ? +1 : -1) * absa ;

      gsl_complex tau;
      GSL_REAL(tau) = (beta_r - GSL_REAL(alpha)) / beta_r ;
      GSL_IMAG(tau) = - GSL_IMAG(alpha) / beta_r ;

      {
        gsl_complex beta = gsl_complex_rect (beta_r, 0.0);
        gsl_vector_complex_set (v, 0, beta) ;
      }

      return tau;
    }
  else
    {
      gsl_complex tau ;
      double beta_r;

      gsl_vector_complex_view x = gsl_vector_complex_subvector (v, 1, n - 1) ;
      gsl_complex alpha = gsl_vector_complex_get (v, 0) ;
      double absa = gsl_complex_abs (alpha);
      double xnorm = gsl_blas_dznrm2 (&x.vector);

      if (xnorm == 0 && GSL_IMAG(alpha) == 0)
        {
          gsl_complex zero = gsl_complex_rect(0.0, 0.0);
          return zero; /* tau = 0 */
        }

      beta_r = - (GSL_REAL(alpha) >= 0 ? +1 : -1) * hypot(absa, xnorm) ;

      GSL_REAL(tau) = (beta_r - GSL_REAL(alpha)) / beta_r ;
      GSL_IMAG(tau) = - GSL_IMAG(alpha) / beta_r ;

      {
        gsl_complex amb = gsl_complex_sub_real(alpha, beta_r);
        gsl_complex s = gsl_complex_inverse(amb);
        gsl_blas_zscal (s, &x.vector);
      }

      {
        gsl_complex beta = gsl_complex_rect (beta_r, 0.0);
        gsl_vector_complex_set (v, 0, beta) ;
      }

      return tau;
    }
}

int
gsl_linalg_complex_householder_hm (gsl_complex tau, const gsl_vector_complex * v, gsl_matrix_complex * A)
{
  /* applies a householder transformation v,tau to matrix m */

  size_t i, j;

  if (GSL_REAL(tau) == 0.0 && GSL_IMAG(tau) == 0.0)
    {
      return GSL_SUCCESS;
    }

  /* w = (v' A)^T */

  for (j = 0; j < A->size2; j++)
    {
      gsl_complex tauwj;
      gsl_complex wj = gsl_matrix_complex_get(A,0,j);

      for (i = 1; i < A->size1; i++)  /* note, computed for v(0) = 1 above */
        {
          gsl_complex Aij = gsl_matrix_complex_get(A,i,j);
          gsl_complex vi = gsl_vector_complex_get(v,i);
          gsl_complex Av = gsl_complex_mul (Aij, gsl_complex_conjugate(vi));
          wj = gsl_complex_add (wj, Av);
        }

      tauwj = gsl_complex_mul (tau, wj);

      /* A = A - v w^T */

      {
        gsl_complex A0j = gsl_matrix_complex_get (A, 0, j);
        gsl_complex Atw = gsl_complex_sub (A0j, tauwj);
        /* store A0j - tau  * wj */
        gsl_matrix_complex_set (A, 0, j, Atw);
      }

      for (i = 1; i < A->size1; i++)
        {
          gsl_complex vi = gsl_vector_complex_get (v, i);
          gsl_complex tauvw = gsl_complex_mul(vi, tauwj);
          gsl_complex Aij = gsl_matrix_complex_get (A, i, j);
          gsl_complex Atwv = gsl_complex_sub (Aij, tauvw);
          /* store Aij - tau * vi * wj */
          gsl_matrix_complex_set (A, i, j, Atwv);
        }
    }

  return GSL_SUCCESS;
}


text:
/*
 *
 * A cross-platform way of providing access to gsl_ieee_env_setup
 * from Fortran on any system MITgcm runs on currently.
 */

#ifdef USE_GSL_IEEE
#include <gsl/gsl_math.h>
#include <gsl/gsl_ieee_utils.h>

void fgsl_ieee_env_setup ()
{
  gsl_ieee_env_setup ();
}

void fgsl_ieee_env_setup_ ()
{
  gsl_ieee_env_setup ();
}

void fgsl_ieee_env_setup__ ()
{
  gsl_ieee_env_setup ();
}

void FGSL_IEEE_ENV_SETUP ()
{
  gsl_ieee_env_setup ();
}
#endif

max_stars_repo_name: ElizabethYankovsky/MITgcm

max_issues_repo_name: ElizabethYankovsky/MITgcm

max_forks_repo_name: ElizabethYankovsky/MITgcm


text:
#ifndef __SphericalHarmonicsInterpolationFilter_h
#define __SphericalHarmonicsInterpolationFilter_h

// this code builds on the ImageToImageFilter.
#include "itkImageToImageFilter.h"
#include "itkImage.h"

// Legendre code from GNU Scientific Library
// #include <gsl/gsl_sf_legendre.h>

#include <string.h>
#include <math.h>
#include <vector>
#include <set>

// itk includes go here
#include <itkImage.h>
#include <itkImageFileReader.h>
#include <itkImageFileWriter.h>
#include <itkExceptionObject.h>
#include <itkMetaDataObject.h>
#include <itkNrrdImageIO.h>
#include <itkNthElementImageAdaptor.h>
#include <itkImageToVectorImageFilter.h>
#include <itkResampleImageFilter.h>
#include <itkVersion.h>

#define DIMENSION 3
#define PI 3.14159265358

namespace itk
{
/** \class SphericalHarmonicsInterpolationFilter
 *
 * \sa Image
 */
template <class TInputImage, class TOutputImage>
class ITK_EXPORT SphericalHarmonicsInterpolationFilter :
  public         ImageToImageFilter<TInputImage, TOutputImage>
{
public:
  /** Extract dimension from input image. */
  itkStaticConstMacro(InputImageDimension, unsigned int,
                      TInputImage::ImageDimension);
  itkStaticConstMacro(OutputImageDimension, unsigned int,
                      TInputImage::ImageDimension);

  /** Convenient typedefs for simplifying declarations. */
  typedef TInputImage                                      InputImageType;
  typedef typename    InputImageType::Pointer              InputImagePointer;
  typedef typename    InputImageType::RegionType           InputImageRegionType;
  typedef typename    InputImageType::PixelType            InputImagePixelType; // variable sized vector of valuetype
  typedef typename    InputImageType::PixelType::ValueType ImageValueType;      // USE THIS TO DETERMINE INDIVIDUAL
                                                                                // VECTOR VALUES
  typedef TOutputImage                                     OutputImageType;
  typedef typename     OutputImageType::Pointer            OutputImagePointer;
  typedef typename     OutputImageType::RegionType         OutputImageRegionType;
  typedef typename     OutputImageType::PixelType          OutputImagePixelType;

  /** Standard class typedefs. */
  typedef SphericalHarmonicsInterpolationFilter               Self;
  typedef ImageToImageFilter<InputImageType, OutputImageType> Superclass;
  typedef SmartPointer<Self>                                  Pointer;
  typedef SmartPointer<const Self>                            ConstPointer;

  /** Typedef to describe the output and input image index and size types. */
  typedef typename TOutputImage::IndexType OutputImageIndexType;
  typedef typename TInputImage::IndexType  InputImageIndexType;
  typedef typename TOutputImage::SizeType  OutputImageSizeType;
  typedef typename TInputImage::SizeType   InputImageSizeType;
  typedef InputImageSizeType               SizeType;

  /** Method for creation through the object factory. */
  itkNewMacro(Self);

  /** Run-time type information (and related methods). */
  itkTypeMacro(SphericalHarmonicsInterpolationFilter, ImageToImageFilter);

  // TODO:
  /** Extra public functions go here */
  void SetOrder(int neworder)
  {
    order = neworder;
  };
  void SetLambda(double lmbda)
  {
    lambda = lmbda;
  };
  void SetNewGradients(vnl_matrix<double> & newgrads)
  {
    newgvectors = newgrads; numPixelComponents = (newgrads.rows() + countBaselines() );
  };
  // void SetNewGradients(vnl_matrix<double> &newgrads) { newgvectors = newgrads; };
  int GetBValue()
  {
    return bvalue;
  };
  int GetNumBaselines()
  {
    return num_baselines;
  };

#ifdef ITK_USE_CONCEPT_CHECKING
  /** Begin concept checking */
  itkConceptMacro(SameDimensionCheck,
                  (Concept::SameDimension<InputImageDimension, OutputImageDimension> ) );
  itkConceptMacro(InputConvertibleToOutputCheck,
                  (Concept::Convertible<InputImagePixelType, OutputImagePixelType> ) );
  /** End concept checking */
#endif
protected:
  SphericalHarmonicsInterpolationFilter()
  {
    // newgvectors = NULL;
    newgvectors.clear();
    order = -1;
    bvalue = -1;
    numPixelComponents = -1;
    lambda = 0;
  }

  virtual ~SphericalHarmonicsInterpolationFilter()
  {
  }

  void PrintSelf(std::ostream& os, Indent indent) const;

  void GenerateOutputInformation();

  void GenerateInputRequestedRegion() throw (InvalidRequestedRegionError);

  void GenerateData();

private:
  SphericalHarmonicsInterpolationFilter(const Self &); // purposely not implemented
  void operator=(const Self &);                        // purposely not implemented

  // extra private functions go here
  double factorial(int n);

  int findMaxOrder(int numgvectors);

  int getNumTerms(int orderval);

  void cart2sph(vnl_matrix<double> & gvectors, vnl_matrix<double> & sgvectors );

  double getSphericalHarmonic(int l, int m, double theta, double phi);

  double getBasisMatrixValue(int curterm, double theta, double phi);

  void generateSHBasisMatrix(int numterms, int numgradients, vnl_matrix<double> & sgradients,
                             vnl_matrix<double> & basismat);

  int countBaselines();

  // extra private variables go here
  vnl_matrix<double> newgvectors;
  int                order;
  double             lambda;
  int                bvalue;
  int                numPixelComponents;
  int                num_baselines;
  std::set<int>      baseline_indices;

};

} // end namespace itk

#ifndef ITK_MANUAL_INSTANTIATION
#include "SphericalHarmonicsInterpolationFilter.txx"
#endif

#endif

max_stars_repo_name: BRAINSia/DTIProcessToolkit

max_issues_repo_name: BRAINSia/DTIProcessToolkit

max_forks_repo_name: BRAINSia/DTIProcessToolkit


text:
/***************************************************************************
 *   Copyright (C) 2007 by Mikhail Zaslavskiy   *
 *      *
 *                                                                         *
 *   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.             *
 ***************************************************************************/
#ifndef GRAPH_H
#define GRAPH_H

#include <iostream>
#include "rpc.h"
#include <math.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_linalg.h>
#include <vector>

/**
Graph class, it implements graph matching algorithms

	@author Mikhail Zaslavskiy
*/
class graph : public rpc
{
public:
    //class constructor
    graph(std::string fconfig="config.txt");
    graph(graph & _graph);
    const graph & operator =(graph &);

    //load graph routine
    int load_graph(std::string fgraph,char ftype='A',char cformat='D',std::string fvertexlst_name="");
    //adjacency matrix control
    const gsl_matrix*  get_adjmatrix(){return gm_A;};
    int set_adjmatrix(const gsl_matrix* _gm_A);
    gsl_matrix* get_descmatrix(char dmt);
    int printout(std::string fname_out);
    void printdot(std::string fname_out, gsl_matrix* =NULL);
    ~graph();
    int add_dummy_nodes(int id);
    long long getN();
protected:
    gsl_matrix* gm_A;//adjacency matrix
    long long N;
};
double gsl_matrix_norm(const gsl_matrix* gm,double p);
int gsl_matrix_abs(gsl_matrix* gm);
double gsl_matrix_min(gsl_matrix* gm_A,double ic);//minimum greater than ic
double gsl_vector_sum(gsl_vector* gv);
int gsl_matrix_printout(const gsl_matrix * gm,std::string,std::string);
int gsl_matrix_printout(gsl_vector * gv,std::string,std::string);
int gsl_matrix_printout(gsl_permutation * gv,std::string,std::string);
int gsl_matrix_printout(std::string sout,std::string );
double abs(double);
void gsl_matrix_sum(gsl_matrix* A,int idim,gsl_vector* gv_res);
double gsl_matrix_sum(gsl_matrix* A);
double gsl_matrix_max_abs(gsl_matrix* A);
double gsl_matrix_min_abs(gsl_matrix* gm);
//global min in [0,1]^2 of a quadratic function
void qglob_min(double a_11,double a_12,double a_22,double a_1,double a_2,double &dalpha_1,double &dalpha_2);
#endif


text:
/**
 * @file coroutine/unix.h
 * @author github.com/luncliff (luncliff@gmail.com)
 * @copyright CC BY 4.0
 */
#ifndef COROUTINE_SYSTEM_WRAPPER_H
#define COROUTINE_SYSTEM_WRAPPER_H
#if !(defined(unix) || defined(__APPLE__) || defined(__FreeBSD__))
#error "expect UNIX platform for this file"
#endif
#include <sys/event.h> // for BSD kqueue

#include <gsl/gsl>

#include <coroutine/return.h>

/**
 * @defgroup BSD
 */

namespace coro {

/**
 * @brief RAII wrapping for kqueue file descriptor
 * @ingroup BSD
 */
class kqueue_owner final {
    int64_t kqfd;

  public:
    /**
     * @brief create a fd with `kqueue`. Throw if the function fails.
     * @see kqeueue
     * @throw system_error
     */
    kqueue_owner() noexcept(false);
    /**
     * @brief close the current kqueue file descriptor
     */
    ~kqueue_owner() noexcept;
    kqueue_owner(const kqueue_owner&) = delete;
    kqueue_owner(kqueue_owner&&) = delete;
    kqueue_owner& operator=(const kqueue_owner&) = delete;
    kqueue_owner& operator=(kqueue_owner&&) = delete;

  public:
    /**
     * @brief bind the event to kqueue
     * @param req
     * @see kevent64
     * @throw system_error
     *
     * The function is named `change` because
     * the given argument is used for 'change list' fo `kqueue64`
     */
    void change(kevent64_s& req) noexcept(false);

    /**
     * @brief fetch all events for the given kqeueue descriptor
     * @param wait_time
     * @param list
     * @return ptrdiff_t
     * @see kevent64
     * @throw system_error
     *
     * The function is named `events` because
     * the given argument is used for 'event list' fo `kqueue64`
     *
     * Timeout is not an error for this function
     */
    ptrdiff_t events(const timespec& wait_time,
                     gsl::span<kevent64_s> list) noexcept(false);

  public:
    /**
     * @brief return temporary awaitable object for given event
     * @param req input for `change` operation
     * @see change
     *
     * There is no guarantee of reusage of returned awaiter object
     * When it is awaited, and `req.udata` is null(0),
     * the value is set to `coroutine_handle<void>`
     *
     * @code
     * auto read_async(kqueue_owner& kq, uint64_t fd) -> frame_t {
     *     kevent64_s req{.ident = fd,
     *                    .filter = EVFILT_READ,
     *                    .flags = EV_ADD | EV_ENABLE | EV_ONESHOT};
     *     co_await kq.submit(req);
     *     // ...
     *     co_await kq.submit(req);
     *     // ...
     * }
     * @endcode
     */
    [[nodiscard]] auto submit(kevent64_s& req) noexcept {
        class awaiter final : public suspend_always {
            kqueue_owner& kq;
            kevent64_s& req;

          public:
            constexpr awaiter(kqueue_owner& _kq, kevent64_s& _req)
                : kq{_kq}, req{_req} {
            }

          public:
            void await_suspend(coroutine_handle<void> coro) noexcept(false) {
                if (req.udata == 0)
                    req.udata = reinterpret_cast<uint64_t>(coro.address());
                return kq.change(req);
            }
        };
        return awaiter{*this, req};
    }
};

} // namespace coro

#endif // COROUTINE_SYSTEM_WRAPPER_H

max_stars_repo_name: lanza/coroutine

max_issues_repo_name: lanza/coroutine

max_forks_repo_name: lanza/coroutine


text:
/* Copyright (c) 2016, 2017 Bradley Worley <geekysuavo@gmail.com>
 * Released under the MIT License
 */

/* ensure once-only inclusion. */
#ifndef __MATTE_OBJECT_H__
#define __MATTE_OBJECT_H__

/* include required standard c library headers. */
#include <stdio.h>
#include <stdlib.h>
#include <stdarg.h>
#include <stddef.h>
#include <string.h>
#include <stdbool.h>

/* include required c math library headers. */
#include <math.h>
#include <complex.h>

/* include the blas/lapack headers. */
#include <cblas.h>
#include <clapack.h>

/* include the required matte headers. */
#include <matte/zone.h>

/* type definitions of blas/lapack enumerations:
 *  MatteTranspose: transpose options.
 *  MatteTriangle: triangle options.
 *  MatteDiagonal: diagonal options.
 */
typedef enum CBLAS_TRANSPOSE MatteTranspose;
typedef enum CBLAS_UPLO MatteTriangle;
typedef enum CBLAS_DIAG MatteDiagonal;

/* MATTE_TYPE: macro to obtain the type structure pointer of an object.
 */
#define MATTE_TYPE(obj) \
  (((Object) obj)->type)

/* MATTE_TYPE_CHECK: macro to check that a matte object has a specific type.
 */
#define MATTE_TYPE_CHECK(obj,typ) \
  ((obj) && MATTE_TYPE(obj) == typ)

/* OBJECT_BASE: macro to place at the beginning of every matte object.
 * this allows all objects to be directly cast into the base object
 * type for memory management, dynamic dispatch, etc.
 */
#define OBJECT_BASE \
  Object base;

/* ObjectType: pointer to a struct _ObjectType. */
typedef struct _ObjectType *ObjectType;

/* Object: pointer to a struct _Object. */
typedef struct _Object *Object;

/* ObjectMethods: pointer to a struct _ObjectMethod. */
typedef struct _ObjectMethod *ObjectMethods;

/* general-purpose function pointer type definition:
 */
typedef Object (*matte_func) (Zone, Object);

/* function pointer type definitions for matte objects:
 */
typedef Object (*obj_constructor) (Zone, Object);
typedef void   (*obj_destructor)  (Zone, Object);
typedef int    (*obj_display)     (Zone, Object);
typedef int    (*obj_assert)      (Object);
typedef Object (*obj_unary)       (Zone, Object);
typedef Object (*obj_binary)      (Zone, Object, Object);
typedef Object (*obj_ternary)     (Zone, Object, Object, Object);
typedef Object (*obj_variadic)    (Zone, int, va_list);
typedef Object (*obj_method)      (Zone, Object, Object);

/* _ObjectMethod: structure that holds information about a matte object
 * method.
 */
struct _ObjectMethod {
  /* @name: string name of the method.
   * @fn: method function pointer.
   */
  const char *name;
  obj_method fn;
};

/* _ObjectType: structure that holds the core set of information required
 * by the matte object system.
 */
struct _ObjectType {
  /* core object information:
   *  @name: string name of the matte type.
   *  @size: base size of the type struct.
   *  @precedence: object precedence for dynamic dispatch.
   */
  const char *name;
  unsigned int size;
  unsigned int precedence;

  /* core object method table:
   */
  obj_constructor fn_new;
  obj_constructor fn_copy;
  obj_destructor  fn_delete;
  obj_display     fn_disp;
  obj_assert      fn_true;

  /* numeric method table:
   */
  obj_binary   fn_plus;        /* a+b      binary addition.             */
  obj_binary   fn_minus;       /* a-b      binary subtraction.          */
  obj_unary    fn_uminus;      /* -a       unary minus.                 */
  obj_binary   fn_times;       /* a.*b     element-wise multiplication. */
  obj_binary   fn_mtimes;      /* a*b      matrix multiplication.       */
  obj_binary   fn_rdivide;     /* a./b     element-wise right division. */
  obj_binary   fn_ldivide;     /* a.\b     element-wise left division.  */
  obj_binary   fn_mrdivide;    /* a/b      matrix right division.       */
  obj_binary   fn_mldivide;    /* a\b      matrix left division.        */
  obj_binary   fn_power;       /* a.^b     element-wise power.          */
  obj_binary   fn_mpower;      /* a^b      matrix power.                */
  obj_binary   fn_lt;          /* a<b      less than.                   */
  obj_binary   fn_gt;          /* a>b      greater than.                */
  obj_binary   fn_le;          /* a<=b     less than or equal to.       */
  obj_binary   fn_ge;          /* a>=b     greater than or equal to.    */
  obj_binary   fn_ne;          /* a!=b     inequality.                  */
  obj_binary   fn_eq;          /* a==b     equality.                    */
  obj_binary   fn_and;         /* a&b      logical and.                 */
  obj_binary   fn_or;          /* a|b      logical or.                  */
  obj_binary   fn_mand;        /* a&&b     matrix logical and.          */
  obj_binary   fn_mor;         /* a||b     matrix logical or.           */
  obj_unary    fn_not;         /* !a       logical negation.            */
  obj_ternary  fn_colon;       /* a:d:b    colon operator.              */
  obj_unary    fn_ctranspose;  /* a'       conjugate transpose.         */
  obj_unary    fn_transpose;   /* a.'      matrix transpose.            */
  obj_variadic fn_horzcat;     /* [a,b]    horizontal concatenation.    */
  obj_variadic fn_vertcat;     /* [a;b]    vertical concatenation.      */
  obj_binary   fn_subsref;     /* a(s)     subscripted reference.       */
  obj_ternary  fn_subsasgn;    /* a(s)=b   subscripted assignment.      */
  obj_unary    fn_subsindex;   /* b(a)     subscript index.             */

  /* general-purpose method table:
   */
  ObjectMethods methods;
};

/* _Object: structure for type-casting between any number of matte objects.
 * holds a pointer to the object type structure as its first element.
 */
struct _Object {
  /* @type: the type of matte object that is located
   * at the current address in memory.
   */
  ObjectType type;

  /* in all other matte objects, object instance variables go here. */
};

/* function declarations (object.c): */

Object object_alloc (Zone z, ObjectType type);

Object object_copy (Zone z, Object obj);

void object_free (Zone z, void *ptr);

void object_free_all (Zone z);

int object_disp (Zone z, Object obj);

int object_display (Zone z, Object obj, const char *var);

int object_true (Object obj);

/* object method declarations (object.c): */

Object object_plus       (Zone z, Object a, Object b);
Object object_minus      (Zone z, Object a, Object b);
Object object_uminus     (Zone z, Object a);
Object object_times      (Zone z, Object a, Object b);
Object object_mtimes     (Zone z, Object a, Object b);
Object object_rdivide    (Zone z, Object a, Object b);
Object object_ldivide    (Zone z, Object a, Object b);
Object object_mrdivide   (Zone z, Object a, Object b);
Object object_mldivide   (Zone z, Object a, Object b);
Object object_power      (Zone z, Object a, Object b);
Object object_mpower     (Zone z, Object a, Object b);
Object object_lt         (Zone z, Object a, Object b);
Object object_gt         (Zone z, Object a, Object b);
Object object_le         (Zone z, Object a, Object b);
Object object_ge         (Zone z, Object a, Object b);
Object object_ne         (Zone z, Object a, Object b);
Object object_eq         (Zone z, Object a, Object b);
Object object_and        (Zone z, Object a, Object b);
Object object_or         (Zone z, Object a, Object b);
Object object_mand       (Zone z, Object a, Object b);
Object object_mor        (Zone z, Object a, Object b);
Object object_not        (Zone z, Object a);
Object object_colon      (Zone z, Object a, Object b, Object c);
Object object_ctranspose (Zone z, Object a);
Object object_transpose  (Zone z, Object a);
Object object_horzcat    (Zone z, int n, ...);
Object object_vertcat    (Zone z, int n, ...);
Object object_subsref    (Zone z, Object a, Object b);
Object object_subsasgn   (Zone z, Object a, Object b, Object c);
Object object_subsindex  (Zone z, Object a);

/* utility function declarations: */

char *strdup (const char *s);

#endif /* !__MATTE_OBJECT_H__ */

max_stars_repo_name: geekysuavo/matte

max_issues_repo_name: geekysuavo/matte

max_forks_repo_name: geekysuavo/matte


text:
#pragma once
#include "utility/scope_guard.h"
#include "utility/types.h"
#include <gsl/gsl>
#include <cfenv>
#include <cassert>

#define scoped_fesetround(mode)                                  \
    int SCOPE_GUARD_PP_UNIQUE(_rounding_mode) = fegetround();    \
    int SCOPE_GUARD_PP_UNIQUE(_rounding_ret) = fesetround(mode); \
    (void)SCOPE_GUARD_PP_UNIQUE(_rounding_ret);                  \
    assert(SCOPE_GUARD_PP_UNIQUE(_rounding_ret) == 0);           \
    SCOPE(exit) { fesetround(SCOPE_GUARD_PP_UNIQUE(_rounding_mode)); }

template <class T> inline T square(T a) { return a * a; }

template <class T> inline T cube(T a) { return a * a * a; }

template <class T> inline T clamp(T x, T min, T max)
{
    x = (x < min) ? min : x;
    x = (x > max) ? max : x;
    return x;
}

// linear interpolator
template <class R> R itp_linear(const R y[], R mu)
{
    return y[0] * (1 - mu) + y[1] * mu;
}

// Catmull-Rom interpolator
template <class R> R itp_catmull(const R y[], R mu)
{
    R mu2 = mu * mu;
    R mu3 = mu2 * mu;
    R a[] = {-R(0.5) * y[0] + R(1.5) * y[1] - R(1.5) * y[2] + R(0.5) * y[3],
             y[0] - R(2.5) * y[1] + R(2) * y[2] - R(0.5) * y[3],
             -R(0.5) * y[0] + R(0.5) * y[2], y[1]};
    return a[0] * mu3 + a[1] * mu2 + a[2] * mu + a[3];
}

//------------------------------------------------------------------------------
// Fixed point routines
//   suffix 'xN', N = bit size of the integer part
//   prefix 'l' = 64-bit operation, 32-bit otherwise
//
//   'fx' number to fixed
//   'ix' integer part, 'rx' fractional part
//   'ffx' fixed to f32, 'dfx' fixed to f64

// number to Q32,32
template <class R> inline constexpr i64 lfx32(R x)
{
    return i64(x * (i64)4294967296);
}
// number to Q16,16
template <class R> inline constexpr i32 fx16(R x)
{
    return i32(x * (i32)65536);
}
// number to Q16,48
template <class R> inline constexpr i64 lfx16(R x)
{
    return i64(x * (i64)281474976710656);
}
// number to Q8,24
template <class R> inline constexpr i32 fx8(R x)
{
    return i32(x * (i32)16777216);
}
// number to Q8,56
template <class R> inline constexpr i64 lfx8(R x)
{
    return i64(x * (i64)72057594037927936);
}

// Q32,32 integer part
inline constexpr i64 lix32(i64 x) { return x / 4294967296; }
inline constexpr u64 lix32(u64 x) { return x / 4294967296; }
template <class T> T lix32(T) = delete;
// Q16,16 integer part
inline constexpr i32 ix16(i32 x) { return x / 65536; }
inline constexpr u32 ix16(u32 x) { return x / 65536; }
template <class T> T ix16(T) = delete;
// Q16,48 integer part
inline constexpr i64 lix16(i64 x) { return x / 281474976710656; }
inline constexpr u64 lix16(u64 x) { return x / 281474976710656; }
template <class T> T lix16(T) = delete;
// Q8,24 integer part
inline constexpr i32 ix8(i32 x) { return x / 16777216; }
inline constexpr u32 ix8(u32 x) { return x / 16777216; }
template <class T> T ix8(T) = delete;
// Q8,56 integer part
inline constexpr i64 lix8(i64 x) { return x / 72057594037927936; }
inline constexpr u64 lix8(u64 x) { return x / 72057594037927936; }
template <class T> T lix8(T) = delete;

// Q32,32 fractional part
inline constexpr u64 lrx32(i64 x) { return x & 4294967295; }
inline constexpr u64 lrx32(u64 x) { return x & 4294967295; }
template <class T> T lrx32(T) = delete;
// Q16,16 fractional part
inline constexpr u32 rx16(i32 x) { return x & 65535; }
inline constexpr u32 rx16(u32 x) { return x & 65535; }
template <class T> T rx16(T) = delete;
// Q16,48 fractional part
inline constexpr u64 lrx16(i64 x) { return x & 281474976710655; }
inline constexpr u64 lrx16(u64 x) { return x & 281474976710655; }
template <class T> T lrx16(T) = delete;
// Q8,24 fractional part
inline constexpr u32 rx8(i32 x) { return x & 16777215; }
inline constexpr u32 rx8(u32 x) { return x & 16777215; }
template <class T> T rx8(T) = delete;
// Q8,56 fractional part
inline constexpr u64 lrx8(i64 x) { return x & 72057594037927935; }
inline constexpr u64 lrx8(u64 x) { return x & 72057594037927935; }
template <class T> T lrx8(T) = delete;

// Q32,32 number to real
inline constexpr f32 lffx32(i64 x) { return x / 4294967296.0f; }
inline constexpr f32 lffx32(u64 x) { return x / 4294967296.0f; }
template <class T> f32 lffx32(T x) = delete;
inline constexpr f64 ldfx32(i64 x) { return x / 4294967296.0; }
inline constexpr f64 ldfx32(u64 x) { return x / 4294967296.0; }
template <class T> f64 ldfx32(T x) = delete;
// Q16,16 number to real
inline constexpr f32 ffx16(i32 x) { return x / 65536.0f; }
inline constexpr f32 ffx16(u32 x) { return x / 65536.0f; }
template <class T> f32 ffx16(T x) = delete;
inline constexpr f64 dfx16(i32 x) { return x / 65536.0; }
inline constexpr f64 dfx16(u32 x) { return x / 65536.0; }
template <class T> f64 dfx16(T x) = delete;
// Q16,48 number to real
inline constexpr f32 lffx16(i64 x) { return x / 281474976710656.0f; }
inline constexpr f32 lffx16(u64 x) { return x / 281474976710656.0f; }
template <class T> f32 lffx16(T x) = delete;
inline constexpr f64 ldfx16(i64 x) { return x / 281474976710656.0; }
inline constexpr f64 ldfx16(u64 x) { return x / 281474976710656.0; }
template <class T> f64 ldfx16(T x) = delete;
// Q8,24 number to real
inline constexpr f32 ffx8(i32 x) { return x / 16777216.0f; }
inline constexpr f32 ffx8(u32 x) { return x / 16777216.0f; }
template <class T> f32 ffx8(T x) = delete;
inline constexpr f64 dfx8(i32 x) { return x / 16777216.0; }
inline constexpr f64 dfx8(u32 x) { return x / 16777216.0; }
template <class T> f64 dfx8(T x) = delete;
// Q8,56 number to real
inline constexpr f32 lffx8(i64 x) { return x / 72057594037927936.0f; }
inline constexpr f32 lffx8(u64 x) { return x / 72057594037927936.0f; }
template <class T> f32 lffx8(T x) = delete;
inline constexpr f64 ldfx8(i64 x) { return x / 72057594037927936.0; }
inline constexpr f64 ldfx8(u64 x) { return x / 72057594037927936.0; }
template <class T> f64 ldfx8(T x) = delete;


text:
/* gsl_histogram_maxval.c
 * Copyright (C) 2000  Simone Piccardi
 *
 * This library 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 library; if not, write to the Free Software Foundation, Inc.,
 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */
/***************************************************************
 *
 * File gsl_histogram_maxval.c:
 * Routine to find maximum and minumum content of a hisogram.
 * Need GSL library and header.
 * Contains the routines:
 * gsl_histogram_max_val find max content values
 * gsl_histogram_min_val find min content values
 * gsl_histogram_bin_max find coordinates of max contents bin
 * gsl_histogram_bin_min find coordinates of min contents bin
 *
 * Author: S. Piccardi
 * Jan. 2000
 *
 ***************************************************************/
#include <config.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_histogram.h>

double
gsl_histogram_max_val (const gsl_histogram * h)
{
  const size_t n = h->n;
  size_t i;
  double max = h->bin[0];
  for (i = 0; i < n; i++)
    {
      if (h->bin[i] > max)
        {
          max = h->bin[i];
        }
    }
  return max;
}

size_t
gsl_histogram_max_bin (const gsl_histogram * h)
{
  size_t i;
  size_t imax = 0;
  double max = h->bin[0];
  for (i = 0; i < h->n; i++)
    {
      if (h->bin[i] > max)
        {
          max = h->bin[i];
          imax = i;
        }
    }
  return imax;
}

double
gsl_histogram_min_val (const gsl_histogram * h)
{
  size_t i;
  double min = h->bin[0];
  for (i = 0; i < h->n; i++)
    {
      if (h->bin[i] < min)
        {
          min = h->bin[i];
        }
    }
  return min;
}

size_t
gsl_histogram_min_bin (const gsl_histogram * h)
{
  size_t i;
  size_t imin = 0;
  double min = h->bin[0];
  for (i = 0; i < h->n; i++)
    {
      if (h->bin[i] < min)
        {
          min = h->bin[i];
          imin = i;
        }
    }
  return imin;
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
/*
** overall concordance correlation coefficient (OCCC)
** Barnhart et al. (2002), Biometrics, 58, 1020-1027.
**
** G.Lohmann, July 2010
*/
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_errno.h>

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

#define SQR(x) ((x) * (x))
#define ABS(x) ((x) > 0 ? (x) : -(x))

double Cov(double *arr1,double *arr2,int n)
{
  int i;
  double nx,sum,ave1,ave2;

  nx = (double)n;

  sum = 0;
  for (i=0; i<n; i++) sum += arr1[i];
  ave1 = sum/nx;

  sum = 0;
  for (i=0; i<n; i++) sum += arr2[i];
  ave2 = sum/nx;

  sum = 0;
  for (i=0; i<n; i++)
    sum += (arr1[i] - ave1)*(arr2[i] - ave2);

  return sum/(nx-1);
}

void AveVar(double *arr,int n,double *ave,double *var)
{
  int i;
  double u,sum1,sum2,nx,mean;

  nx = (double)n;
  sum1 = sum2 = 0;
  for (i=0; i<n; i++) {
    u = arr[i];
    sum1 += u;
    sum2 += u*u;
  }
  mean = sum1/nx;

  *ave = mean;
  *var = (sum2 - nx * mean * mean) / (nx - 1.0);
}

double VOCCC(gsl_matrix *data,int nvoxels,int verbose)
{
  int i,j,k;
  double sum1,sum2,sum3,mean,varx,u,v;
  static gsl_matrix *cov=NULL;
  static gsl_vector *ave=NULL,*var=NULL;

  int nsubjects = data->size1;

  if (ave == NULL) {
    ave = gsl_vector_calloc(nsubjects);
    var = gsl_vector_calloc(nsubjects);
    cov = gsl_matrix_calloc(nsubjects,nsubjects);
  }

  for (i=0; i<nsubjects; i++) {
    double *ptr1 = gsl_matrix_ptr(data,i,0);
    AveVar(ptr1,nvoxels,&mean,&varx);
    gsl_vector_set(ave,i,mean);
    gsl_vector_set(var,i,varx);

    for (j=0; j<nsubjects; j++) {
      double *ptr1 = gsl_matrix_ptr(data,i,0);
      double *ptr2 = gsl_matrix_ptr(data,j,0);
      u = Cov(ptr1,ptr2,nvoxels);
      gsl_matrix_set(cov,i,j,u);
    }
  }
  if (verbose > 0) {
    for (i=0; i<nsubjects; i++) {
      for (j=0; j<nsubjects; j++) {
	fprintf(stderr," %8.5f",gsl_matrix_get(cov,i,j));
      }
      fprintf(stderr,"\n");
    }
    fprintf(stderr,"\n");
    exit(0);
  }

  sum1 = 0;
  for (j=1; j<=nsubjects-1; j++) {
    for (k=j+1; k<=nsubjects; k++) {
      u = gsl_matrix_get(cov,j-1,k-1);
      sum1 += u;
    }
  }
  sum1 *= 2.0;

  sum2 = 0;
  for (j=0; j<nsubjects; j++) {
    sum2 += gsl_vector_get(var,j);
  }
  sum2 *= (double)(nsubjects-1);

  sum3 = 0;
  for (j=1; j<=nsubjects-1; j++) {
    for (k=j+1; k<=nsubjects; k++) {
      u = gsl_vector_get(ave,j-1);
      v = gsl_vector_get(ave,k-1);
      sum3 += (u-v)*(u-v);
    }
  }

  if (sum2 + sum3 < 1.0e-6) return 0;
  return sum1 / (sum2 + sum3);
}

max_stars_repo_name: zrajna/lipsia

max_issues_repo_name: zrajna/lipsia

max_forks_repo_name: zrajna/lipsia


text:
#pragma once

#include "gdx/algo/algorithm.h"
#include "gdx/exception.h"

#include <gsl/span>
#include <numeric>
#include <type_traits>
#include <unordered_map>

namespace gdx {

template <class T>
struct has_sum_member
{
    template <typename U, U u>
    struct reallyHas;

    template <typename C>
    static constexpr std::true_type test(reallyHas<double (C::*)(), &C::sum>* /*unused*/)
    {
        return {};
    }

    template <typename C>
    static constexpr std::true_type test(reallyHas<double (C::*)() const, &C::sum>* /*unused*/)
    {
        return {};
    }

    template <typename>
    static std::false_type test(...)
    {
        return {};
    }

    static constexpr bool value = decltype(test<T>(0))::value;
};

template <typename T>
inline constexpr bool has_sum_member_v = has_sum_member<T>::value;

template <typename RasterType>
typename std::enable_if_t<has_sum_member_v<RasterType>, double> sum(const RasterType& ras)
{
    return ras.template sum<double>();
}

template <typename RasterType>
typename std::enable_if_t<!has_sum_member_v<RasterType>, double> sum(const RasterType& ras)
{
    double sum      = 0.0;
    const auto size = ras.size();

    for (std::size_t i = 0; i < size; ++i) {
        if (!ras.is_nodata(i)) {
            sum += ras[i];
        }
    }

    return sum;
}

template <typename RasterType>
double ssum(const RasterType& ras)
{
    using T = typename RasterType::value_type;
    return static_cast<double>(std::accumulate(value_begin(ras), value_end(ras), T(0), std::plus<T>()));
}

}

max_stars_repo_name: VITObelgium/geodynamix

max_issues_repo_name: VITObelgium/geodynamix

max_forks_repo_name: VITObelgium/geodynamix


text:
#include "../include/paralleltt.h"

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

#define NUM_THREAD 1
#define min(a,b) ((a)>(b)?(b):(a))
#define VERBOSE ((int) 1)

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

#define BUF ((int) 2)

tensor_train* TT_init(const int d, const int* restrict n)
{
    tensor_train* sim = (tensor_train*) malloc(sizeof(tensor_train));
    sim->d = d;
    sim->n = (int*) malloc(sizeof(int)*d);
    sim->r = (int*) malloc(sizeof(int)*(d+1));
    sim->trains = (double**) malloc(sizeof(double*)*d);

    sim->r[0] = 1;
    for (int j = 0; j < d; ++j) {
        sim->n[j] = n[j];

        if (j < d-1) {
            sim->r[j+1] = n[j];
        }
    }
    sim->r[d] = 1;

    for (int j = 0; j < d; ++j) {
        sim->trains[j] = (double*) malloc(sizeof(double)*(sim->r[j])*(sim->n[j])*(sim->r[j+1]));
    }

    return sim;
}

tensor_train* TT_init_rank(const int d, const int* restrict n, const int* restrict r)
{
    tensor_train* sim = (tensor_train*) malloc(sizeof(tensor_train));
    sim->d = d;
    sim->n = (int*) malloc(sizeof(int)*d);
    sim->r = (int*) malloc(sizeof(int)*(d+1));
    sim->trains = (double**) malloc(sizeof(double*)*d);

    for (int j = 0; j < d; ++j) {
        sim->n[j] = n[j];
    }
    for (int j = 0; j < d+1; ++j) {
        sim->r[j] = r[j];
    }

    for (int j = 0; j < d; ++j) {
        sim->trains[j] = (double*) malloc(sizeof(double)*(sim->r[j])*(sim->n[j])*(sim->r[j+1]));
    }

    return sim;
}

tensor_train* TT_copy(tensor_train* X)
{
    int d = X->d;
    int* n = X->n;
    int* r = X->r;
    tensor_train* Y = TT_init_rank(d, n, r);
    for (int j = 0; j < X->d; ++j) {
        int N = r[j] * n[j] * r[j+1];
        memcpy(Y->trains[j], X->trains[j], N * sizeof(double));
    }

    return Y;
}

void TT_free(tensor_train* sim)
{
    for (int j = 0; j < sim->d; ++j) {
        free(sim->trains[j]); sim->trains[j] = NULL;
    }
    free(sim->trains); sim->trains = NULL;
    free(sim->r);      sim->r = NULL;
    free(sim->n);      sim->n = NULL;
    free(sim);
}

void TT_print(tensor_train* tt)
{

    int d = tt->d;
    int* n = tt->n;
    int* r = tt->r;
    double** trains = tt->trains;
    printf("Printing the tensor train with d = %d\n",d);

    for (int ii = 0; ii < d; ++ii){
        printf("\nTensor %d/d (r[%d] = %d, n[%d]=%d, r[%d]=%d):\n",ii,ii,r[ii],ii,n[ii],ii+1,r[ii+1]);
        matrix_tt* A = matrix_tt_wrap(r[ii], n[ii] * r[ii+1], trains[ii]);
        matrix_tt_print(A, 1);
        free(A); A = NULL;
    }
}

void tt_broadcast(MPI_Comm comm, tensor_train* tt){
    for (int ii = 0; ii < tt->d; ++ii){
        matrix_tt* A = matrix_tt_wrap((tt->r[ii]) * (tt->n[ii]) * (tt->r[ii+1]), 1, tt->trains[ii]);
        matrix_tt_broadcast(comm, A);
        free(A);
    }
}

double get_compression(const tensor_train* tt)
{
    long N = 1;
    for (int ii = 0; ii < tt->d; ++ii)
    {
        N = N * tt->n[ii];
    }

    long Ntt = 0;
    for (int ii = 0; ii < tt->d; ++ii)
    {
        Ntt = Ntt + ((tt->r[ii])*(tt->n[ii])*(tt->r[ii+1]));
    }

    double compression = (double) Ntt/N;
    return compression;
}


text:
#ifndef INCLUDED_BLAS_DETAIL_H
#define INCLUDED_BLAS_DETAIL_H

#include <cblas.h>
#include "telepath/blas/traits.h"

namespace blas{

    namespace matrix{

        namespace matrix_detail{

            template< typename T, bool = BlasCast<T>::value >
            struct BlasMatImpl { };

            template< typename T >
            struct BlasMatImpl< T, true >{
                using mat_t = decltype( BlasCast<T>::blas_cast(
                    std::declval<T>() ) );
                mat_t mat;
                constexpr BlasMatImpl( const T& m )
                    : mat( BlasCast<T>::blas_cast( m ) ) { }
                operator mat_t&() { return mat; }
                operator const mat_t&() const{ return mat; }
            };

            template< typename T >
            struct BlasMatImpl< T, false >{
                using mat_t = T;
                mat_t& mat;
                constexpr BlasMatImpl( T& m ) : mat( m ) { }
                operator mat_t&() { return mat; }
                operator const mat_t&() const{ return mat; }
            };

            template< typename T, bool = BlasCast<T>::value >
            struct ConstBlasMatImpl { };

            template< typename T >
            struct ConstBlasMatImpl< T, true >{
                using mat_t = decltype( BlasCast<T>::blas_cast(
                    std::declval<T>() ) );
                mat_t mat;
                constexpr ConstBlasMatImpl( const T& m )
                    : mat( BlasCast<T>::blas_cast( m ) ) { }
                operator const mat_t&() const{ return mat; }
            };

            template< typename T >
            struct ConstBlasMatImpl< T, false >{
                using mat_t = T;
                const mat_t& mat;
                constexpr ConstBlasMatImpl( const T& m ) : mat( m ) { }
                operator const mat_t&() const{ return mat; }
            };

        } //namespace matrix_detail

        template< typename T >
        struct Mat : public matrix_detail::BlasMatImpl<T> {
            using impl = typename matrix_detail::BlasMatImpl<T>;
            using impl::impl;
        };

        template< typename T >
        struct ConstMat : public matrix_detail::ConstBlasMatImpl<T>{
            using impl = typename matrix_detail::ConstBlasMatImpl<T>;
            using impl::impl;
        };

        template< typename T >
        using scalar_t = typename MatrixTraits<T>::scalar_t;

        template< typename T >
        constexpr CBLAS_ORDER layout( const T& x ){
            return MatrixTraits<T>::layout( x );
        }

        template< typename T >
        constexpr CBLAS_TRANSPOSE trans( const T& x ){
            return MatrixTraits<T>::trans( x );
        }

        template< typename T >
        constexpr CBLAS_UPLO uplo( const T& x ){
            return MatrixTraits<T>::uplo( x );
        }

        template< typename T >
        constexpr std::size_t nrows( const T& x ){
            return MatrixTraits<T>::nrows( x );
        }

        template< typename T >
        constexpr std::size_t ncols( const T& x ){
            return MatrixTraits<T>::ncols( x );
        }

        template< typename T >
        constexpr std::size_t ld( const T& x ){
            return MatrixTraits<T>::ld( x );
        }

        template< typename T >
        constexpr auto array( T& x ){
            return MatrixTraits<T>::array( x );
        }

        template< typename T >
        constexpr auto array( const T& x ){
            return MatrixTraits<T>::array( x );
        }

    } //namespace matrix

    template< typename T >
    inline auto blas_matrix( T& x ){
        return matrix::Mat<T>( x );
    }

    template< typename T >
    inline auto blas_matrix( const T& x ){
        return matrix::ConstMat<T>( x );
    }

    template< typename T >
    struct MatrixTraits< matrix::Mat<T> >
        : public MatrixTraits< typename matrix::Mat<T>::mat_t > { };

    template< typename T >
    struct MatrixTraits< matrix::ConstMat<T> >
        : public MatrixTraits< typename matrix::ConstMat<T>::mat_t > { };

    namespace vector{

        namespace vector_detail{

            template< typename T, bool = BlasCast<T>::value >
            struct BlasVecImpl { };

            template< typename T >
            struct BlasVecImpl< T, true >{
                using vec_t = decltype( BlasCast<T>::blas_cast(
                    std::declval<T>() ) );
                vec_t vec;
                constexpr BlasVecImpl( const T& v )
                    : vec( BlasCast<T>::blas_cast( v ) ) { }
                operator vec_t&() { return vec; }
                operator const vec_t&() const{ return vec; }
            };

            template< typename T >
            struct BlasVecImpl< T, false >{
                using vec_t = T;
                vec_t& vec;
                BlasVecImpl( const T& v ) : vec( v ) { }
                operator vec_t&() { return vec; }
                operator const vec_t&() const{ return vec; }
            };

        } //namespace vector_detail

        template< typename T >
        struct Vec : public vector_detail::BlasVecImpl<T> {
            using impl = typename vector_detail::BlasVecImpl<T>;
            using impl::impl;
        };

        template< typename T >
        using scalar_t = typename VectorTraits<T>::scalar_t;

        template< typename T >
        constexpr auto length( const T& x ){
            return VectorTraits<T>::length( x );
        }

        template< typename T >
        constexpr auto inc( const T& x ){
            return VectorTraits<T>::inc( x );
        }

        template< typename T >
        constexpr auto array( T& x ){
            return VectorTraits<T>::array( x );
        }

        template< typename T >
        constexpr auto array( const T& x ){
            return VectorTraits<T>::array( x );
        }

    } //namespace vector

    template< typename T >
    struct VectorTraits< vector::Vec<T> >
        : public VectorTraits< typename vector::Vec<T>::vec_t > { };

} //namespace blas

#endif

max_stars_repo_name: tbepler/telepath

max_issues_repo_name: tbepler/telepath

max_forks_repo_name: tbepler/telepath


text:
//
// Created by bessermt on 6/3/18.
//

#ifndef SRSGRADE_TIMESTAMP_H
#define SRSGRADE_TIMESTAMP_H

#include <iostream>
#include <iomanip>
#include <chrono>
#include <gsl.h>

class Timestamp
{
private:
    std::chrono::system_clock::time_point time_point_{};

    void set(const std::string date_string, const std::string time_string);

    friend std::istream& operator>>(std::istream& is, Timestamp& timestamp)
    {
        std::string date_string;

		is >> date_string;
		Ensures(!date_string.empty());

		char time_cstr[256]{};
		constexpr auto max_time_len{ std::extent<decltype(time_cstr)>::value };
		is >> std::ws;
		is.getline(time_cstr, max_time_len, ',');
		const auto time_len{ strlen(time_cstr) };
		Ensures(time_len > 0 && time_len < max_time_len);

		const auto time_string{ std::string(time_cstr) };

        timestamp.set(date_string, time_string);

        return is;
    }

	friend bool operator<(const Timestamp lhs, const Timestamp rhs)
    {
		return lhs.time_point_ < rhs.time_point_;
    }

public:
    auto is_good() const
    {
		const auto is_past{ time_point_ < std::chrono::system_clock::now() };

        return is_past; // TODO: Any better tests?
    }
};

#endif // SRSGRADE_TIMESTAMP_H

max_stars_repo_name: labermt/SRSGrade

max_issues_repo_name: labermt/SRSGrade

max_forks_repo_name: labermt/SRSGrade