text:
#pragma once

#include "commons.h"
#include "Instancing.h"

#include <renderer/VertexSpecification.h>

#include <gsl/gsl>

namespace ad {

namespace vertex3D
{

struct Data
{
    Vec<3, GLfloat> mPosition;
    Vec<3, GLfloat> mNormal{0.f, 0.f, 1.f};
};

constexpr AttributeDescriptionList gDescription = {
    { 0, 3, offsetof(Data, mPosition), MappedGL<GLfloat>::enumerator},
    { 1, 3, offsetof(Data, mNormal), MappedGL<GLfloat>::enumerator},
};

inline Data operator*(const Data & aLhs, const Matrix<4> & aRhs)
{
    Vec<4> position{aLhs.mPosition.x(), aLhs.mPosition.y(), aLhs.mPosition.z(), 1.0f};
    Vec<4> normal{aLhs.mNormal.x(), aLhs.mNormal.y(), aLhs.mNormal.z(), 0.0f};

    position *= aRhs;
    normal *= aRhs;

    return {
        {position.x(), position.y(), position.z()},
        {normal.x(), normal.y(), normal.z()},
    };
}

inline const gsl::span<Data> & operator*=(const gsl::span<Data> & aLhs, const Matrix<4> & aRhs)
{
    std::transform(aLhs.begin(), aLhs.end(), aLhs.begin(),
                   [&aRhs](const Data & data){ return data * aRhs; });
    return aLhs;
}

} // namespace vertex3D

template <>
class BufferData_trait<vertex3D::Data>
{
public:
    static constexpr AttributeDescriptionList description_value = vertex3D::gDescription;
};

namespace ship3D {

static const std::array<vertex3D::Data, 4> gVertices = {
    vertex3D::Data{
        {-1.25f, -1.5f, 0.0f},
    },
    vertex3D::Data{
        { 3.75f,  0.0f, 0.0f},
    },
    vertex3D::Data{
        {-1.25f,  1.5f, 0.0f},
    },
    vertex3D::Data{
        { 0.0f,  0.0f,  0.0f},
    },
};

} // namespace ship3D

namespace cube {

static constexpr float cpi4 = 0.70710678f;

static const std::array<vertex3D::Data, 8*3> gVertices = {
    // Frontface
    vertex3D::Data{
        {-1.0f, -1.0f, 1.0f},
        {-cpi4,  0.0f, cpi4},
    },
    vertex3D::Data{
        {-1.0f,  1.0f, 1.0f},
        {-cpi4,  0.0f, cpi4},
    },
    vertex3D::Data{
        {-1.0f,  1.0f, 1.0f},
        { 0.0f,  cpi4, cpi4},
    },
    vertex3D::Data{
        { 1.0f,  1.0f, 1.0f},
        { 0.0f,  cpi4, cpi4},
    },
    vertex3D::Data{
        { 1.0f,  1.0f, 1.0f},
        { cpi4,  0.0f, cpi4},
    },
    vertex3D::Data{
        { 1.0f,  -1.0f, 1.0f},
        { cpi4,  0.0f,  cpi4},
    },
    vertex3D::Data{
        { 1.0f, -1.0f, 1.0f},
        { 0.0f, -cpi4, cpi4},
    },
    vertex3D::Data{
        {-1.0f, -1.0f, 1.0f},
        { 0.0f, -cpi4, cpi4},
    },

    // Backface
    vertex3D::Data{
        {-1.0f, -1.0f, -1.0f},
          {-cpi4,  0.0f, -cpi4},
    },
    vertex3D::Data{
        {-1.0f,  1.0f, -1.0f},
        {-cpi4,  0.0f, -cpi4},
    },
    vertex3D::Data{
        {-1.0f,  1.0f, -1.0f},
        { 0.0f,  cpi4, -cpi4},
    },
    vertex3D::Data{
        { 1.0f,  1.0f, -1.0f},
        { 0.0f,  cpi4, -cpi4},
    },
    vertex3D::Data{
        { 1.0f,  1.0f, -1.0f},
        { cpi4,  0.0f, -cpi4},
    },
    vertex3D::Data{
        { 1.0f,  -1.0f, -1.0f},
        { cpi4,   0.0f, -cpi4},
    },
    vertex3D::Data{
        { 1.0f, -1.0f, -1.0f},
        { 0.0f, -cpi4, -cpi4},
    },
    vertex3D::Data{
        {-1.0f, -1.0f, -1.0f},
        { 0.0f, -cpi4, -cpi4},
    },

    // Connect
    vertex3D::Data{
        {-1.0f, -1.0f, -1.0f},
        {-cpi4, -cpi4,  0.0f},
    },
    vertex3D::Data{
        {-1.0f, -1.0f, +1.0f},
        {-cpi4, -cpi4,  0.0f},
    },
    vertex3D::Data{
        {-1.0f,  1.0f, -1.0f},
        {-cpi4,  cpi4,  0.0f},
    },
    vertex3D::Data{
        {-1.0f,  1.0f, +1.0f},
        {-cpi4,  cpi4,  0.0f},
    },
    vertex3D::Data{
        { 1.0f,  1.0f, -1.0f},
        { cpi4,  cpi4,  0.0f},
    },
    vertex3D::Data{
        { 1.0f,  1.0f, +1.0f},
        { cpi4,  cpi4,  0.0f},
    },
    vertex3D::Data{
        { 1.0f,  -1.0f, -1.0f},
        { cpi4,  -cpi4,  0.0f},
    },
    vertex3D::Data{
        { 1.0f,  -1.0f, +1.0f},
        { cpi4,  -cpi4,  0.0f},
    },
};

} // namespace cube

namespace circle3D {

enum Type
{
    Circle,
    Disc,
};

template <GLsizei N_vertices, Type N_type = Circle>
std::vector<vertex3D::Data> makeVertices(GLfloat aRadius)
{
    std::vector<vertex3D::Data> result;
    if (N_type == Disc)
    {
        result.push_back(vertex3D::Data{
            {0.f, 0.f, 0.f},
        });
    }
    //result.reserve(N_vertices+2);
    for(int i=0; i<=N_vertices; ++i)
    {
        result.emplace_back(vertex3D::Data{{
                                                cos(2*pi<float>*i/N_vertices) * aRadius,
                                                sin(2*pi<float>*i/N_vertices) * aRadius,
                                                0.f,
                                           },
                                           {
                                               0.f,
                                               0.f,
                                               N_type == Disc ? -1.f : 1.0f,
                                           }});
    }
    return result;
}

} // namespace circle3D

namespace rectangle3D {

inline std::vector<vertex3D::Data> makeVertices(Size<2> aDimensions)
{
    return {
        { {+aDimensions.width()/2.f, +aDimensions.height()/2.f, 0.f} },
        { {-aDimensions.width()/2.f, +aDimensions.height()/2.f, 0.f} },
        { {-aDimensions.width()/2.f, -aDimensions.height()/2.f, 0.f} },
        { {+aDimensions.width()/2.f, -aDimensions.height()/2.f, 0.f} },
    };
}

} // namespace rectangle3D

namespace instance3D {

struct Data
{
    Vec<2, GLfloat> mPosition;
    Matrix<4, GLfloat> mOrientationScale{Matrix<4, GLfloat>::Identity()};
};

constexpr AttributeDescriptionList gDescription = {
    {2, 2, offsetof(Data, mPosition),      MappedGL<GLfloat>::enumerator},
    {3, 4, offsetof(Data, mOrientationScale),   MappedGL<GLfloat>::enumerator},
    {4, 4, offsetof(Data, mOrientationScale) +  4*sizeof(GLfloat),   MappedGL<GLfloat>::enumerator},
    {5, 4, offsetof(Data, mOrientationScale) +  8*sizeof(GLfloat),   MappedGL<GLfloat>::enumerator},
    {6, 4, offsetof(Data, mOrientationScale) + 12*sizeof(GLfloat),   MappedGL<GLfloat>::enumerator},
};

} // namespace instance3D

template <>
class BufferData_trait<instance3D::Data>
{
public:
    static constexpr AttributeDescriptionList description_value = instance3D::gDescription;
};

} // namespace ad

ext: h

max_stars_repo_name: FranzPoize/shmurp

max_issues_repo_name: FranzPoize/shmurp

max_forks_repo_name: FranzPoize/shmurp


text:
/* This file is auto-generated from countpairs_s_mu_mocks_impl.c.src */
#ifndef DOUBLE_PREC
#define DOUBLE_PREC
#endif
// # -*- mode: c -*-
/* File: countpairs_s_mu_mocks_impl.c.src */
/*
  This file is a part of the Corrfunc package
  Copyright (C) 2015-- Manodeep Sinha (manodeep@gmail.com)
  License: MIT LICENSE. See LICENSE file under the top-level
  directory at https://github.com/manodeep/Corrfunc/
*/

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <signal.h>
#include <unistd.h>
#include <gsl/gsl_interp.h>

#include "countpairs_s_mu_mocks_impl_double.h"
#include "countpairs_s_mu_mocks_kernels_double.c"
#include "cellarray_mocks_double.h"
#include "gridlink_mocks_impl_double.h"

#include "defs.h"
#include "utils.h"
#include "cosmology_params.h"
#include "set_cosmo_dist.h"
#include "cpu_features.h"
#include "progressbar.h"
#include "proj_functions_double.h"

#if defined(_OPENMP)
#include <omp.h>
#endif

int interrupt_status_DDsmu_mocks_double=EXIT_SUCCESS;

void interrupt_handler_countpairs_s_mu_mocks_double(int signo)
{
    fprintf(stderr,"Received signal = `%s' (signo = %d). Aborting \n",strsignal(signo), signo);
    interrupt_status_DDsmu_mocks_double = EXIT_FAILURE;
}

int check_ra_dec_cz_s_mu_double(const int64_t N, double *phi, double *theta, double *cz)
{

    if(N==0) {
        return EXIT_SUCCESS;
    }
    if(phi == NULL || theta == NULL || cz == NULL) {
        fprintf(stderr,"Input arrays can not be NULL. Have RA = %p DEC = %p cz = %p\n", phi, theta, cz);
        return EXIT_FAILURE;
    }

    int fix_cz  = 0;
    int fix_ra  = 0;
    int fix_dec = 0;

    const double max_cz_threshold = 10.0;//if I find that max cz is smaller than this threshold, then I will assume z has been supplied rather than cz
    double max_cz = 0.0;
    //Check input cz -> ensure that cz contains cz and not z
    for(int64_t i=0;i<N;i++) {
        if(cz[i] > max_cz) max_cz = cz[i];
        if(phi[i] < 0.0) {
            fix_ra = 1;
        }
        if(theta[i] > 90.0) {
            fix_dec = 1;
        }
        if(theta[i] > 180) {
            fprintf(stderr,"theta[%"PRId64"] = %"REAL_FORMAT"should be less than 180 deg\n", i, theta[i]);
            return EXIT_FAILURE;
        }
    }
    if(max_cz < max_cz_threshold) fix_cz = 1;

    //Only run the loop if something needs to be fixed
    if(fix_cz==1 || fix_ra == 1 || fix_dec == 1) {
        if(fix_ra == 1) {
            fprintf(stderr,"%s> Out of range values found for ra. Expected ra to be in the range [0.0,360.0]. Found ra values in [-180,180] -- fixing that\n", __FUNCTION__);
        }
        if(fix_dec == 1) {
            fprintf(stderr,"%s> Out of range values found for dec. Expected dec to be in the range [-90.0,90.0]. Found dec values in [0,180] -- fixing that\n", __FUNCTION__);
        }
        if(fix_cz == 1)  {
            fprintf(stderr,"%s> Out of range values found for cz. Expected input to be `cz' but found `z' instead. max_cz (found in input) = %"REAL_FORMAT" threshold "
                    "= %"REAL_FORMAT"\n",__FUNCTION__,max_cz,max_cz_threshold);
        }

        for(int64_t i=0;i<N;i++) {
            if(fix_ra==1) {
                phi[i] += (double) 180.0;
            }
            if(fix_dec==1) {
                theta[i] -= (double) 90.0;
            }
            if(fix_cz == 1) {
                cz[i] *= (double) SPEED_OF_LIGHT;//input was z -> convert to cz
            }
        }
    }

    return EXIT_SUCCESS;
}

countpairs_mocks_func_ptr_double countpairs_s_mu_mocks_driver_double(const struct config_options *options)
{

    static countpairs_mocks_func_ptr_double function = NULL;
    static isa old_isa=-1;
    if(old_isa == options->instruction_set) {
        return function;
    }

    /* Array of function pointers */
    countpairs_mocks_func_ptr_double allfunctions[] = {
#ifdef __AVX__
        countpairs_s_mu_mocks_avx_intrinsics_double,
#endif
#ifdef __SSE4_2__
      countpairs_s_mu_mocks_sse_intrinsics_double,
#endif
      countpairs_s_mu_mocks_fallback_double
    };

    const int num_functions = sizeof(allfunctions)/sizeof(void *);
    const int fallback_offset = num_functions - 1;
#if defined(__AVX__) || defined __SSE4_2__
    const int highest_isa = instrset_detect();
#endif
    int curr_offset = 0;

    /* Now check if AVX is supported by the CPU */
    int avx_offset = fallback_offset;
#ifdef __AVX__
    avx_offset = highest_isa >= 7 ? curr_offset:fallback_offset;
    curr_offset++;
#endif

    /* Is the SSE function supported at runtime and enabled at compile-time?*/
    int sse_offset = fallback_offset;
#ifdef __SSE4_2__
    sse_offset = highest_isa >= 6 ? curr_offset:fallback_offset;
    curr_offset++;
#endif
    if( curr_offset != fallback_offset) {
      fprintf(stderr,"ERROR: Bug in code (current offset = %d *should equal* fallback function offset = %d)\n",
              curr_offset, fallback_offset);
      return NULL;
    }

    int function_dispatch=0;
    /* Check that cpu supports feature */
    if(options->instruction_set >= 0) {
        switch(options->instruction_set) {
        case(AVX512F):
        case(AVX2):
        case(AVX):function_dispatch=avx_offset;break;
        case(SSE42): function_dispatch=sse_offset;break;
        default:function_dispatch=fallback_offset;break;
        }
    }

    if(function_dispatch >= num_functions) {
      fprintf(stderr,"In %s> ERROR: Could not resolve the correct function.\n Function index = %d must lie between [0, %d)\n",
              __FUNCTION__, function_dispatch, num_functions);
      return NULL;
    }
    function = allfunctions[function_dispatch];
    old_isa = options->instruction_set;

    if(options->verbose){
        // This must be first (AVX/SSE may be aliased to fallback)
        if(function_dispatch == fallback_offset){
            fprintf(stderr,"Using fallback kernel\n");
        } else if(function_dispatch == avx_offset){
            fprintf(stderr,"Using AVX kernel\n");
        } else if(function_dispatch == sse_offset){
            fprintf(stderr,"Using SSE kernel\n");
        } else {
            printf("Unknown kernel!\n");
        }
    }

    return function;
}

int countpairs_mocks_s_mu_double(const int64_t ND1, double *ra1, double *dec1, double *czD1,
                                 const int64_t ND2, double *ra2, double *dec2, double *czD2,
                                 const int numthreads,
                                 const int autocorr,
                                 const char *sbinfile,
                                 const double max_mu,
                                 const int nmu_bins,
                                 const int cosmology,
                                 results_countpairs_mocks_s_mu *results,
                                 struct config_options *options, struct extra_options *extra)
{

    if(options->float_type != sizeof(double)) {
        fprintf(stderr,"ERROR: In %s> Can only handle arrays of size=%zu. Got an array of size = %zu\n",
                __FUNCTION__, sizeof(double), options->float_type);
        return EXIT_FAILURE;
    }

    // If no extra options were passed, create dummy options
    // This allows us to pass arguments like "extra->weights0" below;
    // they'll just be NULLs, which is the correct behavior
    struct extra_options dummy_extra;
    if(extra == NULL){
      weight_method_t dummy_method = NONE;
      dummy_extra = get_extra_options(dummy_method);
      extra = &dummy_extra;
    }

    int need_weightavg = extra->weight_method != NONE;

    options->sort_on_z = 1;
    struct timeval t0;
    if(options->c_api_timer) {
        gettimeofday(&t0, NULL);
    }
    if(options->fast_divide_and_NR_steps >= MAX_FAST_DIVIDE_NR_STEPS) {
        fprintf(stderr, ANSI_COLOR_MAGENTA"Warning: The number of requested Newton-Raphson steps = %u is larger than max. allowed steps = %u."
                " Switching to a standard divide"ANSI_COLOR_RESET"\n",
                options->fast_divide_and_NR_steps, MAX_FAST_DIVIDE_NR_STEPS);
        options->fast_divide_and_NR_steps = 0;
    }

    //Check inputs
    if(ND1 == 0 || (autocorr == 0 && ND2 == 0)) {
        return EXIT_SUCCESS;
    }

    //Check inputs
    int status1 = check_ra_dec_cz_s_mu_double(ND1, ra1, dec1, czD1);
    if(status1 != EXIT_SUCCESS) {
        return status1;
    }
    if(autocorr==0) {
        int status2 = check_ra_dec_cz_s_mu_double(ND2, ra2, dec2, czD2);
        if(status2 != EXIT_SUCCESS) {
            return status2;
        }
    }

#if defined(_OPENMP)
    omp_set_num_threads(numthreads);
#else
    (void) numthreads;
#endif

    if(options->max_cells_per_dim == 0) {
        fprintf(stderr,"Warning: Max. cells per dimension is set to 0 - resetting to `NLATMAX' = %d\n", NLATMAX);
        options->max_cells_per_dim = NLATMAX;
    }
    for(int i=0;i<3;i++) {
        if(options->bin_refine_factors[i] < 1) {
            fprintf(stderr,"Warning: bin refine factor along axis = %d *must* be >=1. Instead found bin refine factor =%d\n",
                    i, options->bin_refine_factors[i]);
            reset_bin_refine_factors(options);
            break;/* all factors have been reset -> no point continuing with the loop */
        }
    }

    /* setup interrupt handler -> mostly useful during the python execution.
       Let's Ctrl-C abort the extension  */
    SETUP_INTERRUPT_HANDLERS(interrupt_handler_countpairs_s_mu_mocks_double);

    //Try to initialize cosmology - code will exit if comoslogy is not implemented.
    //Putting in a different scope so I can call the variable status
    {
        int status = init_cosmology(cosmology);
        if(status != EXIT_SUCCESS) {
            return status;
        }
    }

    /***********************
     *initializing the  bins
     ************************/
    double *supp;
    int nsbin;
    double smin,smax;
    setup_bins(sbinfile,&smin,&smax,&nsbin,&supp);
    if( ! (smin > 0.0 && smax > 0.0 && smin < smax && nsbin > 0)) {
        fprintf(stderr,"Error: Could not setup with S bins correctly. (smin = %lf, smax = %lf, with nbins = %d). Expected non-zero smin/smax with smax > smin and nbins >=1 \n",
                smin, smax, nsbin);
        return EXIT_FAILURE;
    }

    if(max_mu <= 0.0 || max_mu > 1.0) {
        fprintf(stderr,"Error: max_mu (max. value for the cosine of the angle with line of sight) must be greater than 0 and at most 1).\n"
                "The passed value is max_mu = %lf. Please change it to be > 0 and <= 1.0\n", max_mu);
        return EXIT_FAILURE;
    }

    if(nmu_bins < 1 ) {
        fprintf(stderr,"Error: Number of mu bins = %d must be at least 1\n", nmu_bins);
        return EXIT_FAILURE;
    }

    //Change cz into co-moving distance
    double *D1 = NULL, *D2 = NULL;
    if(options->is_comoving_dist == 0) {
        D1 = my_malloc(sizeof(*D1),ND1);
        D2 = autocorr == 0 ? my_malloc(sizeof(*D2),ND2):D1;
    } else {
        D1 = czD1;
        D2 = autocorr == 0 ? czD2:czD1;
    }

    if(D1 == NULL || D2 == NULL) {
        free(D1);free(D2);
        return EXIT_FAILURE;
    }

    if(options->is_comoving_dist == 0) {
        //Setup variables to do the cz->comoving distance
        double czmax = 0.0;
        const double inv_speed_of_light = 1.0/SPEED_OF_LIGHT;
        get_max_double(ND1, czD1, &czmax);
        if(autocorr == 0) {
            get_max_double(ND2, czD2, &czmax);
        }
        const double zmax = czmax * inv_speed_of_light + 0.01;

        const int workspace_size = 10000;
        double *interp_redshift  = my_calloc(sizeof(*interp_redshift), workspace_size);//the interpolation is done in 'z' and not in 'cz'
        double *interp_comoving_dist = my_calloc(sizeof(*interp_comoving_dist),workspace_size);
        int Nzdc = set_cosmo_dist(zmax, workspace_size, interp_redshift, interp_comoving_dist, cosmology);
        if(Nzdc < 0) {
            free(interp_redshift);free(interp_comoving_dist);
            return EXIT_FAILURE;
        }

        gsl_interp *interpolation;
        gsl_interp_accel *accelerator;
        accelerator =  gsl_interp_accel_alloc();
        interpolation = gsl_interp_alloc (gsl_interp_linear,Nzdc);
        gsl_interp_init(interpolation, interp_redshift, interp_comoving_dist, Nzdc);
        for(int64_t i=0;i<ND1;i++) {
            D1[i] = gsl_interp_eval(interpolation, interp_redshift, interp_comoving_dist, czD1[i]*inv_speed_of_light, accelerator);
        }

        if(autocorr==0) {
            for(int64_t i=0;i<ND2;i++) {
                D2[i] = gsl_interp_eval(interpolation, interp_redshift, interp_comoving_dist, czD2[i]*inv_speed_of_light, accelerator);
            }
        }
        free(interp_redshift);free(interp_comoving_dist);
        gsl_interp_free(interpolation);
        gsl_interp_accel_free(accelerator);
    }

    double *X1 = my_malloc(sizeof(*X1), ND1);
    double *Y1 = my_malloc(sizeof(*Y1), ND1);
    double *Z1 = my_malloc(sizeof(*Z1), ND1);
    if(X1 == NULL || Y1 == NULL || Z1 == NULL) {
        free(X1);free(Y1);free(Z1);
        return EXIT_FAILURE;
    }
    for(int64_t i=0;i<ND1;i++) {
        X1[i] = D1[i]*COSD(dec1[i])*COSD(ra1[i]);
        Y1[i] = D1[i]*COSD(dec1[i])*SIND(ra1[i]);
        Z1[i] = D1[i]*SIND(dec1[i]);
    }

    double *X2,*Y2,*Z2;
    if(autocorr==0) {
        X2 = my_malloc(sizeof(*X2), ND2);
        Y2 = my_malloc(sizeof(*Y2), ND2);
        Z2 = my_malloc(sizeof(*Z2), ND2);
        for(int64_t i=0;i<ND2;i++) {
            X2[i] = D2[i]*COSD(dec2[i])*COSD(ra2[i]);
            Y2[i] = D2[i]*COSD(dec2[i])*SIND(ra2[i]);
            Z2[i] = D2[i]*SIND(dec2[i]);
        }
    } else {
        X2 = X1;
        Y2 = Y1;
        Z2 = Z1;
    }

    double supp_sqr[nsbin];
    for(int i=0; i < nsbin;i++) {
        supp_sqr[i] = supp[i]*supp[i];
    }
    const double mu_max = (double) max_mu;

    double xmin=1e10,ymin=1e10,zmin=1e10;
    double xmax=-1e10,ymax=-1e10,zmax=-1e10;
    get_max_min_data_double(ND1, X1, Y1, Z1, &xmin, &ymin, &zmin, &xmax, &ymax, &zmax);

    if(autocorr==0) {
        get_max_min_data_double(ND2, X2, Y2, Z2, &xmin, &ymin, &zmin, &xmax, &ymax, &zmax);
    }

    const double xdiff = xmax-xmin;
    const double ydiff = ymax-ymin;
    const double zdiff = zmax-zmin;
    if(get_bin_refine_scheme(options) == BINNING_DFL) {
        if(smax < 0.05*xdiff) {
            options->bin_refine_factors[0] = 1;
      }
        if(smax < 0.05*ydiff) {
            options->bin_refine_factors[1] = 1;
        }
        if(smax < 0.05*zdiff) {
            options->bin_refine_factors[2] = 1;
        }
    }

    /*---Create 3-D lattice--------------------------------------*/
    int nmesh_x=0,nmesh_y=0,nmesh_z=0;
    cellarray_mocks_index_particles_double *lattice1 = gridlink_mocks_index_particles_double(ND1, X1, Y1, Z1, D1, &(extra->weights0),
                                                                                xmin, xmax, ymin, ymax, zmin, zmax,
                                                                                smax, smax, smax,
                                                                                options->bin_refine_factors[0],
                                                                                options->bin_refine_factors[1],
                                                                                options->bin_refine_factors[2],
                                                                                &nmesh_x, &nmesh_y, &nmesh_z,
                                                                                options);
    if(lattice1 == NULL) {
        return EXIT_FAILURE;
    }

    /* If there too few cells (BOOST_CELL_THRESH is ~10), and the number of cells can be increased, then boost bin refine factor by ~1*/
    const double avg_np = ((double)ND1)/(nmesh_x*nmesh_y*nmesh_z);
    const int8_t max_nmesh = fmax(nmesh_x, fmax(nmesh_y, nmesh_z));
    if((max_nmesh <= BOOST_CELL_THRESH || avg_np >= BOOST_NUMPART_THRESH)
        && max_nmesh < options->max_cells_per_dim) {
      fprintf(stderr,"%s> gridlink seems inefficient. nmesh = (%d, %d, %d); avg_np = %.3g. ", __FUNCTION__, nmesh_x, nmesh_y, nmesh_z, avg_np);
      if(get_bin_refine_scheme(options) == BINNING_DFL) {
            fprintf(stderr,"Boosting bin refine factor - should lead to better performance\n");
            // Only boost the first two dimensions.  Prevents excessive refinement.
            for(int i=0;i<2;i++) {
              options->bin_refine_factors[i] += BOOST_BIN_REF;
            }

            free_cellarray_mocks_index_particles_double(lattice1, nmesh_x * (int64_t) nmesh_y * nmesh_z);
            lattice1 = gridlink_mocks_index_particles_double(ND1, X1, Y1, Z1, D1, &(extra->weights0),
                                                             xmin, xmax, ymin, ymax, zmin, zmax,
                                                             smax, smax, smax,
                                                             options->bin_refine_factors[0],
                                                             options->bin_refine_factors[1],
                                                             options->bin_refine_factors[2],
                                                             &nmesh_x, &nmesh_y, &nmesh_z,
                                                             options);
            if(lattice1 == NULL) {
                return EXIT_FAILURE;
            }
        } else {
            fprintf(stderr,"Boosting bin refine factor could have helped. However, since custom bin refine factors "
                  "= (%d, %d, %d) are being used - continuing with inefficient mesh\n", options->bin_refine_factors[0],
                  options->bin_refine_factors[1], options->bin_refine_factors[2]);

        }
    }

    cellarray_mocks_index_particles_double *lattice2 = NULL;
    if(autocorr==0) {
        int ngrid2_x=0,ngrid2_y=0,ngrid2_z=0;
        lattice2 = gridlink_mocks_index_particles_double(ND2, X2, Y2, Z2, D2, &(extra->weights1),
                                                         xmin, xmax,
                                                         ymin, ymax,
                                                         zmin, zmax,
                                                         smax, smax, smax,
                                                         options->bin_refine_factors[0],
                                                         options->bin_refine_factors[1],
                                                         options->bin_refine_factors[2],
                                                         &ngrid2_x, &ngrid2_y, &ngrid2_z, options);
        if(lattice2 == NULL) {
            return EXIT_FAILURE;
        }
        if( ! (nmesh_x == ngrid2_x && nmesh_y == ngrid2_y && nmesh_z == ngrid2_z) ) {
            fprintf(stderr,"Error: The two sets of 3-D lattices do not have identical bins. First has dims (%d, %d, %d) while second has (%d, %d, %d)\n",
                    nmesh_x, nmesh_y, nmesh_z, ngrid2_x, ngrid2_y, ngrid2_z);
            return EXIT_FAILURE;
        }
    } else {
        lattice2 = lattice1;
    }
    free(X1);free(Y1);free(Z1);
    if(autocorr == 0) {
        free(X2);free(Y2);free(Z2);
    }

    if(options->is_comoving_dist == 0) {
        free(D1);
        if(autocorr == 0) {
            free(D2);
        }
    }

    const int64_t totncells = (int64_t) nmesh_x * (int64_t) nmesh_y * (int64_t) nmesh_z;
    {
        int status = assign_ngb_cells_mocks_index_particles_double(lattice1, lattice2, totncells,
                                                                   options->bin_refine_factors[0], options->bin_refine_factors[1], options->bin_refine_factors[2],
                                                                   nmesh_x, nmesh_y, nmesh_z,
                                                                   autocorr);
        if(status != EXIT_SUCCESS) {
            free_cellarray_mocks_index_particles_double(lattice1, totncells);
            if(autocorr == 0) {
                free_cellarray_mocks_index_particles_double(lattice2, totncells);
            }
            free(supp);
            return EXIT_FAILURE;
        }
    }
    /*---Gridlink-variables----------------*/
    const int totnbins = (nmu_bins+1)*(nsbin+1);
    const int nprojbins = nsbin-1;
#if defined(_OPENMP)
    uint64_t **all_npairs = (uint64_t **) matrix_calloc(sizeof(uint64_t), numthreads, totnbins);
    double **all_savg = NULL;
    if(options->need_avg_sep){
        all_savg = (double **) matrix_calloc(sizeof(double),numthreads,totnbins);
    }
    double **all_weightavg = NULL;
    if(need_weightavg) {
      all_weightavg = (double **) matrix_calloc(sizeof(double),numthreads,totnbins);
    }
    double **all_projpairs = (double **) matrix_calloc(sizeof(double),numthreads,nprojbins);
    double **all_projpairs_tensor = (double **) matrix_calloc(sizeof(double),numthreads,nprojbins*nprojbins);

#else //USE_OMP
    uint64_t npairs[totnbins];
    double savg[totnbins], weightavg[totnbins], projpairs[nprojbins];
    double projpairs_tensor[nprojbins*nprojbins];

    for(int i=0; i <totnbins;i++) {
        npairs[i] = 0;
        if(options->need_avg_sep) {
            savg[i] = ZERO;
        }
        if(need_weightavg) {
            weightavg[i] = ZERO;
        }
    }
    for(int i=0;i<nprojbins;i++) {
        projpairs[i] = ZERO;
        for(int j=0;j<nprojbins;j++) {
            projpairs_tensor[i*nprojbins+j] = ZERO;
        }
    }
#endif //USE_OMP

    /* runtime dispatch - get the function pointer */
    countpairs_mocks_func_ptr_double countpairs_s_mu_mocks_function_double = countpairs_s_mu_mocks_driver_double(options);
    if(countpairs_s_mu_mocks_function_double == NULL) {
        return EXIT_FAILURE;
    }

    int interrupted=0,numdone=0, abort_status=EXIT_SUCCESS;
    if(options->verbose) {
        init_my_progressbar(totncells,&interrupted);
    }

#if defined(_OPENMP)
#pragma omp parallel shared(numdone, abort_status, interrupt_status_DDsmu_mocks_double)
    {
        const int tid = omp_get_thread_num();
        uint64_t npairs[totnbins];
        double savg[totnbins], weightavg[totnbins], projpairs[nprojbins];
        double projpairs_tensor[nprojbins*nprojbins];

        for(int i=0;i<totnbins;i++) {
            npairs[i] = 0;
            if(options->need_avg_sep) {
                savg[i] = ZERO;
            }
            if(need_weightavg) {
                weightavg[i] = ZERO;
            }
        }
        for(int i=0;i<nprojbins;i++) {
            projpairs[i] = ZERO;
            for(int j=0;j<nprojbins;j++) {
                projpairs_tensor[i*nprojbins+j] = ZERO;
            }
        }

#pragma omp for  schedule(dynamic)
#endif//USE_OMP

        /*---Loop-over-Data1-particles--------------------*/
        for(int64_t index1=0;index1<totncells;index1++) {

#if defined(_OPENMP)
#pragma omp flush (abort_status, interrupt_status_DDsmu_mocks_double)
#endif
            if(abort_status == EXIT_SUCCESS && interrupt_status_DDsmu_mocks_double == EXIT_SUCCESS) {
                //omp cancel was introduced in omp 4.0 - so this is my way of checking if loop needs to be cancelled
                /* If the verbose option is not enabled, avoid outputting anything unnecessary*/
                if(options->verbose) {
#if defined(_OPENMP)
                    if (omp_get_thread_num() == 0)
#endif
                        my_progressbar(numdone,&interrupted);

#if defined(_OPENMP)
#pragma omp atomic
#endif
                    numdone++;
                }

                const cellarray_mocks_index_particles_double *first  = &(lattice1[index1]);
                if(first->nelements == 0) {
                    continue;
                }
                double *x1 = first->x;
                double *y1 = first->y;
                double *z1 = first->z;
                double *d1 = first->cz;
                const weight_struct_double *weights1 = &(first->weights);
                const int64_t N1 = first->nelements;

                if(autocorr == 1) {
                    int same_cell = 1;
                    double *this_savg = options->need_avg_sep ? &(savg[0]):NULL;
                    double *this_weightavg = need_weightavg ? weightavg:NULL;
                    const int status = countpairs_s_mu_mocks_function_double(N1, x1, y1, z1, d1, weights1,
                                                                             N1, x1, y1, z1, d1, weights1,
                                                                             same_cell,
                                                                             options->fast_divide_and_NR_steps,
                                                                             smax, smin, nsbin,
                                                                             nmu_bins, supp_sqr, mu_max,
                                                                             this_savg, npairs, projpairs,
                                                                             projpairs_tensor,
                                                                             this_weightavg, extra->weight_method);
                    /* This actually causes a race condition under OpenMP - but mostly
                       I care that an error occurred - rather than the exact value of
                       the error status */
                    abort_status |= status;
                }

                for(int64_t ngb=0;ngb<first->num_ngb;ngb++){
                    const cellarray_mocks_index_particles_double *second = first->ngb_cells[ngb];
                    if(second->nelements == 0) {
                        continue;
                    }
                    const int same_cell = 0;
                    double *x2 = second->x;
                    double *y2 = second->y;
                    double *z2 = second->z;
                    double *d2 = second->cz;
                    const weight_struct_double *weights2 = &(second->weights);
                    const int64_t N2 = second->nelements;
                    double *this_savg = options->need_avg_sep ? &(savg[0]):NULL;
                    double *this_weightavg = need_weightavg ? weightavg:NULL;
                    const int status = countpairs_s_mu_mocks_function_double(N1, x1, y1, z1, d1, weights1,
                                                                             N2, x2, y2, z2, d2, weights2,
                                                                             same_cell,
                                                                             options->fast_divide_and_NR_steps,
                                                                             smax, smin, nsbin,
                                                                             nmu_bins, supp_sqr, mu_max,
                                                                             this_savg, npairs, projpairs,
                                                                             projpairs_tensor,
                                                                             this_weightavg, extra->weight_method);
                    /* This actually causes a race condition under OpenMP - but mostly
                       I care that an error occurred - rather than the exact value of
                       the error status */
                    abort_status |= status;
                }//loop over ngb cells
            }//abort_status check
        }//i loop over ND1 particles
#if defined(_OPENMP)
        for(int i=0;i<totnbins;i++) {
            all_npairs[tid][i] = npairs[i];
            if(options->need_avg_sep) {
                all_savg[tid][i] = savg[i];
            }
            if(need_weightavg) {
                all_weightavg[tid][i] = weightavg[i];
            }
        }
        for (int i=0;i<nprojbins;i++) {
            all_projpairs[tid][i] = projpairs[i];
            for(int j=0;j<nprojbins;j++) {
                all_projpairs_tensor[tid][i*nprojbins+j] = projpairs_tensor[i*nprojbins+j];
            }
        }
    }//close the omp parallel region
#endif//USE_OMP

    free_cellarray_mocks_index_particles_double(lattice1,totncells);
    if(autocorr == 0) {
        free_cellarray_mocks_index_particles_double(lattice2,totncells);
    }

    if(abort_status != EXIT_SUCCESS || interrupt_status_DDsmu_mocks_double != EXIT_SUCCESS) {
        /* Cleanup memory here if aborting */
        free(supp);
#if defined(_OPENMP)
        matrix_free((void **) all_npairs, numthreads);
        if(options->need_avg_sep) {
            matrix_free((void **) all_savg, numthreads);
        }
        if(need_weightavg) {
            matrix_free((void **) all_weightavg, numthreads);
        }
        matrix_free((void **) all_projpairs, numthreads);
        matrix_free((void **) all_projpairs_tensor, numthreads);
#endif
        return EXIT_FAILURE;
    }

    if(options->verbose) {
        finish_myprogressbar(&interrupted);
    }

#if defined(_OPENMP)
    uint64_t npairs[totnbins];
    double savg[totnbins], weightavg[totnbins], projpairs[nprojbins];
    double projpairs_tensor[nprojbins*nprojbins];
    for(int i=0;i<totnbins;i++) {
        npairs[i] = 0;
        if(options->need_avg_sep) {
            savg[i] = ZERO;
        }
        if(need_weightavg) {
            weightavg[i] = ZERO;
        }
    }
    for(int i=0;i<nprojbins;i++) {
        projpairs[i] = ZERO;
        for(int j=0;j<nprojbins;j++) {
            projpairs_tensor[i*nprojbins+j] = ZERO;
        }
    }

    for(int i=0;i<numthreads;i++) {
        for(int j=0;j<totnbins;j++) {
            npairs[j] += all_npairs[i][j];
            if(options->need_avg_sep) {
                savg[j] += all_savg[i][j];
            }
            if(need_weightavg) {
                weightavg[j] += all_weightavg[i][j];
            }
        }
        for(int j=0;j<nprojbins;j++) {
            projpairs[j] += all_projpairs[i][j];
            for(int k=0;k<nprojbins;k++) {
                projpairs_tensor[j*nprojbins+k] += all_projpairs_tensor[i][j*nprojbins+k];
            }
        }
    }
    matrix_free((void **) all_npairs, numthreads);
    if(options->need_avg_sep) {
        matrix_free((void **) all_savg, numthreads);
    }
    if(need_weightavg) {
        matrix_free((void **) all_weightavg, numthreads);
    }
    matrix_free((void **) all_projpairs, numthreads);
    matrix_free((void **) all_projpairs_tensor, numthreads);

#endif //USE_OMP

    //The code does not double count for autocorrelations
    //which means the npairs and savg values need to be doubled;
    if(autocorr == 1) {
        const uint64_t int_fac = 2;
        const double dbl_fac = (double) 2.0;
        for(int i=0;i<totnbins;i++) {
            npairs[i] *= int_fac;
            if(options->need_avg_sep) {
                savg[i] *= dbl_fac;
            }
            if(need_weightavg) {
                weightavg[i] *= dbl_fac;
            }
        }
        //TODO: do i also want to double this? think so
        for(int i=0;i<nprojbins;i++) {
            projpairs[i] *= dbl_fac;
            for(int j=0;j<nprojbins;j++) {
                projpairs_tensor[i*nprojbins+j] *= dbl_fac;
            }
        }

    }

    for(int i=0;i<totnbins;i++) {
        if(npairs[i] > 0) {
            if(options->need_avg_sep) {
                savg[i] /= (double) npairs[i] ;
            }
            if(need_weightavg) {
                weightavg[i] /= (double) npairs[i];
            }
        }
    }
    // don't need proj_pairs here, not averaging

    results->nsbin   = nsbin;
    results->nmu_bins = nmu_bins;
    results->mu_max = max_mu;//NOTE max_mu which is double and not mu_max (which might be float)
    results->mu_min = ZERO;
    results->npairs = my_malloc(sizeof(*(results->npairs)), totnbins);
    results->projpairs = my_malloc(sizeof(*(results->npairs)), nprojbins);
    results->projpairs_tensor = my_malloc(sizeof(*(results->npairs)), nprojbins*nprojbins);
    results->supp   = my_malloc(sizeof(*(results->supp))  , nsbin);
    results->savg  = my_malloc(sizeof(*(results->savg)) , totnbins);
    results->weightavg  = my_calloc(sizeof(double)  , totnbins);
    if(results->npairs == NULL || results->supp == NULL || results->savg == NULL || results->weightavg == NULL || results->projpairs == NULL) {
        free_results_mocks_s_mu(results);
        free(supp);
        return EXIT_FAILURE;
    }

    for(int i=0;i<nsbin;i++) {
        results->supp[i] = supp[i];
        for(int j=0;j<nmu_bins;j++) {
            const int index = i*(nmu_bins+1) + j;
            if( index >= totnbins ) {
                fprintf(stderr, "ERROR: In %s> index = %d must be in range [0, %d)\n", __FUNCTION__, index, totnbins);
                free_results_mocks_s_mu(results);
                free(supp);
                return EXIT_FAILURE;
            }
            results->npairs[index] = npairs[index];
            results->savg[index] = ZERO;
            results->weightavg[index] = ZERO;
            if(options->need_avg_sep) {
                results->savg[index] = savg[index];
            }
            if(need_weightavg) {
                results->weightavg[index] = weightavg[index];
            }
        }
    }
    for(int i=0;i<nprojbins;i++) {
        results->projpairs[i] = projpairs[i];
        for(int j=0;j<nprojbins;j++) {
                results->projpairs_tensor[i*nprojbins+j] = projpairs_tensor[i*nprojbins+j];
            }
    }
    free(supp);

    /* reset interrupt handlers to default */
    RESET_INTERRUPT_HANDLERS();
    reset_bin_refine_factors(options);

    if(options->c_api_timer) {
        struct timeval t1;
        gettimeofday(&t1, NULL);
        options->c_api_time = ADD_DIFF_TIME(t0, t1);
    }

    return EXIT_SUCCESS;
}

ext: c

max_stars_repo_name: kstoreyf/foo

max_issues_repo_name: kstoreyf/foo

max_forks_repo_name: kstoreyf/foo


text:
#ifndef __GSL_BLOCK_H__
#define __GSL_BLOCK_H__

#include <gsl/gsl_block_complex_long_double.h>
#include <gsl/gsl_block_complex_double.h>
#include <gsl/gsl_block_complex_float.h>

#include <gsl/gsl_block_long_double.h>
#include <gsl/gsl_block_double.h>
#include <gsl/gsl_block_float.h>

#include <gsl/gsl_block_ulong.h>
#include <gsl/gsl_block_long.h>

#include <gsl/gsl_block_uint.h>
#include <gsl/gsl_block_int.h>

#include <gsl/gsl_block_ushort.h>
#include <gsl/gsl_block_short.h>

#include <gsl/gsl_block_uchar.h>
#include <gsl/gsl_block_char.h>

#endif /* __GSL_BLOCK_H__ */

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: ielomariala/Hex-Game


text:
#include <string.h>
#include <cblas.h>

#include "util.h"
#include "dense.h"

denseLayer denseLayer_init(int M, int N)
{
     denseLayer dl; DENSEL_INIT(dl); dl.M=M; dl.N=N;
     return dl;
}

void denseLayer_free(denseLayer *dl)
{
     ftens_free(&dl->W);   ftens_free(&dl->b);
     ftens_free(&dl->dW);  ftens_free(&dl->db);
     ftens_free(&dl->out); ftens_free(&dl->in);
}

void denseLayer_print_shape(denseLayer *dl)
{
     printf("dense: %d %d\n", dl->M, dl->N);
}

void denseLayer_set(ftens *W, denseLayer *dl)
{
     const int M=W->M, N=W->N;
     ASSERT(W->D==1 && W->L==1, "err: dense shape\n");
     ftens_free(&dl->W);
     dl->M = M; dl->N = N;
     dl->W = ftens_copy(W);
}

void denseLayer_forward(ftens *t, denseLayer *dl, int save)
{
     const int D=t->D, M=dl->M, N=dl->N;
     ASSERT(t->MNL == dl->N,  "err: dense shape\n");

     if (save) {
          int M=t->M, N=t->N, L=t->L;
          if (!dl->in.data) dl->in = ftens_init(D,M,N,L);
          memcpy(dl->in.data, t->data, t->bytes);
     }

     if (!dl->out.data) dl->out = ftens_init(D, 1, M, 1);
     const float *a=dl->W.data;
     const float *b=t->data;
     float       *c=dl->out.data;

     cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasTrans,
                 D, M, N, 1, b, N, a, N, 0, c, M);

}

void denseLayer_backward(ftens *dout, denseLayer *dl)
{
     fprintf(stderr, "not implemented yet\n");
     exit(-2);
}

max_stars_repo_name: pppoe/espresso

max_issues_repo_name: pppoe/espresso

max_forks_repo_name: pppoe/espresso


text:
#pragma once

#include "Cesium3DTiles/BoundingVolume.h"
#include "Cesium3DTiles/Library.h"
#include "Cesium3DTiles/TileContentLoadInput.h"
#include "Cesium3DTiles/TileContentLoadResult.h"
#include "Cesium3DTiles/TileContentLoader.h"
#include "Cesium3DTiles/TileID.h"
#include "Cesium3DTiles/TileRefine.h"

#include <gsl/span>
#include <spdlog/fwd.h>

#include <cstddef>
#include <functional>
#include <memory>
#include <optional>
#include <unordered_map>

namespace Cesium3DTiles {
class TileContent;
class Tileset;

/**
 * @brief Creates {@link TileContentLoadResult} objects from a {@link
 * TileContentLoadInput}.
 *
 * The class offers a lookup functionality for registering {@link
 * TileContentLoader} instances that can create {@link TileContentLoadResult}
 * instances from a {@link TileContentLoadInput}.
 *
 * The loaders are registered based on the magic header or the content type
 * of the input data. The raw data (i.e. the `data` of the {@link
 * TileContentLoadInput}) is usually received as a response to a  network
 * request, and the first four bytes of the raw data form the magic header.
 * Based on this header or the content type of the network response, the loader
 * that will be used for processing the input can be looked up.
 */
class CESIUM3DTILES_API TileContentFactory final {
public:
  TileContentFactory() = delete;

  /**
   * @brief Register the given function for the given magic header.
   *
   * The given magic header is a 4-character string. It will be compared
   * to the first 4 bytes of the raw input data, to decide whether the
   * given factory function should be used to create the
   * {@link TileContentLoadResult} from the input data.
   *
   * @param magic The string describing the magic header.
   * @param pLoader The loader that will be used to create the tile content.
   */
  static void registerMagic(
      const std::string& magic,
      const std::shared_ptr<TileContentLoader>& pLoader);

  /**
   * @brief Register the given function for the given content type.
   *
   * The given string describes the content type of a network response.
   * It is used for deciding whether the given loader should be
   * used to create the {@link TileContentLoadResult} from a
   * {@link TileContentLoadInput} with the same `contentType`.
   *
   * @param contentType The string describing the content type.
   * @param pLoader The loader that will be used to create the tile content
   */
  static void registerContentType(
      const std::string& contentType,
      const std::shared_ptr<TileContentLoader>& pLoader);

  /**
   * @brief Creates the {@link TileContentLoadResult} from the given {@link
   * TileContentLoadInput}.
   *
   * This will look up the {@link TileContentLoader} that can be used to
   * process the given input data, based on all loaders that
   * have been registered with {@link TileContentFactory::registerMagic}
   * or {@link TileContentFactory::registerContentType}.
   *
   * It will first try to find a loader based on the magic header
   * of the `data` in the given input. If no matching loader is found, then
   * it will look up a loader based on the `contentType` of the given input.
   * (This will ignore any parameters that may appear after a `;` in the
   * `contentType` string).
   *
   * If no such loader is found then `nullptr` is returned.
   *
   * If a matching loader is found, it will be applied to the given
   * input, and the result will be returned.
   *
   * @param input The {@link TileContentLoadInput}.
   * @return The {@link TileContentLoadResult}, or `nullptr` if there is
   * no loader registered for the magic header of the given
   * input, and no loader for the content type of the input.
   */
  static std::unique_ptr<TileContentLoadResult>
  createContent(const TileContentLoadInput& input);

private:
  static std::optional<std::string>
  getMagic(const gsl::span<const std::byte>& data);

  static std::unordered_map<std::string, std::shared_ptr<TileContentLoader>>
      _loadersByMagic;
  static std::unordered_map<std::string, std::shared_ptr<TileContentLoader>>
      _loadersByContentType;
};

} // namespace Cesium3DTiles

max_stars_repo_name: TJKoury/cesium-native

max_issues_repo_name: MAPSWorks/cesium-native

max_forks_repo_name: MAPSWorks/cesium-native


text:
/* Starting from version 7.8, MATLAB BLAS expects ptrdiff_t arguments for integers */
#if MATLAB_VERSION >= 0x0708
#include <stddef.h>
#include <stdlib.h>
#endif
#include <string.h>

/* Define MX_HAS_INTERLEAVED_COMPLEX for version <9.4 */
#ifndef MX_HAS_INTERLEAVED_COMPLEX
#define MX_HAS_INTERLEAVED_COMPLEX 0
#endif

/* Starting from version 7.6, MATLAB BLAS is seperated */
#if MATLAB_VERSION >= 0x0705
#include <blas.h>
#endif
#include <lapack.h>
#include "f2c.h"

// Conversion of optimal problems with coupling weighting terms to standard problems
#define sb02mt FORTRAN_WRAPPER(sb02mt)
extern void sb02mt(
        const char *jobg,
        const char *jobl,
        const char *fact,
        const char *uplo,
        const ptrdiff_t *n,
        const ptrdiff_t *m,
        double *a,
        const ptrdiff_t *lda,
        double *b,
        const ptrdiff_t *ldb,
        double *q,
        const ptrdiff_t *ldq,
        double *r,
        const ptrdiff_t *ldr,
        double *l,
        const ptrdiff_t *ldl,
        ptrdiff_t *ipiv,
        ptrdiff_t *oufact,
        double *g,
        const ptrdiff_t *ldg,
        ptrdiff_t *iwork,
        double *dwork,
        const ptrdiff_t *ldwork,
        ptrdiff_t *info
        );

// Conversion of optimal problems with coupling weighting terms to standard problems (more flexibility)
#define sb02mx FORTRAN_WRAPPER(sb02mx)
extern void sb02mx(
        const char *jobg,
        const char *jobl,
        const char *fact,
        const char *uplo,
        const char *trans,
        const char *flag,
        const char *def,
        const ptrdiff_t *n,
        const ptrdiff_t *m,
        double *a,
        const ptrdiff_t *lda,
        double *b,
        const ptrdiff_t *ldb,
        double *q,
        const ptrdiff_t *ldq,
        double *r,
        const ptrdiff_t *ldr,
        double *l,
        const ptrdiff_t *ldl,
        ptrdiff_t *ipiv,
        ptrdiff_t *oufact,
        double *g,
        const ptrdiff_t *ldg,
        ptrdiff_t *iwork,
        double *dwork,
        const ptrdiff_t *ldwork,
        ptrdiff_t *info
        );

// Constructing the 2n-by-2n Hamiltonian or symplectic matrix for linear-quadratic optimization problems
#define sb02mu FORTRAN_WRAPPER(sb02mu)
extern void sb02mu(
        const char *dico,
        const char *hinv,
        const char *uplo,
        const ptrdiff_t *n,
        double *a,
        const ptrdiff_t *lda,
        const double *g,
        const ptrdiff_t *ldg,
        const double *q,
        const ptrdiff_t *ldq,
        double *s,
        const ptrdiff_t *lds,
        ptrdiff_t *iwork,
        double *dwork,
        const ptrdiff_t *ldwork,
        ptrdiff_t *info
        );

// Constructing the 2n-by-2n Hamiltonian or symplectic matrix for linear-quadratic optimization problems (improved)
#define sb02ru FORTRAN_WRAPPER(sb02ru)
extern void sb02ru(
        const char *dico,
        const char *hinv,
        const char *trana,
        const char *uplo,
        const ptrdiff_t *n,
        const double *a,
        const ptrdiff_t *lda,
        double *g,
        const ptrdiff_t *ldg,
        double *q,
        const ptrdiff_t *ldq,
        double *s,
        const ptrdiff_t *lds,
        ptrdiff_t *iwork,
        double *dwork,
        const ptrdiff_t *ldwork,
        ptrdiff_t *info
        );

// Constructing the extended Hamiltonian or symplectic matrix pairs for linear-quadratic optimization problems, and compressing them to 2N-by-2N matrices
#define sb02oy FORTRAN_WRAPPER(sb02oy)
extern void sb02oy(
        const char *type,
        const char *dico,
        const char *jobb,
        const char *fact,
        const char *uplo,
        const char *jobl,
        const char *jobe,
        const ptrdiff_t *n,
        const ptrdiff_t *m,
        const ptrdiff_t *p,
        const double *a,
        const ptrdiff_t *lda,
        const double *b,
        const ptrdiff_t *ldb,
        const double *q,
        const ptrdiff_t *ldq,
        const double *r,
        const ptrdiff_t *ldr,
        const double *l,
        const ptrdiff_t *ldl,
        const double *e,
        const ptrdiff_t *lde,
        double *af,
        const ptrdiff_t *ldaf,
        double *bf,
        const ptrdiff_t *ldbf,
        double *tol,
        ptrdiff_t *iwork,
        double *dwork,
        const ptrdiff_t *ldwork,
        ptrdiff_t *info
        );

// Optimal state feedback matrix for an optimal control problem
#define sb02nd FORTRAN_WRAPPER(sb02nd)
extern void sb02nd(
        const char *dico,
        const char *fact,
        const char *uplo,
        const char *jobl,
        const ptrdiff_t *n,
        const ptrdiff_t *m,
        const ptrdiff_t *p,
        const double *a,
        const ptrdiff_t *lda,
        double *b,
        const ptrdiff_t *ldb,
        double *r,
        const ptrdiff_t *ldr,
        ptrdiff_t *ipiv,
        const double *l,
        const ptrdiff_t *ldl,
        double *x,
        const ptrdiff_t *ldx,
        const double *rnorm,
        double *f,
        const ptrdiff_t *ldf,
        ptrdiff_t *oufact,
        ptrdiff_t *iwork,
        double *dwork,
        const ptrdiff_t *ldwork,
        ptrdiff_t *info
        );

// Solution of continuous- or discrete-time algebraic Riccati equations for descriptor systems
#define sg02nd FORTRAN_WRAPPER(sg02nd)
extern void sg02nd(
        const char *dico,
        const char *jobe,
        const char *job,
        const char *jobx,
        const char *fact,
        const char *uplo,
        const char *jobl,
        const char *trans,
        const ptrdiff_t *n,
        const ptrdiff_t *m,
        const ptrdiff_t *p,
        const double *a,
        const ptrdiff_t *lda,
        const double *e,
        const ptrdiff_t *lde,
        double *b,
        const ptrdiff_t *ldb,
        double *r,
        const ptrdiff_t *ldr,
        ptrdiff_t *ipiv,
        const double *l,
        const ptrdiff_t *ldl,
        double *x,
        const ptrdiff_t *ldx,
        const double *rnorm,
        double *k,
        const ptrdiff_t *ldk,
        double *h,
        const ptrdiff_t *ldh,
        double *xe,
        const ptrdiff_t *ldxe,
        ptrdiff_t *oufact,
        ptrdiff_t *iwork,
        double *dwork,
        const ptrdiff_t *ldwork,
        ptrdiff_t *info
        );

// Column interchanges in a complex matrix
#define ma02gz FORTRAN_WRAPPER(ma02gz)
extern void ma02gz(
        const ptrdiff_t *n,
        double *a,
        const ptrdiff_t *lda,
        const ptrdiff_t *k1,
        const ptrdiff_t *k2,
        const ptrdiff_t *ipiv,
        const ptrdiff_t *incx
        );

// Solution of linear equations X op(A) = B
#define mb02vd FORTRAN_WRAPPER(mb02vd)
extern void mb02vd(
        const char *trans,
        const ptrdiff_t *m,
        const ptrdiff_t *n,
        double *a,
        const ptrdiff_t *lda,
        ptrdiff_t *ipiv,
        double *b,
        const ptrdiff_t *ldb,
        ptrdiff_t *info
        );

// Periodic Hessenberg form of a product of p matrices using orthogonal similarity transformations
#define mb03vd FORTRAN_WRAPPER(mb03vd)
extern void mb03vd(
        const ptrdiff_t *n,
        const ptrdiff_t *p,
        const ptrdiff_t *ilo,
        const ptrdiff_t *ihi,
        double *a,
        const ptrdiff_t *lda1,
        const ptrdiff_t *lda2,
        double *tau,
        const ptrdiff_t *ldtau,
        double *dwork,
        ptrdiff_t *info
        );

// Orthogonal matrices for reduction to periodic Hessenberg form of a product of matrices
#define mb03vy FORTRAN_WRAPPER(mb03vy)
extern void mb03vy(
        const ptrdiff_t *n,
        const ptrdiff_t *p,
        const ptrdiff_t *ilo,
        const ptrdiff_t *ihi,
        double *a,
        const ptrdiff_t *lda1,
        const ptrdiff_t *lda2,
        double *tau,
        const ptrdiff_t *ldtau,
        double *dwork,
        const ptrdiff_t *ldwork,
        ptrdiff_t *info
        );

// Schur decomposition and eigenvalues of a product of matrices in periodic Hessenberg form
#define mb03wd FORTRAN_WRAPPER(mb03wd)
extern void mb03wd(
        const char *job,
        const char *compz,
        const ptrdiff_t *n,
        const ptrdiff_t *p,
        const ptrdiff_t *ilo,
        const ptrdiff_t *ihi,
        const ptrdiff_t *iloz,
        const ptrdiff_t *ihiz,
        double *h,
        const ptrdiff_t *ldh1,
        const ptrdiff_t *ldh2,
        double *z,
        const ptrdiff_t *ldz1,
        const ptrdiff_t *ldz2,
        double *wr,
        double *wi,
        double *dwork,
        const ptrdiff_t *ldwork,
        ptrdiff_t *info
        );

// Computes the general product of K complex scalars trying to avoid over- and underflow.
#define zlapr1 FORTRAN_WRAPPER(zlapr1)
extern int zlapr1(
        doublereal *base,
        integer *k,
        integer *s,
        doublecomplex *a,
        integer *inca,
        doublecomplex *alpha,
        doublecomplex *beta,
        integer *scal
        );

// Finding the eigenvalues of the complex generalized matrix product.
#define zpgeqz FORTRAN_WRAPPER(zpgeqz)
extern int zpgeqz(
        char *job,
        char *compq,
        integer *k,
        integer *n,
        integer *ilo,
        integer *ihi,
        integer *s,
        doublecomplex *a,
        integer *lda1,
        integer *lda2,
        doublecomplex *alpha,
        doublecomplex *beta,
        integer *scal,
        doublecomplex *q,
        integer *ldq1,
        integer *ldq2,
        doublereal *dwork,
        integer *ldwork,
        doublecomplex *zwork,
        integer *lzwork,
        integer *info
        );

// Swaps adjacent diagonal 1-by-1 blocks in a complex generalized matrix product.
#define zpgex2 FORTRAN_WRAPPER(zpgex2)
extern int zpgex2(
        logical *wantq,
        integer *k,
        integer *n,
        integer *j,
        integer *s,
        doublecomplex *a,
        integer *lda1,
        integer *lda2,
        doublecomplex *q,
        integer *ldq1,
        integer *ldq2,
        doublecomplex *zwork,
        integer *info
        );

// Swaps adjacent diagonal 1-by-1 blocks in a complex generalized matrix product.
#define zpghrd FORTRAN_WRAPPER(zpghrd)
extern int zpghrd(
        char *compq,
        integer *k,
        integer *n,
        integer *ilo,
        integer *ihi,
        integer *s,
        doublecomplex *a,
        integer *lda1,
        integer *lda2,
        doublecomplex *q,
        integer *ldq1,
        integer *ldq2,
        doublereal *dwork,
        integer *ldwork,
        doublecomplex *zwork,
        integer *lzwork,
        integer *info
        );

// Reorders the periodic Schur decomposition of a complex generalized matrix product.
#define zpgord FORTRAN_WRAPPER(zpgord)
extern int zpgord(
        logical *wantq,
        integer *k,
        integer *n,
        integer *s,
        logical *select,
        doublecomplex *a,
        integer *lda1,
        integer *lda2,
        doublecomplex *alpha,
        doublecomplex *beta,
        integer *scal,
        doublecomplex *q,
        integer *ldq1,
        integer *ldq2,
        integer *m,
        doublecomplex *zwork,
        integer *lzwork,
        integer *info
        );

max_stars_repo_name: iwoodsawyer/dpre

max_issues_repo_name: iwoodsawyer/dpre

max_forks_repo_name: iwoodsawyer/dpre


text:
/*This class introduces a method to work with experimental fitness landscapes
given in a data file in the format
sequence          fitness
1   0   1   ...   1.987
and brings the adjacency matrix as well as other helpful properties

(cc) Johannes Neidhart, 2012*/

#ifndef LANDSCAPE_H
#define LANDSCAPE_H

#include <iostream>
#include <vector>
#include <fstream>
#include <sstream>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_multifit.h>
#include <gsl/gsl_multifit_nlin.h>
#include <vector>
#include <algorithm>
#include <cmath>
#include <map>
#include <numeric>

using namespace std;

class Clandscape {
    private:
        vector < vector < int> > linekey;
        int n, L;
        string name;
        int delta(int, int);
        int hamming( vector < int >, vector < int > );
    public:
        /*Clandscape gets a vector containing the fitness values, in the order of the
        linenumbers of the data file, a vector of vectors, containing the correspondance
        of the line number and the sequence. Sequence length and number of sequences are
        also given.
        */

        double binomial (int, int);
        Clandscape (string, int , int);
        int distance (int, int);
        vector < double > fitness;
        void fitit( vector < double >, vector < double >, vector < double > &);
        void fit(vector <double>, vector < double > &, vector < double > &, double);
        vector < vector < int > > adjacency;
        vector < double > spectrum();
        vector < vector < int > > hypercube;

};

#endif

max_stars_repo_name: PidZero/Amplitudemaker

max_issues_repo_name: PidZero/Amplitudemaker

max_forks_repo_name: PidZero/Amplitudemaker


text:
#ifndef setup_h
#define setup_h

#include <ceed.h>
#include <petsc.h>
#include <petscdmplex.h>
#include <petscfe.h>
#include <petscksp.h>
#include <stdbool.h>
#include <string.h>
#include "problems/problems.h"
#include "include/cl-options.h"
#include "include/matops.h"
#include "include/misc.h"
#include "include/structs.h"
#include "include/setup-dm.h"
#include "include/setup-libceed.h"
#include "include/utils.h"

#endif // setup_h

max_stars_repo_name: AdelekeBankole/libCEED

max_issues_repo_name: AdelekeBankole/libCEED

max_forks_repo_name: AdelekeBankole/libCEED


text:
#include <cblas.h>

const char* dgemm_desc = "Reference dgemm.";

/*
 * This routine performs a dgemm operation
 *  C := C + A * B
 * where A, B, and C are lda-by-lda matrices stored in column-major format.
 * On exit, A and B maintain their input values.
 * This function wraps a call to the BLAS-3 routine DGEMM,
 * via the standard FORTRAN interface - hence the reference semantics.
 */
void square_dgemm(int n, double* A, double* B, double* C) {
    cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, n, n, n, 1., A, n, B, n, 1., C, n);
}


text:
#ifndef UTILS_STRINGS_ONE_AWAY_H_INCLUDED
#define UTILS_STRINGS_ONE_AWAY_H_INCLUDED

#include <string> // std::string
#include <gsl/gsl>

bool isOneAway(const std::string &, const std::string &);

#endif


text:
/*==========================================================
 * pca.c - Principal component analysis of diffraction data
 *
 * Calculates the two main components of the data
 * by Principal Component Analysis, as detailed in
 * Bernhardt et.al (2016), Biophys. J.
 *
 * The calling syntax is:
 *
 *		pca(data, mask, qx, qy, w, phi, m, n)
 *
 * Copyright 2017 Institute for X-ray Physics (University of Göttingen)
 *
 * Permission is hereby granted, free of charge, to any person obtaining
 * a copy of this software and associated documentation files (the "Software"),
 * to deal in the Software without restriction, including without limitation
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
 * and/or sell copies of the Software, and to permit persons to whom the
 * Software is furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included
 * in all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
 * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
 * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
 * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
 * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE
 * OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
 *
 *========================================================*/

#include <math.h>
#include <gsl/gsl_fft_real.h>
#include <gsl/gsl_fft_halfcomplex.h>

#define PI 3.14159265

/* The computational routine */
void symmetry(double *x, double *bgr, double sigma, double noise_level, int m, int n)
{
    int delta_theta = 20; /* hard coded theta window size */

    int col;
    int row;

    /* calculate snr */
    for (row=0; row<m; row++) {
        for (col=0; col<n; col++) {
            if (bgr[row*m + col] == 0){
                x[row*m + col] = 0;
            } else {
                x[row*m + col] = (x[row*m + col] - bgr[row*m + col]) / sqrt(bgr[row*m + col]);
        }
    }

    /* next we implement
% Gaussian filter (remove noise)
h = fspecial('gaussian', [10 10], sigma);
h = h./sum(h(:));
snr = conv2(snr,h,'same');

% binary image (snr threshold)
binary_image = (snr >= noise_level);
        */

    /* prepare fft */
    work = gsl_fft_real_workspace_alloc (n);
    real = gsl_fft_real_wavetable_alloc (n);

    /* fourier transform the rows */
    for (int row = 0; row < m; row++){
        double* data = x[row*m];
        gsl_fft_real_transform (data, 1, n, real, work);
        gsl_fft_real_wavetable_free (real);
        /* apply filter HERE */
        hc = gsl_fft_halfcomplex_wavetable_alloc (n);
        gsl_fft_halfcomplex_inverse (data, 1, n, hc, work);
        gsl_fft_halfcomplex_wavetable_free (hc);
        gsl_fft_real_workspace_free (work);
    }
    /* fourier transform the rows */

	/* apply gauss filter */

	/* fourier-backtransform the rows */
    /* fourier-backtransform the cols */

	/* apply threshold */
    for (row=0; row<m; row++) {
        for (col=0; col<n; col++) {
            if (x[row*m + col] > noise_level){
                x[row*m + col] = 1;
            } else {
                x[row*m + col] = 0;
        }
    }

}

max_stars_repo_name: irpgoe/nanodiffraction

max_issues_repo_name: irpgoe/nanodiffraction

max_forks_repo_name: irpgoe/nanodiffraction


text:
#include "redisgears.h"
#include "redisai.h"
#include "minmax_heap.h"
#include <math.h>
#include "redisgears_memory.h"
#include <cblas.h>
#include <sys/time.h>

#define STR1(a) #a
#define STR(e) STR1(e)

#define MIN(a,b) (((a)<(b))?(a):(b))
#define MAX(a,b) (((a)>(b))?(a):(b))

static RedisModuleCtx* staticCtx;

static RecordType* ScoreRecordType = NULL;
static RecordType* HeapRecordType = NULL;

#define VEC_SIZE 128

#define VEC_HOLDER_SIZE 1024 * 1024

typedef struct VecsHolder VecsHolder;

typedef struct VecDT{
    size_t index;
    VecsHolder* holder;
    RedisModuleString* keyName;
}VecDT;

typedef struct VecsHolder{
    size_t size;
    VecDT* vecDT[VEC_HOLDER_SIZE];
    float vecs[VEC_HOLDER_SIZE * VEC_SIZE];
}VecsHolder;

#define HOLDER_VECDT(h, i) (h->vecDT[i])
#define HOLDER_VEC(h, i) (h->vecs[i * VEC_SIZE])

VecsHolder** vecList = NULL;
RedisModuleType *vecRedisDT;

typedef struct VecReaderCtx{
    size_t index;
    Record** pendings;
    float vec[VEC_SIZE];
    size_t topK;
}VecReaderCtx;

typedef struct TopKArg{
    size_t topK;
}TopKArg;

typedef struct ScoreRecord{
    Record baseRecord;
    RedisModuleString* key;
    float score;
}ScoreRecord;

typedef struct HeapRecord{
    Record baseRecord;
    heap_t* heap;
}HeapRecord;

static VecReaderCtx* VecReaderCtx_Create(float* data, size_t topK){
    VecReaderCtx* ctx = RG_ALLOC(sizeof(*ctx));
    ctx->index = 0;
    ctx->pendings = array_new(Record*, 10);
    ctx->topK = topK;
    if(data){
        memcpy(ctx->vec, data, VEC_SIZE * sizeof(*data));
        float denom_vec = cblas_snrm2(VEC_SIZE, ctx->vec, 1);
        for(size_t i = 0 ; i < VEC_SIZE ; ++i){
            ctx->vec[i] /= denom_vec;
        }
    }
    return ctx;
}

static void VecReaderCtx_Free(VecReaderCtx* ctx){
    for(size_t i = 0 ; i < array_len(ctx->pendings) ; ++i){
        Record* r = ctx->pendings[i];
        RedisGears_FreeRecord(r);
    }

    array_free(ctx->pendings);

    RG_FREE(ctx);
}

static Record* to_score_records(ExecutionCtx* rctx, Record *data, void* arg){
    HeapRecord* hr = (HeapRecord*)data;

    Record* lr = RedisGears_ListRecordCreate(hr->heap->count);

    size_t count = hr->heap->count;
    for(size_t i = 0 ; i < count ; ++i){
        ScoreRecord* sr = mmh_pop_min(hr->heap);
        RedisGears_ListRecordAdd(lr, &sr->baseRecord);
    }

    RedisGears_FreeRecord(data);

    return lr;
}

static int heap_cmp(const void *a, const void *b, const void *udata){
    ScoreRecord* s1 = (ScoreRecord*)a;
    ScoreRecord* s2 = (ScoreRecord*)b;
    if(s1->score < s2->score){
        return -1;
    }else if(s1->score > s2->score){
        return 1;
    }else {
        return 0;
    }
}

static Record* top_k(ExecutionCtx* rctx, Record *accumulate, Record *r, void* arg){
    TopKArg* topKArg = arg;

    HeapRecord* heap = (HeapRecord*)accumulate;

    if(!heap){
        heap = (HeapRecord*)RedisGears_RecordCreate(HeapRecordType);
        heap->heap = mmh_init_with_size(topKArg->topK, heap_cmp, NULL, (mmh_free_func)RedisGears_FreeRecord);
    }

    heap_t* h = heap->heap;

    ScoreRecord* currSr = (ScoreRecord*)r;

    // If the queue is not full - we just push the result into it
    // If the pool size is 0 we always do that, letting the heap grow dynamically
    if (h->count < topKArg->topK) {
        mmh_insert(h, r);
    } else {
        // find the min result
        ScoreRecord* sr = mmh_peek_min(h);
        if(sr->score < currSr->score){
            void* temp = mmh_pop_min(h);
            RedisModule_Assert(temp == sr);
            RedisGears_FreeRecord(&(sr->baseRecord));
            mmh_insert(h, r);
        }else{
            RedisGears_FreeRecord(r);
        }
    }

    return &(heap->baseRecord);
}

static void on_done(ExecutionPlan* ctx, void* privateData){
    RedisModuleBlockedClient *bc = privateData;
    RedisModuleCtx *rctx = RedisModule_GetThreadSafeContext(bc);
    RedisGears_ReturnResultsAndErrors(ctx, rctx);
    RedisModule_UnblockClient(bc, NULL);
    RedisGears_DropExecution(ctx);
    RedisModule_FreeThreadSafeContext(rctx);
}

VecDT* vec_insert(RedisModuleString *keyName, const float* data){
    VecsHolder* holder = NULL;
    if(!vecList){
        vecList = array_new(VecsHolder*, 1);
    }
    if(array_len(vecList) == 0){
        holder = RG_CALLOC(1, sizeof(VecsHolder));
        vecList = array_append(vecList, holder);
    }else{
        holder = vecList[array_len(vecList) - 1];
    }

    if(holder->size >= VEC_HOLDER_SIZE){
        // we need to create a new holder
        holder = RG_CALLOC(1, sizeof(VecsHolder));
        vecList = array_append(vecList, holder);
    }

    float* v = &HOLDER_VEC(holder, holder->size);
    memcpy(v, data, sizeof(float) * VEC_SIZE);

    float demon = cblas_snrm2(VEC_SIZE, v, 1);

    for(size_t i = 0 ; i < VEC_SIZE ; ++i){
        v[i] /= demon;
    }

    VecDT* vDT = RG_CALLOC(1, sizeof(*vDT));
    vDT->holder = holder;
    vDT->index = holder->size;
    vDT->keyName = keyName;
    RedisModule_RetainString(NULL, vDT->keyName);
    HOLDER_VECDT(holder, holder->size) = vDT;

    ++holder->size;

    return vDT;
}

/*
 * rg.vec_add <k> <blob>
 */
int vec_add_command(RedisModuleCtx *ctx, RedisModuleString **argv, int argc){
    if(argc != 3){
        return RedisModule_WrongArity(ctx);
    }

    size_t dataLen;
    float* data = (float*)RedisModule_StringPtrLen(argv[2], &dataLen);
    if(dataLen != (VEC_SIZE * sizeof(float))){
        RedisModule_ReplyWithError(ctx, "Given blob is not float vector of size " STR(VEC_SIZE));
        return REDISMODULE_OK;
    }

    RedisModuleKey *kp = RedisModule_OpenKey(ctx, argv[1], REDISMODULE_WRITE);
    if(RedisModule_KeyType(kp) != REDISMODULE_KEYTYPE_EMPTY){
        RedisModule_ReplyWithError(ctx, "Key is not empty");
        RedisModule_CloseKey(kp);
        return REDISMODULE_OK;
    }

    VecDT* vDT = vec_insert(argv[1], data);

    RedisModule_ModuleTypeSetValue(kp, vecRedisDT, vDT);

    RedisModule_CloseKey(kp);

    RedisModule_ReplicateVerbatim(ctx);

    RedisModule_ReplyWithSimpleString(ctx, "OK");

    return REDISMODULE_OK;
}

/*
 * rg.vec_sim <k> <blob>
 */
int vec_sim_command(RedisModuleCtx *ctx, RedisModuleString **argv, int argc){

    if(argc != 3){
        return RedisModule_WrongArity(ctx);
    }

    char* err = NULL;

    long long topK;
    if(RedisModule_StringToLongLong(argv[1], &topK) != REDISMODULE_OK){
        RedisModule_ReplyWithError(ctx, "Failed extracting <k>");
        return REDISMODULE_OK;
    }

    size_t dataSize;
    float* data = (float*)RedisModule_StringPtrLen(argv[2], &dataSize);
    if(dataSize != (VEC_SIZE * sizeof(float))){
        RedisModule_ReplyWithError(ctx, "Given blob is not at the right size");
        return REDISMODULE_OK;
    }

    TopKArg* topKArg2 = RG_ALLOC(sizeof(*topKArg2));
    topKArg2->topK = topK;

    FlatExecutionPlan* fep = RGM_CreateCtx(VecReader, &err);

    VecReaderCtx* rCtx = VecReaderCtx_Create(data, topK);

    RGM_Collect(fep);

    RGM_Accumulate(fep, top_k, topKArg2);

    RGM_FlatMap(fep, to_score_records, NULL);

    ExecutionPlan* ep = RGM_Run(fep, ExecutionModeAsync, rCtx, NULL, NULL, &err);

    if(!ep){
        RedisGears_FreeFlatExecution(fep);
        RedisModule_ReplyWithError(ctx, err);
        return REDISMODULE_OK;
    }

    RedisModuleBlockedClient *bc = RedisModule_BlockClient(ctx, NULL, NULL, NULL, 0);
    RedisGears_AddOnDoneCallback(ep, on_done, bc);
    RedisGears_FreeFlatExecution(fep);

    return REDISMODULE_OK;
}

static int ScoreRecord_SendReply(Record* base, RedisModuleCtx* rctx){
    ScoreRecord* sr = (ScoreRecord*)base;
    RedisModule_ReplyWithArray(rctx, 2);
    RedisModule_ReplyWithString(rctx, sr->key);
    RedisModule_ReplyWithDouble(rctx, sr->score);
    return REDISMODULE_OK;
}

static int ScoreRecord_RecordSerialize(ExecutionCtx* ctx, Gears_BufferWriter* bw, Record* base){
    ScoreRecord* sr = (ScoreRecord*)base;
    const char* keyStr = RedisModule_StringPtrLen(sr->key, NULL);
    RedisGears_BWWriteString(bw, keyStr);
    RedisGears_BWWriteBuffer(bw, (char*)(&(sr->score)), sizeof(sr->score));
    return REDISMODULE_OK;
}

static Record* ScoreRecord_RecordDeserialize(ExecutionCtx* ctx, Gears_BufferReader* br){
    ScoreRecord* sr = (ScoreRecord*)RedisGears_RecordCreate(ScoreRecordType);
    const char* keyStr = RedisGears_BRReadString(br);
    sr->key = RedisModule_CreateString(NULL, keyStr, strlen(keyStr));
    size_t len;
    char* data = RedisGears_BRReadBuffer(br, &len);
    RedisModule_Assert(len == sizeof(float));
    sr->score = *((float*)(data));

    return &sr->baseRecord;
}

static void ScoreRecord_RecordFree(Record* base){
    ScoreRecord* sr = (ScoreRecord*)base;
    if(sr->key){
        RedisModule_FreeString(NULL, sr->key);
    }

}

static int HeapRecord_SendReply(Record* base, RedisModuleCtx* rctx){
    RedisModule_Assert(false);
    return REDISMODULE_OK;
}

static int HeapRecord_Serialize(ExecutionCtx* ctx, Gears_BufferWriter* bw, Record* base){
    RedisModule_Assert(false);
    return REDISMODULE_OK;
}

static Record* HeapRecord_Deserialize(ExecutionCtx* ctx, Gears_BufferReader* br){
    RedisModule_Assert(false);
    return REDISMODULE_OK;
}

static void HeapRecord_Free(Record* base){
    HeapRecord* hr = (HeapRecord*)base;
    mmh_free(hr->heap);
}

static void TopKArg_ObjectFree(void* arg){
    RG_FREE(arg);
}

static void* TopKArg_ArgDuplicate(void* arg){
    TopKArg* topK = arg;
    TopKArg* topKDup = RG_ALLOC(sizeof(*topKDup));
    topKDup->topK = topK->topK;
    return topKDup;
}

static int TopKArg_ArgSerialize(FlatExecutionPlan* fep, void* arg, Gears_BufferWriter* bw, char** err){
    TopKArg* topK = arg;
    RedisGears_BWWriteLong(bw, topK->topK);
    return REDISMODULE_OK;
}

static void* TopKArg_ArgDeserialize(FlatExecutionPlan* fep, Gears_BufferReader* br, int version, char** err){
    TopKArg* topKDup = RG_ALLOC(sizeof(*topKDup));
    topKDup->topK = RedisGears_BRReadLong(br);
    return topKDup;
}

static char* TopKArg_ArgToString(void* arg){
    return RG_STRDUP("I am topk argument :)");
}

#define TopKTypeVersion 1

#define VS_PLUGIN_NAME "VECTOR_SIM"
#define REDISGEARSJVM_PLUGIN_VERSION 1

#define VEC_TYPE_VERSION 1

static void* VecDT_Load(RedisModuleIO *rdb, int encver){
    RedisModuleString *keyName = RedisModule_LoadString(rdb);
    size_t dataLen;
    float* data = (float*)RedisModule_LoadStringBuffer(rdb, &dataLen);
    RedisModule_Assert(dataLen == sizeof(float) * VEC_SIZE);

    VecDT* vDT = vec_insert(keyName, data);

    RedisModule_FreeString(NULL, keyName);
    RedisModule_Free(data);

    return vDT;
}

static void VecDT_Save(RedisModuleIO *rdb, void *value){
    VecDT* vDT = value;

    RedisModule_SaveString(rdb, vDT->keyName);
    RedisModule_SaveStringBuffer(rdb, (char*)&HOLDER_VEC(vDT->holder, vDT->index), sizeof(float) * VEC_SIZE);
}

static void VecDT_Free(void *value){
    VecDT* vDT = value;
    VecsHolder* holder = vDT->holder;
    size_t index = vDT->index;

    RedisModule_FreeString(NULL, vDT->keyName);
    RG_FREE(vDT);

    if(!holder){
        // we probably inside flush, the vector DT was detached and we can just return.
        return;
    }

    // get the last vector
    VecsHolder* lastVH = vecList[array_len(vecList) - 1];
    --lastVH->size;
    VecDT* lastVDT = HOLDER_VECDT(lastVH, lastVH->size);

    if(lastVDT != vDT){
        // swap last with current
        memmove(&HOLDER_VEC(holder, index), &HOLDER_VEC(lastVH, lastVH->size), VEC_SIZE * sizeof(float));

        HOLDER_VECDT(holder, index) = lastVDT;
        lastVDT->holder = holder;
        lastVDT->index = index;
    }

    if(lastVH->size == 0){
        // free the holder, it has no more data.
        RG_FREE(lastVH);
        if(array_len(vecList) > 1){
            vecList = array_trimm_cap(vecList, array_len(vecList) - 1);
        }else{
            array_free(vecList);
            vecList = NULL;
        }
    }
}

static float scores[VEC_HOLDER_SIZE];

static Record* VecReader_Next(ExecutionCtx* rctx, void* ctx){
//    struct timeval stop, start;
    if(!vecList){
        return NULL;
    }

    VecReaderCtx* readerCtx = ctx;
    RedisModuleCtx* redisCtx = RedisGears_GetRedisModuleCtx(rctx);
    if(array_len(readerCtx->pendings) > 0){
        return array_pop(readerCtx->pendings);
    }

    RedisGears_LockHanlderAcquire(redisCtx);

    const float* b1 = readerCtx->vec;

    while(readerCtx->index < array_len(vecList)){
        VecsHolder* holder = vecList[readerCtx->index++];

        cblas_sgemv(CblasRowMajor, CblasNoTrans, holder->size, VEC_SIZE, 1, holder->vecs, VEC_SIZE, b1, 1, 0, scores, 1);

        for(size_t i = 0 ; i < MIN(holder->size, readerCtx->topK) ; ++i){
            size_t index = cblas_isamax(holder->size, scores, 1);
            ScoreRecord* s = (ScoreRecord*)RedisGears_RecordCreate(ScoreRecordType);
            s->key = HOLDER_VECDT(holder, index)->keyName;
            RedisModule_RetainString(NULL, s->key);
            s->score = scores[index];
            readerCtx->pendings = array_append(readerCtx->pendings, &s->baseRecord);
            scores[index] = 0;
        }

        if(array_len(readerCtx->pendings) > 0){
            RedisGears_LockHanlderRelease(redisCtx);
            return array_pop(readerCtx->pendings);
        }
    }

    RedisGears_LockHanlderRelease(redisCtx);

    return NULL;
}

static void VecReader_Free(void* ctx){
    VecReaderCtx_Free(ctx);
}

static int VecReader_Serialize(ExecutionCtx* ectx, void* ctx, Gears_BufferWriter* bw){
    VecReaderCtx* readerCtx = ctx;
    RedisGears_BWWriteBuffer(bw, (char*)readerCtx->vec, VEC_SIZE * sizeof(float));
    RedisGears_BWWriteLong(bw, readerCtx->topK);
    return REDISMODULE_OK;
}

static int VecReader_Deserialize(ExecutionCtx* ectx, void* ctx, Gears_BufferReader* br){
    VecReaderCtx* readerCtx = ctx;
    size_t dataLen;
    float* data = (float*)RedisGears_BRReadBuffer(br, &dataLen);
    RedisModule_Assert(dataLen == VEC_SIZE * sizeof(*data));

    readerCtx->topK = RedisGears_BRReadLong(br);

    memcpy(readerCtx->vec, data, VEC_SIZE * sizeof(*data));

    return REDISMODULE_OK;
}

static Reader* VecReader_CreateReaderCallback(void* arg){
    VecReaderCtx* ctx = arg;
    if(!ctx){
        ctx = VecReaderCtx_Create(NULL, 0);
    }
    Reader* r = RG_ALLOC(sizeof(*r));
    *r = (Reader){
        .ctx = ctx,
        .next = VecReader_Next,
        .free = VecReader_Free,
        .serialize = VecReader_Serialize,
        .deserialize = VecReader_Deserialize,
    };
    return r;
}

RedisGears_ReaderCallbacks VecReader = {
        .create = VecReader_CreateReaderCallback,
};

static void OnFlush(struct RedisModuleCtx *ctx, RedisModuleEvent eid, uint64_t subevent, void *data){
    if(subevent != REDISMODULE_SUBEVENT_FLUSHDB_START){
        return;
    }

    if(!vecList){
        return;
    }

    // before flush we need to clean all the Vector Holders and disconnect the keys
    for(size_t i = 0 ; i < array_len(vecList) ; ++i){
        VecsHolder* holder = vecList[i];
        for(size_t j = 0 ; j < holder->size ; ++j){
            VecDT* vDT = HOLDER_VECDT(holder, j);
            vDT->holder = NULL;
        }
        RG_FREE(holder);
    }

    array_free(vecList);

    vecList = NULL;
}

int RedisGears_OnLoad(RedisModuleCtx *ctx) {
    openblas_set_num_threads(1);

    if(RedisGears_InitAsGearPlugin(ctx, VS_PLUGIN_NAME, REDISGEARSJVM_PLUGIN_VERSION) != REDISMODULE_OK){
        RedisModule_Log(ctx, "warning", "Failed initialize RedisGears API");
        return REDISMODULE_ERR;
    }

    RedisModule_Log(ctx, "warning", "OpenBlac num of threads: %d", openblas_get_num_threads());

    staticCtx = RedisModule_GetThreadSafeContext(NULL);

    RedisModuleTypeMethods vecDT = {
        .version = REDISMODULE_TYPE_METHOD_VERSION,
        .rdb_load = VecDT_Load,
        .rdb_save = VecDT_Save,
        .free = VecDT_Free,
    };

    vecRedisDT = RedisModule_CreateDataType(ctx, "vec_index", VEC_TYPE_VERSION, &vecDT);
    if (vecRedisDT == NULL) {
        RedisModule_Log(ctx, "error", "Could not create vector type");
        return REDISMODULE_ERR;
    }

    RGM_RegisterReader(VecReader);

    ScoreRecordType = RedisGears_RecordTypeCreate("ScoreRecord",
                                                   sizeof(ScoreRecord),
                                                   ScoreRecord_SendReply,
                                                   ScoreRecord_RecordSerialize,
                                                   ScoreRecord_RecordDeserialize,
                                                   ScoreRecord_RecordFree);

    HeapRecordType = RedisGears_RecordTypeCreate("HeapRecord",
                                                 sizeof(HeapRecord),
                                                 HeapRecord_SendReply,
                                                 HeapRecord_Serialize,
                                                 HeapRecord_Deserialize,
                                                 HeapRecord_Free);

    ArgType* TopKType = RedisGears_CreateType("TopKType",
                                              TopKTypeVersion,
                                              TopKArg_ObjectFree,
                                              TopKArg_ArgDuplicate,
                                              TopKArg_ArgSerialize,
                                              TopKArg_ArgDeserialize,
                                              TopKArg_ArgToString);

    RGM_RegisterMap(to_score_records, NULL);
    RGM_RegisterAccumulator(top_k, TopKType);

    if (RedisModule_CreateCommand(ctx, "rg.vec_sim", vec_sim_command, "readonly", 0, 0, 0) != REDISMODULE_OK) {
        RedisModule_Log(ctx, "warning", "could not register command rg.vec_sim");
        return REDISMODULE_ERR;
    }

    if (RedisModule_CreateCommand(ctx, "rg.vec_add", vec_add_command, "write deny-oom", 1, 1, 1) != REDISMODULE_OK) {
        RedisModule_Log(ctx, "warning", "could not register command rg.vec_add");
        return REDISMODULE_ERR;
    }

    RedisModule_SubscribeToServerEvent(ctx, RedisModuleEvent_FlushDB, OnFlush);

    return REDISMODULE_OK;
}

max_stars_repo_name: rafie/VecSim

max_issues_repo_name: rafie/VecSim

max_forks_repo_name: rafie/VecSim


text:
#include <gsl/gsl_rng.h>
#include <easy_rng.h>
#include <time.h>

#ifndef NRUNS
#define NRUNS 10000000
#endif

void run_test(const gsl_rng_type *gtype, const easy_rng_type *etype) {
	gsl_rng *grng = gsl_rng_alloc(gtype);
	easy_rng *erng = easy_rng_alloc(etype);

	clock_t gstart, gend;
	clock_t estart, eend;
	int i;

	gstart = clock();

	for (i = 0 ; i < NRUNS ; i++)
		gsl_rng_get(grng);

	gend = clock();

	estart = clock();

	for (i = 0 ; i < NRUNS ; i++)
		easy_rng_get(erng);

	eend = clock();

	gsl_rng_free(grng);
	easy_rng_free(erng);

	fprintf(stdout, "Comparing easyRNG's %s with GSL's %s\n", etype->name, gtype->name);
	fprintf(stdout, "easyRNG %g\n", ((double) (eend - estart))/CLOCKS_PER_SEC);
	fprintf(stdout, "GSL %g\n", ((double) (gend - gstart))/CLOCKS_PER_SEC);
}

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

	run_test(gsl_rng_mt19937, easy_rng_mt19937);
	run_test(gsl_rng_ranlux, easy_rng_ranlux24);
	run_test(gsl_rng_ranlux389, easy_rng_ranlux48);
	run_test(gsl_rng_minstd, easy_rng_minstd_rand0);
	run_test(gsl_rng_fishman20, easy_rng_minstd_rand);

	return 0;
}

max_stars_repo_name: tschoonj/easyRNG

max_issues_repo_name: tschoonj/easyRNG

max_forks_repo_name: tschoonj/easyRNG


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

#include <stdio.h>
#include <gsl/gsl_poly.h>

int
main ()
{
  int i;
  double a[6] = { -1, 0, 0, 0, 0, 1 };  /*  P(x) =  x^5 - 1  */
  double z[10];

  gsl_poly_complex_workspace * w = gsl_poly_complex_workspace_alloc (6) ;

  gsl_poly_complex_solve (a, 6, w, z) ;

  gsl_poly_complex_workspace_free (w) ;

  for (i = 0; i < 5 ; i++)
    {
      printf("z%d = %+.18f %+.18f\n", i, z[2*i], z[2*i+1]) ;
    }
}

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim


text:
/*
 * Copyright 2016 Maikel Nadolski
 *
 * 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 HIDDEN_MARKOV_MODEL_H_
#define HIDDEN_MARKOV_MODEL_H_

#include <Eigen/Dense>
#include <gsl_util.h>

#include "maikel/hmm/stochastical_conditions.h"

namespace maikel { namespace hmm {

  struct hmm_errors: public std::runtime_error {
    hmm_errors(std::string s): std::runtime_error(s) {}
  };

  struct dimensions_not_consistent: public hmm_errors {
    dimensions_not_consistent(const std::string& a): hmm_errors(a) {}
  };

  template <class T>
    class hidden_markov_model
      {
        public:
          using matrix     = typename Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>;
          using row_vector = typename Eigen::Matrix<T, 1, Eigen::Dynamic>;
          using size_type  = typename matrix::Index;

          struct arguments_not_probability_arrays: public hmm_errors {
              matrix A;
              matrix B;
              row_vector pi;
              arguments_not_probability_arrays(
                  const matrix& A_,
                  const matrix& B_,
                  const row_vector& pi_,
                  const std::string& a): hmm_errors(a), A{A_}, B{B_}, pi{pi_} {}
          };

          hidden_markov_model(
              const matrix&     transition_matrix,
              const matrix&     symbol_matrix,
              const row_vector& initial_dist )
           : A { transition_matrix },
             B { symbol_matrix     },
             pi{ initial_dist      },
             num_states  { B.rows() },
             num_symbols { B.cols() }
          {
            if (!rows_are_probability_arrays(A) || !rows_are_probability_arrays(B)
                || !is_probability_array(pi.array()))
              throw arguments_not_probability_arrays
                { A, B, pi, "Some inputs in constructor do not have the stochastical property." };
            if (A.rows() != A.cols() || A.rows() != B.rows() || A.rows() != pi.cols())
              throw dimensions_not_consistent
                { "Dimensions of input matrices are not consistent with each other." };
          }

          inline size_type states() const noexcept  { return num_states;  }
          inline size_type symbols() const noexcept { return num_symbols; }

          inline const matrix&     transition_matrix()    const noexcept { return A; }
          inline const matrix&     symbol_probabilities() const noexcept { return B; }
          inline const row_vector& initial_distribution() const noexcept { return pi; }

        private:
          matrix A;
          matrix B;
          row_vector pi;
          size_type num_states;
          size_type num_symbols;
      };

} // namespace hmm
} // namespace maikel

#endif /* HIDDEN_MARKOV_MODEL_H_ */

max_stars_repo_name: maikel/hidden-markov-model

max_issues_repo_name: maikel/Hidden-Markov-Model

max_forks_repo_name: maikel/Hidden-Markov-Model


text:
#include <gsl/gsl_interp.h>
#include <gsl/gsl_bspline.h>
#include <gsl/gsl_version.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_multifit.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_statistics_double.h>


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

#include <config.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_bst.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_test.h>

enum array_order
{
  ORD_RANDOM = 0,               /* random order */
  ORD_ASCENDING,                /* ascending order */
  ORD_DESCENDING,               /* descending order */
  ORD_BALANCED,                 /* balanced tree order */
  ORD_ZIGZAG,                   /* zig-zag order */
  ORD_ASCENDING_SHIFTED,        /* ascending from middle, then beginning */
  ORD_END_NODUP,                /* end of no-duplicate ordering */

  ORD_RANDOM_DUP                /* random order with duplicates */
};

/* fill array[] with random integers in [lower,upper] with duplicates allowed */
static void
random_integers(const size_t n, const int lower, const int upper,
                int array[], gsl_rng * r)
{
  size_t i;

  for (i = 0; i < n; ++i)
    array[i] = (int) ((upper - lower) * gsl_rng_uniform(r) + lower);
}

/* fills array[] with a random permutation of the integers between 0 and n - 1 */
static void
random_permuted_integers (const size_t n, int array[], gsl_rng * r)
{
  size_t i;

  for (i = 0; i < n; i++)
    array[i] = i;

  for (i = 0; i < n; i++)
    {
      size_t j = i + (unsigned) (gsl_rng_uniform(r) * (n - i));
      int t = array[j];
      array[j] = array[i];
      array[i] = t;
    }
}

static int
compare_ints(const void *pa, const void *pb, void *params)
{
  const int *a = pa;
  const int *b = pb;

  (void) params;

  return (*a < *b) ? -1 : (*a > *b);
}

/* Generates a list of integers that produce a balanced tree when
   inserted in order into a binary tree in the usual way.
   |min| and |max| inclusively bound the values to be inserted.
   Output is deposited starting at |*array|. */
static void
gen_balanced_tree (const int min, const int max, int **array)
{
  int i;

  if (min > max)
    return;

  i = (min + max + 1) / 2;
  *(*array)++ = i;
  gen_balanced_tree (min, i - 1, array);
  gen_balanced_tree (i + 1, max, array);
}

/* generates a permutation of the integers |0| to |n - 1| */
static void
gen_int_array (const size_t n, const enum array_order order, int array[], gsl_rng * r)
{
  size_t i;

  switch (order)
    {
      case ORD_RANDOM:
        random_permuted_integers (n, array, r);
        break;

      case ORD_ASCENDING:
        for (i = 0; i < n; i++)
          array[i] = i;
        break;

      case ORD_DESCENDING:
        for (i = 0; i < n; i++)
          array[i] = n - i - 1;
        break;

      case ORD_BALANCED:
        gen_balanced_tree (0, n - 1, &array);
        break;

      case ORD_ZIGZAG:
        for (i = 0; i < n; i++)
          {
            if (i % 2 == 0)
              array[i] = i / 2;
            else
              array[i] = n - i / 2 - 1;
          }
        break;

      case ORD_ASCENDING_SHIFTED:
        for (i = 0; i < n; i++)
          {
            array[i] = i + n / 2;
            if ((size_t) array[i] >= n)
              array[i] -= n;
          }
        break;

      case ORD_RANDOM_DUP:
        random_integers(n, -10, 10, array, r);
        break;

      default:
        assert (0);
    }
}

static void
check_traverser(const size_t n, const enum array_order order, gsl_bst_trav * trav, int data,
                const char *desc, const gsl_bst_workspace * w)
{
  int *prev, *cur, *next;

  prev = gsl_bst_trav_prev(trav);
  if (prev != NULL)
    {
      gsl_test(*prev > data, "bst %s[n=%zu,order=%d] %s traverser ahead of %d, but should be ahead of %d",
               gsl_bst_name(w), n, order, desc, *prev, data);
    }
  gsl_bst_trav_next(trav);

  cur = gsl_bst_trav_cur(trav);
  gsl_test(*cur != data, "bst %s[n=%zu,order=%d] %s traverser at %d, but should be at %d",
           gsl_bst_name(w), n, order, desc, *cur, data);

  next = gsl_bst_trav_next(trav);
  if (next != NULL)
    {
      gsl_test(*next < data, "bst %s[n=%zu,order=%d] %s traverser behind %d, but should be behind %d",
               gsl_bst_name(w), n, order, desc, *next, data);
    }
  gsl_bst_trav_prev(trav);
}

static void
test_bst_int(const size_t n, const gsl_bst_type * T, const enum array_order order, gsl_rng * r)
{
  int *data = malloc(n * sizeof(int));
  int *data_delete = malloc(n * sizeof(int));
  int *sorted_data = malloc(n * sizeof(int));
  gsl_bst_workspace * w = gsl_bst_alloc(T, NULL, compare_ints, NULL);
  gsl_bst_trav trav;
  int *p;
  int i;
  size_t nodes;

  /* generate data to be inserted in tree */
  gen_int_array(n, order, data, r);

  for (i = 0; i < (int) n; ++i)
    sorted_data[i] = data[i];

  gsl_sort_int(sorted_data, 1, n);

  if (order != ORD_RANDOM_DUP)
    {
      /* generate random order to delete data from tree */
      gen_int_array(n, ORD_RANDOM, data_delete, r);
    }
  else
    {
      for (i = 0; i < (int) n; ++i)
        data_delete[i] = sorted_data[i];
    }

  /* insert data */
  for (i = 0; i < (int) n; ++i)
    {
      p = gsl_bst_insert(&data[i], w);
      gsl_test(p != NULL, "bst_int %s[n=%zu,order=%d] insert i=%d", gsl_bst_name(w), n, order, i);
    }

  if (order != ORD_RANDOM_DUP)
    {
      nodes = gsl_bst_nodes(w);
      gsl_test(nodes != n, "bst_int %s[n=%zu,order=%d] after insertion count = %zu/%zu",
               gsl_bst_name(w), n, order, nodes, n);
    }

  /* test data was inserted and can be found */
  for (i = 0; i < (int) n; ++i)
    {
      p = gsl_bst_find(&data[i], w);
      gsl_test(*p != data[i], "bst_int %s[n=%zu,order=%d] find [%d,%d]",
               gsl_bst_name(w), n, order, *p, data[i]);

      p = gsl_bst_trav_find(&data[i], &trav, w);
      gsl_test(p == NULL, "bst_int %s[n=%zu,order=%d] trav_find unable to find item %d",
               gsl_bst_name(w), n, order, data[i]);

      check_traverser(n, order, &trav, data[i], "post-insertion", w);
    }

  /* traverse tree in-order */
  p = gsl_bst_trav_first(&trav, w);
  i = 0;
  while (p != NULL)
    {
      int *q = gsl_bst_trav_cur(&trav);

      gsl_test(*p != sorted_data[i], "bst_int %s[n=%zu,order=%d] traverse i=%d [%d,%d]",
               gsl_bst_name(w), n, order, i, *p, sorted_data[i]);

      gsl_test(*p != *q, "bst_int %s[n=%zu,order=%d] traverse cur i=%d [%d,%d]",
               gsl_bst_name(w), n, order, i, *p, *q);

      p = gsl_bst_trav_next(&trav);
      ++i;
    }

  gsl_test(i != (int) n, "bst_int %s[n=%zu,order=%d] traverse number=%d",
           gsl_bst_name(w), n, order, i);

  /* traverse tree in reverse order */
  p = gsl_bst_trav_last(&trav, w);
  i = n - 1;
  while (p != NULL)
    {
      int *q = gsl_bst_trav_cur(&trav);

      gsl_test(*p != sorted_data[i], "bst_int %s[n=%zu,order=%d] traverse reverse i=%d [%d,%d]",
               gsl_bst_name(w), n, order, i, *p, sorted_data[i]);

      gsl_test(*p != *q, "bst_int %s[n=%zu,order=%d] traverse reverse cur i=%d [%d,%d]",
               gsl_bst_name(w), n, order, i, *p, *q);

      p = gsl_bst_trav_prev(&trav);
      --i;
    }

  gsl_test(i != -1, "bst_int %s[n=%zu,order=%d] traverse reverse number=%d",
           gsl_bst_name(w), n, order, i);

  /* test traversal during tree modifications */
  for (i = 0; i < (int) n; ++i)
    {
      gsl_bst_trav x, y, z;

      gsl_bst_trav_find(&data[i], &x, w);
      check_traverser(n, order, &x, data[i], "pre-deletion", w);

      if (data[i] == data_delete[i])
        continue;

      p = gsl_bst_remove(&data_delete[i], w);
      gsl_test(*p != data_delete[i], "bst_int %s[n=%zu,order=%d] remove i=%d [%d,%d]",
               gsl_bst_name(w), n, order, i, *p, data_delete[i]);

      p = gsl_bst_trav_copy(&y, &x);
      gsl_test(*p != data[i], "bst_int %s[n=%zu,order=%d] copy i=%d [%d,%d]",
               gsl_bst_name(w), n, order, i, *p, data[i]);

      /* re-insert item */
      p = gsl_bst_trav_insert(&data_delete[i], &z, w);

      check_traverser(n, order, &x, data[i], "post-deletion", w);
      check_traverser(n, order, &y, data[i], "copied", w);
      check_traverser(n, order, &z, data_delete[i], "insertion", w);

#if 0
      /* delete again */
      gsl_bst_remove(&data[i], w);
#endif
    }

  /* emmpty tree */
  gsl_bst_empty(w);

  nodes = gsl_bst_nodes(w);
  gsl_test(nodes != 0, "bst_int %s[n=%zu,order=%d] empty count = %zu",
           gsl_bst_name(w), n, order, nodes);

  gsl_bst_free(w);
  free(data);
  free(data_delete);
  free(sorted_data);
}

static void
test_bst(const gsl_bst_type * T, gsl_rng * r)
{
  enum array_order order;

  for (order = 0; order < ORD_END_NODUP; ++order)
    {
      test_bst_int(50, T, order, r);
      test_bst_int(100, T, order, r);
      test_bst_int(500, T, order, r);
    }
}

int
main(void)
{
  gsl_rng * r = gsl_rng_alloc(gsl_rng_default);

  test_bst(gsl_bst_avl, r);
  test_bst(gsl_bst_rb, r);

  gsl_rng_free(r);

  exit (gsl_test_summary());
}

max_forks_repo_name: karanbirsandhu/nu-sense