text:
/* popGenTools.h
/
/
/ Andrew Kern
/
/
*/
#ifndef PGTOOLS_INC
#define PGTOOLS_INC

//#define MAXSNPS 1000000
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_min.h>
#include "mpi.h"

struct my_F_params{
	int i, n, ascSize, ascFreq, derivedFlag;
	gsl_integration_workspace * w;
};

struct trans_params{
	int xi, i, n, ascSize, ascFreq, derivedFlag;
	double t, q, p, tau,v;
	gsl_integration_workspace * w;
};

struct prf_params{
	int xi, i, n, ascSize, ascFreq, derivedFlag;
};

struct my_snpProb_params{
  int i, n, ascSize;
  double denom;
	gsl_integration_workspace * w;
};

struct my_f_params{
  int i,n, ascSize;
  double beta;
};

struct prfLik_params{
  int snpNumber, maxSampleSize, ascSize, runMode;
  gsl_vector *weights, *sampleSizeVector, *betas;
  struct snp *data;
  gsl_matrix *sfsCounts;
};

struct my_lik_params{
  int snpNumber, maxSampleSize, ascSize, runMode;
  gsl_vector *weights, *sampleSizeVector, *betas;
  gsl_integration_workspace * w;
  struct snp *data;
  gsl_matrix *sfsBools;
  gsl_min_fminimizer *min;

};
struct my_likCI_params{
  int snpNumber, maxSampleSize;
  gsl_vector *weights, *sampleSizeVector;
  struct snp *data;
  gsl_matrix *sfsBools;
  double lMax, logUnits, beta_hat;
};

struct my_jointLik_params{
  int snpNumber, maxSampleSize, ascSize, nstates;
  gsl_vector *weights, *sampleSizeVector, *betas, *startingPoint;
  struct snp *data;
  gsl_matrix *sfsBools, *posts;

};

double my_f(double q, void * params);
double my_F(double beta, void * params);
double my_F2(double beta, void * p);
double my_FAdapt(double beta, void * p);
double my_lik(double beta, void * params);
double bustProb(double beta, double theta, void *p);
double sfsLikBetaVector(gsl_vector *beta, void * p);
double prfLik(double beta, void * p);
double weightedLik(double beta, void * params);
double weightedLikLook(double beta, void * params);
double weightedLikLookCI(double beta, void * params);
double likWrap(double beta);
double snpProb(double beta, void * params);
double snpProbDenomLookup(double beta, void * params);
double snpProbDenom(double beta, int sampleSize);
gsl_vector *makeSnpProbDenomVector(double beta, int maxSampleSize, gsl_vector *sampleSizeVec);
double ml_est(double * lik_beta_hat);
double ml_prf(double * lik_beta_hat, void * p);
double weighted_ml_est(double * lik_beta_hat, void * p);
double weighted_ml_est_lookup(double * lik_beta_hat, void * p);
double weighted_ml_est_lookup_errReport(double * lik_beta_hat, void * p);
double weighted_ml_CILower_lookup(double * lik_beta_hat, void * p);
double weighted_ml_CIUpper_lookup(double * lik_beta_hat, void * p);
double wlikWrap(double beta, void * p);
double  wlikWrapLook(double beta, void * p);
double  wlikCIWrapLook(double beta, void * p);
double  prfLikWrap(double beta, void * p);
gsl_vector *sampleSizeVector(struct snp *data, int snpNumber, int maxSampleSize);
void sampleSizeVectorFill(struct snp *data, int snpNumber, int maxSampleSize, gsl_vector *bools);
gsl_matrix *summarizeSFS(int maxSampleSize, struct snp *data, int snpNumber);
gsl_matrix *summarizeSFSBool(int max, struct snp *data, int snpNumber);
void summarizeSFSBoolFill(gsl_matrix *bool, struct snp *data, int snpNumber);
void summarizeSFSFill(gsl_matrix *counts, struct snp *data, int snpNumber);
int summarizeSFSFillFromTo(gsl_matrix *counts, struct snp *data, int snpNumber, int from, int to, int startPos);
gsl_matrix *snpProbMatrix(double beta,  int maxSampleSize, gsl_vector *sampleSizeVector, gsl_matrix *sfsBools);
gsl_matrix *snpProbMatrixNotLog(double beta,  int maxSampleSize, gsl_vector *sampleSizeVector, gsl_matrix *sfsBools);
gsl_matrix *snpProbMatrixNotLogFull(double beta,  int maxSampleSize,gsl_vector *sampleSizeVector);
gsl_vector *snpProbVectorNotLog(double beta,  int sampleSize);
int simulateSiteFreq(int sampleSize, double alpha, void *r);
double simpleAscertain(double sampleFreq, void *p);
double probAscertainmentGivenModel(double beta, void *p);
gsl_matrix *snpAscMatrix(int maxSampleSize, gsl_vector *sampleSizeVector, gsl_matrix *sfsBools, int ascSize);
gsl_matrix *estimateAscSFS(int maxSampleSize, gsl_vector *sampleSizeVector,  gsl_matrix *sfsBools, int ascSize, gsl_matrix *sfsSummary);
double weightedLikLookAsc(double beta, void * p);
double  wlikWrapLookAsc(double beta, void * p);
double weighted_ml_est_lookup_asc(double * lik_beta_hat, void * p);
double probAscertainmentGivenModelLookup(double beta, int sampSize, gsl_matrix *snpProbs, gsl_matrix *ascProbs);
double probAscertainmentGivenModelHemiLookup(double beta, int sampSize, gsl_vector *snpProbs, gsl_matrix *ascProbs);
double outgroupAscertain(double sampleFreq, void * p);
gsl_matrix *snpOutgroupAscMatrix(int maxSampleSize, gsl_vector *sampleSizeVector, int ascSize, int derived, int runMode);
double weightedLikLookOutgroupAsc(double beta, void * p);
double  wlikWrapLookOutgroupAsc(double beta, void * p);
double weighted_ml_est_lookup_outgroup_asc(double * lik_beta_hat, void * p);
double probOutgroupAscertainmentGivenModel(double beta, void *p);
double sfsLikBetaVectorOutgroupAsc(gsl_vector *betas, void * p);
double sfsLikBetaVectorOutgroupAscPre(gsl_vector *betas, gsl_matrix *ascProbsAncestral,gsl_matrix *ascProbsDerived, gsl_vector *snpProbs, void * p);
//pre-alloced versions
gsl_matrix *snpOutgroupAscMatrixPre(int maxSampleSize, gsl_vector *sampleSizeVector, int ascSize, int derived, \
	int runMode, gsl_matrix *probs);
gsl_vector *snpProbVectorNotLogPre(double beta,  int sampleSize, gsl_vector *probs);
//MPI stuff
double sfsLikBetaVectorMPI(gsl_vector *betas, void * p, int myrank, int mprocs, MPI_Comm comm);
double sfsLikBetaVectorOutgroupAscPreMPI(gsl_vector *betas, gsl_matrix *ascProbsAncestral,gsl_matrix *ascProbsDerived, \
	gsl_vector *snpProbs, void * p,int myrank, int mprocs, MPI_Comm comm);

double weightedLikLookFolded(double beta, void * p);
double weighted_ml_est_lookup_folded(double * lik_beta_hat, void * p);
double weighted_ml_CILower_lookupFolded(double * lik_beta_hat, void * p);
double weighted_ml_CIUpper_lookupFolded(double * lik_beta_hat, void * p);
double  wlikWrapLookFolded(double beta, void * p);

int segSites(struct snp *data, int sites);
double theta_pi(struct snp *data, int sites);
double a1f(int n);
double theta_w(int segSites, int n);
double a2f(int);
double b1f(int);
double b2f(int);
double c1f(double, double);
double c2f(int, double, double, double);
double e1f(double, double);
double e2f(double, double, double);
double tajdSNP(int, int, double) ;

#endif

max_stars_repo_name: andrewkern/pgLib

max_issues_repo_name: andrewkern/pgLib

max_forks_repo_name: andrewkern/pgLib


text:
#pragma once

#include <thrust/iterator/counting_iterator.h>
#include <thrust/iterator/transform_iterator.h>

#include <cuda/runtime_api.hpp>

#include <gsl-lite/gsl-lite.hpp>

#include <cub/cub.cuh>

#include <thrustshift/math.h>
#include <thrustshift/not-a-vector.h>

namespace cub {

template <typename _UnsignedBits, typename T>
struct BaseTraits<SIGNED_INTEGER,
                  false,
                  false,
                  _UnsignedBits,
                  thrustshift::AbsView<T>> {
	typedef _UnsignedBits UnsignedBits;

	static const Category CATEGORY = SIGNED_INTEGER;
	static const UnsignedBits HIGH_BIT = UnsignedBits(1)
	                                     << ((sizeof(UnsignedBits) * 8) - 1);
	static const UnsignedBits LOWEST_KEY = HIGH_BIT;
	static const UnsignedBits MAX_KEY = UnsignedBits(-1) ^ HIGH_BIT;

	enum {
		PRIMITIVE = true,
		NULL_TYPE = false,
	};

	static __device__ __forceinline__ UnsignedBits TwiddleIn(UnsignedBits key) {
		T val = key;
		if (val >= 0) {
			return val;
		}
		UnsignedBits new_val = std::abs(val);
		return new_val | HIGH_BIT;
	};

	static __device__ __forceinline__ UnsignedBits
	TwiddleOut(UnsignedBits key) {
		if (key & HIGH_BIT) {
			auto del_high_bit = key ^ HIGH_BIT;
			T x = -T(del_high_bit);
			return x;
		}
		return key;
	};

	static __host__ __device__ __forceinline__ T Max() {
		return std::numeric_limits<T>::max();
	}

	static __host__ __device__ __forceinline__ T Lowest() {
		return std::numeric_limits<T>::lowest();
	}
};

template <typename _UnsignedBits, typename T>
struct BaseTraits<FLOATING_POINT,
                  false,
                  false,
                  _UnsignedBits,
                  thrustshift::AbsView<T>> {
	typedef _UnsignedBits UnsignedBits;

	static const Category CATEGORY = FLOATING_POINT;
	static const UnsignedBits HIGH_BIT = UnsignedBits(1)
	                                     << ((sizeof(UnsignedBits) * 8) - 1);
	static const UnsignedBits LOWEST_KEY = UnsignedBits(-1);
	static const UnsignedBits MAX_KEY = UnsignedBits(-1) ^ HIGH_BIT;

	enum {
		PRIMITIVE = true,
		NULL_TYPE = false,
	};

	static __device__ __forceinline__ UnsignedBits TwiddleIn(UnsignedBits key) {
		return key;
	};

	static __device__ __forceinline__ UnsignedBits
	TwiddleOut(UnsignedBits key) {
		return key;
	};

	static __host__ __device__ __forceinline__ T Max() {
		return FpLimits<T>::Max();
	}

	static __host__ __device__ __forceinline__ T Lowest() {
		return FpLimits<T>::Lowest();
	}
};

template <>
struct NumericTraits<thrustshift::AbsView<int>>
    : BaseTraits<SIGNED_INTEGER,
                 false,
                 false,
                 unsigned int,
                 thrustshift::AbsView<int>> {};
template <>
struct NumericTraits<thrustshift::AbsView<long int>>
    : BaseTraits<SIGNED_INTEGER,
                 false,
                 false,
                 long unsigned int,
                 thrustshift::AbsView<long int>> {};
template <>
struct NumericTraits<thrustshift::AbsView<long long int>>
    : BaseTraits<SIGNED_INTEGER,
                 false,
                 false,
                 long long unsigned int,
                 thrustshift::AbsView<long long int>> {};
template <>
struct NumericTraits<thrustshift::AbsView<float>>
    : BaseTraits<FLOATING_POINT,
                 false,
                 false,
                 unsigned int,
                 thrustshift::AbsView<float>> {};
template <>
struct NumericTraits<thrustshift::AbsView<double>>
    : BaseTraits<FLOATING_POINT,
                 false,
                 false,
                 long long unsigned int,
                 thrustshift::AbsView<double>> {};

} // namespace cub

namespace thrustshift {

namespace async {

/*! \brief Batched sort of keys with values.
 *
 *  Uses CUB's radix sort.
 */
template <class KeyInRange,
          class KeyOutRange,
          class ValueInRange,
          class ValueOutRange,
          class MemoryResource>
void sort_batched_descending(cuda::stream_t& stream,
                             KeyInRange&& keys_in,
                             KeyOutRange&& keys_out,
                             ValueInRange&& values_in,
                             ValueOutRange&& values_out,
                             std::size_t batch_len,
                             MemoryResource& delayed_memory_resource) {

	const std::size_t N = keys_in.size();

	gsl_Expects(N % batch_len == 0);
	gsl_Expects(batch_len > 0);
	gsl_Expects(keys_out.size() == N);
	gsl_Expects(values_in.size() == N);
	gsl_Expects(values_out.size() == N);

	auto cit = thrust::make_counting_iterator(0);
	auto tit = thrust::make_transform_iterator(
	    cit, [batch_len] __device__(int i) { return i * batch_len; });

	const std::size_t num_batches = N / batch_len;

	size_t tmp_bytes_size = 0;
	void* tmp_ptr = nullptr;

	using KeyT = typename std::remove_reference<KeyOutRange>::type::value_type;

	auto exec = [&] {
		cuda::throw_if_error(cub::DeviceSegmentedRadixSort::SortPairsDescending(
		    tmp_ptr,
		    tmp_bytes_size,
		    keys_in.data(),
		    keys_out.data(),
		    values_in.data(),
		    values_out.data(),
		    gsl_lite::narrow<int>(N),
		    gsl_lite::narrow<int>(num_batches),
		    tit,
		    tit + 1,
		    0, // default value by CUB
		    sizeof(KeyT) * 8, // default value by CUB
		    stream.handle()));
	};
	exec();
	auto tmp =
	    make_not_a_vector<uint8_t>(tmp_bytes_size, delayed_memory_resource);
	tmp_ptr = tmp.to_span().data();
	exec();
}

/*! \brief Batched sort of keys with respect to their absolute values.
 *
 *  Example:
 *
 *  ```
 *  batch_len = 5
 *  keys_in = {-8, 7, 10, -6, 5}
 *  // sort_batched_abs
 *  keys_out = {5, -6, 7, -8, 10}
 *  ```
 *
 *  ```
 *  batch_len = 3
 *  keys_in = {-8, 7, 10,   -6, 5, 4}
 *  // sort_batched_abs
 *  keys_out = {7, -8, 10,   4, 5, -6}
 *  ```
 */
template <class KeyInRange,
          class KeyOutRange,
          class ValueInRange,
          class ValueOutRange,
          class MemoryResource>
void sort_batched_abs(cuda::stream_t& stream,
                      KeyInRange&& keys_in,
                      KeyOutRange&& keys_out,
                      ValueInRange&& values_in,
                      ValueOutRange&& values_out,
                      std::size_t batch_len,
                      MemoryResource& delayed_memory_resource) {

	const std::size_t N = keys_in.size();

	gsl_Expects(N % batch_len == 0);
	gsl_Expects(batch_len > 0);
	gsl_Expects(keys_out.size() == N);
	gsl_Expects(values_in.size() == N);
	gsl_Expects(values_out.size() == N);

	using KeyT = typename std::remove_reference<KeyOutRange>::type::value_type;
	using AbsT = AbsView<KeyT>;

	auto cit = thrust::make_counting_iterator(0);
	auto tit = thrust::make_transform_iterator(
	    cit, [batch_len] __device__(int i) { return i * batch_len; });

	const std::size_t num_batches = N / batch_len;

	size_t tmp_bytes_size = 0;
	void* tmp_ptr = nullptr;

	auto exec = [&] {
		cuda::throw_if_error(cub::DeviceSegmentedRadixSort::SortPairs(
		    tmp_ptr,
		    tmp_bytes_size,
		    reinterpret_cast<const AbsT*>(keys_in.data()),
		    reinterpret_cast<AbsT*>(keys_out.data()),
		    values_in.data(),
		    values_out.data(),
		    gsl_lite::narrow<int>(N),
		    gsl_lite::narrow<int>(num_batches),
		    tit,
		    tit + 1,
		    0, // first key bit for comparison
		    sizeof(KeyT) * 8 - 1, // highest bit is used to save sign
		    stream.handle()));
	};
	exec();
	auto tmp =
	    make_not_a_vector<uint8_t>(tmp_bytes_size, delayed_memory_resource);
	tmp_ptr = tmp.to_span().data();
	exec();
}

} // namespace async

} // namespace thrustshift

ext: h

max_stars_repo_name: pauleonix/thrustshift

max_issues_repo_name: pauleonix/thrustshift

max_forks_repo_name: pauleonix/thrustshift


text:
#ifndef CTETRA_WEIGHTS_H
#define CTETRA_WEIGHTS_H

#include <stdbool.h>
#include <gsl/gsl_vector.h>
#include "ecache.h"
#include "dos.h"

void WeightsAtK(double E_Fermi, int i, int j, int k, EnergyCache *Ecache, double *result);

void sortEsKs(double *Es, int **ks);

void WeightContrib(double E_Fermi, double E1, double E2, double E3, double E4, double num_tetra, double *ws);

void addCurvatureCorrection(double E_Fermi, double E1, double E2, double E3, double E4, double num_tetra, double *ws);

#endif // CTETRA_WEIGHTS_H

max_stars_repo_name: tflovorn/ctetra

max_issues_repo_name: tflovorn/ctetra

max_forks_repo_name: tflovorn/ctetra


text:
#include <algorithm>
#include <numeric>

#ifdef HAS_MKL
#include <mkl.h>
#else
#include <cblas.h>
#endif

#ifdef __SSE4_2__
#include <emmintrin.h>
#endif

typedef double Real;

inline void GetTotalTraction_(Real *TotalTraction,
    const Real *ElectricDisplacementx, int ndim) {
    // if (ndim==3) {
    //     TotalTraction[0] = ElectricDisplacementx[0];
    //     TotalTraction[1] = ElectricDisplacementx[1];
    //     TotalTraction[2] = ElectricDisplacementx[2];
    // }
    // else if (ndim == 2) {
    //     TotalTraction[0] = ElectricDisplacementx[0];
    //     TotalTraction[1] = ElectricDisplacementx[1];
    // }
    std::copy(ElectricDisplacementx,ElectricDisplacementx+ndim,TotalTraction);
}

inline void FillConstitutiveB_(Real *B, const Real* SpatialGradient,
                     int ndim, int rows, int cols) {
    int i = 0;

    if (ndim == 3) {

        for (; i<rows; ++i) {
            // Store in registers
            const Real a0 = SpatialGradient[i*ndim];
            const Real a1 = SpatialGradient[i*ndim+1];
            const Real a2 = SpatialGradient[i*ndim+2];

            // ELECTROSTATIC TERMS
            B[i*cols+0]     = a0;
            B[i*cols+1]     = a1;
            B[i*cols+2]     = a2;
        }
    }

    else if (ndim == 2) {

        for (; i<rows; ++i) {
            // Store in registers
            const Real a0 = SpatialGradient[i*ndim];
            const Real a1 = SpatialGradient[i*ndim+1];

            // ELECTROSTATIC TERMS
            B[i*cols+0]     = a0;
            B[i*cols+1]     = a1;
        }
    }
}
#include <iostream>
inline void _ConstitutiveStiffnessLaplacian_Filler_(Real *stiffness, Real *traction,
    const Real* SpatialGradient,
    const Real* ElectricDisplacementx,
    const Real* H_Voigt,
    const Real* detJ,
    int ngauss,
    int noderpelem,
    int ndim,
    int nvar,
    int H_VoigtSize,
    int requires_geometry_update) {

    int local_size = nvar*noderpelem;

    Real *t;

#ifdef __SSE4_2__
    t = (Real*)_mm_malloc(ndim*sizeof(Real),32);
    Real *B = (Real*)_mm_malloc(H_VoigtSize*local_size*sizeof(Real),32);
    Real *HBT = (Real*)_mm_malloc(H_VoigtSize*local_size*sizeof(Real),32);
    Real *BDB_1 = (Real*)_mm_malloc(local_size*local_size*sizeof(Real),32);
#else
    t = (Real*)malloc(ndim*sizeof(Real));
    Real *B = (Real*)malloc(H_VoigtSize*local_size*sizeof(Real));
    Real *HBT = (Real*)malloc(H_VoigtSize*local_size*sizeof(Real));
    Real *BDB_1 = (Real*)malloc(local_size*local_size*sizeof(Real));
#endif

    std::fill(B,B+H_VoigtSize*local_size,0.);

    for (int igauss = 0; igauss < ngauss; ++igauss) {

        FillConstitutiveB_(B,&SpatialGradient[igauss*ndim*noderpelem],ndim,noderpelem,H_VoigtSize);

        cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasTrans,
            H_VoigtSize, local_size, H_VoigtSize, 1.0, &H_Voigt[igauss*H_VoigtSize*H_VoigtSize], H_VoigtSize, B, H_VoigtSize, 0.0, HBT, local_size);

        cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
            local_size, local_size, H_VoigtSize, 1.0, B, H_VoigtSize, HBT, local_size, 0.0, BDB_1, local_size);

        // Multiply stiffness with detJ
        const Real detJ_igauss = detJ[igauss];
        for (int i=0; i<local_size*local_size; ++i) {
            stiffness[i] += BDB_1[i]*detJ_igauss;
        }

        // Compute tractions
        GetTotalTraction_(t,
            &ElectricDisplacementx[igauss*ndim], ndim);

        // Multiply B with traction - for loop is okay
        for (int i=0; i<local_size; ++i) {
            Real tmp = 0;
            for (int j=0; j<H_VoigtSize; ++j) {
                tmp += B[i*H_VoigtSize+j]*t[j];
            }
            traction[i] += tmp*detJ_igauss;
        }
    }

#ifdef __SSE4_2__
    _mm_free(t);
    _mm_free(B);
    _mm_free(HBT);
    _mm_free(BDB_1);
#else
    free(t);
    free(B);
    free(HBT);
    free(BDB_1);
#endif
}

max_stars_repo_name: jdlaubrie/florence

max_issues_repo_name: jdlaubrie/florence

max_forks_repo_name: jdlaubrie/florence


text:
#pragma once

#include <glm/glm.hpp>
#include <gsl/span>
#include <memory>
#include <string>
#include <string_view>
#include <vector>

namespace rev {
class IndexedModel;

class ObjFile {
public:
    ObjFile(const std::string& path);

    struct WavefrontObject {
        std::string materialName;

        using IndexedTriangle = std::array<glm::uvec3, 3>;
        std::vector<IndexedTriangle> triangles;
    };

    gsl::span<const WavefrontObject> getWavefrontObjects() const;
    const glm::vec3& positionAtIndex(size_t index) const;
    const glm::vec2& textureCoordinateAtIndex(size_t index) const;
    const glm::vec3& normalAtIndex(size_t index) const;

private:
    void processLine(const std::string& line);

    std::vector<WavefrontObject> _wfObjects;
    std::vector<glm::vec3> _positions;
    std::vector<glm::vec2> _textureCoordinates;
    std::vector<glm::vec3> _normals;
};
} // namespace rev

max_stars_repo_name: eyebrowsoffire/rev

max_issues_repo_name: eyebrowsoffire/rev

max_forks_repo_name: eyebrowsoffire/rev


text:
/*
 * BRAINS
 * (B)LR (R)everberation-mapping (A)nalysis (I)n AGNs with (N)ested (S)ampling
 * Yan-Rong Li, liyanrong@ihep.ac.cn
 * Thu, Aug 4, 2016
 */

/*! \file dnest_line1d.c
 *  \brief run dnest sampling for 1d line analysis.
 */

#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <math.h>
#include <float.h>
#include <string.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_interp.h>
#include <mpi.h>

#include "brains.h"

DNestFptrSet *fptrset_line1d;

/*!
 * This function setup functions for dnest and run dnest.
 */
int dnest_line1d(int argc, char **argv)
{
  int i;

  set_blr_model1d();

  num_params_blr = num_params_blr_model + 1; /* include line sys err */
  num_params_blr_tot = num_params_blr;
  num_params = parset.n_con_recon + num_params_var + num_params_blr_tot;
  idx_resp = num_params_blr_tot + num_params_drw + num_params_trend;
  idx_difftrend = idx_resp + num_params_resp;

  par_fix = (int *) malloc(num_params * sizeof(int));
  par_fix_val = (double *) malloc(num_params * sizeof(double));

  par_range_model = malloc( num_params * sizeof(double *));
  par_prior_gaussian = malloc(num_params * sizeof(double *));
  for(i=0; i<num_params; i++)
  {
    par_range_model[i] = malloc(2*sizeof(double));
    par_prior_gaussian[i] = malloc(2*sizeof(double));
  }
  par_prior_model = malloc( num_params * sizeof(int));

  fptrset_line1d = dnest_malloc_fptrset();

  /* setup functions used for dnest */
  fptrset_line1d->from_prior = from_prior_line1d;
  fptrset_line1d->print_particle = print_particle_line1d;
  fptrset_line1d->restart_action = restart_action_1d;
  fptrset_line1d->accept_action = accept_action_1d;
  fptrset_line1d->kill_action = kill_action_1d;
  fptrset_line1d->read_particle = read_particle_line1d;
  fptrset_line1d->perturb = perturb_line1d;

  if(parset.flag_exam_prior != 1)
  {
    fptrset_line1d->log_likelihoods_cal_initial = log_likelihoods_cal_initial_line1d;
    fptrset_line1d->log_likelihoods_cal_restart = log_likelihoods_cal_restart_line1d;
    fptrset_line1d->log_likelihoods_cal = log_likelihoods_cal_line1d;
  }
  else
  {
    fptrset_line1d->log_likelihoods_cal_initial = log_likelihoods_cal_line1d_exam;
    fptrset_line1d->log_likelihoods_cal_restart = log_likelihoods_cal_line1d_exam;
    fptrset_line1d->log_likelihoods_cal = log_likelihoods_cal_line1d_exam;
  }

  set_par_range_model1d();
  set_par_fix_blrmodel();

  for(i=num_params_blr_model; i<num_params; i++)
  {
    par_fix[i] = 0;
    par_fix_val[i] = -DBL_MAX;
  }

  /* fix continuum variation parameter sigma and tau if flag_fixvar is true */
  if(parset.flag_fixvar == 1)
  {
    par_fix[num_params_blr + 1] = 1;
    par_fix_val[num_params_blr + 1] = var_param[1];
    par_fix[num_params_blr + 2] = 1;
    par_fix_val[num_params_blr + 2] = var_param[2];
  }

  /* fix systematic error of line */
  if(parset.flag_line_sys_err != 1)
  {
    par_fix[num_params_blr-1] = 1;
    par_fix_val[num_params_blr-1] = log(1.0);
  }

  /* fix systematic error of continuum */
  if(parset.flag_con_sys_err != 1)
  {
    par_fix[num_params_blr] = 1;
    par_fix_val[num_params_blr] = log(1.0);
  }

  /* fix non-linear response */
  if(parset.flag_nonlinear !=1)
  {
    par_fix[idx_resp + 1] = 1;
    par_fix_val[idx_resp + 1] = 0.0;
  }

  print_par_names_model1d();

  force_update = parset.flag_force_update;
  if(parset.flag_para_name != 1)
    logz_line = dnest(argc, argv, fptrset_line1d, num_params, dnest_options_file);

  dnest_free_fptrset(fptrset_line1d);
  return 0;
}

/*!
 * this function setups parameter ranges.
 *
 * The order of parameters is:                              \n
 *   I.   blr model.............()                          \n
 *   II.  systematic error......()                          \n
 *   III. variability...........()                          \n
 *   IV.  long-term trend.......()                          \n
 *   V.   response A and Ag.....()                          \n
 *   VI.  different trend.......(if flag_difftend is ture)  \n
 *   VII. continuum light curve.()                          \n
 */
void set_par_range_model1d()
{
  int i, j, i1, i2;

  /* BLR parameters first */
  for(i=0; i<num_params_blr_model; i++)
  {
    par_range_model[i][0] = blr_range_model[i][0];
    par_range_model[i][1] = blr_range_model[i][1];

    par_prior_model[i] = UNIFORM;
    par_prior_gaussian[i][0] = 0.0;
    par_prior_gaussian[i][1] = 0.0;
  }

  /* note that the last BLR parameters is the systematic error (1d) */
  i = num_params_blr -1;
  par_range_model[i][0] = sys_err_line_range[0];
  par_range_model[i][1] = sys_err_line_range[1];

  par_prior_model[i] = UNIFORM;
  par_prior_gaussian[i][0] = 0.0;
  par_prior_gaussian[i][1] = 0.0;

  /* variability parameters */
  /* first systematic error */
  i = num_params_blr;
  par_range_model[i][0] = var_range_model[0][0];
  par_range_model[i][1] = var_range_model[0][1];

  par_prior_model[i] = UNIFORM;
  par_prior_gaussian[i][0] = 0.0;
  par_prior_gaussian[i][1] = 0.0;

  for(i=num_params_blr+1; i<num_params_drw + num_params_blr; i++)
  {
    if(var_param_std[i-num_params_blr] > 0.0)
    {
      par_range_model[i][0] = var_param[i-num_params_blr] - 5.0 * var_param_std[i-num_params_blr];
      par_range_model[i][1] = var_param[i-num_params_blr] + 5.0 * var_param_std[i-num_params_blr];

      /* make sure that the range lies within the initial range */
      par_range_model[i][0] = fmax(par_range_model[i][0], var_range_model[i-num_params_blr][0]);
      par_range_model[i][1] = fmin(par_range_model[i][1], var_range_model[i-num_params_blr][1]);

      par_prior_model[i] = GAUSSIAN;
      par_prior_gaussian[i][0] = var_param[i-num_params_blr];
      par_prior_gaussian[i][1] = var_param_std[i-num_params_blr];
    }
    else
    {
      par_range_model[i][0] = var_range_model[i-num_params_blr][0];
      par_range_model[i][1] = var_range_model[i-num_params_blr][1];

      par_prior_model[i] = UNIFORM;
      par_prior_gaussian[i][0] = 0.0;
      par_prior_gaussian[i][1] = 0.0;
    }
  }
  /* long-term trend */
  for(i=num_params_drw + num_params_blr; i< num_params_drw + num_params_trend + num_params_blr; i++)
  {
    par_range_model[i][0] = var_range_model[3][0];
    par_range_model[i][1] = var_range_model[3][1];

    par_prior_model[i] = GAUSSIAN;
    par_prior_gaussian[i][0] = 0.0;
    par_prior_gaussian[i][1] = 1.0;
  }
  /* response A and Ag */
  j = 0;
  i1 = idx_resp;
  i2 = idx_resp + num_params_resp;
  for(i=i1; i<i2; i++)
  {
    par_range_model[i][0] = resp_range[j][0];
    par_range_model[i][1] = resp_range[j][1];

    par_prior_model[i] = UNIFORM;
    par_prior_gaussian[i][0] = 0.0;
    par_prior_gaussian[i][1] = 0.0;

    j++;
  }
  /* different trend */
  j = 0;
  i1 = idx_difftrend;
  i2 = idx_difftrend + num_params_difftrend;
  for(i=i1; i<i2; i++)
  {
    par_range_model[i][0] = var_range_model[4 + j][0];
    par_range_model[i][1] = var_range_model[4 + j][1];

    par_prior_model[i] = UNIFORM;
    par_prior_gaussian[i][0] = 0.0;
    par_prior_gaussian[i][1] = 0.0;

    j++;
  }

  /* continuum light curve parameters */
  for(i=num_params_blr+num_params_var; i<num_params; i++)
  {
    par_range_model[i][0] = var_range_model[4+num_params_difftrend][0];
    par_range_model[i][1] = var_range_model[4+num_params_difftrend][1];

    par_prior_model[i] = GAUSSIAN;
    par_prior_gaussian[i][0] = 0.0;
    par_prior_gaussian[i][1] = 1.0;
  }
  return;
}

/*!
 *  print names and prior ranges for parameters
 *
 */
void print_par_names_model1d()
{
  if(thistask!= roottask)
    return;

  int i, j;
  FILE *fp;
  char fname[BRAINS_MAX_STR_LENGTH], str_fmt[BRAINS_MAX_STR_LENGTH];

  sprintf(fname, "%s/%s", parset.file_dir, "data/para_names_model1d.txt");
  fp = fopen(fname, "w");
  if(fp == NULL)
  {
    fprintf(stderr, "# Error: Cannot open file %s.\n", fname);
    exit(0);
  }

  strcpy(str_fmt, "%4d %-15s %10.6f %10.6f %4d %4d %15.6e\n");

  printf("# Print parameter name in %s\n", fname);

  fprintf(fp, "#*************************************************\n");
  fprint_version(fp);
  fprintf(fp, "#*************************************************\n");

  fprintf(fp, "%4s %-15s %10s %10s %4s %4s %15s\n", "#", "Par", "Min", "Max", "Prior", "Fix", "Val");

  i=-1;
  for(j=0; j<num_params_blr_model; j++)
  {
    i++;
    fprintf(fp, str_fmt, i, "BLR model", par_range_model[i][0], par_range_model[i][1], par_prior_model[i],
                            par_fix[i], par_fix_val[i]);
  }

  i++;
  fprintf(fp, str_fmt, i, "sys_err_line", par_range_model[i][0], par_range_model[i][1], par_prior_model[i],
                            par_fix[i], par_fix_val[i]);

  i++;
  fprintf(fp, str_fmt, i, "sys_err_con", par_range_model[i][0], par_range_model[i][1], par_prior_model[i],
                            par_fix[i], par_fix_val[i]);
  i++;
  fprintf(fp, str_fmt, i, "sigmad", par_range_model[i][0], par_range_model[i][1], par_prior_model[i],
                            par_fix[i], par_fix_val[i]);
  i++;
  fprintf(fp, str_fmt, i, "taud", par_range_model[i][0], par_range_model[i][1], par_prior_model[i],
                            par_fix[i], par_fix_val[i]);

  for(j=0; j<num_params_trend; j++)
  {
    i++;
    fprintf(fp, str_fmt, i, "trend", par_range_model[i][0], par_range_model[i][1], par_prior_model[i],
                            par_fix[i], par_fix_val[i]);
  }

  i++;
  fprintf(fp, str_fmt, i, "A", par_range_model[i][0], par_range_model[i][1], par_prior_model[i],
                            par_fix[i], par_fix_val[i]);

  i++;
  fprintf(fp, str_fmt, i, "Ag", par_range_model[i][0], par_range_model[i][1], par_prior_model[i],
                            par_fix[i], par_fix_val[i]);

  for(j=0; j<num_params_difftrend; j++)
  {
    i++;
    fprintf(fp, str_fmt, i, "diff trend", par_range_model[i][0], par_range_model[i][1], par_prior_model[i],
                            par_fix[i], par_fix_val[i]);
  }

  for(j=0; j<parset.n_con_recon; j++)
  {
    i++;
    fprintf(fp, str_fmt, i, "time series", par_range_model[i][0], par_range_model[i][1], par_prior_model[i],
                            par_fix[i], par_fix_val[i]);
  }

  fclose(fp);
}

/*!
 * This function generates a sample from the prior.
 */
void from_prior_line1d(void *model)
{
  int i;
  double *pm = (double *)model;

  for(i=0; i<num_params; i++)
  {
    if(par_prior_model[i] == GAUSSIAN)
    {
      pm[i] = dnest_randn()*par_prior_gaussian[i][1] + par_prior_gaussian[i][0];
      dnest_wrap(&pm[i], par_range_model[i][0], par_range_model[i][1]);
    }
    else
    {
      pm[i] = par_range_model[i][0] + dnest_rand() * (par_range_model[i][1] - par_range_model[i][0]);
    }
  }
  /* cope with fixed parameters. */
  for(i=0; i<num_params_blr + num_params_var; i++)
  {
    if(par_fix[i] == 1)
      pm[i] = par_fix_val[i];
  }

  /* all parameters need to update at the initial step */
  which_parameter_update = -1;
  return;
}

/*!
 * This function calculate log likelihood probability at the initial step.
 */
double log_likelihoods_cal_initial_line1d(const void *model)
{
  double logL;
  logL = prob_initial_line1d(model);
  return logL;
}

/*!
 * This function calculate log likelihood probability at the initial step.
 */
double log_likelihoods_cal_restart_line1d(const void *model)
{
  double logL;
  logL = prob_restart_line1d(model);
  return logL;
}

/*!
 * This function calculate log likelihood probability.
 */
double log_likelihoods_cal_line1d(const void *model)
{
  double logL;
  logL = prob_line1d(model);
  return logL;
}

/*!
 * This function print out the particle into the file.
 */
void print_particle_line1d(FILE *fp, const void *model)
{
  int i;
  double *pm = (double *)model;

  for(i=0; i<num_params; i++)
  {
    fprintf(fp, "%e ", pm[i] );
  }
  fprintf(fp, "\n");
  return;
}

/*!
 * This function read the particle from the file.
 */
void read_particle_line1d(FILE *fp, void *model)
{
  int j;
  double *psample = (double *)model;

  for(j=0; j < dnest_num_params; j++)
  {
    if(fscanf(fp, "%lf", psample+j) < 1)
    {
      printf("%f\n", *psample);
      fprintf(stderr, "#Error: Cannot read file %s.\n", options.sample_file);
      exit(0);
    }
  }
  return;
}

/*!
 * this function perturbs the parameters.
 */
double perturb_line1d(void *model)
{
  double *pm = (double *)model;
  double logH = 0.0, limit1, limit2, width, rnd;
  int which, which_level;

  /*
   * sample BLR and variability parameters more frequently;
   * fixed parameter needs not to update.
   */
  do
  {
    rnd = dnest_rand();
    if(rnd < fmax(0.2, 1.0*(num_params_blr+num_params_var)/num_params))
      which = dnest_rand_int(num_params_blr + num_params_var);
    else
      which = dnest_rand_int(parset.n_con_recon) + num_params_blr + num_params_var;

  }while(par_fix[which]==1);

  which_parameter_update = which;

  /* level-dependent width */
  which_level_update = dnest_get_which_level_update();
  which_level = which_level_update > (size_levels - 10)?(size_levels -10):which_level_update;

  if( which_level > 0)
  {
    limit1 = limits[(which_level-1) * num_params *2 + which *2];
    limit2 = limits[(which_level-1) * num_params *2 + which *2 + 1];
    width = limit2 - limit1;
  }
  else
  {
    width = ( par_range_model[which][1] - par_range_model[which][0] );
  }

  width /= (2.35);

  if(par_prior_model[which] == GAUSSIAN)
  {
    logH -= (-0.5*pow((pm[which] - par_prior_gaussian[which][0])/par_prior_gaussian[which][1], 2.0) );
    pm[which] += dnest_randh() * width;
    dnest_wrap(&pm[which], par_range_model[which][0], par_range_model[which][1]);
    logH += (-0.5*pow((pm[which] - par_prior_gaussian[which][0])/par_prior_gaussian[which][1], 2.0) );
  }
  else
  {
    pm[which] += dnest_randh() * width;
    dnest_wrap(&(pm[which]), par_range_model[which][0], par_range_model[which][1]);
  }

  return logH;
}

/*!
 * This function calculate log likelihood probability.
 */
double log_likelihoods_cal_line1d_exam(const void *model)
{
  return 0.0;
}

void accept_action_1d()
{
  int param;
  double *ptemp;

  // the parameter previously updated
  param = which_parameter_update;
  /* continuum parameter is updated */
  if( param >= num_params_blr )
  {
    /*
     *note that (response) Fline is also changed as long as Fcon is changed.
     *num_params_blr-the parameter is the systematic error of continuum.
     *the change of this parameter also changes continuum reconstruction.
     */

    ptemp = Fcon_rm_particles[which_particle_update];
    Fcon_rm_particles[which_particle_update] = Fcon_rm_particles_perturb[which_particle_update];
    Fcon_rm_particles_perturb[which_particle_update] = ptemp;

    if(force_update != 1)
    {
      ptemp = Fline_at_data_particles[which_particle_update];
      Fline_at_data_particles[which_particle_update] = Fline_at_data_particles_perturb[which_particle_update];
      Fline_at_data_particles_perturb[which_particle_update] = ptemp;
    }
  }
  else if( param < num_params_blr-1 && force_update != 1)
  {
    /* BLR parameter is updated
     * Note a) that the (num_par_blr-1)-th parameter is systematic error of line.
     * when this parameter is updated, Trans1D and Fline are unchanged.
     *      b) Fline is always changed, except param = num_params_blr-1.
     */

    {

      ptemp = TransTau_particles[which_particle_update];
      TransTau_particles[which_particle_update] = TransTau_particles_perturb[which_particle_update];
      TransTau_particles_perturb[which_particle_update] = ptemp;

      ptemp = Trans1D_particles[which_particle_update];
      Trans1D_particles[which_particle_update] = Trans1D_particles_perturb[which_particle_update];
      Trans1D_particles_perturb[which_particle_update] = ptemp;

      ptemp = Fline_at_data_particles[which_particle_update];
      Fline_at_data_particles[which_particle_update] = Fline_at_data_particles_perturb[which_particle_update];
      Fline_at_data_particles_perturb[which_particle_update] = ptemp;
    }
  }

  return;
}

/*
 * action when particle i is killed in cdnest sampling.
 * particle i_copy's properties is copyed to particle i.
 */
void kill_action_1d(int i, int i_copy)
{
  memcpy(Fcon_rm_particles[i], Fcon_rm_particles[i_copy], parset.n_con_recon * sizeof(double));
  memcpy(Fline_at_data_particles[i], Fline_at_data_particles[i_copy], n_line_data * sizeof(double));
  memcpy(TransTau_particles[i], TransTau_particles[i_copy], parset.n_tau*sizeof(double));
  memcpy(Trans1D_particles[i], Trans1D_particles[i_copy], parset.n_tau * sizeof(double));
  return;
}

ext: c

max_stars_repo_name: yzxamos/BRAINS

max_issues_repo_name: yzxamos/BRAINS

max_forks_repo_name: yzxamos/BRAINS


text:
/*
  @copyright Russell Standish 2000-2013
  @author Russell Standish
  This file is part of EcoLab

  Open source licensed under the MIT license. See LICENSE for details.
*/

/**\file
\brief random generator based on GNUSL implementation
*/
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

namespace ecolab
{

  class urand: public random_gen
  {
    void operator=(urand&);
    CLASSDESC_ACCESS(urand);
  public:
    gsl_rng *gen;
    urand(const gsl_rng_type *descr=gsl_rng_mt19937) {gen=gsl_rng_alloc(descr);}
    ~urand() {gsl_rng_free(gen);}
    void Seed(int s) {gsl_rng_set(gen,s);}
    void seed(TCL_args args) {Seed(args);}
    double rand();
    /* select a different uniform random generator according the GSL's rng
       string interface */
    void set_gen(TCL_args args);
    void Set_gen(const string& descr)
    {
      static const gsl_rng_type ** rngTypes=gsl_rng_types_setup();
      const gsl_rng_type **g=rngTypes;
      for (; *g; ++g)
        if (descr==(*g)->name)
          {
            gsl_rng_free(gen);
            gen=gsl_rng_alloc(*g);
            break;
          }
      if (!*g) throw error("Cannot create generator %s",descr.c_str());
    }
  };

  class gaussrand: public random_gen
  {
    void operator=(gaussrand&);
    CLASSDESC_ACCESS(gaussrand);
  public:
    urand uni;
    double rand();
  };
}

 /* don't attempt to send these objects over the wire ! */
#ifdef CLASSDESC
#pragma omit pack  ecolab::urand
#pragma omit unpack  ecolab::urand
#endif
namespace classdesc_access
{
  template <>
  struct access_pack<ecolab::urand>:
    public classdesc::NullDescriptor<classdesc::pack_t> {};
  template <>
  struct access_unpack<ecolab::urand>:
    public classdesc::NullDescriptor<classdesc::unpack_t> {};
  template <>
  struct access_pack<ecolab::urand*>:
    public classdesc::NullDescriptor<classdesc::pack_t> {};
  template <>
  struct access_unpack<ecolab::urand*>:
    public classdesc::NullDescriptor<classdesc::unpack_t> {};
}

max_stars_repo_name: digiperfect/ecolab

max_issues_repo_name: digiperfect/ecolab

max_forks_repo_name: digiperfect/ecolab


text:
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_statistics_double.h>
#include <utilities.h>
#include "tz_error.h"
#include "tz_constant.h"
#include "tz_image_lib.h"
#include "tz_objdetect.h"
#include "tz_imatrix.h"
#include "tz_stack_math.h"
#include "tz_stack_bwdist.h"
#include "tz_pixel_array.h"
#include "tz_stack_bwmorph.h"
#include "tz_dimage_lib.h"
#include "tz_darray.h"
#include "tz_u8array.h"
#include "tz_stack_objlabel.h"
#include "tz_voxel_linked_list.h"
#include "tz_voxel_graphics.h"
#include "tz_dmatrix.h"
#include "tz_stack_sampling.h"
#include "tz_stack_draw.h"

INIT_EXCEPTION_MAIN(e)

int main(int argc, char* argv[])
{
  char neuron_name[100];
  if (argc == 1) {
    strcpy(neuron_name, "fly_neuron");
  } else {
    strcpy(neuron_name, argv[1]);
  }

  char file_path[100];
  sprintf(file_path, "../data/%s/mask.tif", neuron_name);

  Stack *mask = Read_Stack(file_path);

  //sprintf(file_path, "../data/%s/seeds.pa", neuron_name);
  //Pixel_Array *pa = Pixel_Array_Read(file_path);
  //double *pa_array = (double *) pa->array;
  //double max_dist = gsl_stats_max(pa_array, 1, pa->size);

  sprintf(file_path, "../data/%s/seeds.ar", neuron_name);
  int length;
  double *pa_array = darray_read(file_path, &length);
  gsl_sort(pa_array, 1, length);

  //double mean_dist = gsl_stats_median_from_sorted_data(pa_array, 1, length);
  double mean_dist = darray_mean(pa_array, length);
  double std_dist = sqrt(gsl_stats_variance(pa_array, 1, length));

  /* for fly neuron data, a good eroding size is around 6*/
  int erode_size = (int) (mean_dist + std_dist + 0.5);
  //int erode_size = 6;

  //int erode_size = (int)
  //  round(gsl_stats_quantile_from_sorted_data(pa_array, 1, length, 0.85));
  //int erode_size = 6;

  //darray_write("../data/test.bn",pa_array, length);

  printf("%d\n", erode_size);

  int dilate_size = erode_size * 2;

  Struct_Element *se = Make_Ball_Se(erode_size);
  Stack *erode_mask = Stack_Erode_Fast(mask, NULL, se);

  Free_Se(se);

  se = Make_Ball_Se(dilate_size);
  Stack *dilate_mask = Stack_Dilate(erode_mask, NULL, se);

  /*
  Stack *erode_region = Copy_Stack(erode_mask);
  int region_start = 2;
  int nobj = Stack_Label_Objects_N(erode_region, NULL, 1, region_start, 26);
  double *dist = darray_malloc(nobj);
  int i;
  for (i = 0; i < nobj; i++) {
    dist[i] = 0.0;
  }

  sprintf(file_path, "../data/%s/erode_region.tif", neuron_name);
  Write_Stack(file_path, erode_region);

  sprintf(file_path, "../data/%s/seeds.tif", neuron_name);
  Stack *seed_stack = Read_Stack(file_path);
  Voxel_List *list = Stack_To_Voxel_List(seed_stack);
  Kill_Stack(seed_stack);
  Pixel_Array *region_pa = Voxel_List_Sampling(erode_region, list);
  uint8 *region_array = (uint8 *) region_pa->array;
  //  u8array_print2(region_array, region_pa->size, 1);

  dist -= region_start;
  for (i = 0; i < region_pa->size; i++) {
    if (region_array[i] >= region_start) {
      if (dist[region_array[i]] <  pa_array[i]) {
	dist[region_array[i]] = pa_array[i];
      }
    }
  }
  dist += region_start;

  darray_print2(dist, nobj, 1);

  Stack *dilate_mask = Copy_Stack(erode_mask);
  Stack *tmp_dilate = Copy_Stack(erode_mask);
  for (i = 0; i < nobj; i++) {
    se = Make_Ball_Se((int) (dist[i] + 0.5));
    Stack_Dilate_Rm(erode_mask, tmp_dilate, se,
		    erode_region, region_start + i);
    Stack_Or(dilate_mask, tmp_dilate, dilate_mask);
    Free_Se(se);
  }

  Kill_Stack(tmp_dilate);
  */

  Stack_And(mask, dilate_mask, dilate_mask);

  sprintf(file_path, "../data/%s/blobmask.tif", neuron_name);
  Write_Stack(file_path, dilate_mask);
  printf("%s created.\n", file_path);

  sprintf(file_path, "../data/%s/large_seeds.tif", neuron_name);
  Stack *seeds = Read_Stack(file_path);
  Stack *blob = dilate_mask;
  Stack_Not(blob, blob);
  Stack_And(seeds, blob, seeds);

  Kill_Stack(mask);
  Kill_Stack(dilate_mask);
  Kill_Stack(erode_mask);

  sprintf(file_path, "../data/%s.tif", neuron_name);
  Stack *signal = Read_Stack(file_path);
  Stack *canvas = Translate_Stack(signal, COLOR, 0);
  Stack_Draw_Voxel_Mc(canvas, seeds, 4, 255, 0, 0);
  sprintf(file_path, "../data/%s/seedimg.tif", neuron_name);
  Write_Stack(file_path, canvas);
  printf("%s (show where tracing seeds are located) created.\n", file_path);

  return 0;
}


text:
//#####################################################################
// Header blas
//#####################################################################
#ifdef GEODE_BLAS

// On Mac, we can't include Accelerate.h unless __SSE__ is available
#if !defined(__APPLE__) || (defined(__SSE__) && !defined(__COVERITY__))

#ifndef __blas_wrap_iterating__
#ifndef __blas_wrap_h__
#define __blas_wrap_h__

#ifdef GEODE_MKL
#include <mkl_lapack.h>
#include <mkl_cblas.h>
#include <mkl_blas.h>
#include <mkl_trans.h>
#else
#ifdef __APPLE__
#include <Accelerate/Accelerate.h>
#else
extern "C" {
#include <cblas.h>
}
#endif
#endif
#include <geode/array/RawArray.h>
#include <geode/array/Array2d.h>
#include <geode/geometry/Box.h>
namespace geode {

// Problem-dependent parameter inspection
template<class T> GEODE_CORE_EXPORT int ilaenv(int ispec,const char* name,const char* opts,int m,int n);

template<> struct FromPython<CBLAS_TRANSPOSE>{GEODE_CORE_EXPORT static CBLAS_TRANSPOSE convert(PyObject* object);};

#define WRAP(declaration) \
  static declaration GEODE_UNUSED; \
  static declaration

#define __blas_wrap_iterating__

#ifdef GEODE_MKL
#  define BC(name) s##name
#  define MKL_BC(name) mkl_s##name
#  define DECLARE(name,...)
#else
#  define BC(name) s##name##_
#  define DECLARE(name,...) extern "C" { int s##name##_(__VA_ARGS__); }
#endif
#define CBC(name) cblas_s##name
#define SA(name) ("a" #name)
#define T float
#include "blas.h"
#undef SA
#undef BC
#undef CBC
#undef MKL_BC
#undef DECLARE
#undef T

#ifdef GEODE_MKL
#  define BC(name) d##name
#  define MKL_BC(name) mkl_d##name
#  define DECLARE(name,...)
#else
#  define BC(name) d##name##_
#  define DECLARE(name,...) extern "C" { int d##name##_(__VA_ARGS__); }
#endif
#define CBC(name) cblas_d##name
#define SA(name) ("d" #name)
#define T double
#include "blas.h"
#undef SA
#undef BC
#undef CBC
#undef MKL_BC
#undef DECLARE
#undef T

#undef __blas_wrap_iterating__
#undef WRAP

}
#endif
#else // __blas_wrap_iterating__

// Scale a vector: x = ax
WRAP(void scal(int n, T alpha, T* x)) {
  if(alpha) CBC(scal)(n,alpha,x,1);
  else memset(x,0,n*sizeof(T)); // mkl scal doesn't work for alpha = 0
}

// Scale a vector: x = ax
WRAP(void scal(T alpha, RawArray<T> x)) {
  scal(x.size(),alpha,x.data());
}

// Scale a vector: x = ax
WRAP(void scal(T alpha, Subarray<T> x)) {
  if(alpha) CBC(scal)(x.m,alpha,x.data(),x.stride);
  else x.fill(0);
}

DECLARE(geqrf,int*,int*,T*,int*,T*,T*,int*,int*)
DECLARE(gelqf,int*,int*,T*,int*,T*,T*,int*,int*)

// QR factorize a dense matrix without pivoting.  Call geqrf_work to define the required work size
WRAP(void geqrf(CBLAS_UPLO uplo, RawArray<T,2> A, RawArray<T> tau, Array<T>& work)) {
  int info,lda=max(1,A.n),r=min(A.m,A.n),lwork=work.max_size();
  GEODE_ASSERT(tau.size()==r);if(!r) return; // lapack requires tau.size()=1 here, so do nothing instead
  (uplo==CblasUpper?BC(gelqf):BC(geqrf))(const_cast<int*>(&A.n),const_cast<int*>(&A.m),A.data(),&lda,tau.data(),work.data(),&lwork,&info);
  GEODE_ASSERT(!info);
}

// Compute work.size() required for geqrf
WRAP(int geqrf_work(CBLAS_UPLO uplo, RawArray<T,2> A)) {
  if(!A.m || !A.n) return 0;
  T work;int info,lda=max(1,A.n),lwork=-1;
  (uplo==CblasUpper?BC(gelqf):BC(geqrf))(const_cast<int*>(&A.n),const_cast<int*>(&A.m),0,&lda,0,&work,&lwork,&info);
  return max((int)work,A.m,A.n); // Mkl returns the wrong work size, so grow it if necessary
}

DECLARE(geqp3,int*,int*,T*,int*,int*,T*,T*,int*,int*)

// QR factorize a dense matrix with column pivoting
WRAP(void geqp3(RawArray<T,2> A, RawArray<int> p, RawArray<T> tau, Array<T>& work)) {
  GEODE_ASSERT(tau.size()==max(1,min(A.m,A.n)) && p.size()==A.m);
  int info,lda=max(1,A.n),lwork=work.max_size();p.fill(0);
  BC(geqp3)(const_cast<int*>(&A.n),const_cast<int*>(&A.m),A.data(),&lda,p.data(),tau.data(),work.data(),&lwork,&info);
}

// Compute work.size() required for geqp3
WRAP(int geqp3_work(RawArray<T,2> A)) {
  T work;int info,lda=max(1,A.n),lwork=-1;
  BC(geqp3)(const_cast<int*>(&A.n),const_cast<int*>(&A.m),0,&lda,0,0,&work,&lwork,&info);
  return (int)work;
}

DECLARE(ormqr,char*,char*,int*,int*,int*,T*,int*,T*,T*,int*,T*,int*,int*)
DECLARE(ormlq,char*,char*,int*,int*,int*,T*,int*,T*,T*,int*,T*,int*,int*)

// Multiply a matrix C by Q or Q^T, where Q is computed by geqrf
WRAP(void ormqr(CBLAS_UPLO uplo, CBLAS_SIDE s, CBLAS_TRANSPOSE t, Subarray<const T,2> A, RawArray<const T> tau, Subarray<T,2> C, Array<T>& work)) {
  const char *side=s==CblasLeft?"R":"L",*trans=t==CblasTrans?"T":"N"; // side flipped due to row-major vs. column major
  int k=min(A.m,A.n),c=s==CblasLeft?C.m:C.n,lda=max(1,A.stride),ldc=max(1,C.stride),lwork=work.max_size(),info;
  GEODE_ASSERT(k<=c && c<=(uplo==CblasUpper?A.m:A.n) && tau.size()==k);
  if(!k) return; // lapack requires tau.size()=1 here, so do nothing instead
  (uplo==CblasUpper?BC(ormlq):BC(ormqr))(const_cast<char*>(side),const_cast<char*>(trans),const_cast<int*>(&C.n),const_cast<int*>(&C.m),&k,
      const_cast<T*>(A.data()),&lda,const_cast<T*>(tau.data()),C.data(),&ldc,work.data(),&lwork,&info);
  GEODE_ASSERT(!info);
}

// Compute work.size() required by ormqr
WRAP(int ormqr_work(CBLAS_UPLO uplo, CBLAS_SIDE s, CBLAS_TRANSPOSE t, Subarray<const T,2> A, Subarray<T,2> C)) {
  const char *side=s==CblasLeft?"R":"L",*trans=t==CblasTrans?"T":"N"; // side flipped due to row-major vs. column major
  int k=min(A.m,A.n),lda=max(1,A.stride),ldc=max(1,C.stride),lwork=-1,info;T work;
  GEODE_ASSERT((uplo==CblasUpper?A.m:A.n)>=(s==CblasLeft?C.m:C.n));
  if(!k) return 0; // lapack requires tau.size()=1 here, so do nothing instead
  (uplo==CblasUpper?BC(ormlq):BC(ormqr))(const_cast<char*>(side),const_cast<char*>(trans),const_cast<int*>(&C.n),const_cast<int*>(&C.m),&k,0,&lda,0,0,&ldc,&work,&lwork,&info);
  GEODE_ASSERT(!info);
  return (int)work;
}

// Multiply a vector C by Q or Q^T, where Q is computed by geqrf (the single vector case of ormqr above)
WRAP(void ormqr(CBLAS_UPLO uplo, CBLAS_SIDE s, CBLAS_TRANSPOSE t, Subarray<const T,2> A, RawArray<const T> tau, RawArray<T> C, Array<T>& work)) {
  int r=C.size(),c=1;if(s==CblasRight) swap(r,c);
  ormqr(uplo,s,t,A,tau,C.reshape(r,c),work);
}

// Compute work.size() required by the vector case of ormqr
WRAP(int ormqr_work(CBLAS_UPLO uplo,CBLAS_SIDE s,CBLAS_TRANSPOSE t,Subarray<const T,2> A,RawArray<T> C)) {
  int r=C.size(),c=1;if(s==CblasRight) swap(r,c);
  return ormqr_work(uplo,s,t,A,C.reshape(r,c));
}

DECLARE(orgqr,int*,int*,int*,T*,int*,T*,T*,int*,int*)
DECLARE(orglq,int*,int*,int*,T*,int*,T*,T*,int*,int*)

// Expand the QR factorization computed by geqrf into a dense matrix in standard format
WRAP(void orgqr(CBLAS_UPLO uplo, int r, Subarray<T,2> A, RawArray<const T> tau, Array<T>& work)) {
  int info,m=A.m,n=A.n,lda=max(1,A.stride),lwork=work.max_size();
  GEODE_ASSERT(r<=min(m,n) && (uplo==CblasUpper?n<=m:m<=n) && tau.size()==r);
  (uplo==CblasUpper?BC(orglq):BC(orgqr))(&n,&m,&r,A.data(),&lda,const_cast<T*>(tau.data()),work.data(),&lwork,&info);
}

// Compute work.size() required by orgqr
WRAP(int orgqr_work(CBLAS_UPLO uplo, int r, Subarray<T,2> A)) {
  T work;int info,m=A.m,n=A.n,lda=max(1,A.stride),lwork=-1;
  GEODE_ASSERT(r<=min(m,n) && (uplo==CblasUpper?n<=m:m<=n));
  (uplo==CblasUpper?BC(orglq):BC(orgqr))(&n,&m,&r,0,&lda,0,&work,&lwork,&info);
  return (int)work;
}

// y = alpha*A^op x + beta*y
WRAP(void gemv(CBLAS_TRANSPOSE t, T alpha, Subarray<const T,2> A, const T* x, T beta, T* y)) {
  int m=A.m,n=A.n;if(t==CblasTrans) swap(m,n);
  if(!n) scal(m,beta,y); // work around bug in mkl
  else CBC(gemv)(CblasRowMajor,t,A.m,A.n,alpha,A.data(),max(1,A.stride),x,1,beta,y,1);
}

// y = alpha*A^op x + beta*y
WRAP(void gemv(CBLAS_TRANSPOSE t, T alpha, Subarray<const T,2> A, RawArray<const T> x, T beta, RawArray<T> y)) {
  int m=A.m,n=A.n;if(t==CblasTrans) swap(m,n);
  GEODE_ASSERT(x.size()==n && m==y.size());
  gemv(t,alpha,A,x.data(),beta,y.data());
}

// y = alpha*A x + beta*y
WRAP(void gemv(T alpha, Subarray<const T,2> A, const T* x, T beta, T* y)) {
  gemv(CblasNoTrans,alpha,A,x,beta,y);
}

// y = alpha*A x + beta*y
WRAP(void gemv(T alpha, Subarray<const T,2> A, RawArray<const T> x, T beta, RawArray<T> y)) {
  gemv(CblasNoTrans,alpha,A,x,beta,y);
}

// y = alpha*A x + beta*y for symmetric A
WRAP(void symv(CBLAS_UPLO uplo, T alpha, Subarray<const T,2> A, RawArray<const T> x, T beta, RawArray<T> y)) {
  GEODE_ASSERT(A.m==A.n && x.size()==A.m && y.size()==A.m);
  CBC(symv)(CblasRowMajor,uplo,A.m,alpha,A.data(),max(1,A.stride),x.data(),1,beta,y.data(),1);
}

// // y = alpha*A^op+beta*y for banded A with kl subdiagonals and ku superdiagonals.
WRAP(void gbmv(CBLAS_TRANSPOSE t, int m, int n, int kl, int ku, T alpha, RawArray<const T> A, const T* x, T beta, T* y)) {
  int lda=kl+ku+1;GEODE_ASSERT(kl>=0 && ku>=0 && A.size()==n*lda);
  CBC(gbmv)(CblasColMajor,t,m,n,kl,ku,alpha,A.data(),lda,x,1,beta,y,1);
}

DECLARE(gesdd,char*,int*,int*,T*,int*,T*,T*,int*,T*,int*,T*,int*,int*,int*)

// Compute the singular value decomposition of a dense matrix.  A is destroyed.
WRAP(void gesdd(Array<T,2>& A, Array<T>& s, Array<T,2>& U, Array<T,2>& Vt, Array<T>& work, bool all)) {
  char job=all?'A':'S';int r=min(A.m,A.n),liwork=8*r,lwork=(work.max_size()*sizeof(T)-liwork*sizeof(int))/sizeof(T),info;
  U.resize(A.m,all?A.m:r);Vt.resize(all?A.n:r,A.n);s.resize(r);
  BC(gesdd)(&job,&A.n,&A.m,A.data(),&A.n,s.data(),Vt.data(),&Vt.n,U.data(),&U.n,work.data(),&lwork,reinterpret_cast<int*>(work.data()+lwork),&info);
}

// Compute work.size() required by gesdd with singular vectors
WRAP(int gesdd_work(Array<T,2>& A, Array<T,2>& U, Array<T,2>& Vt, bool all)) {
  T work;char job=all?'A':'S';int r=min(A.m,A.n),liwork=max(1,8*r),lwork=-1,iwork,info;
  BC(gesdd)(&job,&A.n,&A.m,0,&A.n,0,0,&A.n,0,&A.m,&work,&lwork,&iwork,&info);
  return int(work+(liwork*sizeof(int)+sizeof(T)-1)/sizeof(T));
}

// Compute the singular values of a dense matrix without the singular vectors.  A is destroyed.
WRAP(void gesdd(Array<T,2>& A, Array<T>& s, Array<T>& work)) {
  char job='N';int r=min(A.m,A.n),one=1,info,liwork=8*r,lwork=(work.max_size()*sizeof(T)-liwork*sizeof(int))/sizeof(T);
  s.resize(r);
  BC(gesdd)(&job,&A.n,&A.m,A.data(),&A.n,s.data(),0,&one,0,&one,work.data(),&lwork,reinterpret_cast<int*>(work.data()+lwork),&info);
}

// Compute work.size() required by gesdd without singular vectors
WRAP(int gesdd_work(Array<T,2>& A)) {
  char job='N';T work;int r=min(A.m,A.n),one=1,info,liwork=8*r,iwork,lwork=-1;
  BC(gesdd)(&job,&A.n,&A.m,0,&A.n,0,0,&one,0,&one,&work,&lwork,&iwork,&info);
  return int(work+(liwork*sizeof(int)+sizeof(T)-1)/sizeof(T));
}

DECLARE(syevd,char*,char*,int*,T*,int*,T*,T*,int*,int*,int*,int*)

// Compute the eigenvalues of a symmetric matrix.  A is destroyed.
WRAP(void syevd(CBLAS_UPLO uplo, Subarray<T,2> A, Array<T>& w, Array<T>& work)) {
  int n=A.n;
  GEODE_ASSERT(A.m==n);
  T lwork_;int liwork,query=-1,lda=A.stride,info;
  char *job=const_cast<char*>("N"),*uplo_=const_cast<char*>(uplo==CblasLower?"U":"L");
  BC(syevd)(job,uplo_,&n,A.data(),&lda,0,&lwork_,&query,&liwork,&query,&info);
  int lwork=(int)lwork_;
  GEODE_ASSERT(lwork*sizeof(T)+liwork*sizeof(int)>=work.max_size()*sizeof(T));
  w.resize(n);
  BC(syevd)(job,uplo_,&n,A.data(),&lda,w.data(),work.data(),&lwork,reinterpret_cast<int*>(work.data()+lwork),&liwork,&info);
}

// Compute work.size() required by syevd
WRAP(int syevd_work(CBLAS_UPLO uplo, Subarray<const T,2> A)) {
  int n=A.n;
  GEODE_ASSERT(A.m==n);
  T work;int iwork,lda=A.stride,info,lwork=-1;
  char *job=const_cast<char*>("N"),*uplo_=const_cast<char*>(uplo==CblasLower?"U":"L");
  BC(syevd)(job,uplo_,&n,const_cast<T*>(A.data()),&lda,0,&work,&lwork,&iwork,&lwork,&info);
  return int(work+(iwork*sizeof(int)+sizeof(T)-1)/sizeof(T));
}

DECLARE(stedc,char*,int*,T*,T*,T*,int*,T*,int*,int*,int*,int*)

// Compute all eigenvalues and eigenvectors of a symmetric tridiagonal matrix.  work is automatically resized.
WRAP(void stedc(RawArray<T> d, RawArray<T> e, Array<T,2>& z, Array<char>& work)) {
  int n=d.size();
  GEODE_ASSERT(n-1<=e.size() && e.size()<=n);
  z.resize(n,n,uninit);
  if(!n) return;
  char *compz=const_cast<char*>("I");
  T lwork_T;int liwork,query=-1,ldz=z.n,info;
  BC(stedc)(compz,&n,d.data(),e.data(),z.data(),&ldz,&lwork_T,&query,&liwork,&query,&info);
  GEODE_ASSERT(!info);
  int lwork=(int)lwork_T;
  static_assert(sizeof(T)/sizeof(int)*sizeof(int)==sizeof(T),"");
  work.preallocate(lwork*sizeof(T)+liwork*sizeof(int));
  BC(stedc)(compz,&n,d.data(),e.data(),z.data(),&ldz,reinterpret_cast<T*>(work.data()),&lwork,
    reinterpret_cast<int*>(work.data()+lwork*sizeof(T)),&liwork,&info);
  GEODE_ASSERT(!info);
}

DECLARE(stegr,char*,char*,int*,T*,T*,T*,T*,int*,int*,T*,int*,T*,T*,int*,int*,T*,int*,int*,int*,int*)

// Compute selected eigenvalues and eigenvectors of a symmetric tridiagonal matrix.  WARNING: Does not work in Mkl 10.2.3.
WRAP(void stegr(RawArray<T> d,RawArray<T> e,Box<T> range,Array<T>& w,Array<T,2>& z,Array<int>& isuppz,Array<char>& work)) {
  int n=d.size();
  GEODE_ASSERT(n==e.size());
  w.resize(n,uninit);z.resize(n,n,uninit);isuppz.resize(2*n,uninit);
  if(!n) return;
  char *jobz=const_cast<char*>("V"),*ran=const_cast<char*>(range.contains(Box<T>::full_box())?"A":"V");
  T lwork_T,unused_T;int liwork,m,query=-1,unused,ldz=z.n,info;
  BC(stegr)(jobz,ran,&n,d.data(),e.data(),&range.min,&range.max,&unused,&unused,&unused_T,&m,w.data(),z.data(),&ldz,
      isuppz.data(),&lwork_T,&query,&liwork,&query,&info);
  GEODE_ASSERT(!info);
  int lwork=(int)lwork_T;
  static_assert(sizeof(T)/sizeof(int)*sizeof(int)==sizeof(T),"");
  work.preallocate(lwork*sizeof(T)+liwork*sizeof(int));
  BC(stegr)(jobz,ran,&n,d.data(),e.data(),&range.min,&range.max,&unused,&unused,&unused_T,&m,w.data(),z.data(),&ldz,
      isuppz.data(),reinterpret_cast<T*>(work.data()),&lwork,reinterpret_cast<int*>(work.data()+lwork*sizeof(T)),&liwork,&info);
  GEODE_ASSERT(!info);
  w.resize(m);z.resize(m,n);isuppz.resize(m);
}

// Matrix-matrix multiply: C = alpha*A^ta B^tb+beta*C
WRAP(void gemm(CBLAS_TRANSPOSE ta, CBLAS_TRANSPOSE tb, T alpha, Subarray<const T,2> A, Subarray<const T,2> B, T beta, Subarray<T,2> C)) {
  int k=ta==CblasTrans?A.m:A.n;
  GEODE_ASSERT(C.m==A.m+A.n-k && k==(tb==CblasTrans?B.n:B.m) && B.m+B.n-k==C.n);
  int lda=max(1,A.stride),ldb=max(1,B.stride),ldc=max(1,C.stride);
  CBC(gemm)(CblasRowMajor,ta,tb,C.m,C.n,k,alpha,A.data(),lda,B.data(),ldb,beta,C.data(),ldc);
}

// Matrix-matrix multiply: C = alpha*AB+beta*C
WRAP(void gemm(T alpha, Subarray<const T,2> A, Subarray<const T,2> B, T beta, Subarray<T,2> C)) {
  gemm(CblasNoTrans,CblasNoTrans,alpha,A,B,beta,C);
}

// Perform a symmetric rank-k update of a matrix: C = alpha*A^t A^!t + beta*C
WRAP(void syrk(CBLAS_UPLO uplo, CBLAS_TRANSPOSE t, T alpha, Subarray<const T,2> A, T beta, Subarray<T,2> C)) {
  int k=t==CblasTrans?A.m:A.n,lda=max(A.stride,1),ldc=max(C.stride,1);
  GEODE_ASSERT(C.m==C.n && C.m==A.m+A.n-k);
  CBC(syrk)(CblasRowMajor,uplo,t,C.m,k,alpha,A.data(),lda,beta,C.data(),ldc);
}

DECLARE(gels,char*,int*,int*,int*,T*,int*,T*,int*,T*,int*,int*)

// Solve a full rank least squares problem.  work is automatically resized.
WRAP(void gels(char t,Subarray<T,2> A,RawArray<T> b,Array<T>& work)) {
  GEODE_ASSERT((t=='N' || t=='T') && max(A.m,A.n)==b.size() && (t=='N'?A.m:A.n)==b.size());
  char ot=t=='N'?'T':'N';int one=1,ldb=max(1,b.size()),lwork=work.size(),info;
  if(!lwork){work.resize(1);lwork=-1;}
  BC(gels)(&ot,const_cast<int*>(&A.n),const_cast<int*>(&A.m),&one,A.data(),const_cast<int*>(&A.stride),b.data(),&ldb,work.data(),&lwork,&info);
  if(lwork<0) work.resize(max(1,(int)work[0]));
}

DECLARE(gelsy,int*,int*,int*,T*,int*,T*,int*,int*,T*,int*,T*,int*,int*)

// Compute a minimum-norm solution to the linear least squares problem A^T x = b.  work is automatically resized.
WRAP(void gelsy(Subarray<T,2> A,RawArray<T> b,Array<T>& work)) {
  GEODE_ASSERT(max(A.m,A.n)==b.size() && A.n==b.size());
  int one=1,ldb=max(1,b.size()),lwork=work.size(),info;
  if(!lwork){work.resize(1);lwork=-1;}
  Array<int> jpvt(A.m);T rcond;int rank;
  BC(gelsy)(const_cast<int*>(&A.n),const_cast<int*>(&A.m),&one,A.data(),const_cast<int*>(&A.stride),b.data(),&ldb,jpvt.data(),&rcond,&rank,work.data(),&lwork,&info);
  if(lwork<0) work.resize(max(1,(int)work[0]));
}

DECLARE(gelsd,int*,int*,int*,T*,int*,T*,int*,T*,T*,int*,T*,int*,int*,int*)

// If iwork is zero size, we will call gelsd in workspace query mode.
// This function will then resize work to be the proper size.
// Nothing else will have been done, so you have to call this function again.
// b.n must be equal to max(A.m,A.n). If A.m > A.n, only the first A.n rows
// contain the solution, the rest contain residuals.
WRAP(void gelsd(Subarray<T,2> A,RawArray<T> b,Array<T>& work,Array<int>& iwork)) { // solves A^T x = b
  int m=A.m,n=A.n;GEODE_ASSERT(max(m,n)==b.size());
  int one=1,ldb=max(1,b.size()),lwork=work.size(),info;
  if(!lwork){work.resize(1);iwork.resize(1);lwork=-1;}
  Array<T> s(max(1,min(m,n)),uninit);T rcond;int rank;
  BC(gelsd)(&n,&m,&one,A.data(),const_cast<int*>(&A.stride),b.data(),&ldb,s.data(),&rcond,&rank,work.data(),&lwork,iwork.data(),&info);
  if(lwork<0){work.resize(max(1,(int)work[0]));iwork.resize(max(1,(int)iwork[0]));}
  if (info)
    std::cout << "error in gelsd: " << info << std::endl;
}

// If iwork is zero size, we will call gelsd in workspace query mode.
// This function will then resize work to be the proper size.
// Nothing else will have been done, so you have to call this function again.
WRAP(void gelsd(Subarray<T,2> A,Array<T,2>b,Array<T>& work,Array<int>& iwork)) { // solves A^T X^T = B^T for all columns of X and B
  int m=A.m,n=A.n;
  GEODE_ASSERT(max(m,n)==b.n);
  int nrhs=b.m, ldb=b.n, lwork=work.size(), info;
  if(!lwork) {
    work.resize(1);
    iwork.resize(1);
    lwork=-1;
  }
  Array<T> s(min(m,n),uninit);
  T rcond;
  int rank;
  BC(gelsd)(&n,&m,&nrhs,A.data(),const_cast<int*>(&A.stride),b.data(),&ldb,s.data(),&rcond,&rank,work.data(),&lwork,iwork.data(),&info);
  if (lwork<0) {
    work.resize((int)work[0]);
    iwork.resize((int)iwork[0]);
  }
  if (info)
    std::cout << "error in gelsd: " << info << std::endl;
}

// y = ax + y
WRAP(void axpy(int n, T alpha, const T* x, T* y)) {
  CBC(axpy)(n,alpha,x,1,y,1);
}

// y = ax + y
WRAP(void axpy(T alpha, RawArray<const T> x, RawArray<T> y)) {
  int n=x.size();GEODE_ASSERT(n==y.size());
  CBC(axpy)(n,alpha,x.data(),1,y.data(),1);
}

// Squared norm: |x|^2 = dot(x,x)
WRAP(T nrm2(RawArray<const T> x)) {
  return CBC(nrm2)(x.size(),x.data(),1);
}

// sum_i x_i
WRAP(T asum(RawArray<const T> x)) {
  return CBC(asum)(x.size(),x.data(),1);
}

// dot(x,y)
WRAP(T dot(RawArray<const T> x, RawArray<const T> y)) {
  int n=x.size();GEODE_ASSERT(n==y.size());
  return CBC(dot)(n,x.data(),1,y.data(),1);
}

// Swap the contents of two subarrays
WRAP(void blas_swap(Subarray<T> x, Subarray<T> y)) {
  int n=x.size();GEODE_ASSERT(n==y.size());CBC(swap)(n,x.data(),x.stride,y.data(),y.stride);
}

// Rank one update: A = A + alpha*x y^T
WRAP(void ger(T alpha, Subarray<const T> x, Subarray<const T> y, Subarray<T,2> A)) {
  GEODE_ASSERT(x.size()==A.m && A.n==y.size());
  int lda=max(1,A.stride),xinc=max(1,x.stride),yinc=max(1,y.stride);
  CBC(ger)(CblasRowMajor,A.m,A.n,alpha,x.data(),xinc,y.data(),yinc,A.data(),lda);
}

DECLARE(laset,char*,int*,int*,T*,T*,T*,int*)

// Initialize the diagonal and offdiagonal of a matrix to the given values
WRAP(void laset(char uplo, T offdiagonal, T diagonal, Subarray<T,2> A)) {
  char uplo_s[2]={uplo,0};int lda=max(1,A.stride);
  BC(laset)(uplo_s,const_cast<int*>(&A.n),const_cast<int*>(&A.m),&offdiagonal,&diagonal,A.data(),&lda);
}

// Compute a scalar matrix-matrix product with one matrix triangular.  TODO: Give details.
WRAP(void trmm(CBLAS_SIDE side, CBLAS_UPLO uplo, CBLAS_TRANSPOSE transA, CBLAS_DIAG diag, T alpha, Subarray<const T,2> A, Subarray<T,2> B)) {
  int k=side==CblasLeft?B.m:B.n;
  GEODE_ASSERT(A.m>=k && A.n>=k && (uplo==CblasLower?A.m:A.n)==k);
  CBC(trmm)(CblasRowMajor,side,uplo,transA,diag,B.m,B.n,alpha,A.data(),max(1,A.stride),B.data(),max(1,B.stride));
}

// Solve a triangular matrix equation.  TODO: Give details.
WRAP(void trsm(CBLAS_SIDE side, CBLAS_UPLO uplo, CBLAS_TRANSPOSE transA, CBLAS_DIAG diag, T alpha, Subarray<const T,2> A, Subarray<T,2> B)) {
  int k=side==CblasLeft?B.m:B.n;
  GEODE_ASSERT(A.m>=k && A.n>=k && (uplo==CblasLower?A.m:A.n)==k);
  CBC(trsm)(CblasRowMajor,side,uplo,transA,diag,B.m,B.n,alpha,A.data(),max(1,A.stride),B.data(),max(1,B.stride));
}

// Solve a triangular matrix equation.  TODO: Give details.
WRAP(void trsm(CBLAS_SIDE side, CBLAS_UPLO uplo, CBLAS_TRANSPOSE transA, CBLAS_DIAG diag, T alpha, Subarray<const T,2> A, RawArray<T> B))
{
  int r=B.size(),c=1;if(side==CblasRight) swap(r,c);
  trsm(side,uplo,transA,diag,alpha,A,B.reshape(r,c));
}

DECLARE(gesv,int*,int*,T*,int*,int*,T*,int*,int*)

// Solve a linear system: b = A^-T b.  WARNING: Note that the matrix is transposed.
WRAP(void gesv(RawArray<T,2> A,RawArray<T> b)) {
  GEODE_ASSERT(A.m==A.n && A.m==b.size());
  int one=1,lda=max(1,A.n),info;
  Array<int> ipiv(A.m,uninit);
  BC(gesv)(const_cast<int*>(&A.m),&one,A.data(),&lda,ipiv.data(),b.data(),const_cast<int*>(&A.n),&info);
}

DECLARE(getrf,int*,int*,T*,int*,int*,int*)

// Compute the pivoted LU factorization of a square matrix.
WRAP(void getrf(RawArray<T,2> A, RawArray<int> ipiv)) {
  int r=min(A.m,A.n),info;GEODE_ASSERT(ipiv.size()==r);
  if(!r) return;
  BC(getrf)(&A.n,&A.m,A.data(),&A.n,ipiv.data(),&info);
}

DECLARE(gecon,char*,int*,T*,int*,T*,T*,T*,int*,int*)

// Estimate the reciprocal of the condition number of a general matrix in the 1-norm or the infinity-norm.
WRAP(T gecon(char norm, RawArray<T,2> A, T anorm)) {
  GEODE_ASSERT(A.m==A.n && (norm=='1' || norm=='0' || norm=='I'));
  T rcond;int info;
  Array<T> work(max(1,4*A.m),uninit);
  Array<int> iwork(max(1,A.m),uninit);
  BC(gecon)(&norm,&A.m,A.data(),&A.n,&anorm,&rcond,work.data(),iwork.data(),&info);
  return rcond;
}

DECLARE(sytrf,char*,int*,T*,int*,int*,T*,int*,int*)

// Compute the Bunch-Kaufman factorization of a symmetric matrix.
WRAP(void sytrf(CBLAS_UPLO uplo, Subarray<T,2> A, RawArray<int> ipiv, Array<T>& work)) {
  GEODE_ASSERT(A.m==A.n && A.m==ipiv.size());
  const char *u=uplo==CblasLower?"L":"U";
  int n=A.m,lda=max(1,A.stride),lwork=work.max_size(),info;
  BC(sytrf)(const_cast<char*>(u),&n,A.data(),&lda,ipiv.data(),work.data(),&lwork,&info);
}

// Compute work.size() required by sytrf
WRAP(int sytrf_work(CBLAS_UPLO uplo,Subarray<T,2> A)) {
  // Mkl appears to have a bug in the sytrf work query, so we do it ourselves
  static const int nb=max(ilaenv<T>(1,SA(sytrf),"L",1024,-1),ilaenv<T>(1,SA(sytrf),"U",1024,-1));
  return max(1,nb*A.m);
}

DECLARE(sytrs,char*,int*,int*,T*,int*,int*,T*,int*,int*)

// Solve a linear system given the Bunch-Kaufman factorization computed by sytrf.  WARNING: Computes B = B A^{-1} since Lapack is column-major.
WRAP(void sytrs(CBLAS_UPLO uplo,Subarray<const T,2> A,RawArray<const int> ipiv,Subarray<T,2> B)) {
  GEODE_ASSERT(A.m==A.n && A.m==ipiv.size() && B.n==A.m);
  const char *u=uplo==CblasLower?"L":"U";
  int n=A.n,nrhs=B.m,lda=max(1,A.stride),ldb=max(1,B.stride),info;
  BC(sytrs)(const_cast<char*>(u),&n,&nrhs,const_cast<T*>(A.data()),&lda,const_cast<int*>(ipiv.data()),B.data(),&ldb,&info);
}

// Solve a linear system given the Bunch-Kaufman factorization computed by sytrf.
WRAP(void sytrs(CBLAS_UPLO uplo,Subarray<const T,2> A,RawArray<const int> ipiv,RawArray<T> B)) {
  sytrs(uplo,A,ipiv,B.reshape(1,B.size()));
}

#ifdef GEODE_MKL

// Perform scaling and in-place transposition/copying of matrices.
WRAP(void imatcopy(T alpha,Subarray<T,2> A)) {
  MKL_BC(imatcopy)('r','t',A.m,A.n,alpha,A.data(),A.n,A.m);
}

// Perform scaling and out-place transposition/copying of matrices.
WRAP(void omatcopy(CBLAS_TRANSPOSE trans,T alpha,Subarray<const T,2> A,Subarray<T,2> B)) {
  int k=trans==CblasTrans?A.n:A.m,lda=max(1,A.stride),ldb=max(1,B.stride);
  GEODE_ASSERT(k==B.m && A.m+A.n-k==B.n);
  MKL_BC(omatcopy)('r',trans==CblasTrans?'t':'n',A.m,A.n,alpha,A.data(),lda,B.data(),ldb);
}

#endif

#endif
#endif
#endif

max_stars_repo_name: jjqcat/geode

max_issues_repo_name: bantamtools/geode

max_forks_repo_name: bantamtools/geode


text:
/**
 * @copyright (c) 2012- King Abdullah University of Science and
 *                      Technology (KAUST). All rights reserved.
 **/

/**
 * @file testing/testing_helper.h

 * KBLAS is a high performance CUDA library for subset of BLAS
 *    and LAPACK routines optimized for NVIDIA GPUs.
 * KBLAS is provided by KAUST.
 *
 * @version 3.0.0
 * @author Wajih Halim Boukaram
 * @author Ali Charara
 * @date 2018-11-14
 **/

#ifndef __TESTING_HELPER_H__
#define __TESTING_HELPER_H__

#include <cuda.h>
#include <cuda_runtime.h>
#include <cuda_runtime_api.h>
#include <cublas_v2.h>
#include <kblas.h>
#include <cusolverDn.h>

#ifdef USE_MAGMA
#include <magma_v2.h>
#endif
#ifdef USE_MKL
#include <mkl.h>
#else
#include <cblas.h>		//TODO: if MKL not set we need to use other libs
#include <lapacke.h>
#endif
#ifdef USE_OPENMP
#include <omp.h>
#endif
#include <kblas.h>

#ifdef __cplusplus
extern "C" {
#endif

inline int iDivUp( int a, int b ) { return (a % b != 0) ? (a / b + 1) : (a / b); }

#define kmin(a,b) ((a)>(b)?(b):(a))
#define kmax(a,b) ((a)<(b)?(b):(a))
int kblas_roundup(int x, int y);

////////////////////////////////////////////////////////////
// Error checking
////////////////////////////////////////////////////////////
#define check_error(ans) { gpuAssert((ans), __FILE__, __LINE__); }
void gpuAssert(cudaError_t code, const char *file, int line);

#define check_cublas_error(ans) { gpuCublasAssert((ans), __FILE__, __LINE__); }
void gpuCublasAssert(cublasStatus_t code, const char *file, int line);

#define check_kblas_error(ans) { gpuKblasAssert((ans), __FILE__, __LINE__); }
void gpuKblasAssert(int code, const char *file, int line);

#define check_cusolver_error(ans) { gpuCusolverAssert((ans), __FILE__, __LINE__); }
void gpuCusolverAssert(cusolverStatus_t code, const char *file, int line);

////////////////////////////////////////////////////////////
// Timers and related stuff
////////////////////////////////////////////////////////////
double gettime(void);

struct GPU_Timer;
typedef GPU_Timer* GPU_Timer_t;

GPU_Timer_t newGPU_Timer(cudaStream_t stream);
void deleteGPU_Timer(GPU_Timer_t timer);
void gpuTimerTic(GPU_Timer_t timer);
void gpuTimerRecordEnd(GPU_Timer_t timer);
double gpuTimerToc(GPU_Timer_t timer);

void avg_and_stdev(double* values, int num_vals, double* avg, double* std_dev, int warmup);

////////////////////////////////////////////////////////////
// Generate array of pointers from a strided array
////////////////////////////////////////////////////////////
void generateDArrayOfPointers(double* original_array, double** array_of_arrays, int stride, int num_arrays, cudaStream_t stream);
void generateSArrayOfPointers(float* original_array, float** array_of_arrays, int stride, int num_arrays, cudaStream_t stream);
void generateDArrayOfPointersHost(double* original_array, double** array_of_arrays, int stride, int num_arrays);
void generateSArrayOfPointersHost(float* original_array, float** array_of_arrays, int stride, int num_arrays);

// Reductions
int getMaxElement(int* a, size_t elements, cudaStream_t stream);
int getMinElement(int* a, size_t elements, cudaStream_t stream);

////////////////////////////////////////////////////////////
// Allocations
////////////////////////////////////////////////////////////
#define TESTING_MALLOC_CPU( ptr, T, size)                       \
{			\
  if ( (ptr = (T*) malloc( (size)*sizeof( T ) ) ) == NULL) {    \
    fprintf( stderr, "!!!! malloc_cpu failed for: %s\n", #ptr ); \
    exit(-1);                                                   \
  } \
}
#define TESTING_MALLOC_DEV( ptr, T, size) check_error( cudaMalloc( (void**)&ptr, (size)*sizeof(T) ) )
#define TESTING_MALLOC_PIN( ptr, T, size) check_error( cudaHostAlloc ( (void**)&ptr, (size)*sizeof( T ), cudaHostAllocPortable  ))

#define TESTING_FREE_CPU(ptr)	{ if( (ptr) ) free( (ptr) ); }
#define TESTING_FREE_DEV(ptr)	check_error( cudaFree( (ptr) ) );

////////////////////////////////////////////////////////////
// Data generation
////////////////////////////////////////////////////////////
void generateDrandom(double* random_data, long num_elements, int num_ops);
void generateSrandom(float* random_data, long num_elements, int num_ops);
void srand_matrix(long rows, long cols, float* A, long LDA);
void drand_matrix(long rows, long cols, double* A, long LDA);
void crand_matrix(long rows, long cols, cuFloatComplex* A, long LDA);
void zrand_matrix(long rows, long cols, cuDoubleComplex* A, long LDA);

void smatrix_make_hpd(int N, float* A, int lda);
void dmatrix_make_hpd(int N, double* A, int lda);
void cmatrix_make_hpd(int N, cuFloatComplex* A, int lda);
void zmatrix_make_hpd(int N, cuDoubleComplex* A, int lda);

void generateRandDimensions(int minDim, int maxDim, int* randDims, int seed, int num_ops);
void fillIntArray(int* array_vals, int value, int num_elements);

void fillGPUIntArray(int* array_vals, int value, int num_elements, cudaStream_t stream);
void copyGPUPointerArray(void** originalPtrs, void** copyPtrs, int num_ptrs, cudaStream_t stream);

// set cond = 0 to use exp decay
void generateDrandomMatrices(
	double* M_strided, int stride_M, double* svals_strided, int stride_S, int rows, int cols,
	double cond, double exp_decay, int seed, int num_ops, int threads
);
void generateSrandomMatrices(
	float* M_strided, int stride_M, float* svals_strided, int stride_S, int rows, int cols,
	float cond, float exp_decay, int seed, int num_ops, int threads
);

// Mode = 0: use provided singualr values
// Mode = 1: sets singular values to random numbers in the range (1/cond, 1)
// Mode = 2: generate singular values such that S(i) = exp(-exp_decay * i);
//	     	 where exp_decay is a random value between exp_decay_min and exp_decay_max
void generateDrandomMatricesArray(
	double** M_ptrs, int* ldm_array, double** svals_ptrs, int *rows_array, int *cols_array,
	int mode, double cond, double exp_decay_min, double exp_decay_max, int seed, int num_ops, int threads
);
void generateSrandomMatricesArray(
	float** M_ptrs, int* ldm_array, float** svals_ptrs, int *rows_array, int *cols_array,
	int mode, float cond, float exp_decay_min, float exp_decay_max, int seed, int num_ops, int threads
);

void generateSsingular_values(float** svals_ptrs, int rank, float min_sval, float max_sval, int seed, int batchCount);
void generateDsingular_values(double** svals_ptrs, int rank, double min_sval, double max_sval, int seed, int batchCount);

////////////////////////////////////////////////////////////
// Result checking
////////////////////////////////////////////////////////////
// Vectors
float sget_max_error(float* ref, float *res, int n, int inc);
double dget_max_error(double* ref, double *res, int n, int inc);
float cget_max_error(cuFloatComplex* ref, cuFloatComplex *res, int n, int inc);
double zget_max_error(cuDoubleComplex* ref, cuDoubleComplex *res, int n, int inc);
// Matrices
float sget_max_error_matrix(float* ref, float *res, long m, long n, long lda);
double dget_max_error_matrix(double* ref, double *res, long m, long n, long lda);
float cget_max_error_matrix(cuFloatComplex* ref, cuFloatComplex *res, long m, long n, long lda);
double zget_max_error_matrix(cuDoubleComplex* ref, cuDoubleComplex *res, long m, long n, long lda);

double dget_max_error_matrix_uplo(double* ref, double *res, long m, long n, long lda, char uplo);
// float sget_max_error_matrix_uplo(float* ref, float *res, long m, long n, long lda, char uplo);
// double dget_max_error_matrix_symm(double* ref, double *res, long m, long n, long lda, char uplo);

#define printMatrix(m, n, A, lda, out) { \
  for(int r = 0; r < (m); r++){ \
    for(int c = 0; c < (n); c++){ \
      fprintf((out), "%.7e  ", (A)[r + c * (lda)]); \
    } \
    fprintf((out), "\n"); \
  } \
  fprintf((out), "\n"); \
}

////////////////////////////////////////////////////////////
// Command line parser
////////////////////////////////////////////////////////////
#define MAX_NTEST 1000

typedef struct kblas_opts
{
	// matrix size
	int ntest;
	int msize[ MAX_NTEST ];
	int nsize[ MAX_NTEST ];
	int ksize[ MAX_NTEST ];
	int rsize[ MAX_NTEST ];

	// scalars
	int devices[MAX_NGPUS];
	int nstream;
	int ngpu;
	int niter;
	int nruns;
	double      tolerance;
	int check;
	int verbose;
	int nb;
	int db;
	int custom;
	int warmup;
	int time;
	int lapack;
	int magma;
	int svd;
	int cuda;
	int nonUniform;
	//int bd[KBLAS_BACKDOORS];
	int batchCount;
	int strided;
	int btest, batch[MAX_NTEST];
	int rtest, rank[MAX_NTEST];
	int omp_numthreads;
  char LR;//Low rank format
  int version;

	// lapack flags
	char uplo;
	char transA;
	char transB;
	char side;
	char diag;
} kblas_opts;

int parse_opts(int argc, char** argv, kblas_opts *opts);

#ifdef __cplusplus
}
#endif

#ifdef __cplusplus

#define kmin(a,b) ((a)>(b)?(b):(a))
#define kmax(a,b) ((a)<(b)?(b):(a))

#include "kblas_operators.h"
template<class T, class T_complex>
T get_magnitude(T_complex a){ return sqrt(a.x * a.x + a.y * a.y); }
template<class T>
T get_magnitude(T ax, T ay){ return sqrt(ax * ax + ay * ay); }

template<typename T>
bool kblas_laisnan(T val1, T val2){
  return val1 != val2;
}

template<typename T>
bool kblas_isnan(T val){
  return kblas_laisnan(val,val);
}

inline float Xabs(float a){return fabs(a);}
inline double Xabs(double a){return fabs(a);}
inline float Xabs(cuFloatComplex a){return get_magnitude<float,cuFloatComplex>(a);}
inline double Xabs(cuDoubleComplex a){return get_magnitude<double,cuDoubleComplex>(a);}

template<class T>
void matrix_make_hpd(int N, T* A, int lda, T diag)
{
  ssize_t ldas = ssize_t(lda);
  for(ssize_t i = 0; i < ssize_t(N); i++)
  {
    A[i + i * ldas] = A[i + i * ldas] + diag;
    for(ssize_t j = 0; j < i; j++){
      A[j + i*ldas] = A[i + j*ldas];
    }
  }
}

template<typename T, typename R>
R kblas_lange(char type, int M, int N, T* arr, int lda){
  R value = make_zero<R>();
  R temp;
  for(int j = 0; j < N; j++){
    for(int i = 0; i < M; i++){
      temp = Xabs(arr[i + j * lda]);
      if( kblas_isnan(temp) ){
        value = temp;
        printf("NAN encountered (%d,%d)\n", i, j);
        return value;
      }
      if( value < temp)
        value = temp;
    }
  }
  return value;
}

template<typename T>
void kblasXaxpy (int n, T alpha, const T *x, int incx, T *y, int incy){
  int ix = 0, iy = 0;
  if(incx < 0) ix = 1 - n * incx;
  if(incy < 0) iy = 1 - n * incy;
  for(int i = 0; i < n; i++, ix+=incx, iy+=incy){
    y[iy] += alpha * x[ix];
  }
}

template<class T>
void hilbertMatrix(int m, int n, T* A, int lda, T scal){
  for(int r = 0; r < m; r++){
    for(int c = 0; c < n; c++){
      A[r + c*lda] = T(scal)/T(r+c+1);
    }
  }
}

#endif //__cplusplus

#ifdef DBG_MSG
#define ECHO_I(_val) printf("%s(%d) ", #_val, (_val));fflush( stdout )
#define ECHO_f(_val) printf("%s(%e) ", #_val, (_val));fflush( stdout )
#define ECHO_p(_val) printf("%s(%p) ", #_val, (_val));fflush( stdout )
#define ECHO_LN printf("\n");fflush( stdout )
#else
#define ECHO_I(_val)
#define ECHO_f(_val)
#define ECHO_p(_val)
#define ECHO_LN
#endif

#endif // __TESTING_HELPER_H__

max_stars_repo_name: ecrc/kblas-gpu

max_issues_repo_name: ecrc/kblas-gpu

max_forks_repo_name: ecrc/kblas-gpu


text:
/*
Copyright (c) 2017, LibMars Developers.
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

* Redistributions of source code must retain the above copyright notice, this
  list of conditions and the following disclaimer.

* 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.

* Neither the name of the copyright holder nor the names of its
  contributors may be used to endorse or promote products derived from
  this software without specific prior written permission.

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.
*/

#ifndef MARS_HEADER_H
#define MARS_HEADER_H 1

#include <iostream>
#include <cstring>
#include <math.h>
#include <inttypes.h>
#include <fstream>
#include <sstream>
#include <iostream>
#include <map>
#include <iomanip>
#include <iterator>
#include <chrono>
#include <typeinfo>
#include <type_traits>
#include <complex>
#include <vector>
#include <algorithm>
#include <assert.h>
#include <cassert>
#include <csignal>
#include <ctime>
#include <stdarg.h>
#include <stdexcept>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string>
#include <string.h>
#include <list>
#include <set>
#include <utility>
#include <ctype.h>
#include <deque>
#include <initializer_list>
#include <unordered_map>
#include <unordered_set>
#include <cstdint>
#include <limits.h>

#include <new>
#include <locale>

#define UNW_LOCAL_ONLY

//#include <thread>
#include <cxxabi.h>

#include <dirent.h>
#include <errno.h>
#include <exception>
#include <execinfo.h>
#include <fcntl.h>

#include <libgen.h>
#include <libunwind.h>

#include <regex>
#include <semaphore.h>

#include <signal.h>

#include <sys/mman.h>
#include <sys/shm.h>
#include <sys/stat.h>
#include <sys/time.h>
#include <sys/types.h>
#include <sys/wait.h>
#include <sys/resource.h>

#include <unistd.h>

#include <gsl/gsl_cdf.h>
#include <gsl/gsl_histogram.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_sf.h>
#include <gsl/gsl_statistics.h>

#include <fftw3.h>

#include <boost/range/irange.hpp>
#include <boost/functional/hash.hpp> /* hash std::pair<long,long> */

#define MARS_COLOR_MODE   1
#ifdef  MARS_COLOR_MODE

#define MARS_COLOR_RED     "\033[22;31m"
#define MARS_COLOR_GREEN   "\033[22;32m"
#define MARS_COLOR_YELLOW  "\033[22;33m"
#define MARS_COLOR_BLUE    "\033[22;34m"
#define MARS_COLOR_MAGENTA "\033[22;35m"
//#define MARS_COLOR_CYAN    "\033[22;36m"
#define MARS_COLOR_CYAN    "\e[1;36m" // bright cyan
#define MARS_COLOR_WHITE   "\e[1;37m"

#define codef          "\033[0;22m"
#else
#define MARS_COLOR_RED     ""
#define MARS_COLOR_GREEN   ""
#define MARS_COLOR_YELLOW  ""
#define MARS_COLOR_BLUE    ""
#define MARS_COLOR_MAGENTA ""
#define MARS_COLOR_CYAN    ""

#define CODEF         ""
#endif

#define COFUNC   MARS_COLOR_CYAN
#define COFILE   MARS_COLOR_CYAN
#define COFOLDER MARS_COLOR_BLUE
#define COVAR    MARS_COLOR_RED
#define COVAL    MARS_COLOR_WHITE
#define COERR    MARS_COLOR_RED
#define COWARN   MARS_COLOR_YELLOW

#if __STDC_VERSION__ < 199901L
#if __GNUC__ >= 2
#define __FCN__ __PRETTY_FUNCTION__
#else
#define __FCN__ __func__
#endif
#endif

#endif

max_stars_repo_name: marsrobot/Tensor

max_issues_repo_name: marsrobot/Tensor

max_forks_repo_name: marsrobot/Tensor