text:
#include "cwc.h"
#include "../ccv_internal.h"
#ifdef HAVE_GSL
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#endif
#include "cwc_ext.h"
#ifdef USE_DISPATCH
#include <dispatch/dispatch.h>
#endif

#ifdef HAVE_GSL
static void _cwc_convnet_random_image_manipulation(gsl_rng* rng, ccv_dense_matrix_t* image, float image_manipulation)
{
	assert(rng && CCV_GET_CHANNEL(image->type) == CCV_C3 && image_manipulation > 0 && image_manipulation <= 1);
	int ord[3] = {0, 1, 2};
	gsl_ran_shuffle(rng, ord, 3, sizeof(int));
	int i;
	for (i = 0; i < 3; i++)
		// change the applying order
		switch (ord[i])
		{
			case 0:
				// introduce some brightness changes to the original image
				ccv_scale(image, (ccv_matrix_t**)&image, 0, gsl_rng_uniform_pos(rng) * image_manipulation * 2 + (1 - image_manipulation));
				break;
			case 1:
				// introduce some saturation changes to the original image
				ccv_saturation(image, &image, 0, gsl_rng_uniform_pos(rng) * image_manipulation * 2 + (1 - image_manipulation));
				break;
			case 2:
				// introduce some contrast changes to the original image
				ccv_contrast(image, &image, 0, gsl_rng_uniform_pos(rng) * image_manipulation * 2 + (1 - image_manipulation));
				break;
		}
}

void cwc_convnet_batch_formation(gsl_rng* rng, ccv_array_t* categorizeds, ccv_dense_matrix_t* mean_activity, ccv_dense_matrix_t* eigenvectors, ccv_dense_matrix_t* eigenvalues, float image_manipulation, float color_gain, int* idx, ccv_size_t dim, int min_dim, int max_dim, int rows, int cols, int channels, int category_count, int symmetric, int batch, int offset, int size, float* b, int* c)
{
	assert(size > 0 && size <= batch);
	assert(min_dim >= rows && min_dim >= cols);
	assert(max_dim >= min_dim);
	float* channel_gains = (float*)alloca(sizeof(float) * channels * size);
	memset(channel_gains, 0, sizeof(float) * channels * size);
	int i;
	gsl_rng** rngs = (gsl_rng**)alloca(sizeof(gsl_rng*) * size);
	memset(rngs, 0, sizeof(gsl_rng*) * size);
	if (rng)
		for (i = 0; i < size; i++)
		{
			rngs[i] = gsl_rng_alloc(gsl_rng_default);
			gsl_rng_set(rngs[i], gsl_rng_get(rng));
		}
	parallel_for(i, size) {
		int j, k;
		assert(offset + i < categorizeds->rnum);
		ccv_categorized_t* categorized = (ccv_categorized_t*)ccv_array_get(categorizeds, idx ? idx[offset + i] : offset + i);
		assert(categorized->c < category_count && categorized->c >= 0); // now only accept classes listed
		if (c)
			c[i] = categorized->c;
		ccv_dense_matrix_t* image = 0;
		switch (categorized->type)
		{
			case CCV_CATEGORIZED_DENSE_MATRIX:
				image = categorized->matrix;
				break;
			case CCV_CATEGORIZED_FILE:
				image = 0;
				ccv_read(categorized->file.filename, &image, CCV_IO_ANY_FILE | CCV_IO_RGB_COLOR);
				break;
		}
		if (image)
		{
			// first resize to between min_dim and max_dim
			ccv_dense_matrix_t* input = 0;
			ccv_size_t resize = dim;
			if (rngs[i]) // randomize the resized dimensions
			{
				int d = gsl_rng_uniform_int(rngs[i], max_dim - min_dim + 1) + min_dim;
				resize = ccv_size(d, d);
			}
			// if neither side is the same as expected, we have to resize first
			if (resize.height != image->rows && resize.width != image->cols)
				ccv_convnet_input_formation(resize, image, &input);
			else
				input = image;
			if (rngs[i] && image_manipulation > 0)
				_cwc_convnet_random_image_manipulation(rngs[i], input, image_manipulation);
			ccv_dense_matrix_t* patch = 0;
			if (input->cols != cols || input->rows != rows)
			{
				int x = rngs[i] ? gsl_rng_uniform_int(rngs[i], input->cols - cols + 1) : (input->cols - cols + 1) / 2;
				int y = rngs[i] ? gsl_rng_uniform_int(rngs[i], input->rows - rows + 1) : (input->rows - rows + 1) / 2;
				ccv_slice(input, (ccv_matrix_t**)&patch, CCV_32F, y, x, rows, cols);
			} else
				ccv_shift(input, (ccv_matrix_t**)&patch, CCV_32F, 0, 0); // converting to 32f
			if (input != image) // only unload if we created new input
				ccv_matrix_free(input);
			// we loaded image in, deallocate it now
			if (categorized->type != CCV_CATEGORIZED_DENSE_MATRIX)
				ccv_matrix_free(image);
			// random horizontal reflection
			if (symmetric && rngs[i] && gsl_rng_uniform_int(rngs[i], 2) == 0)
				ccv_flip(patch, &patch, 0, CCV_FLIP_X);
			int x = rngs[i] ? gsl_rng_uniform_int(rngs[i], mean_activity->cols - cols + 1) : (mean_activity->cols - cols + 1) / 2;
			int y = rngs[i] ? gsl_rng_uniform_int(rngs[i], mean_activity->rows - rows + 1) : (mean_activity->rows - rows + 1) / 2;
			ccv_dense_matrix_t mean_patch = ccv_reshape(mean_activity, y, x, rows, cols);
			ccv_subtract(patch, &mean_patch, (ccv_matrix_t**)&patch, 0);
			assert(channels == CCV_GET_CHANNEL(patch->type));
			if (color_gain > 0 && rngs[i] && eigenvectors && eigenvalues)
			{
				assert(channels == 3); // only support RGB color gain
				memset(channel_gains + channels * i, 0, sizeof(float) * channels);
				for (j = 0; j < channels; j++)
				{
					float alpha = gsl_ran_gaussian(rngs[i], color_gain) * eigenvalues->data.f64[j];
					for (k = 0; k < channels; k++)
						channel_gains[k + i * channels] += eigenvectors->data.f64[j * channels + k] * alpha;
				}
			}
			for (j = 0; j < channels; j++)
				for (k = 0; k < rows * cols; k++)
					b[(j * rows * cols + k) * batch + i] = patch->data.f32[k * channels + j] + channel_gains[j + i * channels];
			ccv_matrix_free(patch);
		} else
			PRINT(CCV_CLI_ERROR, "cannot load %s.\n", categorized->file.filename);
	} parallel_endfor
	if (rng)
		for (i = 0; i < size; i++)
			gsl_rng_free(rngs[i]);
}
#endif

void cwc_convnet_mean_formation(ccv_array_t* categorizeds, ccv_size_t dim, int channels, int symmetric, ccv_dense_matrix_t** b)
{
	int i, count = 0;
	ccv_dense_matrix_t* c = ccv_dense_matrix_new(dim.height, dim.width, channels | CCV_64F, 0, 0);
	ccv_zero(c);
	ccv_dense_matrix_t* db = *b = ccv_dense_matrix_renew(*b, dim.height, dim.width, channels | CCV_32F, channels | CCV_32F, 0);
	for (i = 0; i < categorizeds->rnum; i++)
	{
		if (i % 23 == 0 || i == categorizeds->rnum - 1)
			FLUSH(CCV_CLI_INFO, " - compute mean activity %d / %d", i + 1, categorizeds->rnum);
		ccv_categorized_t* categorized = (ccv_categorized_t*)ccv_array_get(categorizeds, i);
		ccv_dense_matrix_t* image = 0;
		switch (categorized->type)
		{
			case CCV_CATEGORIZED_DENSE_MATRIX:
				image = categorized->matrix;
				break;
			case CCV_CATEGORIZED_FILE:
				ccv_read(categorized->file.filename, &image, CCV_IO_ANY_FILE | CCV_IO_RGB_COLOR);
				break;
		}
		if (!image)
		{
			PRINT(CCV_CLI_ERROR, "cannot load %s.\n", categorized->file.filename);
			continue;
		}
		ccv_dense_matrix_t* patch = 0;
		if (image->cols != dim.width || image->rows != dim.height)
		{
			int x = (image->cols - dim.width + 1) / 2;
			int y = (image->rows - dim.height + 1) / 2;
			assert(x == 0 || y == 0);
			ccv_slice(image, (ccv_matrix_t**)&patch, CCV_32F, y, x, dim.height, dim.width);
		} else
			ccv_shift(image, (ccv_matrix_t**)&patch, CCV_32F, 0, 0); // converting to 32f
		if (categorized->type != CCV_CATEGORIZED_DENSE_MATRIX)
			ccv_matrix_free(image);
		ccv_add(patch, c, (ccv_matrix_t**)&c, CCV_64F);
		++count;
		ccv_matrix_free(patch);
	}
	if (symmetric)
	{
		int j, k;
		double p = 0.5 / count;
		double* cptr = c->data.f64;
		float* dbptr = db->data.f32;
		for (i = 0; i < db->rows; i++)
		{
			for (j = 0; j < db->cols; j++)
				for (k = 0; k < channels; k++)
					dbptr[j * channels + k] = p * (cptr[j * channels + k] + cptr[(c->cols - j - 1) * channels + k]);
			dbptr += db->cols * channels;
			cptr += c->cols * channels;
		}
	} else {
		double p = 1.0 / count;
		for (i = 0; i < dim.height * dim.width * channels; i++)
			db->data.f32[i] = p * c->data.f64[i];
	}
	ccv_matrix_free(c);
	PRINT(CCV_CLI_INFO, "\n");
}

void cwc_convnet_channel_eigen(ccv_array_t* categorizeds, ccv_dense_matrix_t* mean_activity, ccv_size_t dim, int channels, ccv_dense_matrix_t** eigenvectors, ccv_dense_matrix_t** eigenvalues)
{
	assert(channels == 3); // this function cannot handle anything other than 3x3 covariance matrix
	double* mean_value = (double*)alloca(sizeof(double) * channels);
	memset(mean_value, 0, sizeof(double) * channels);
	assert(CCV_GET_CHANNEL(mean_activity->type) == channels);
	assert(mean_activity->rows == dim.height);
	assert(mean_activity->cols == dim.width);
	int i, j, k, c, count = 0;
	for (i = 0; i < dim.height * dim.width; i++)
		for (k = 0; k < channels; k++)
			mean_value[k] += mean_activity->data.f32[i * channels + k];
	for (i = 0; i < channels; i++)
		mean_value[i] = mean_value[i] / (dim.height * dim.width);
	double* covariance = (double*)alloca(sizeof(double) * channels * channels);
	memset(covariance, 0, sizeof(double) * channels * channels);
	for (c = 0; c < categorizeds->rnum; c++)
	{
		if (c % 23 == 0 || c == categorizeds->rnum - 1)
			FLUSH(CCV_CLI_INFO, " - compute covariance matrix for data augmentation (color gain) %d / %d", c + 1, categorizeds->rnum);
		ccv_categorized_t* categorized = (ccv_categorized_t*)ccv_array_get(categorizeds, c);
		ccv_dense_matrix_t* image = 0;
		switch (categorized->type)
		{
			case CCV_CATEGORIZED_DENSE_MATRIX:
				image = categorized->matrix;
				break;
			case CCV_CATEGORIZED_FILE:
				ccv_read(categorized->file.filename, &image, CCV_IO_ANY_FILE | CCV_IO_RGB_COLOR);
				break;
		}
		if (!image)
		{
			PRINT(CCV_CLI_ERROR, "cannot load %s.\n", categorized->file.filename);
			continue;
		}
		ccv_dense_matrix_t* patch = 0;
		if (image->cols != dim.width || image->rows != dim.height)
		{
			int x = (image->cols - dim.width + 1) / 2;
			int y = (image->rows - dim.height + 1) / 2;
			assert(x == 0 || y == 0);
			ccv_slice(image, (ccv_matrix_t**)&patch, CCV_32F, y, x, dim.height, dim.width);
		} else
			ccv_shift(image, (ccv_matrix_t**)&patch, CCV_32F, 0, 0); // converting to 32f
		if (categorized->type != CCV_CATEGORIZED_DENSE_MATRIX)
			ccv_matrix_free(image);
		for (i = 0; i < dim.width * dim.height; i++)
			for (j = 0; j < channels; j++)
				for (k = j; k < channels; k++)
					covariance[j * channels + k] += (patch->data.f32[i * channels + j] - mean_value[j]) * (patch->data.f32[i * channels + k] - mean_value[k]);
		++count;
		ccv_matrix_free(patch);
	}
	for (i = 0; i < channels; i++)
		for (j = 0; j < i; j++)
			covariance[i * channels + j] = covariance[j * channels + i];
	double p = 1.0 / ((double)count * dim.height * dim.width);
	for (i = 0; i < channels; i++)
		for (j = 0; j < channels; j++)
			covariance[i * channels + j] *= p; // scale down
	ccv_dense_matrix_t covm = ccv_dense_matrix(3, 3, CCV_64F | CCV_C1, covariance, 0);
	ccv_eigen(&covm, eigenvectors, eigenvalues, CCV_64F, 1e-8);
	PRINT(CCV_CLI_INFO, "\n");
}

max_stars_repo_name: ChiahungTai/ccv

max_issues_repo_name: ChiahungTai/ccv

max_forks_repo_name: ChiahungTai/ccv


text:
// clang -DHAVE_INLINE -O3 mussel.c -o mussel -lgsl -lgslcblas

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

int main(int argc, char *argv[])
{
    // Setup the random number generator:
    gsl_rng_env_setup();
    gsl_rng *rng = gsl_rng_alloc(gsl_rng_taus2);
    const unsigned int seed = (argc == 2)? atof(argv[1]) : time(NULL);
    gsl_rng_set(rng, seed);

    // Parameters:
    const int T = 1000;
    const double A0 = 0.5;
    const double A1 = 0.025;
    const double A2 = 0.2;
    const int N = 1000;
    const int TWO_N = N * N;

    // Allocate memory for the matrix:
    int **m = (int**)malloc(N * sizeof(int*));
    int i;
    for (i = 0; i < N; ++i)
    {
        m[i] = (int*)malloc(N * sizeof(int));
    }

    // Fill the matrix:
    int x0, y0;
    for (x0 = 0; x0 < N; ++x0)
    {
        for (y0 = 0; y0 < N; ++y0)
        {
            m[x0][y0] = (gsl_rng_uniform(rng) < 0.5) ? 0 : 2;
        }
    }

    // Main loop
    int n, x, y, x1, y1, t, cell;
    for (t = 0; t < T; ++t)
    {
        // Loop around the lattice at random to update cells;
        for (cell = 0; cell < TWO_N; ++cell)
        {
            // Select at random the cell to update
            x0 = (int)(gsl_rng_uniform(rng) * N);
            y0 = (int)(gsl_rng_uniform(rng) * N);

            // Update the mussels
            // Disturbed cells automatically transformed into empty cells
            if (m[x0][y0] == 0)
            {
                m[x0][y0] = 1;
            }
            else if (m[x0][y0] == 1) // If the cell is empty, test if there is colonization
            {
                // First calculate the number of neighbours
                n = 0;
                for (x = -1; x <= 1; ++x)
                {
                    for (y = -1; y <= 1; ++y)
                    {
                        if (x0 + x == -1)
                        {
                            x1 = N - 1;
                        }
                        else if (x0 + x == N)
                        {
                            x1 = 0;
                        }
                        else
                        {
                            x1 = x0 + x;
                        }

                        if (y0 + y == -1)
                        {
                            y1 = N - 1;
                        }
                        else if (y0 + y == N)
                        {
                            y1 = 0;
                        }
                        else
                        {
                            y1 = y0 + y;
                        }

                        if (x != 0 && y != 0 && m[x1][y1] == 2)
                        {
                            ++n;
                        }
                    }
                }
                // Calculate if the status if the cell is changed
                if (gsl_rng_uniform(rng) < A2 * n * 0.125)
                {
                    m[x0][y0] = 2;
                }
            }
            else if (m[x0][y0] == 2)
            { // If the cell is occupied, test if there is disturbance
                n = 0;
                for (x = -1; x <= 1; x++)
                {
                    for (y = -1; y <= 1; y++)
                    {
                        if (x0 + x == -1)
                        {
                            x1 = N - 1;
                        }
                        else if (x0 + x == N)
                        {
                            x1 = 0;
                        }
                        else
                        {
                            x1 = x0 + x;
                        }

                        if (y0 + y == -1)
                        {
                            y1 = N - 1;
                        }
                        else if (y0 + y == N)
                        {
                            y1 = 0;
                        }
                        else
                        {
                            y1 = y0 + y;
                        }

                        if (x != 0 && y != 0 && m[x1][y1] == 0)
                        {
                            n = 1;
                        }
                    }
                }
                // Calculate if the status of the cell is changed
                if (gsl_rng_uniform(rng) < A0 * n * 0.125 + A1)
                {
                    m[x0][y0] = 0;
                }
            }
        }
    }

    char buffer[100];
    sprintf(buffer, "mussel-%u.txt", seed);
    FILE *out = fopen(buffer, "w");

    for (x0 = 0; x0 < N; ++x0)
    {
        for (y0 = 0; y0 < N; ++y0)
        {
            fprintf(out, "%d ", m[x0][y0]);
        }
        fprintf(out, "\n");
    }
    fclose(out);

    gsl_rng_free(rng);
    for (i = 0; i < N; ++i)
    {
        free(m[i]);
    }
    free(m);
    return EXIT_SUCCESS;
}

ext: c


text:
#ifndef UTILITIES_H
#define UTILITIES_H
#include <vector>
#include <math.h>
#include <random>
#include <iostream>
#include <fstream>
#include <eigen3/Eigen/Dense>
#include <boost/random/mersenne_twister.hpp>
#include <boost/random/normal_distribution.hpp>
#include <boost/math/distributions/normal.hpp>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <opencv2/opencv.hpp>
#include <opencv2/ml/ml.hpp>
#include  "Structures.h"
#include <map>

using namespace std;
using namespace boost::math;
using namespace Eigen;
using namespace cv;

namespace Eigen {
    namespace internal {
        template<typename Scalar>
        struct scalar_normal_dist_op
        {
          static boost::mt19937 rng;    // The uniform pseudo-random algorithm
          mutable boost::normal_distribution<Scalar> norm;  // The gaussian combinator
          EIGEN_EMPTY_STRUCT_CTOR(scalar_normal_dist_op)
          template<typename Index>
          inline const Scalar operator() (Index, Index = 0) const { return norm(rng); }
        };
        template<typename Scalar> boost::mt19937 scalar_normal_dist_op<Scalar>::rng;
        template<typename Scalar>
        struct functor_traits<scalar_normal_dist_op<Scalar> >{
            enum{
                Cost = 50 * NumTraits<Scalar>::MulCost, PacketAccess = false, IsRepeatable = false
            };
        };
    }
}
class Utilities{
    public:
        Utilities();
        virtual ~Utilities();
        //Read a file that has a specific number of dimensions for each attribute
        virtual vector < vector < vector<double> > >  readFile(string CloudSeperator, string filepath);
        double multivariateNormalPDF(Vector3d instance, Vector3d mu, Matrix3d covar, int dimensionality);
        Eigen::MatrixXd sampleMultivariateNormal( Eigen::VectorXd mean, Eigen::MatrixXd covar, int numOfSamples, int dimensionality);
        Eigen::MatrixXd sampleMultinomial(vector<double> probabilities, int samples);
        int randcat( vector<double> * vec);
        Eigen::Matrix3d iwishrnd( Matrix3d tau, double nu, int dimensionality, int df);
        Eigen::VectorXd exprnd(double rate, int samples);
        double exppdf(double x , double lambda);
        double gammarnd(double alpha , double beta);
        RowVectorXd dirrnd(RowVectorXd q0);
        double catpdf(int index , vector<double>  probabilities);
        // Exponential distribution distances
        double Expsquaredhellinger(double lambda1, double lambda2);
        double ExpKLDivergence(double lambda1, double lambda2);
        // Gaussian distribution distances
        double GaussKLDivergence(std::vector<double> mean1, Matrix3d covar1, std::vector<double> mean2, Matrix3d covar2 );
        double Wasserstein(std::vector<double> mean1, Matrix3d covar1, std::vector<double> mean2, Matrix3d covar2 );

        vector<float>  categoricalhistogramCompare( float histA[], float histB[] , int N);
        float categoricalKLDivergence( cv::Mat * mat1, cv::Mat * mat2);
    protected:
    private:
};
#endif // UTILITIES_H

ext: h

max_stars_repo_name: hadjichristslave/SMC

max_issues_repo_name: hadjichristslave/SMC

max_forks_repo_name: hadjichristslave/SMC


text:
/**
 *
 * @file core_ssetvar.c
 *
 *  PLASMA core_blas kernel
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Mark Gates
 * @date 2010-11-15
 * @generated s Tue Jan  7 11:44:49 2014
 *
 **/
#include <lapacke.h>
#include "common.h"

/***************************************************************************//**
 *
 * @ingroup CORE_float
 *
 * CORE_ssetvar sets a single variable, x := alpha.
 *
 *******************************************************************************
 *
 * @param[in] alpha
 *         Scalar to set x to, passed by pointer so it can depend on runtime value.
 *
 * @param[out] x
 *         On exit, x = alpha.
 *
 ******************************************************************************/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_ssetvar = PCORE_ssetvar
#define CORE_ssetvar PCORE_ssetvar
#endif
void CORE_ssetvar(const float *alpha, float *x)
{
    *x = *alpha;
}


text:
#include "../globals.h"
#include <gsl/gsl_sort.h>
#include <gsl/gsl_sort_vector.h>
#include <gsl/gsl_vector.h>

void setup_Sedov_Blast()
{
    Problem.Boxsize[0] = 3;
    Problem.Boxsize[1] = 3;
    Problem.Boxsize[2] = 3;

    sprintf ( Problem.Name, "IC_SedovBlast" );

    const double rho = 1.24E7;

    Problem.Rho_Max = rho;

    Density_Func_Ptr = &Sedov_Blast_Density;
    U_Func_Ptr = &Sedov_Blast_U;
    PostProcessing_Func_Ptr = &Sedov_Blast_PostProcessing;
}

float Sedov_Blast_Density ( const int ipart , const double bias )
{
    return 1.24E7;
}

// This function calculates the distance from each particle to the center of the Box. It returns the distance to the NNpart times furthest neighbour of zero.

float Sedov_Blast_abs ()
{

    int NNpart = 296;

    gsl_vector *abs_of_zero = gsl_vector_alloc ( Param.Npart );

    for ( int i = 0; i < Param.Npart; i++ ) {

        double r = sqrt ( ( P[i].Pos[0] - 0.5 * Problem.Boxsize[0] ) * ( P[i].Pos[0] - 0.5 * Problem.Boxsize[0] ) + ( P[i].Pos[1] * 0.5 * Problem.Boxsize[1] ) * ( P[i].Pos[1] - 0.5 * Problem.Boxsize[1] )
                          + ( P[i].Pos[2] - 0.5 * Problem.Boxsize[2] ) * ( P[i].Pos[2] - 0.5 * Problem.Boxsize[2] ) );
        gsl_vector_set ( abs_of_zero, i, r );
    }

    gsl_sort_vector ( abs_of_zero );

    double dist = gsl_vector_get ( abs_of_zero, NNpart - 1 );

    gsl_vector_free ( abs_of_zero );

    return dist;

}

#ifdef KINETIC_SEDOV
float Sedov_Blast_kinetic ()
{

    int NNpart = 32;
    Radius =
        double partpos[3][32];

    float maxDistance = Sedov_Blast_abs ();

    if ( Radius <= maxDistance ) {

        for ( int i =; i < 3; i++ ) {

            for ( int j = 0; j < NNpart, j++ ) {

                double pos = P[i].Pos[j];
                partpos[i][j] = pos;

            }

        }

    }

}
#endif // KINETIC_SEDOV

float Sedov_Blast_U ( const int ipart )
{
    return 0.0;
}

//! @todo improvement: use gsl_sort_vector_index
void Sedov_Blast_PostProcessing ()
{
    /*const double u_sn = 4.18971E5;
    int sn_count = pow ( Param.Npart / 3200., 3.0 );
    sn_count = min ( sn_count, 1 );
    //! @todo assign u_sn to sn_count innermost particles*/

    float maxDistance = Sedov_Blast_abs ();

    for ( int ipart = 0; ipart < Param.Npart; ++ipart ) {

        const double x = P[ipart].Pos[0] - Problem.Boxsize[0] * 0.5;
        const double y = P[ipart].Pos[1] - Problem.Boxsize[1] * 0.5;
        const double z = P[ipart].Pos[2] - Problem.Boxsize[2] * 0.5;
        double Radius = sqrt ( x * x + y * y + z * z );
        // int NNpart = 32;
        const double  Supernova_Energy = 0.00502765; // This is the value of one supernova energy per unit mass in Gadget units 1e51 erg
        // const double hsml = SphP[i].Hsml;

        if ( Radius <=  maxDistance ) {
            SphP[ipart].U = Supernova_Energy;
        }
    }

}

max_stars_repo_name: elehcim/WVTICs

max_issues_repo_name: elehcim/WVTICs

max_forks_repo_name: elehcim/WVTICs


text:
/*
** random swap functions
**
** M. Kuhlmann, MPI-KYB, Sept 2015
**/

#include <viaio/Vlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <stdlib.h>

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_combination.h>

/* binomial coefficient, n choose k */
double binom (int n,int k)
{
  double prod,ix,nx,kx;
  int i;
  nx = (double)n;
  kx = (double)k;

  prod = ix = 1;
  for (i=1; i<=k; i++) {
    prod *= (nx-kx+ix)/ix;
    ix++;
  }
  return prod;
}

double factorial(int n) {
  int i;
  double ix = 1;
  double prod = 1;
  for (i=1; i<=n; i++) {
    prod *= ix;
    ix++;
  }
  return prod;
}

char **enumerate_swaps(int n)
{
  size_t i,j,num,b;

  long maxiter = (size_t) gsl_pow_int(2.0,(int)n);

  char **swaps = (char **)VCalloc(maxiter,sizeof(char *));
  for(i = 0; i < maxiter; i++)
    swaps[i] = (char *)VCalloc(n,sizeof(char));

  /* generate swap table */
  for (i=0; i<maxiter; i++) {
    /* decimal to binary */
    j = 0;
    num = i;
    while(num > 0)
    {
      b=num%2;
      swaps[i][j] = b;
      num=num/2;
      j++;
    }
  }

  return swaps;
}

char **randomSwaps(long maxiter,int n,unsigned long int seed)
{
  size_t i,j;

  char **swaps = (char **)VCalloc(maxiter,sizeof(char *));
  for(i = 0; i < maxiter; i++)
    swaps[i] = (char *)VCalloc(n,sizeof(char));

  gsl_rng_env_setup();
  const gsl_rng_type *T = gsl_rng_default;
  gsl_rng *rx = gsl_rng_alloc(T);
  gsl_rng_set(rx,(unsigned long int)seed);

  for(i = 0; i < maxiter; i++)
    for(j = 0; j < n; j++)
      swaps[i][j] = gsl_ran_bernoulli(rx, 0.5);

  return swaps;
}

char **randomSwaps_noduplicates(long maxiter,int n,unsigned long int seed)
{
  size_t iter,i,j;

  /* if number of maxiter exceeds number of all possible combinations */
  if(maxiter > (size_t) gsl_pow_int(2.0,(int)n))
    VError("Error: Number of iterations larger than possible permutations!");

  gsl_rng_env_setup();
  const gsl_rng_type *T = gsl_rng_default;
  gsl_rng *rx = gsl_rng_alloc(T);
  gsl_rng_set(rx,(unsigned long int)seed);

  char **swaps = (char **)VCalloc(maxiter,sizeof(char *));
  for(i = 0; i < maxiter; i++)
    swaps[i] = (char *)VCalloc(n,sizeof(char));

  /* generate swap table */
  iter = 0;
  while (iter < maxiter) {
    for (j=0; j<n; j++)
      swaps[iter][j] = gsl_ran_bernoulli(rx,0.5);
    /* check if already there */
    /* decrementing iter will result in overwriting the swaps */
    /* since there are no duplicates in can only be decremented once */
    for (j=0; j<iter; j++)
      if (memcmp(swaps[iter],swaps[j],n*sizeof(char)) == 0)
	iter--;
    iter++;
  }

  return swaps;
}

/*
 * shuffle blocks; needed for permutation_test
 */
int **randomBlocks(long int maxiter, int nimages, int block_size, unsigned long int seed)
{
  int nblocks = nimages/block_size;
  int i,j,iter;

  /* ini with real order */
  int **table = (int **)VCalloc(maxiter,sizeof(int*));
  for(i = 0; i < maxiter; i++)
  {
    table[i] = (int *)VCalloc(2*nimages,sizeof(int));
    for(j = 0; j < nimages; j++)
    {
      table[i][2*j] = j;
      table[i][2*j+1] = j+nimages;
    }
  }

  /* each row of table has order: 0, nimages, 1, nimages+1, ... */
  /* shuffle blocks of size 2*block_size */
  gsl_rng_env_setup();
  const gsl_rng_type *T = gsl_rng_default;
  gsl_rng *rx = gsl_rng_alloc(T);
  gsl_rng_set(rx,(unsigned long int)seed);

  for (iter=0; iter<maxiter; iter++)
    for (j=0; j<nblocks; j++)
      gsl_ran_shuffle(rx, (void *)&table[iter][j*2*block_size], 2*block_size, sizeof(int));
  return table;
}
/*
 * shuffle blocks; needed for permutation_test
 * maxiter is the number of random permutations
 * nimages is the number of images in one condition (2*nimages in total)
 * block_size is the size of the permutation blocks
 */
int **randomBlocks_noduplicates(long int maxiter, int nimages, int block_size, unsigned long int seed)
{
  int nblocks = nimages/block_size;
  int i,j,iter;

  /* ini with real order */
  int **table = (int **)VCalloc(maxiter,sizeof(int*));
  for(i = 0; i < maxiter; i++)
  {
    table[i] = (int *)VCalloc(2*nimages,sizeof(int));
    for(j = 0; j < nimages; j++)
    {
      table[i][2*j] = j;
      table[i][2*j+1] = j+nimages;
    }
  }

  /* each row of table has order: 0, nimages, 1, nimages+1, ... */
  /* shuffle blocks of size 2*block_size */
  gsl_rng_env_setup();
  const gsl_rng_type *T = gsl_rng_default;
  gsl_rng *rx = gsl_rng_alloc(T);
  gsl_rng_set(rx,(unsigned long int)seed);

  iter = 0;
  while (iter < maxiter) {
    /* shuffle */
    for (j = 0; j<nblocks; j++)
      gsl_ran_shuffle(rx, (void *)&table[iter][j*2*block_size], 2*block_size, sizeof(int));

    /* sort for one side */
    int k, l, m;
    int index = 0;
    int temp = 0;
    for (m = 0; m < nblocks; m++) {
      for (k = 0; k < block_size-1; k++) {
	index = k;
	for (l = k+1; l < block_size; l++) {
	  if (table[iter][m*2*block_size+2*index] < table[iter][m*2*block_size+2*l])
	    index = l;
	}
	/* swap pairs 'k' and 'index' */
	temp = table[iter][m*2*block_size+2*k];
	table[iter][m*2*block_size+2*k] = table[iter][m*2*block_size+2*index];
	table[iter][m*2*block_size+2*index] = temp;
	temp = table[iter][m*2*block_size+2*k+1];
	table[iter][m*2*block_size+2*k+1] = table[iter][m*2*block_size+2*index+1];
	table[iter][m*2*block_size+2*index+1] = temp;
      }
    }

    /*
     * decrementing iter will result in overwriting the last swaps row
     * since there are no duplicates it can only be decremented once
     */
    for(j=0; j<iter; j++)
      if (memcmp(table[iter],table[j],2*nimages*sizeof(int)) == 0)
	iter--;
    iter++;
  }
  return table;
}

max_stars_repo_name: zrajna/lipsia

max_issues_repo_name: zrajna/lipsia

max_forks_repo_name: zrajna/lipsia


text:
#if !defined(RANDOM_H_INCLUDED)
#define RANDOM_H_INCLUDED

#include <gsl/gsl>
#include <iosfwd>

class Random
{
public:
	static int usage(::std::ostream& strm, const ::std::string& progName,
		const char* pMsg);

	Random(::gsl::span<const char*const> args);
	int run() const;

	Random(const Random&) = delete;
	Random& operator=(const Random&) = delete;
	Random(Random&&) = delete;
	Random& operator=(Random&&) = delete;

private:
	static long getRandomInteger(long low, long high);

	long	m_lowerBound;
	long	m_upperBound;
	long	m_count;
};

#endif // RANDOM_H_INCLUDED

max_stars_repo_name: IanEmmons/CmdLineUtil

max_issues_repo_name: IanEmmons/CmdLineUtil

max_forks_repo_name: IanEmmons/CmdLineUtil


text:
/*
    This section is aimed to calculate the gaussian integrals that occurs in everywhere

    This part is greatly inspired by May, Andrew James. Density fitting in explicitly correlated electronic structure theory. Diss. University of Bristol, 2006.
 */

/*

Here is the table of the relation between the names and the definition:

S   <a|b>
J   (a|b) = <a|1/r_{12}|b>
G   (a|e^{-\gamma r_{12}^2}|b)
F   (a|f_{12}|b) = (a|c_i exp(-\gamma_i r_{12}^2)|b)
GJ  (a|e^{-\gamma r_{12}^2 r_12^{-1}}|b)
FJ  (a|f_{12} r_{12}^{-1}|b)
FF  (a|f_{12}^2|b)
FT  (ab|[\hat{t}_1, f_{12}]|c) = (ab|[- \frac{1}{2} \nabla ^2, f_{12}]|c)
FTF (a|\frac{1}{2}[f_12,[\hat{t}_1 + \hat{t}_2, f_{12}]]|b)
F-F (a|f_{12}|b|f_{23}|c)
J-F (a|r_{12}^{-1}|b|f_{23}|c)
X   (ab|[\hat{t}_1,r_{12}]|c)
Y   (ab|[\hat{t}_1,r_{12}^{-1}])
Z   (a|\frac{1}{r_{1Z}}|b)

Currently only S integrals and J integrals are supported.

 */

#ifndef __BASIS_H__
#define __BASIS_H__
#include "basis.h"
#endif

#include <math.h>

#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_sf.h>

typedef struct gaussian_chain{
    double R[3];
    int a[3];

    double exponent;

    double coefficient;

    gaussian_chain * NEXT;

}gaussian_chain;

gaussian_chain * gaussian_chain_calloc();

void gaussian_chain_free(gaussian_chain * HEAD);

double gaussian_chain_get(gaussian_chain * HEAD,double x, double y, double z);
double orbital_get(orbital * orbital,double x, double y, double z);

double Gamma(double z);
double Boys(double x, int n);
double Binomials(int n, int k);

double f(int k, int a, int b, double PA, double PB);
double tranformation_coefficient(int a[3], int b[3], int p[3], double PA[3], double PB[3], double xi, double AB);

double SIntegral(double ra[3], double rb[3], int ax, int ay, int az, int bx, int by, int bz, double alpha,double beta);
double JIntegral(double ra[3], double rb[3], int ax, int ay, int az, int bx, int by, int bz, double alpha,double beta, int m);
double ZIntegral(double ra[3], double rb[3], double rz[3], int ax, int ay, int az, int bx, int by, int bz, double alpha, double beta, int m);

double gaussian_chain_SIntegral(gaussian_chain * a, gaussian_chain * b);
double gaussian_chain_JIntegral(gaussian_chain * a, gaussian_chain * b);
double gaussian_chain_ZIntegral(gaussian_chain * a, gaussian_chain * b, double rz[3]);
double gaussian_chain_full_SIntegral(gaussian_chain * a_HEAD, gaussian_chain * b_HEAD);
double gaussian_chain_full_JIntegral(gaussian_chain * a_HEAD, gaussian_chain * b_HEAD);
double gaussian_chain_full_ZIntegral(gaussian_chain * a_HEAD, gaussian_chain * b_HEAD, double rz[3]);

void gaussian_chain_derivative(gaussian_chain * dest, gaussian_chain * src, int key);
void gaussian_chain_second_derivative(gaussian_chain * dest, gaussian_chain * src, int key);
void gaussian_chain_laplacian(gaussian_chain * dest, gaussian_chain * src);

double gaussian_chain_kinetic_energy(gaussian_chain * a_HEAD, gaussian_chain * b_HEAD);

void single_electron_transform(gaussian_chain * HEAD, orbital * a);
void two_electron_transform(gaussian_chain * HEAD, orbital * a, orbital * b);

double orbital_SIntegral(orbital * a,orbital * b);
double orbital_JIntegral(orbital * a,orbital * b);
double orbital_ZIntegral(orbital * a,orbital * b, double rz[3]);
//double orbital_GIntegral(orbital *,orbital *);
//double orbital_FIntegral(orbital *,orbital *);
//double orbital_GJIntegral(orbital *,orbital *);
//double orbital_FJIntegral(orbital *,orbital *);
//double orbital_FFIntegral(orbital *,orbital *);

double two_electron_JIntegral(orbital * a_1, orbital * b_1, orbital * c_2, orbital * d_2);

double orbital_kinetic_energy(orbital * a, orbital * b);

void orbital_S_matrix(gsl_matrix * dest, orbital * HEAD);

max_stars_repo_name: Walter-Feng/Hartree-Fock

max_issues_repo_name: Walter-Feng/Hartree-Fock

max_forks_repo_name: Walter-Feng/Hartree-Fock


text:
/*
 pbrt source code is Copyright(c) 1998-2016
 Matt Pharr, Greg Humphreys, and Wenzel Jakob.

 This file is part of pbrt.

 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.

 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.

 */

#if defined(_MSC_VER)
#define NOMINMAX
#pragma once
#endif

#ifndef PBRT_CAMERAS_REALISTICEYE_H
#define PBRT_CAMERAS_REALISTICEYE_H

// cameras/realistic.h*
#include "pbrt.h"
#include "camera.h"
#include "film.h"
#include <gsl/gsl_randist.h>
#include "spectrum.h" // This is necessary to declare Spectrum class in this header file.

namespace pbrt {

    // This is slightly different than our normal lens class. We have more variables to be able to describe the elements of the eye including biconic surfaces and unique ocular mediums.
    struct LensElementEye{
        float radiusX;
        float radiusY;
        float thickness;
        float mediumIndex; // Corresponds to the mediumElement. Describes media directly "behind" (-z direction, toward the retina) the surface.
        float semiDiameter;
        float conicConstantX;
        float conicConstantY;
    };

    // RealisticCamera Declarations
    class RealisticEye : public Camera {
    public:
        // RealisticCamera Public Methods
        RealisticEye(const AnimatedTransform &CameraToWorld,
                     Float shutterOpen,
                     Float shutterClose,
                     bool simpleWeighting,
                     bool noWeighting,
                     Film *film,
                     const Medium *medium,
                     std::string specfile,
                     Float pupilDiameter,
                     Float retinaDistance,
                     Float retinaRadius,
                     Float retinaSemiDiam,
                     std::vector<Spectrum> iorSpectra,
                     bool flipRad,
                     bool mmUnits,
                     bool diffractionEnabled);

        Float GenerateRay(const CameraSample &sample, Ray *) const;

    private:

        const bool simpleWeighting;
        const bool noWeighting;

        // Lens information
        std::vector<LensElementEye> lensEls;
        Float effectiveFocalLength;

        // Specific parameters for the human eye
        Float pupilDiameter;
        Float retinaDistance;
        Float retinaRadius;
        Float retinaSemiDiam;
        Float retinaDiag; // This will take the place of "film->diag"
        Float frontThickness; // The distance from the back of the lens to the front of the eye.
        std::vector<Spectrum> iorSpectra;

        // Flags for conventions
        bool diffractionEnabled;
        float lensScaling;

        // Private methods for tracing through lens
        bool IntersectLensElAspheric(const Ray &r, Float *tHit, LensElementEye currElement, Float zShift, Vector3f *n) const;
        void applySnellsLaw(Float n1, Float n2, Float lensRadius, Vector3f &normalVec, Ray * ray ) const;
        Float lookUpIOR(int mediumIndex, const Ray &ray) const;
        void diffractHURB(Point3f intersect, Float apertureRadius, const Float wavelength, const Vector3f oldDirection, Vector3f *newDirection) const;

        // Handy method to explicity solve for the z(x,y) at a given point (x,y), for the biconic SAG
        Float BiconicZ(Float x, Float y, LensElementEye currElement) const;

        // GSL seed(?) for random number generation
        gsl_rng * r;

    };

    RealisticEye *CreateRealisticEye(const ParamSet &params,
                                     const AnimatedTransform &cam2world,
                                     Film *film, const Medium *medium);

}  // namespace pbrt

#endif  // PBRT_CAMERAS_REALISTICEYE_H


text:
#ifndef _CCS_RNG_H
#define _CCS_RNG_H

#ifdef __cplusplus
extern "C" {
#endif

#include <gsl/gsl_rng.h>

/**
 * @file rng.h
 * CCS rng define random number generators. For now they are wrappers over gsl
 * random number generators.
 */

/**
 * Create a new random number generator using the gsl default type (see
 * gsl_rng_default).
 * @param [out] rng_ret a pointer to the variable that will contain the returned
 *              random number generator
 * @return #CCS_SUCCESS on success
 * @return -#CCS_INVALID_VALUE if \p rng_ret is NULL
 * @return -#CCS_OUT_OF_MEMORY if there was not enough memory to allocate the
 *                             new random number generator
 */
extern ccs_result_t
ccs_rng_create(ccs_rng_t *rng_ret);

/**
 * Create a new random number generator using the provided gsl type (see
 * gsl_rng_type).
 * @param [in] rng_type a pointer to the type of gsl random number generator to
 *                      use.
 * @param [out] rng_ret a pointer to the variable that will contain the returned
 *              random number generator
 * @return #CCS_SUCCESS on success
 * @return -#CCS_INVALID_VALUE if \p rng_ret or \p rng_type are NULL
 * @return -#CCS_OUT_OF_MEMORY if there was not enough memory to allocate the
 *                             new random number generator
 */
extern ccs_result_t
ccs_rng_create_with_type(const gsl_rng_type *rng_type,
                         ccs_rng_t          *rng_ret);

/**
 * Get the gsl type of a random number generator.
 * @param [in] rng
 * @param [out] rng_type_ret a pointer that will contained a pointer to the
 *                           returned gsl random number generator type.
 * @return #CCS_SUCCESS on success
 * @return -#CCS_INVALID_OBJECT if \p rng is not a valid CCS random number
 *                              generator
 * @return -#CCS_INVALID_VALUE if \p rng_type_ret is NULL
 */
extern ccs_result_t
ccs_rng_get_type(ccs_rng_t            rng,
                 const gsl_rng_type **rng_type_ret);

/**
 * Set the seed of a random number generator.
 * @param [in] rng
 * @param [in] seed the seed to use with the random number generator
 * @return #CCS_SUCCESS on success
 * @return -#CCS_INVALID_OBJECT if \p rng is not a valid CCS random number
 *                              generator
 */
extern ccs_result_t
ccs_rng_set_seed(ccs_rng_t         rng,
                 unsigned long int seed);

/**
 * Get a random integer from a random number generator. Integer is contained
 * between the value returned by #ccs_rng_min and #ccs_rng_max, both included.
 * @param [in] rng
 * @param [out] value_ret a pointer to the variable that will contain the
 *                        returned value
 * @return #CCS_SUCCESS on success
 * @return -#CCS_INVALID_OBJECT if \p rng is not a valid CCS random number
 *                              generator
 * @return -#CCS_INVALID_VALUE if \p value_ret is NULL
 */
extern ccs_result_t
ccs_rng_get(ccs_rng_t          rng,
            unsigned long int *value_ret);

/**
 * Get a random floating point value uniformly sampled in the interval [0.0,
 * 1.0).
 * @param [in] rng
 * @param [out] value_ret a pointer to the variable that will contain the
 *                        returned value
 * @return #CCS_SUCCESS on success
 * @return -#CCS_INVALID_OBJECT if \p rng is not a valid CCS random number
 *                              generator
 * @return -#CCS_INVALID_VALUE if \p value_ret is NULL
 */
extern ccs_result_t
ccs_rng_uniform(ccs_rng_t    rng,
                ccs_float_t *value_ret);

/**
 * Get the underlying gsl random number generator.
 * @param [in] rng
 * @param [out] gsl_rng_ret a pointer to the variable that will contain a
 *                          pointer to the underlying random number generator
 * @return #CCS_SUCCESS on success
 * @return -#CCS_INVALID_OBJECT if \p rng is not a valid CCS random number
 *                              generator
 * @return -#CCS_INVALID_VALUE if \p gsl_rng_ret is NULL
 */
extern ccs_result_t
ccs_rng_get_gsl_rng(ccs_rng_t   rng,
                    gsl_rng   **gsl_rng_ret);

/**
 * Get the minimum value that can be returned by #ccs_rng_get.
 * @param [in] rng
 * @param [out] value_ret a pointer to the variable that will contain the
 *                        returned value
 * @return #CCS_SUCCESS on success
 * @return -#CCS_INVALID_OBJECT if \p rng is not a valid CCS random number
 *                              generator
 * @return -#CCS_INVALID_VALUE if \p value_ret is NULL
 */
extern ccs_result_t
ccs_rng_min(ccs_rng_t          rng,
            unsigned long int *value_ret);

/**
 * Get the maximum value that can be returned by #ccs_rng_get.
 * @param [in] rng
 * @param [out] value_ret a pointer to the variable that will contain the
 *                        returned value
 * @return #CCS_SUCCESS on success
 * @return -#CCS_INVALID_OBJECT if \p rng is not a valid CCS random number
 *                              generator
 * @return -#CCS_INVALID_VALUE if \p value_ret is NULL
 */
extern ccs_result_t
ccs_rng_max(ccs_rng_t          rng,
            unsigned long int *value_ret);

#ifdef __cplusplus
}
#endif

#endif //_CCS_RNG_H

max_stars_repo_name: deephyper/CCS

max_issues_repo_name: deephyper/CCS

max_forks_repo_name: deephyper/CCS


text:
/*
** Kendall's coefficient of concordance - Kendall's W
**
** G.Lohmann, July 2010
*/
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_sort_double.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_cdf.h>

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

extern void gsl_sort_vector_index (gsl_permutation *,gsl_vector *);

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

double kendall(size_t **xrank,double *r,int n,int m,int otype)
{
  int i,j;
  double S,W,R,chi_square=0;
  double nx=0,mx=0;

  nx = (double)n;
  mx = (double)m;

  for (i=0; i<n; i++) {
    r[i] = 0;
    for (j=0; j<m; j++) r[i] += xrank[j][i];
  }

  R = 0;
  for (i=0; i<n; i++) R += r[i];
  R /= nx;

  S = 0;
  for (i=0; i<n; i++) S += SQR(r[i] - R);

  W = 12.0*S/(mx*mx*(nx*nx-1)*nx);
  if (otype == 0) return W;

  /* by central limit theorem, chi_square is approx normal for large n */
  /* see en.wikipedia.org/wiki/Chi-square_distribution */
  chi_square = (nx-1)*mx*W;
  return (chi_square - nx)/(sqrt(2.0*nx));
}

double VKendall_W(gsl_matrix *data,int nvoxels,int otype)
{
  int i,j,n,m;
  double W0;
  static double *r=NULL;
  static gsl_vector *vec = NULL;
  static size_t **xrank=NULL;
  static gsl_permutation *perm=NULL,*rank=NULL;

  m = data->size1;
  n = nvoxels;

  if (vec == NULL) {
    r    = (double *) calloc(data->size2,sizeof(double));
    vec  = gsl_vector_calloc(data->size2);
    perm = gsl_permutation_alloc(data->size2);
    rank = gsl_permutation_alloc(data->size2);

    xrank = (size_t **) calloc(m,sizeof(size_t *));
    for (i=0; i<m; i++) xrank[i] = (size_t *) calloc(data->size2,sizeof(size_t));
  }

  gsl_vector_set_all (vec,9999);
  for (i=0; i<m; i++) {
    for (j=0; j<n; j++) {
      gsl_vector_set(vec,j,gsl_matrix_get(data,i,j));
    }
    gsl_sort_vector_index (perm, vec);
    gsl_permutation_inverse (rank, perm);
    for (j=0; j<n; j++) xrank[i][j] = rank->data[j];
  }

  /* Kendall's W of measured data */
  W0 = kendall(xrank,r,n,m,otype);
  return W0;
}


text:
/* $Header$ */

/* Purpose: GSL functions that handle missing values - stubs only
   Extra NCO code derived from GSL removed 2020-04-02 to avoid license conflict between GSL GPL and NCO BSD code */

/* Copyright (C) 1995--present Charlie Zender
   This file is part of NCO, the netCDF Operators. NCO is free software.
   You may redistribute and/or modify NCO under the terms of the
   3-Clause BSD License with exceptions described in the LICENSE file */

#include <gsl/gsl_statistics.h>
#include <gsl/gsl_fit.h>
#include "nco_gsl.h" /* Missing value-aware GSL functions */

/* Fit the data (x_i, y_i) to the linear relationship

   Y = c0 + c1 x

   returning,

   c0, c1  --  coefficients
   cov00, cov01, cov11  --  variance-covariance matrix of c0 and c1,
   sumsq   --   sum of squares of residuals

   This fit can be used in the case where the errors for the data are
   uknown, but assumed equal for all points. The resulting
   variance-covariance matrix estimates the error in the coefficients
   from the observed variance of the points around the best fit line.
*/

int
nco_gsl_fit_linear
(const double *x,
 const size_t xstride,
 const double *y,
 const size_t ystride,
 const size_t n,
 double *c0,
 double *c1,
 double *cov_00,
 double *cov_01,
 double *cov_11,
 double *sumsq,
 const double *mss_val)
{

  return gsl_fit_linear(x, xstride, y, ystride, n, c0, c1, cov_00, cov_01, cov_11, sumsq);

}

static double
compute_covariance
(const double data1[],
 const size_t stride1,
 const double data2[],
 const size_t stride2,
 const size_t n,
 const double mean1,
 const double mean2,
 const double *mss_val) /* Missing value */
{
  /* takes a dataset and finds the covariance */

  long double covariance = 0 ;

  size_t i;

  /* find the sum of the squares */
  if (mss_val==NULL)
  {
    for (i = 0; i < n; i++)
    {
      const long double delta1 = (data1[i * stride1] - mean1);
      const long double delta2 = (data2[i * stride2] - mean2);
      covariance += (delta1 * delta2 - covariance) / (i + 1);
    }

  }
  else
  {
    for (i = 0; i < n; i++)
    {
      if (data2[i * stride2]!=*mss_val)
      {
        const long double delta1 = (data1[i * stride1] - mean1);
        const long double delta2 = (data2[i * stride2] - mean2);
        covariance += (delta1 * delta2 - covariance) / (i + 1);
      }
    }
  }

  return covariance ;
}

double
nco_gsl_stats_covariance_m
(const double data1[],
 const size_t stride1,
 const double data2[],
 const size_t stride2,
 const size_t n,
 const double mean1,
 const double mean2,
 const double *mss_val) /* Missing value */
{
  return gsl_stats_covariance_m(data1, stride1, data2,stride2,n,mean1,mean2);

}

double
nco_gsl_stats_covariance
(const double data1[],
 const size_t stride1,
 const double data2[],
 const size_t stride2,
 const size_t n,
 const double *mss_val) /* Missing value */
{
  return gsl_stats_covariance(data1, stride1, data2,stride2,n);

}

double
nco_gsl_stats_mean
(const double data[],
 const size_t stride,
 const size_t size,
 const double *mss_val) /* Missing value */
{

  return gsl_stats_mean(data, stride, size);

}

max_stars_repo_name: acchen/nco

max_issues_repo_name: acchen/nco

max_forks_repo_name: acchen/nco


text:
#pragma once

#include "Core/Delegate.h"
#include "Core/Pointers.h"
#include "Scene/Entity.h"

#include <gsl/span>

#include <string>
#include <unordered_map>
#include <vector>

class CameraComponent;
class DirectionalLightComponent;
class ModelComponent;
class PointLightComponent;
class SpotLightComponent;

class Scene
{
public:
   Scene();
   ~Scene();

   void tick(float dt);

   float getTime() const
   {
      return time;
   }

   float getDeltaTime() const
   {
      return deltaTime;
   }

   template<typename... ComponentTypes>
   Entity* createEntity()
   {
      entities.push_back(Entity::create<ComponentTypes...>(*this));
      return entities.back().get();
   }

   Entity* createEntity(gsl::span<std::string> componentClassNames)
   {
      entities.push_back(Entity::create(componentClassNames, *this));
      return entities.back().get();
   }

   bool destroyEntity(Entity* entityToDestroy);

   const std::vector<UPtr<Entity>>& getEntities() const
   {
      return entities;
   }

   CameraComponent* getActiveCameraComponent()
   {
      return activeCameraComponent;
   }

   const CameraComponent* getActiveCameraComponent() const
   {
      return activeCameraComponent;
   }

   void setActiveCameraComponent(CameraComponent* newActiveCameraComponent);

   const std::vector<CameraComponent*>& getCameraComponents() const
   {
      return cameraComponents;
   }

   void registerCameraComponent(CameraComponent* cameraComponent);
   void unregisterCameraComponent(CameraComponent* cameraComponent);

   const std::vector<ModelComponent*>& getModelComponents() const
   {
      return modelComponents;
   }

   void registerModelComponent(ModelComponent* modelComponent);
   void unregisterModelComponent(ModelComponent* modelComponent);

   const std::vector<DirectionalLightComponent*>& getDirectionalLightComponents() const
   {
      return directionalLightComponents;
   }

   void registerDirectionalLightComponent(DirectionalLightComponent* directionalLightComponent);
   void unregisterDirectionalLightComponent(DirectionalLightComponent* directionalLightComponent);

   const std::vector<PointLightComponent*>& getPointLightComponents() const
   {
      return pointLightComponents;
   }

   void registerPointLightComponent(PointLightComponent* pointLightComponent);
   void unregisterPointLightComponent(PointLightComponent* pointLightComponent);

   const std::vector<SpotLightComponent*>& getSpotLightComponents() const
   {
      return spotLightComponents;
   }

   void registerSpotLightComponent(SpotLightComponent* spotLightComponent);
   void unregisterSpotLightComponent(SpotLightComponent* spotLightComponent);

private:
   float time;
   float deltaTime;

   std::vector<UPtr<Entity>> entities;

   std::vector<CameraComponent*> cameraComponents;
   CameraComponent* activeCameraComponent;

   std::vector<ModelComponent*> modelComponents;

   std::vector<DirectionalLightComponent*> directionalLightComponents;
   std::vector<PointLightComponent*> pointLightComponents;
   std::vector<SpotLightComponent*> spotLightComponents;
};

max_stars_repo_name: aaronmjacobs/Swap

max_issues_repo_name: aaronmjacobs/Swap

max_forks_repo_name: aaronmjacobs/Swap


text:
#include <cblas.h>
int main() { double *x, *y; cblas_ddot(0, x, 0, y, 0); }

max_stars_repo_name: t-sakashita/rokko

max_issues_repo_name: t-sakashita/rokko

max_forks_repo_name: t-sakashita/rokko


text:
/* specfunc/gsl_sf_psi.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2004 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_PSI_H__
#define __GSL_SF_PSI_H__

#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

/* Poly-Gamma Functions
 *
 * psi(m,x) := (d/dx)^m psi(0,x) = (d/dx)^{m+1} log(gamma(x))
 */

/* Di-Gamma Function  psi(n) = psi(0,n)
 *
 * n > 0
 * exceptions: GSL_EDOM
 */
int     gsl_sf_psi_int_e(const int n, gsl_sf_result * result);
double  gsl_sf_psi_int(const int n);

/* Di-Gamma Function psi(x) = psi(0, x)
 *
 * x != 0.0, -1.0, -2.0, ...
 * exceptions: GSL_EDOM, GSL_ELOSS
 */
int     gsl_sf_psi_e(const double x, gsl_sf_result * result);
double  gsl_sf_psi(const double x);

/* Di-Gamma Function Re[psi(1 + I y)]
 *
 * exceptions: none
 */
int     gsl_sf_psi_1piy_e(const double y, gsl_sf_result * result);
double  gsl_sf_psi_1piy(const double y);

/* Di-Gamma Function psi(z) for general complex argument z = x + iy
 *
 * exceptions: GSL_EDOM
 */
int gsl_sf_complex_psi_e(
  const double x,
  const double y,
  gsl_sf_result * result_re,
  gsl_sf_result * result_im
  );

/* Tri-Gamma Function psi^(1)(n)
 *
 * n > 0
 * exceptions: GSL_EDOM
 */
int     gsl_sf_psi_1_int_e(const int n, gsl_sf_result * result);
double  gsl_sf_psi_1_int(const int n);

/* Tri-Gamma Function psi^(1)(x)
 *
 * x != 0.0, -1.0, -2.0, ...
 * exceptions: GSL_EDOM, GSL_ELOSS
 */
int     gsl_sf_psi_1_e(const double x, gsl_sf_result * result);
double  gsl_sf_psi_1(const double x);

/* Poly-Gamma Function psi^(n)(x)
 *
 * n >= 0, x > 0.0
 * exceptions: GSL_EDOM
 */
int     gsl_sf_psi_n_e(const int n, const double x, gsl_sf_result * result);
double  gsl_sf_psi_n(const int n, const double x);

__END_DECLS

#endif /* __GSL_SF_PSI_H__ */

max_stars_repo_name: shi-bash-cmd/qtTest

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <pthread.h>
#include <semaphore.h>
#include <apop.h>

#include <gsl/gsl_vector.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_multifit.h>

/********************************HOPE: Sum of norm of residuals fitting for Hough Mode Extension Data*****************************************/
//Copyright (C) 2012 Daniel Strano

/*'Hope' is an open-source project based on Oberheide and Forbes GEOPHYSICAL RESEARCH LETTERS, VOL. 35, L04816, doi:10.1029/2007GL032397, 2008. The program was originally designed for the purpose of handling the challenge of empirically fitting multi-component complex wave data, in situations where the goodness of fit of a least squares model is questionable. The program was originally designed to work with three self-consistent variables to fit simultaneously. The variances of the variables are first scaled to satisify the assumptions of constant variance and independence of errors, necessary for a consistent least squares model. Then, the program takes advantage of the finite phase domain to indirectly search for a maximization of the sum of absolute (norms of) residuals.*/

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

#define depVar 3
const int rawblock = 29;
const int block = 33;
const int heightsOut = 8;
const int heightsIn = 8;
const int HME_N = 2;
const int HMEBlockCount = 8;
const int satrows = 264;

const double PI = 3.141592653589793;
const double PHASE_INTERVAL = 0.001;
const int PHASE_INT_INV = 1000;
const double P_I_x_PI = 0.001 * 3.141592653589793;

int outrows=0;
double totalSumSqr;
double TotalSN[depVar];

int global_phase_iter1 = 0;
int global_phase_iter2 = 0;

int threadtracker[4];

pthread_mutex_t iter_update;
pthread_mutex_t compare;

sem_t threadsAvail;

gsl_matrix *Full_Descrip;
gsl_vector *dependent;
gsl_vector *BestModel;

int minIndepOffset;
double *threadIndep;

double BestPhase1;
double BestPhase2;

double BestResid = 0;
double BestRNorm = 0;

int itemscount;

void *threadParallel(void *arg) {

	int *threadNum = (int *)arg; //Argument is used for thread scheduling (used at bottom of thread function)

	int phase_iter1, phase_iter2; //Local phase iteration holder variables

	//Lock the current global phase iteration counters, pull the next value for the local iteration, and update global counters:
	pthread_mutex_lock(&iter_update);
		phase_iter1 = global_phase_iter1;
		phase_iter2 = global_phase_iter2;
		++global_phase_iter1;
		if (global_phase_iter1 >= PHASE_INT_INV) {
			global_phase_iter1 = 0;
			++global_phase_iter2;
		}
	pthread_mutex_unlock(&iter_update);

	int loop1, loop2, loop3, loop4, loop5;
	int returncode;
	float phase_shift[2]= {P_I_x_PI * phase_iter1, P_I_x_PI * phase_iter2}; //Phase shifts calculated from phase iteration counters
	float Rotation[8]; //Used to rotate complex vectors

	int holder, holder2, holder3, holder4, holder5, holder6; //For loop optimization

	double chiSqr; //Sum of squares of residuals of each iterative model
	float coeff_of_determ; //Normalized coefficient of determination ""

	//Calculate rotation factors to apply phase shift corresponding to the current iteration:
	for (loop1 = 0; loop1 < HME_N; loop1++) {
		holder = loop1 * 4;
		holder2 = holder + 2;
		for (loop2 = 0; loop2 < 2; loop2++) {
			Rotation[holder + loop2] = cos(phase_shift[loop1]);
			if (loop2) {
				Rotation[holder2 + loop2] = -1 * sin(phase_shift[loop1]);
			}
			else {
				Rotation[holder2 + loop2] = sin(phase_shift[loop1]);
			}
		}
	}

	//Declare objects for GSL least squares model for fixed phases:
	gsl_matrix *Local_Descrip = gsl_matrix_alloc(itemscount, HME_N); //Local description of amplitude-only model
	gsl_matrix *covariance = gsl_matrix_alloc(HME_N, HME_N);
	gsl_multifit_linear_workspace *modelspace = gsl_multifit_linear_alloc(itemscount, HME_N);
	gsl_vector *beta = gsl_vector_alloc(HME_N); //fitting parameters

	//Apply rotations and initialize 'Local_Descrip'
	for (loop5 = 0; loop5 < HME_N; loop5++) {
		holder = loop5 * 4;
		holder2 = holder + 1;
		holder3 = holder2 + 1;
		holder4 = holder3 + 1;
		holder6 = loop5 * minIndepOffset;
		for (loop1 = 0; loop1 < itemscount; loop1 += 2) {
			holder5 = loop1 + 1;
			gsl_matrix_set(Local_Descrip, loop1, loop5,
				Rotation[holder] * threadIndep[holder6 + loop1] + Rotation[holder4] * threadIndep[holder6 + holder5]);
			gsl_matrix_set(Local_Descrip, holder5, loop5,
				Rotation[holder2] * threadIndep[holder6 + holder5] + Rotation[holder3] * threadIndep[holder6 + loop1]);
		}
	}

	//Find least squares fit:
	returncode = gsl_multifit_linear(Local_Descrip, dependent, beta, covariance, &chiSqr, modelspace);
	coeff_of_determ = 1 - chiSqr/totalSumSqr;

	//Free model space for fit (covariance is also not used):
	gsl_matrix_free(covariance);
	gsl_multifit_linear_free(modelspace);

	//From partial model beta parameters, calculate equivalent full model parameters:
	double sumNormResid[3] = {0, 0, 0};
	double R_norm = 0;
	double R_resid, I_resid;
	double Full_Beta[2 * HME_N];
	for (loop1 = 0; loop1 < HME_N; loop1++) {
		Full_Beta[2 * loop1] = gsl_vector_get(beta, loop1) * cos(phase_shift[loop1]);
		Full_Beta[2 * loop1 + 1] = gsl_vector_get(beta, loop1) * sin(phase_shift[loop1]);
	}

	//Re-weight variables to equal total weight for sum of norms of residuals model,
	//calculate sum of norms of residuals for full model from least squares beta coefficients:
	int maxiter1 = 2 * outrows;
	int maxiter2 = maxiter1 * depVar;
	int maxiter3 = HME_N * 2;
	for (loop2 = 0; loop2 < maxiter2; loop2 += maxiter1) {
		holder3 = (int)(loop2/maxiter1);
		for (loop1 = 0; loop1 < maxiter1; loop1 += 2) {
			R_resid = gsl_vector_get(dependent, loop1 + loop2);
			holder = loop2 + loop1;
			holder2 = holder + 1;
			for (loop3 = 0; loop3 < maxiter3; loop3++) R_resid -= Full_Beta[loop3] * gsl_matrix_get(Full_Descrip, holder, loop3);
			I_resid = gsl_vector_get(dependent, holder2);
			for (loop3 = 0; loop3 < maxiter3; loop3++) I_resid -= Full_Beta[loop3] * gsl_matrix_get(Full_Descrip, holder2, loop3);
			sumNormResid[holder3] += sqrt((R_resid * R_resid) + (I_resid * I_resid));
		}
	}

	//Calculated normalized coefficient of determination for sum of norms of residuals:
	for (loop1 = 0; loop1 < depVar; loop1++) R_norm -= sumNormResid[loop1]/TotalSN[loop1];
	R_norm = R_norm/3;
	R_norm++;

	//Lock global best model comparisons and updates,
	//compare this thread's "R-norm" and update global best model if a better model has been found:
	pthread_mutex_lock(&compare);
		if (BestRNorm < R_norm) {
			BestRNorm = R_norm;
			BestResid = coeff_of_determ;
			BestPhase1 = phase_shift[0];
			BestPhase2 = phase_shift[1];
			for (loop1 = 0; loop1 < HME_N; loop1++) {
				gsl_vector_set(BestModel, loop1, gsl_vector_get(beta, loop1));
			}
		}
	pthread_mutex_unlock(&compare);

	//Free remaining manually allocated memory objects:
	gsl_matrix_free(Local_Descrip);
	gsl_vector_free(beta);

	//Prepare to signal to loop in main() that this thread identity is becoming free, then post to semaphore:
	threadtracker[*threadNum] = 0;
	sem_post(&threadsAvail);

	return NULL;
}

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

	//Thread expects input files as specified below:
	if (argc != 4) {
		printf("usage: hope [saber_t] [tidi_e] [tidi_n]\n");
		printf("Note: Satellite data should be in *.csv form, in a single column, aligned by row, with no column or row names. Dummy value is -999. \n");
	}
	else {

		gsl_matrix *satellite[3];

		double *minIndep[2];

		int lev_rows;
		double height[heightsIn];
		double latitudes[rawblock];

		double *HME_ptr[8][HME_N];

		double **uninterp[HME_N][depVar*2];
		double **interp_lat[HME_N][depVar*2];
		double **interp_lev[HME_N][depVar*2];

		BestModel = gsl_vector_alloc(2);

		pthread_t zero;
		pthread_t one;
		pthread_t two;
		pthread_t three;

		sem_init(&threadsAvail, 0, 4);

		long long maxpage;
		long long pagein = 1;
		int comp, comp_int, tenth_percent_complete = 0, percent_complete = 0;

		int loop, loop1, loop2, loop3, loop4, loop5, loop6;

		int threadID[4];

		double TotalSumSqr[depVar];
		double Variance[depVar];
		double AbsDeviance[depVar];
		double Scale[depVar], DevScale[depVar];

		double scaledRow;

		int holder;
		int holder2;

		for (loop = 0; loop < depVar; loop++) {
			TotalSumSqr[loop] = 0;
			TotalSN[loop] = 0;
			Variance[loop] = 0;  //optional
			Scale[loop] = 0; //optional
		}

		for (loop = 0; loop < 4; loop++) {
			threadtracker[loop] = 0;
			threadID[loop] = loop;
		}

			//Read in raw data:

				//(Independent) descriptor data:
		apop_text_to_db(.text_file="HME1.csv", .tabname="X1", .has_row_names = 0, .has_col_names = 1);
		HME_ptr[0][0] = apop_vector_to_array(apop_query_to_vector("SELECT lat FROM X1"));
		HME_ptr[1][0] = apop_vector_to_array(apop_query_to_vector("SELECT lev FROM X1"));
		HME_ptr[2][0] = apop_vector_to_array(apop_query_to_vector("SELECT amp_t FROM X1"));
		HME_ptr[3][0] = apop_vector_to_array(apop_query_to_vector("SELECT phase_t FROM X1"));
		HME_ptr[4][0] = apop_vector_to_array(apop_query_to_vector("SELECT amp_e FROM X1"));
		HME_ptr[5][0] = apop_vector_to_array(apop_query_to_vector("SELECT phase_e FROM X1"));
		HME_ptr[6][0] = apop_vector_to_array(apop_query_to_vector("SELECT amp_n FROM X1"));
		HME_ptr[7][0] = apop_vector_to_array(apop_query_to_vector("SELECT phase_n FROM X1"));

		apop_text_to_db(.text_file="HME2.csv", .tabname="X2", .has_row_names = 0, .has_col_names = 1);
		HME_ptr[0][1] = apop_vector_to_array(apop_query_to_vector("SELECT lat FROM X2"));
		HME_ptr[1][1] = apop_vector_to_array(apop_query_to_vector("SELECT lev FROM X2"));
		HME_ptr[2][1] = apop_vector_to_array(apop_query_to_vector("SELECT amp_t FROM X2"));
		HME_ptr[3][1] = apop_vector_to_array(apop_query_to_vector("SELECT phase_t FROM X2"));
		HME_ptr[4][1] = apop_vector_to_array(apop_query_to_vector("SELECT amp_e FROM X2"));
		HME_ptr[5][1] = apop_vector_to_array(apop_query_to_vector("SELECT phase_e FROM X2"));
		HME_ptr[6][1] = apop_vector_to_array(apop_query_to_vector("SELECT amp_n FROM X2"));
		HME_ptr[7][1] = apop_vector_to_array(apop_query_to_vector("SELECT phase_n FROM X2"));

				//(Dependent) variable data:
		apop_text_to_db(.text_file=argv[1], .tabname="Y_t", .has_row_names = 0, .has_col_names = 1);
		satellite[0] = apop_query_to_matrix("SELECT * FROM Y_t");

		apop_text_to_db(.text_file=argv[2], .tabname="Y_e", .has_row_names = 0, .has_col_names = 1);
		satellite[1] = apop_query_to_matrix("SELECT * FROM Y_e");

		apop_text_to_db(.text_file=argv[3], .tabname="Y_n", .has_row_names = 0, .has_col_names = 1);
		satellite[2] = apop_query_to_matrix("SELECT * FROM Y_n");

		//Need to know how many non-dummy rows are in satellite data
		//Calculate total sum of squares and total sum of norms of satellite data while tallying non-dummy rows in 'outrows':
		outrows = 0;
		for (loop1 = 0; loop1 < satrows; loop1++) {
			if (gsl_matrix_get(satellite[0], loop1, 0) > -990) {
				if (gsl_matrix_get(satellite[1], loop1, 0) > -990) {
					if (gsl_matrix_get(satellite[2], loop1, 0) > -990) {
						for (loop2 = 0; loop2 < depVar; loop2++) {
							for (loop3 = 0; loop3 < 2; loop3++) {
								TotalSumSqr[loop2] += (gsl_matrix_get(satellite[loop2], loop1, loop3) * gsl_matrix_get(satellite[loop2], loop1, loop3));
							}
							TotalSN[loop2] += sqrt((gsl_matrix_get(satellite[loop2], loop1, 0) * gsl_matrix_get(satellite[loop2], loop1, 0)) + (gsl_matrix_get(satellite[loop2], loop1, 1) * gsl_matrix_get(satellite[loop2], loop1, 1)));
						}
						++outrows;
					}
				}
			}
		}

		itemscount = outrows * depVar * 2; //Total number of data items (rows * columns)

		//Allocate arrays for uninterpolated data, lattitude interpolated data, and data interpolated by both lattitude and elevation:
		for (loop1 = 0; loop1 < HME_N; loop1++) {

			for (loop2 = 0; loop2 < depVar * 2; loop2++) {
				uninterp[loop1][loop2]   = malloc(HMEBlockCount * sizeof(double *));
				interp_lat[loop1][loop2] = malloc(block * sizeof(double *));
				interp_lev[loop1][loop2] = malloc(block * sizeof(double *));
				for (loop3 = 0; loop3 < HMEBlockCount; loop3++) {
					uninterp[loop1][loop2][loop3]   = malloc(rawblock * sizeof(double));
				}
				for (loop3 = 0; loop3 < block; loop3++) {
					interp_lat[loop1][loop2][loop3] = malloc(HMEBlockCount * sizeof(double));
					interp_lev[loop1][loop2][loop3] = malloc(heightsOut * sizeof(double));
				}
			}
		}

		//Initialize block of uninterpolated data:
		for (loop1 = 0; loop1 < HMEBlockCount; loop1++) {
			holder = rawblock * loop1;
			for (loop2 = 0; loop2 < 2 * depVar; loop2 += 2) {
				for (loop3 = 0; loop3 < HME_N; loop3++) {
					for (loop4 = 0; loop4 < rawblock; loop4++) {
						uninterp[loop3][0 + loop2][loop1][(rawblock - 1) - loop4] = *(HME_ptr[2 + loop2][loop3] + (holder + loop4)) * cos(*(HME_ptr[3 + loop2][loop3] + holder + loop4) / 12 * PI);
						uninterp[loop3][1 + loop2][loop1][(rawblock - 1) - loop4] = *(HME_ptr[2 + loop2][loop3] + (holder + loop4)) * sin(*(HME_ptr[3 + loop2][loop3] + holder + loop4) / 12 * PI);
					}
				}
			}
		}

		//Read lattitudes from input:
		for (loop1 = 0; loop1 < rawblock; loop1++) {
			latitudes[(rawblock - 1) - loop1] = *(HME_ptr[0][0] + loop1);
		}

		//Read elevations from input:
		for (loop1 = 0; loop1 < heightsIn; loop1++) {
			height[loop1] = *(HME_ptr[1][0] + (loop1 * rawblock));
		}

	//GSL Spline Interpolation:

		//Lattitude interpolation:
		for (loop = 0; loop < HME_N; loop++) {
			for (loop3 = 0; loop3 < 2; loop3++) {
				for (loop2 = 0; loop2 < HMEBlockCount; loop2++) {
					for (loop4 = 0; loop4 < depVar * 2; loop4 += 2) {
						gsl_interp_accel *lat_acc = gsl_interp_accel_alloc();
						gsl_spline *lat_spline = gsl_spline_alloc(gsl_interp_cspline, rawblock);
						gsl_spline_init(lat_spline, latitudes, uninterp[loop][loop4 + loop3][loop2], rawblock);

						for (loop1 = 0; loop1 < block; loop1++) {
							interp_lat[loop][loop4 + loop3][loop1][loop2] = gsl_spline_eval(lat_spline, (loop1 * 5) - 80, lat_acc);
						}

						gsl_spline_free (lat_spline);
						gsl_interp_accel_free (lat_acc);
					}
				}
			}
		}

		//Altitude ("lev") interpolation:
		holder = 0;
		for (loop = 0; loop < HME_N; loop++) {
			for (loop3 = 0; loop3 < 2; loop3++) {
				for (loop1 = 0; loop1 < block; loop1++) {
					for (loop4 = 0; loop4 < depVar * 2; loop4 += 2) {
						gsl_interp_accel *lev_acc = gsl_interp_accel_alloc();
						gsl_spline *lev_spline = gsl_spline_alloc(gsl_interp_cspline, heightsIn);
						gsl_spline_init(lev_spline, height, interp_lat[loop][loop4 + loop3][loop1], heightsIn);

						for (loop2 = 0; loop2 < heightsOut; loop2++) {
							interp_lev[loop][loop4 + loop3][loop1][loop2] = gsl_spline_eval(lev_spline, 87.5 + (loop2 * 2.5), lev_acc);

						}

						gsl_spline_free (lev_spline);
						gsl_interp_accel_free (lev_acc);
					}
				}
			}
		}

		//Variance and deviance calculation (for normalizing different variables):
		for (loop1 = 0; loop1 < depVar; loop1++) {
			Variance[loop1] = TotalSumSqr[loop1] / (outrows/2);
			AbsDeviance[loop1] = TotalSN[loop1] / (outrows/2);
		}

		for (loop1 = 0; loop1 < depVar; loop1++) {
			Scale[loop1] =  sqrt(Variance[0]/ Variance[loop1]);
			DevScale[loop1] = AbsDeviance[0] / AbsDeviance[loop1];
		}

		//In interative fit, phases are held fixed while amplitudes are fit
		//Phases are incremented after each amplitude fit
		//Create an array with the minimum information from the independent matrix necessary to perform iterative fit:
		for (loop1 = 0; loop1 < HME_N; loop1++) {
			minIndep[loop1] = malloc(satrows * depVar * 2 * sizeof(double));
		}

		loop4 = 0;
		for (loop1 = 0; loop1 < heightsOut; loop1++) {
			for (loop3 = 0; loop3 < depVar; loop3++) {
				for (loop4 = 0; loop4 < block; loop4++) {
					for (loop5 = 0; loop5 < HME_N; loop5++) {
						for (loop = 0; loop < 2; loop++) {
							minIndep[loop5][(loop3 * 2 * satrows) + (block * 2 * loop1) + (loop4 * 2) + loop] = interp_lev[loop5][2 * loop3 + loop][block - 1 - loop4][loop1];
						}
					}
				}
			}
		}

		//Garbage collection after pre-processing:
		for (loop1 = 0; loop1 < HME_N; loop1++) {
			for (loop2 = 0; loop2 < depVar * 2; loop2++) {
				for (loop3 = 0; loop3 < HMEBlockCount; loop3++) {
					free(uninterp[loop1][loop2][loop3]);
				}
				for (loop3 = 0; loop3 < block; loop3++) {
					free(interp_lat[loop1][loop2][loop3]);
					free(interp_lev[loop1][loop2][loop3]);
				}
				free(uninterp[loop1][loop2]);
				free(interp_lat[loop1][loop2]);
				free(interp_lev[loop1][loop2]);
			}
		}

		//Dummy values in satellite data can not be passed to GSL fitting functions
		//Real and imaginary components alternate.
		//The dependent variables counted by "depVar" are concatenated as contigious blocks
		//Remove dummy values and concatenate satellite vectors;
		dependent = gsl_vector_alloc(2 * depVar * outrows);
		totalSumSqr = 0;
		for (loop2 = 0; loop2 < depVar; loop2++) {
			loop5 = 0;
			for (loop1 = 0; loop1 < satrows; loop1++) {
				if (gsl_matrix_get(satellite[0], loop1, 0) > -990) {
					if (gsl_matrix_get(satellite[1], loop1, 0) > -990) {
						if (gsl_matrix_get(satellite[2], loop1, 0) > -990) {
							for (loop3 = 0; loop3 < 2; loop3++) {
								scaledRow = Scale[loop2] * gsl_matrix_get(satellite[loop2], loop1, loop3);
								gsl_vector_set(dependent, (2 * outrows * loop2) + loop5 + loop3, scaledRow);
								totalSumSqr += scaledRow * scaledRow;
								for (loop4 = 0; loop4 < HME_N; loop4++){
									minIndep[loop4][(2 * outrows * loop2) + loop5 + loop3] = Scale[loop2] * minIndep[loop4][(2 * satrows * loop2) + 2 * loop1 + loop3];
								}
							}
							loop5+=2;
						}
					}
				}
			}
		}

		//Reallocate shrink minIndep once dummy values are removed:
		int *error_ptr;
		minIndepOffset = outrows * depVar * 2;
		for (loop1 = 0; loop1 < HME_N; loop1++) {
			 error_ptr = realloc(minIndep[loop1], sizeof(double) * minIndepOffset);
		}

		//Create a 1D array of the minimum data that needs to be passed to threads for iterative fit:
		threadIndep = malloc(sizeof(double) * minIndepOffset * HME_N);
		for (loop1 = 0; loop1 < HME_N; loop1++) {
			for (loop2 = 0; loop2 < minIndepOffset; loop2++) {
				threadIndep[minIndepOffset * loop1 + loop2] = minIndep[loop1][loop2];
			}
		}

		//Create an array with the full 'clean' independent data for an (analytical) R^2 fit to compare against:
		Full_Descrip = gsl_matrix_alloc((outrows * depVar) * 2, HME_N * 2);
		gsl_matrix *covariance = gsl_matrix_alloc(HME_N * 2, HME_N * 2);
		gsl_multifit_linear_workspace *modelspace = gsl_multifit_linear_alloc((outrows * depVar) * 2, HME_N * 2);
		gsl_vector *beta = gsl_vector_alloc(HME_N * 2);

		for (loop1 = 0; loop1 < 2 * outrows * depVar; loop1++) {
			for (loop2 = 0; loop2 < HME_N; loop2++) {
				for (loop3 = 0; loop3 < 2; loop3++) {
					if (loop3 && !(loop1 % 2)) {
						gsl_matrix_set(Full_Descrip, loop1, 2 * loop2 + loop3, -minIndep[loop2][loop1 + loop3]);
					}
					else {
						gsl_matrix_set(Full_Descrip, loop1, 2 * loop2 + loop3, minIndep[loop2][loop1 + loop3]);
					}
				}
			}
		}

		//'minIndep' is no longer needed
		//Release 'minIndep':
		for (loop1 = 0; loop1 < 2; loop1++) {
			free(minIndep[loop1]);
		}

		//Find R^2 fit for full independent array:
		double chiSqr;
		gsl_multifit_linear(Full_Descrip, dependent, beta, covariance, &chiSqr, modelspace);
		double coeff_of_determ = 1 - chiSqr/totalSumSqr;

		//Find sum of absolute norms of residuals for the returned least squares model (normalizing variables weights to new model moment):
		double sumNormResid[3] = {0, 0, 0};
		double R_norm = 0;
		double R_resid, I_resid;
		double Full_Beta[2 * HME_N];
		double amp_phase[2 * HME_N];
		for (loop1 = 0; loop1 < HME_N * 2; loop1++) {
			Full_Beta[loop1] = gsl_vector_get(beta, loop1);
		}

		for (loop1 = 0; loop1 < HME_N * 2; loop1+=2) {
			amp_phase[loop1] = sqrt(Full_Beta[loop1] * Full_Beta[loop1] + Full_Beta[loop1 + 1] * Full_Beta[loop1 + 1]);
			amp_phase[loop1+1] = atan(Full_Beta[loop1 + 1] / Full_Beta[loop1]);
		}

		int maxiter1 = 2 * outrows;
		int maxiter2 = maxiter1 * depVar;
		for (loop2 = 0; loop2 < maxiter2; loop2 += maxiter1) {
			for (loop1 = 0; loop1 < maxiter1; loop1 += 2) {
				R_resid = gsl_vector_get(dependent, loop1 + loop2);
				for (loop3 = 0; loop3 < HME_N * 2; loop3++) R_resid -= Full_Beta[loop3] * gsl_matrix_get(Full_Descrip, loop2 + loop1, loop3);
				I_resid = gsl_vector_get(dependent, loop1 + loop2 + 1);
				for (loop3 = 0; loop3 < HME_N * 2; loop3++) I_resid -= Full_Beta[loop3] * gsl_matrix_get(Full_Descrip, loop2 + loop1 + 1, loop3);
				sumNormResid[(int)(loop2/maxiter1)] += sqrt((R_resid * R_resid) + (I_resid * I_resid));
			}
		}

		for (loop1 = 0; loop1 < depVar; loop1++) R_norm -= sumNormResid[loop1]/TotalSN[loop1];
		R_norm = R_norm/3;
		R_norm += 1;

		//Pre-processing has completed

		//Start calling iterative fitting threads for R-norm:
		sem_wait(&threadsAvail);
		threadtracker[0] = 1;
		pthread_create(&zero,NULL,threadParallel, &threadID[0]);
		++pagein;

		sem_wait(&threadsAvail);
		threadtracker[1] = 1;
		pthread_create(&one,NULL,threadParallel, &threadID[1]);
		++pagein;

		sem_wait(&threadsAvail);
		threadtracker[2] = 1;
		pthread_create(&two,NULL,threadParallel, &threadID[2]);
		++pagein;

		sem_wait(&threadsAvail);
		threadtracker[3] = 1;
		pthread_create(&three,NULL,threadParallel, &threadID[3]);
		++pagein;

		//With every thread initially called, begin thread scheduling:

		maxpage = (PHASE_INT_INV) * (PHASE_INT_INV);

		comp = floor(maxpage/1000);
		comp_int = comp;

		while (pagein < maxpage) {

			sem_wait(&threadsAvail);
			if (!threadtracker[0]) {
				threadtracker[0] = 1;
				pthread_join(zero, NULL);
				pthread_create(&zero,NULL,threadParallel, &threadID[0]);
			}
			else if (!threadtracker[1]) {
				threadtracker[1] = 1;
				pthread_join(one, NULL);
				pthread_create(&one,NULL,threadParallel, &threadID[1]);
			}
			else if (!threadtracker[2]) {
				threadtracker[2] = 1;
				pthread_join(two, NULL);
				pthread_create(&two,NULL,threadParallel, &threadID[2]);
			}
			else {
				threadtracker[3] = 1;
				pthread_join(three, NULL);
				pthread_create(&three, NULL, threadParallel, &threadID[3]);
			}
			++pagein;

			//Output progress on percent or 1/10 of a percent of completion:
			if (pagein > comp) {
				comp = comp + comp_int;
				++tenth_percent_complete;
				if (tenth_percent_complete > 9) {
					++percent_complete;
					tenth_percent_complete = 0;
					//printf("Phase intercept 1: %f \n", BestPhase1);
					//printf("Phase intercept 2: %f \n", BestPhase2);
					//printf("R Squared average: %f \n", BestResid);
					//printf("%d.%d percent complete \n \n", percent_complete, tenth_percent_complete);
				}

				printf("Amplitude 1: %f \n", gsl_vector_get(BestModel, 0));
				printf("Phase intercept 1: %f \n", BestPhase1);
				printf("Amplitude 2: %f \n", gsl_vector_get(BestModel, 1));
				printf("Phase intercept 2: %f \n", BestPhase2);
				printf("Best R norm: %f \n", BestRNorm);
				printf("Corresponding R^2: %f \n", BestResid);
				printf("Best R^2: %f \n", coeff_of_determ);
				printf("Best R^2 Amplitude 1: %f \n", amp_phase[0]);
				printf("Best R^2 Phase intercept 1: %f \n", amp_phase[1]);
				printf("Best R^2 Amplitude 2: %f \n", amp_phase[2]);
				printf("Best R^2 Phase intercept 2: %f \n", amp_phase[3]);
				printf("R norm corresponding to analytical best R^2: %f \n", R_norm);
				printf("%d.%d percent complete \n \n", percent_complete, tenth_percent_complete);

			}
		}

		//All iterations have been called
		//Join last calls:
		pthread_join(zero, NULL);
		pthread_join(one, NULL);
		pthread_join(two, NULL);
		pthread_join(three, NULL);

		//Print final result:
		printf("Best phase 1 intercept: %f \n", BestPhase1);
		printf("Best phase 2 intercept: %f \n", BestPhase2);
		printf("Best R norm: %f \n", BestRNorm);
		printf("Best corresponding R squared: %f \n", BestResid);
		printf("Corresponding R^2: %f \n", BestResid);
		printf("Best R^2: %f \n", coeff_of_determ);
		printf("Best R^2 Amplitude 1: %f \n", amp_phase[0]);
		printf("Best R^2 Phase intercept 1: %f \n", amp_phase[1]);
		printf("Best R^2 Amplitude 2: %f \n", amp_phase[2]);
		printf("Best R^2 Phase intercept 2: %f \n", amp_phase[3]);
		printf("R norm corresponding to analytical best R^2: %f \n", R_norm);

		//Free remaining manually allocated memory objects:
		for (loop1 = 0; loop1 < depVar; loop1++) {
			gsl_matrix_free(satellite[loop1]);
		}
		gsl_vector_free(dependent);
		gsl_vector_free(BestModel);
		gsl_matrix_free(Full_Descrip);

	}
	return 0;
}

max_stars_repo_name: WrathfulSpatula/Hope

max_issues_repo_name: WrathfulSpatula/Hope

max_forks_repo_name: WrathfulSpatula/Hope


text:
/* specfunc/gsl_sf_expint.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 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.
 */

/* Author: G. Jungman */

#ifndef __GSL_SF_EXPINT_H__
#define __GSL_SF_EXPINT_H__

#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

/* E_1(x) := Re[ Integrate[ Exp[-xt]/t, {t,1,Infinity}] ]
 *
 * x != 0.0
 * exceptions: GSL_EDOM, GSL_EOVRFLW, GSL_EUNDRFLW
 */
int     gsl_sf_expint_E1_e(const double x, gsl_sf_result * result);
double     gsl_sf_expint_E1(const double x);

/* E_2(x) := Re[ Integrate[ Exp[-xt]/t^2, {t,1,Infinity}] ]
 *
 * x != 0.0
 * exceptions: GSL_EDOM, GSL_EOVRFLW, GSL_EUNDRFLW
 */
int     gsl_sf_expint_E2_e(const double x, gsl_sf_result * result);
double     gsl_sf_expint_E2(const double x);

/* Ei(x) := PV Integrate[ Exp[-t]/t, {t,-x,Infinity}]
 *
 * x != 0.0
 * exceptions: GSL_EDOM, GSL_EOVRFLW, GSL_EUNDRFLW
 */
int     gsl_sf_expint_Ei_e(const double x, gsl_sf_result * result);
double     gsl_sf_expint_Ei(const double x);

/* Shi(x) := Integrate[ Sinh[t]/t, {t,0,x}]
 *
 * exceptions: GSL_EOVRFLW, GSL_EUNDRFLW
 */
int     gsl_sf_Shi_e(const double x, gsl_sf_result * result);
double     gsl_sf_Shi(const double x);

/* Chi(x) := Re[ M_EULER + log(x) + Integrate[(Cosh[t]-1)/t, {t,0,x}] ]
 *
 * x != 0.0
 * exceptions: GSL_EDOM, GSL_EOVRFLW, GSL_EUNDRFLW
 */
int     gsl_sf_Chi_e(const double x, gsl_sf_result * result);
double     gsl_sf_Chi(const double x);

/* Ei_3(x) := Integral[ Exp[-t^3], {t,0,x}]
 *
 * x >= 0.0
 * exceptions: GSL_EDOM
 */
int     gsl_sf_expint_3_e(const double x, gsl_sf_result * result);
double     gsl_sf_expint_3(double x);

/* Si(x) := Integrate[ Sin[t]/t, {t,0,x}]
 *
 * exceptions: none
 */
int     gsl_sf_Si_e(const double x, gsl_sf_result * result);
double     gsl_sf_Si(const double x);

/* Ci(x) := -Integrate[ Cos[t]/t, {t,x,Infinity}]
 *
 * x > 0.0
 * exceptions: GSL_EDOM
 */
int     gsl_sf_Ci_e(const double x, gsl_sf_result * result);
double     gsl_sf_Ci(const double x);

/* AtanInt(x) := Integral[ Arctan[t]/t, {t,0,x}]
 *
 *
 * exceptions:
 */
int gsl_sf_atanint_e(const double x, gsl_sf_result * result);
double gsl_sf_atanint(const double x);

__END_DECLS

#endif /* __GSL_SF_EXPINT_H__ */


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

#pragma once

#include "core/common/inlined_containers.h"
#include "core/codegen/common/common.h"
#include "core/graph/graph.h"
#include "core/providers/nuphar/common/analysis/analysis.h"
#include <gsl/gsl>

namespace onnxruntime {
namespace nuphar {

class OutputAliasAnalysis : public NupharAnalysis {
 public:
  OutputAliasAnalysis()
      : NupharAnalysis("OutputAliasAnalysis") {}

  ~OutputAliasAnalysis() = default;

  void Evaluate(const onnxruntime::nuphar::NupharSubgraphUnit& graph) override;

  bool IsOutputNode(const onnxruntime::Node* node) const;

  bool IsOutputAlias(const onnxruntime::Node* node) const;

  const onnxruntime::NodeArg* SourceDefOfOutputAlias(const onnxruntime::NodeArg* node) const;

 private:
  // a set for output nodes
  std::set<NodeKey> output_nodes_;
  // a map from an output alias to its input
  std::map<NodeKey, const onnxruntime::NodeArg*> alias_use_defs_;

  void Traverse(gsl::span<const Node* const> nodes,
                const InlinedHashSet<std::string_view>& graph_inputs,
                const InlinedHashSet<std::string_view>& graph_outputs);

 private:
  ORT_DISALLOW_COPY_ASSIGNMENT_AND_MOVE(OutputAliasAnalysis);
};

}  // namespace nuphar
}  // namespace onnxruntime

max_stars_repo_name: SiriusKY/onnxruntime

max_issues_repo_name: SiriusKY/onnxruntime

max_forks_repo_name: SiriusKY/onnxruntime


text:
//See http://www.gnu.org/software/gsl/doc/html/usage.html#an-example-program
//T. Golfinopoulos, 29 Dec. 2017
#include <stdio.h>
#include <gsl/gsl_sf_bessel.h>

int
main (void)
{
  double x = 5.0;
  double y = gsl_sf_bessel_J0 (x);
  printf ("J0(%g) = %.18e\n", x, y);
  return 0;
}

max_stars_repo_name: golfit/C-Mod

max_issues_repo_name: golfit/C-Mod

max_forks_repo_name: golfit/C-Mod


text:
/* ode-initval/evolve.c
 *
 * 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
 */
#include <config.h>
#include <string.h>
#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_odeiv.h>

#include "odeiv_util.h"

gsl_odeiv_evolve *
gsl_odeiv_evolve_alloc (size_t dim)
{
  gsl_odeiv_evolve *e =
    (gsl_odeiv_evolve *) malloc (sizeof (gsl_odeiv_evolve));

  if (e == 0)
    {
      GSL_ERROR_NULL ("failed to allocate space for evolve struct",
                      GSL_ENOMEM);
    }

  e->y0 = (double *) malloc (dim * sizeof (double));

  if (e->y0 == 0)
    {
      free (e);
      GSL_ERROR_NULL ("failed to allocate space for y0", GSL_ENOMEM);
    }

  e->yerr = (double *) malloc (dim * sizeof (double));

  if (e->yerr == 0)
    {
      free (e->y0);
      free (e);
      GSL_ERROR_NULL ("failed to allocate space for yerr", GSL_ENOMEM);
    }

  e->dydt_in = (double *) malloc (dim * sizeof (double));

  if (e->dydt_in == 0)
    {
      free (e->yerr);
      free (e->y0);
      free (e);
      GSL_ERROR_NULL ("failed to allocate space for dydt_in", GSL_ENOMEM);
    }

  e->dydt_out = (double *) malloc (dim * sizeof (double));

  if (e->dydt_out == 0)
    {
      free (e->dydt_in);
      free (e->yerr);
      free (e->y0);
      free (e);
      GSL_ERROR_NULL ("failed to allocate space for dydt_out", GSL_ENOMEM);
    }

  e->dimension = dim;
  e->count = 0;
  e->failed_steps = 0;
  e->last_step = 0.0;

  return e;
}

int
gsl_odeiv_evolve_reset (gsl_odeiv_evolve * e)
{
  e->count = 0;
  e->failed_steps = 0;
  e->last_step = 0.0;
  return GSL_SUCCESS;
}

void
gsl_odeiv_evolve_free (gsl_odeiv_evolve * e)
{
  RETURN_IF_NULL (e);
  free (e->dydt_out);
  free (e->dydt_in);
  free (e->yerr);
  free (e->y0);
  free (e);
}

/* Evolution framework method.
 *
 * Uses an adaptive step control object
 */
int
gsl_odeiv_evolve_apply (gsl_odeiv_evolve * e,
                        gsl_odeiv_control * con,
                        gsl_odeiv_step * step,
                        const gsl_odeiv_system * dydt,
                        double *t, double t1, double *h, double y[])
{
  const double t0 = *t;
  double h0 = *h;
  int step_status;
  int final_step = 0;
  double dt = t1 - t0;  /* remaining time, possibly less than h */

  if (e->dimension != step->dimension)
    {
      GSL_ERROR ("step dimension must match evolution size", GSL_EINVAL);
    }

  if ((dt < 0.0 && h0 > 0.0) || (dt > 0.0 && h0 < 0.0))
    {
      GSL_ERROR ("step direction must match interval direction", GSL_EINVAL);
    }

  /* No need to copy if we cannot control the step size. */

  if (con != NULL)
    {
      DBL_MEMCPY (e->y0, y, e->dimension);
    }

  /* Calculate initial dydt once if the method can benefit. */

  if (step->type->can_use_dydt_in)
    {
      int status = GSL_ODEIV_FN_EVAL (dydt, t0, y, e->dydt_in);

      if (status)
        {
          return status;
        }
    }

try_step:

  if ((dt >= 0.0 && h0 > dt) || (dt < 0.0 && h0 < dt))
    {
      h0 = dt;
      final_step = 1;
    }
  else
    {
      final_step = 0;
    }

  if (step->type->can_use_dydt_in)
    {
      step_status =
        gsl_odeiv_step_apply (step, t0, h0, y, e->yerr, e->dydt_in,
                              e->dydt_out, dydt);
    }
  else
    {
      step_status =
        gsl_odeiv_step_apply (step, t0, h0, y, e->yerr, NULL, e->dydt_out,
                              dydt);
    }

  /* Check for stepper internal failure */

  if (step_status != GSL_SUCCESS)
    {
      *h = h0;  /* notify user of step-size which caused the failure */
      *t = t0;  /* restore original t value */
      return step_status;
    }

  e->count++;
  e->last_step = h0;

  if (final_step)
    {
      *t = t1;
    }
  else
    {
      *t = t0 + h0;
    }

  if (con != NULL)
    {
      /* Check error and attempt to adjust the step. */

      double h_old = h0;

      const int hadjust_status
        = gsl_odeiv_control_hadjust (con, step, y, e->yerr, e->dydt_out, &h0);

      if (hadjust_status == GSL_ODEIV_HADJ_DEC)
        {
          /* Check that the reported status is correct (i.e. an actual
             decrease in h0 occured) and the suggested h0 will change
             the time by at least 1 ulp */

          double t_curr = GSL_COERCE_DBL(*t);
          double t_next = GSL_COERCE_DBL((*t) + h0);

          if (fabs(h0) < fabs(h_old) && t_next != t_curr)
            {
              /* Step was decreased. Undo step, and try again with new h0. */
              DBL_MEMCPY (y, e->y0, dydt->dimension);
              e->failed_steps++;
              goto try_step;
            }
          else
            {
              h0 = h_old; /* keep current step size */
            }
        }
    }

  *h = h0;  /* suggest step size for next time-step */

  return step_status;
}

max_stars_repo_name: ielomariala/Hex-Game


text:
/*
 * Copyright 2021 The DAPHNE Consortium
 *
 * 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.
 */

#ifndef SRC_RUNTIME_LOCAL_KERNELS_GEMV_H
#define SRC_RUNTIME_LOCAL_KERNELS_GEMV_H

#include <runtime/local/context/DaphneContext.h>
#include <runtime/local/datastructures/DataObjectFactory.h>
#include <runtime/local/datastructures/DenseMatrix.h>

#include <cblas.h>

// ****************************************************************************
// Struct for partial template specialization
// ****************************************************************************

template<class DTRes, class DTMat, class DTVec>
struct Gemv {
    static void apply(DTRes *& res, const DTMat * mat, const DTVec * vec, DCTX(ctx)) = delete;
};

// ****************************************************************************
// Convenience function
// ****************************************************************************

template<class DTRes, class DTMat, class DTVec>
void gemv(DTRes *& res, const DTMat * mat, const DTVec * vec, DCTX(ctx)) {
    Gemv<DTRes, DTMat, DTVec>::apply(res, mat, vec, ctx);
}

// ****************************************************************************
// (Partial) template specializations for different data/value types
// ****************************************************************************

// ----------------------------------------------------------------------------
// DenseMatrix <- DenseMatrix
// ----------------------------------------------------------------------------

template<>
struct Gemv<DenseMatrix<double>, DenseMatrix<double>, DenseMatrix<double>> {
    static void apply(DenseMatrix<double> *& res, const DenseMatrix<double> * mat, const DenseMatrix<double> * vec, DCTX(ctx)) {
        const size_t numRows = mat->getNumRows();
        const size_t numCols = mat->getNumCols();

        if(res == nullptr)
            res = DataObjectFactory::create<DenseMatrix<double>>(numCols, 1, false);

        cblas_dgemv(CblasRowMajor,
            CblasTrans,
            numRows,
            numCols,
            1.0,
            mat->getValues(),
            mat->getRowSkip(),
            vec->getValues(),
            vec->getRowSkip(),
            0.0,
            res->getValues(),
            res->getRowSkip()
        );
    }
};

template<>
struct Gemv<DenseMatrix<float>, DenseMatrix<float>, DenseMatrix<float>> {
    static void apply(DenseMatrix<float> *& res, const DenseMatrix<float> * mat, const DenseMatrix<float> * vec, DCTX(ctx)) {
        const size_t numRows = mat->getNumRows();
        const size_t numCols = mat->getNumCols();

        if(res == nullptr)
            res = DataObjectFactory::create<DenseMatrix<float>>(numCols, 1, false);

        cblas_sgemv(CblasRowMajor,
            CblasTrans,
            numRows,
            numCols,
            1.0,
            mat->getValues(),
            mat->getRowSkip(),
            vec->getValues(),
            vec->getRowSkip(),
            0.0,
            res->getValues(),
            res->getRowSkip()
        );
    }
};

#endif //SRC_RUNTIME_LOCAL_KERNELS_GEMV_H

max_stars_repo_name: daphne-eu/daphne

max_issues_repo_name: daphne-eu/daphne

max_forks_repo_name: daphne-eu/daphne