text:
#pragma once

#include <GLES2/gl2.h>
#include <GLES2/gl2ext.h>
#include <GLES3/gl3.h>
#include <GLES3/gl3ext.h>
#include <GLES3/gl3platform.h>
#include <EGL/egl.h>
#include <EGL/eglext.h>
#include <gsl/gsl>

namespace android::OpenGLHelpers
{
    constexpr GLint GetTextureUnit(GLenum texture)
    {
        return texture - GL_TEXTURE0;
    }

    GLuint CreateShaderProgram(const char* vertShaderSource, const char* fragShaderSource);

    namespace GLTransactions
    {
        inline auto SetCapability(GLenum capability, bool isEnabled)
        {
            const auto setCapability{ [capability](bool isEnabled)
            {
                if (isEnabled)
                {
                    glEnable(capability);
                }
                else
                {
                    glDisable(capability);
                }
            }};

            const auto wasEnabled{ glIsEnabled(capability) };
            setCapability(isEnabled);
            return gsl::finally([wasEnabled, setCapability]() { setCapability(wasEnabled); });
        }

        inline auto BindFrameBuffer(GLuint frameBufferId)
        {
            GLint previousFrameBufferId;
            glGetIntegerv(GL_FRAMEBUFFER_BINDING, &previousFrameBufferId);
            glBindFramebuffer(GL_FRAMEBUFFER, frameBufferId);
            return gsl::finally([previousFrameBufferId]() { glBindFramebuffer(GL_FRAMEBUFFER, static_cast<GLuint>(previousFrameBufferId)); });
        }

        inline auto DepthMask(GLboolean depthMask)
        {
            GLboolean previousDepthMask;
            glGetBooleanv(GL_DEPTH_WRITEMASK, &previousDepthMask);
            glDepthMask(depthMask);
            return gsl::finally([previousDepthMask]() { glDepthMask(previousDepthMask); });
        }

        inline auto BindSampler(GLenum unit, GLuint id)
        {
            glActiveTexture(unit);
            GLint previousId;
            glGetIntegerv(GL_SAMPLER_BINDING, &previousId);
            glBindSampler(unit - GL_TEXTURE0, id);
            return gsl::finally([unit, id{ previousId }]() { glActiveTexture(unit); glBindSampler(unit - GL_TEXTURE0, id); });
        }

        inline auto MakeCurrent(EGLDisplay display, EGLSurface drawSurface, EGLSurface readSurface, EGLContext context)
        {
            EGLDisplay previousDisplay{ eglGetDisplay(EGL_DEFAULT_DISPLAY) };
            EGLSurface previousDrawSurface{ eglGetCurrentSurface(EGL_DRAW) };
            EGLSurface previousReadSurface{ eglGetCurrentSurface(EGL_READ) };
            EGLContext previousContext{ eglGetCurrentContext() };
            eglMakeCurrent(display, drawSurface, readSurface, context);
            return gsl::finally([previousDisplay, previousDrawSurface, previousReadSurface, previousContext]() { eglMakeCurrent(previousDisplay, previousDrawSurface, previousReadSurface, previousContext); });
        }
    }
}


text:
/***************************************************************************
    File                 : muParserScript.h
	Project              : QtiPlot
    --------------------------------------------------------------------

    Copyright            : (C) 2006 by Ion Vasilief, Knut Franke
    Email (use @ for *)  : ion_vasilief*yahoo.fr, knut.franke*gmx.de
    Description          : Evaluate mathematical expressions using muParser

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

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

#include "ScriptingEnv.h"
#include "Script.h"

#include <muParser.h>
#include "math.h"
#include <gsl/gsl_sf.h>
#include <q3asciidict.h>

//! TODO
class muParserScript: public Script
{
  Q_OBJECT

  public:
    muParserScript(ScriptingEnv *env, const QString &code, QObject *context=0, const QString &name="<input>", bool checkMultilineCode = true);

  public slots:
    bool compile(bool asFunction=true);
    QVariant eval();
    double evalSingleLine();
    QString evalSingleLineToString(const QLocale& locale, char f, int prec);
    bool exec();
    bool setQObject(QObject *val, const char *name);
    bool setInt(int val, const char* name);
    bool setDouble(double val, const char* name);
    double* defineVariable(const char *name, double val = 0.0);
    int codeLines(){return muCode.size();};

  private:
    double col(const QString &arg);
    double tablecol(const QString &arg);
    double cell(int row, int col);
    double tableCell(int col, int row);
    double *addVariable(const char *name);
    double *addVariableR(const char *name);
    static double *mu_addVariableR(const char *name) { return current->addVariableR(name); }
    static double mu_col(const char *arg) { return current->col(arg); }
    static double mu_cell(double row, double col) { return current->cell(qRound(row), qRound(col)); }
    static double mu_tableCell(double col, double row) { return current->tableCell(qRound(col), qRound(row)); }
    static double mu_tablecol(const char *arg) { return current->tablecol(arg); }
    static double *mu_addVariable(const char *name, void *){ return current->addVariable(name); }
    static double *mu_addVariableR(const char *name, void *) { return current->addVariableR(name); }
    static QString compileColArg(const QString& in);

    mu::Parser parser, rparser;
    Q3AsciiDict<double> variables, rvariables;
    QStringList muCode;
	//! Flag telling is the parser should warn users on multiline code input
	bool d_warn_multiline_code;

  public:
    static muParserScript *current;
};

#endif

ext: h

max_stars_repo_name: hoehnp/SpaceDesignTool

max_issues_repo_name: hoehnp/SpaceDesignTool

max_forks_repo_name: hoehnp/SpaceDesignTool


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

#ifdef _MSC_VER
# define inline __inline
#endif

/*
 * General Cholesky Delete.
 * Remove an element from the cholesky factorization
 * m = columns
 * n = rows
 *
 * TODO: put transpose as an option
 */
static inline void cholesky_delete_dbl(int m, int n, double *L, int go_out)
{
    double c, s;

    /* delete row go_out */
    double *L1 = L + (go_out * m);
    int i;
    for (i = go_out; i < n - 1; ++i) {
        cblas_dcopy (i + 2, L1 + m , 1, L1,  1);
        L1 += m;
    }

    L1 = L + (go_out * m);
    for (i=go_out; i < n - 1; ++i) {

        cblas_drotg(L1 + i, L1 + i + 1, &c, &s);
        if (L1[i] < 0) {
            /* Diagonals cannot be negative */
            L1[i] = fabs(L1[i]);
            c = -c;
            s = -s;
        }
        L1[i+1] = 0.; /* just for cleanup */
        L1 += m;

        cblas_drot(n - (i + 2), L1 + i, m, L1 + i + 1,
                   m, c, s);
    }
}

static inline void cholesky_delete_flt(int m, int n, float *L, int go_out)
{
    float c, s;

    /* delete row go_out */
    float *L1 = L + (go_out * m);
    int i;
    for (i = go_out; i < n - 1; ++i) {
        cblas_scopy (i + 2, L1 + m , 1, L1,  1);
        L1 += m;
    }

    L1 = L + (go_out * m);
    for (i=go_out; i < n - 1; ++i) {

        cblas_srotg(L1 + i, L1 + i + 1, &c, &s);
        if (L1[i] < 0) {
            /* Diagonals cannot be negative */
            L1[i] = fabsf(L1[i]);
            c = -c;
            s = -s;
        }
        L1[i+1] = 0.; /* just for cleanup */
        L1 += m;

        cblas_srot(n - (i + 2), L1 + i, m, L1 + i + 1,
                   m, c, s);
    }
}


text:
/**
 *
 * @file core_zgetrf_reclap.c
 *
 *  PLASMA core_blas kernel
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.8.0
 * @author Hatem Ltaief
 * @author Mathieu Faverge
 * @author Piotr Luszczek
 * @date 2009-11-15
 *
 * @precisions normal z -> c d s
 *
 **/

#include <math.h>
#include <cblas.h>
#include <lapacke.h>
#include "common.h"

struct CORE_zgetrf_data_s {
    volatile PLASMA_Complex64_t *CORE_zamax;
    volatile int                *CORE_zstep;
};

static inline void
CORE_zgetrf_reclap_update(CORE_zgetrf_data_t *data,
                          int M, int column, int n1, int n2,
                          PLASMA_Complex64_t *A, int LDA, int *IPIV,
                          int thidx, int thcnt);
static inline void
CORE_zgetrf_reclap_rec(CORE_zgetrf_data_t *data,
                       int M, int N,
                       PLASMA_Complex64_t *A, int LDA,
                       int *IPIV, int *info,
                       int thidx, int thcnt, int column);

/***************************************************************************//**
 *
 * @ingroup dplasma_cores_complex64
 *
 *  CORE_zgetrf_reclap computes a LU factorization of a general M-by-N
 *  matrix A stored in CCRB layout using partial pivoting with row
 *  interchanges.
 *
 *  The factorization has the form
 *
 *    A = P * L * U
 *
 *  where P is a permutation matrix, L is lower triangular with unit
 *  diagonal elements (lower trapezoidal if m > n), and U is upper
 *  triangular (upper trapezoidal if m < n).
 *
 *  This is the recursive version of the algorithm applied on column
 *  major layout.
 *
 *  WARNINGS:
 *     - The function CORE_zgetrf_reclap_init has to be called prior
 *     to any call to this function.
 *     - You cannot call this kernel on different matrices at the same
 *     time.
 *     - The matrix A cannot be more than one tile wide.
 *     - The number of threads calling this function has to be excatly
 *     the number defined by info[2] with each one of them a different
 *     index between 0 included and info[2] excluded.
 *
 *******************************************************************************
 *
 *  @param[in] data
 *          Common data structure to all threads initialized by
 *          CORE_zgetrf_reclap_init() that contains information for thread
 *          synchronization and maximum reductions. All threads working on a
 *          given matrix A must provide the same data structure.
 *
 *  @param[in] M
 *          The number of rows of the matrix A.  M >= 0.
 *
 *  @param[in] N
 *          The number of columns of the matrix A.  N >= 0.
 *
 *  @param[in,out] A
 *          On entry, the M-by-N matrix to be factorized.
 *          On exit, the factors L and U from the factorization
 *          A = P*L*U; the unit diagonal elements of L are not stored.
 *
 *  @param[in] LDA
 *          The leading dimension of the array A.  LDA >= max(1,M).
 *
 *  @param[out] IPIV
 *          The pivot indices; for 1 <= i <= min(M,N), row i of the
 *          matrix was interchanged with row IPIV(i).
 *          1 <= IPIV[i] <= M.
 *
 *  @param[in,out] info
 *          Array of 3 integers
 *          - info[0], see returned value
 *          - info[1], is the thread index 0 <= info[0] < info[2]
 *          - info[2], on entry is the number of threads trying to
 *                     participate to the factorization,
 *                     on exit is the real number of threads used to
 *                     perform the factorization.
 *    Info[2] threads, and exactly info[2], have to call this function
 *    to avoid dead lock.
 *
 *******************************************************************************
 *
 * @return
 *          \retval PLASMA_SUCCESS successful exit
 *          \retval -k, the k-th argument had an illegal value
 *          \retval k if U(k,k) is exactly zero. The factorization
 *                  has been completed, but the factor U is exactly
 *                  singular, and division by zero will occur if it is used
 *                  to solve a system of equations.
 *
 */
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_zgetrf_reclap = PCORE_zgetrf_reclap
#define CORE_zgetrf_reclap PCORE_zgetrf_reclap
#endif
int CORE_zgetrf_reclap(CORE_zgetrf_data_t *data,
                       int M, int N,
                       PLASMA_Complex64_t *A, int LDA,
                       int *IPIV, int *info)
{
    int thidx = info[1];
    int divmn = M / N;
    int thcnt = min( info[2], ((divmn == 0) && (M != 0))? 1 : divmn );
    int minMN = min(M, N);

    info[0] = 0;
    info[2] = thcnt;

    if( M < 0 ) {
        coreblas_error(1, "illegal value of M");
        return -1;
    }
    if( N < 0 ) {
        coreblas_error(2, "illegal value of N");
        return -2;
    }
    if( LDA < max(1, M) ) {
        coreblas_error(5, "illegal value of LDA");
        return -5;
    }

    /*
     * Quick return
     */
    if ( (M == 0) || (N == 0) || (thidx >= thcnt) ){
      return PLASMA_SUCCESS;
    }

    *info = 0;
    CORE_zgetrf_reclap_rec( data, M, minMN, A, LDA, IPIV, info,
                            thidx, thcnt, 0 );

    if ( N > minMN ) {
        CORE_zgetrf_reclap_update(data, M, 0, minMN, N-minMN,
                                  A, LDA, IPIV,
                                  thidx, thcnt);
    }

    return info[0];
}

/*******************************************************************
 *   Additional routines
 */
static double sfmin = 0.;

CORE_zgetrf_data_t *
CORE_zgetrf_reclap_init(int nbthrd)
{
    int i;
    CORE_zgetrf_data_t *data;

    data = (CORE_zgetrf_data_t*)malloc( nbthrd * (sizeof(PLASMA_Complex64_t)+sizeof(int))
                                        + 2 * sizeof(void*) );
    data->CORE_zamax = (PLASMA_Complex64_t*)((char*)data + 2 * sizeof(void*));
    data->CORE_zstep = (int*)((char*)data + 2 * sizeof(void*) + nbthrd * sizeof(PLASMA_Complex64_t));

    for (i = 0; i < nbthrd; ++i) {
        data->CORE_zamax[i] = 0.;
        data->CORE_zstep[i] = -1;
    }
    if ( sfmin == 0. ) {
        sfmin = LAPACKE_dlamch_work('S');
    }
    return data;
}

static inline void
psplit(int n, int pidx, int pcnt, int *poff_p, int *psiz_p)
{
    int q = n / pcnt, r = n % pcnt;

    if (pidx < r) {
        q++;
        *psiz_p = q;
        *poff_p = pidx * q;
    } else {
        *psiz_p = q;
        *poff_p = r * (q + 1) + (pidx - r) * q;
    }
}

static inline void
CORE_zamax1_thread(CORE_zgetrf_data_t *data,
                   PLASMA_Complex64_t localamx,
                   int thidx, int thcnt, int *thwinner,
                   PLASMA_Complex64_t *globalamx,
                   int pividx, int *ipiv)
{
    volatile PLASMA_Complex64_t *CORE_zamax = data->CORE_zamax;
    volatile int                *CORE_zstep = data->CORE_zstep;

    if (thidx == 0) {
        int i, j = 0;
        PLASMA_Complex64_t curval = localamx, tmp;
        double curamx = cabs(localamx);

        /* make sure everybody filled in their value */
        for (i = 1; i < thcnt; ++i) {
            while (CORE_zstep[i] == -1) { /* wait for thread i to store its value */
            }
        }

        /* better not fuse the loop above and below to make sure data is sync'd */

        for (i = 1; i < thcnt; ++i) {
            tmp = CORE_zamax[i];
            if (cabs(tmp) > curamx) {
                curamx = cabs(tmp);
                curval = tmp;
                j = i;
            }
        }

        if (0 == j)
            ipiv[0] = pividx;

        /* make sure everybody knows the amax value */
        for (i = 1; i < thcnt; ++i)
            CORE_zamax[ i ] = curval;

        CORE_zstep[0] = -j - 2; /* set the index of the winning thread */

        *thwinner = j;
        *globalamx = curval;

        for (i = 1; i < thcnt; ++i)
            CORE_zstep[i] = -3;

        /* make sure everybody read the max value */
        for (i = 1; i < thcnt; ++i) {
            while (CORE_zstep[i] != -1) {
            }
        }

        CORE_zstep[0] = -1;
    } else {
        CORE_zamax[thidx] = localamx;
        CORE_zstep[thidx] = -2;  /* announce to thread 0 that local amax was stored */
        while (CORE_zstep[0] == -1) { /* wait for thread 0 to finish calculating the global amax */
        }
        while (CORE_zstep[thidx] != -3) { /* wait for thread 0 to store amax */
        }
        *thwinner  = -CORE_zstep[0] - 2;
        *globalamx = CORE_zamax[thidx];     /* read the amax from the location adjacent to the one in the above loop */
        CORE_zstep[thidx] = -1;  /* signal thread 0 that this thread is done reading */

        if (thidx == *thwinner)
            ipiv[0] = pividx;

        while (CORE_zstep[0] != -1) { /* wait for thread 0 to finish */
        }
    }
}

static inline void
CORE_zbarrier_thread(CORE_zgetrf_data_t *data,
                     int thidx, int thcnt)
{
    int idum1, idum2;
    PLASMA_Complex64_t ddum2 = 0.;
    /* it's probably faster to implement a dedicated barrier */
    CORE_zamax1_thread( data, 1.0, thidx, thcnt, &idum1, &ddum2, 0, &idum2 );
}

static inline void
CORE_zlaswap1(int ncol, PLASMA_Complex64_t *a, int lda,
              int idxStart, int idxMax, const int *piv)
{
    int i, j;
    PLASMA_Complex64_t tmp;

    for (j = 0; j < ncol; ++j) {
        for (i = idxStart; i < idxMax; ++i) {
            tmp = a[j*lda + piv[i] - 1];
            a[j*lda + piv[i] - 1] = a[i + j*lda];
            a[i + j*lda] = tmp;
        }
    }
}

static inline void
CORE_zgetrf_reclap_update(CORE_zgetrf_data_t *data,
                          int M, int column, int n1, int n2,
                          PLASMA_Complex64_t *A, int LDA, int *IPIV,
                          int thidx, int thcnt)
{
    static PLASMA_Complex64_t posone =  1.0;
    static PLASMA_Complex64_t negone = -1.0;
    PLASMA_Complex64_t *Atop  = A    + column*LDA;
    PLASMA_Complex64_t *Atop2 = Atop + n1    *LDA;
    int coff, ccnt, lm, loff;

    CORE_zbarrier_thread( data, thidx, thcnt );

    psplit( n2, thidx, thcnt, &coff, &ccnt );

    if (ccnt > 0) {
        CORE_zlaswap1( ccnt, Atop2 + coff*LDA, LDA, column, n1 + column, IPIV ); /* swap to the right */

        cblas_ztrsm( CblasColMajor, CblasLeft, CblasLower, CblasNoTrans, CblasUnit,
                     n1, ccnt, CBLAS_SADDR(posone), Atop + column, LDA, Atop2 + coff*LDA + column, LDA );
    }

    /* __sync_synchronize(); */ /* hopefully we will not need memory fences */

    /* need to wait for pivoting and triangular solve to finish */
    CORE_zbarrier_thread( data, thidx, thcnt );

    psplit( M, thidx, thcnt, &loff, &lm );
    if (thidx == 0) {
        loff = column + n1;
        lm  -= column + n1;
    };

    cblas_zgemm( CblasColMajor, CblasNoTrans, CblasNoTrans, lm, n2, n1,
                 CBLAS_SADDR(negone), Atop+loff, LDA, Atop2 + column, LDA, CBLAS_SADDR(posone), Atop2+loff, LDA );
}

static void
CORE_zgetrf_reclap_rec(CORE_zgetrf_data_t *data, int M, int N,
                       PLASMA_Complex64_t *A, int LDA,
                       int *IPIV, int *info,
                       int thidx, int thcnt, int column)
{
    int jp, n1, n2, lm, loff;
    PLASMA_Complex64_t tmp1, tmp2, tmp3;
    PLASMA_Complex64_t *Atop = A + column*LDA;

    /* Assumption: N = min( M, N ); */
    if (N > 1) {
        int coff, ccnt;

        n1 = N / 2;
        n2 = N - n1;

        CORE_zgetrf_reclap_rec( data, M, n1, A, LDA, IPIV, info,
                                thidx, thcnt, column );
        if ( *info != 0 )
            return;

        CORE_zgetrf_reclap_update(data, M, column, n1, n2,
                                  A, LDA, IPIV,
                                  thidx, thcnt);

        CORE_zgetrf_reclap_rec( data, M, n2, A, LDA, IPIV, info,
                                thidx, thcnt, column + n1 );
        if ( *info != 0 )
            return;

        psplit( n1, thidx, thcnt, &coff, &ccnt );

        if (ccnt > 0) {
            CORE_zlaswap1( ccnt, Atop+coff*LDA, LDA, n1 + column, N + column, IPIV ); /* swap to the left */
        }

    } else {
        int thrd;

        CORE_zbarrier_thread( data, thidx, thcnt );

        psplit( M, thidx, thcnt, &loff, &lm );

        if (thidx == 0) {
            loff = column;
            lm -= column;
        }

        tmp2 = Atop[column]; /* all threads read the pivot element in case they need it */

        jp = cblas_izamax( lm, Atop + loff, 1 );
        tmp1 = Atop[loff + jp];

        CORE_zamax1_thread( data, tmp1, thidx, thcnt, &thrd,
                            &tmp3, loff + jp + 1, IPIV + column );

        Atop[column] = tmp3; /* all threads set the pivot element: no need for synchronization */

        if ( tmp3 != 0.0 ) {
            if ( cabs(tmp3) >= sfmin ) {
                PLASMA_Complex64_t tmp = (PLASMA_Complex64_t)1.0 / tmp3;
                n1 = (thidx == 0) ? 1 : 0;
                cblas_zscal( lm - n1, CBLAS_SADDR(tmp), Atop + loff + n1, 1 );
            } else {
                int i;
                PLASMA_Complex64_t *Atop2;
                n1 = (thidx == 0) ? 1 : 0;
                Atop2 = Atop + loff + n1;

                for( i=0; i < lm-n1; i++, Atop2++)
                    *Atop2 = *Atop2 / tmp3;
            }

            if (thrd == thidx) { /* the thread that owns the best pivot */
              if (loff + jp != column) /* if there is a need to exchange the pivot */
                Atop[loff + jp] = tmp2 / tmp3;
            }

        } else {
            *info = column + 1;
            return;
        }

        CORE_zbarrier_thread( data, thidx, thcnt );
    }
}

ext: c

max_stars_repo_name: therault/dplasma

max_issues_repo_name: therault/dplasma

max_forks_repo_name: therault/dplasma


text:
#ifndef COMPUTE_SIGNAL_H
#define COMPUTE_SIGNAL_H

#include <gsl/gsl_rng.h>
#include <gsl/gsl_histogram.h>

void compute_signal(double * signal, const size_t length, const double delta,
    gsl_rng * r, gsl_histogram * hist);

#endif


text:
// Data Structure and Routines for peak fiting

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_multifit_nlin.h>

#include <errno.h>

#include "grid.h"
#include "grid_operations.h"
#include "minmax.h"
#include "point.h"

#ifndef _FITTOFUNCTION_H_
#define _FITTOFUNCTION_H_

typedef  struct {
  size_t n;
  double * y;
} ObservedValues ;

void fitToFunctionLorentz(Grid *image, double *fitx, double *fity,
		   double *background, double *intens,
		   double *widthx, double *widthy, double *tilt,double *chisq);

void fitToFunctionGauss(Grid *image, double *fitx, double *fity,
		   double *background, double *intens,
		   double *widthx, double *widthy, double *tilt,double *chisq);

#endif

max_stars_repo_name: carterbox/cold

max_issues_repo_name: carterbox/cold

max_forks_repo_name: carterbox/cold


text:
/*
    learn.h
*/
#ifndef LEARN_H
#define LEARN_H
#include <stdlib.h>
#include <gsl/gsl_rng.h>
#include "feature.h"

//#define RANDOM ((double)rand()/(double)RAND_MAX)

int sampling_multinomial(gsl_rng *r, double *p, double *cum_sum_p, int len_p);

extern void mvslda_learn(document *data,
        double **resp,
        double *alpha,
        double beta,
        double nu2, double sigma2,
        int nclass, int nlex, int dlenmax, int nresp,
        int maxiter,
        double **phi, double **theta, double **eta,
        int **n_mz, int **n_zw,
        FILE *likp, FILE *hyperp,
        unsigned long int random_seed);

#endif


text:
#include <stdio.h>
#include <math.h>
#include <unistd.h> /* for getopt */
#include <string.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv2.h>

#include "defaults.h"
#include "h2ocalc.h"
//#include"h2ocalc_test.h"

/* gcc h2ocalc.c -o h2ocalc -lgsl -lgslcblas -lm -Wall  */

/* struct ODE_param{ */
/*   int a[NEQ];  /\* dont worry let just run with global vars for now *\/ */
/*   float b[NEQ]; */
/* }; */

int nbprint(double x, double y[], int n) {
    int i;
    double sum = 0;
    printf("%6.4f   ", x);
    for (i = 0; i < n; i++) {
        printf("%6.2f ", y[i]);
        sum += y[i];
    }
    printf("  %6.2f\n",sum);
    return 0;
}

int nbeprint(double x, double y[], int n) {
    int i;
    double sum = 0;

    printf("%6.4e   ",x);
    for (i = 0; i < n; i++) {
        printf("%6.2e ", y[i]);
        sum += y[i];
    }
    printf("  %6.2e\n", sum);
    return 0;
}

int func (double t, const double y[],double f[], void *params) {
    // struct ODE_param *p = (struct ODE_params *) params;
    int i,j,k;
    double rate[NEQ];

    //printf("func\n");

    /* for each species ... */
    for (k = 0; k < NSPECIES; k++) {
        //printf("SPECIES %i\n",k);

        f[k] = 0;

        /* ...find all contributing and removing things */

        for (j = 0; j < NEQ; j++) {

            if (nmatrix[j][k] != 0) {

                rate[j] = nmatrix [j][k]; /* build the final equation */
                //printf("%3i ",nmatrix[j][k]);

                /* add generate I equations */
                rate[j] *= rconst[j]; /* get rate for reaction I_j */

                //printf(" * k%i ",j);
                for (i = 0; i < NSPECIES; i++) {
                    if (nmatrix[j][i] == -1 ) {
                        rate[j] *= y[i]; /* first order kinetics */
                        //printf("* y[%i] ",i);
                    }
                    if (nmatrix[j][i] == -2 ) {
                        rate[j] *= y[i]*y[i]; // second order kinetics
                        //printf("* y[%i]^2 ",i);
                    }
                }

                f[k] += rate[j];
                //printf(" + \n");
            }
        }

        //printf(" ----  k f[k]: %i     %e ---- \n",k,f[k]);

        //    nbeprint(t,y,NSPECIES);
        //    if (k ==8)
        //      exit(0);
    }

    //nbeprint(t,y,NSPECIES);
    //nbeprint(t,f,NSPECIES);

    /* add function for radiation, assumed in last bin */
    //f[NEQ-1] = 0;

    return GSL_SUCCESS;
}

int main(int argc, char **argv) {
    //  struct ODE_param *p = (struct ODE_params *) params;

    double dummy_param = 0;
    gsl_odeiv2_system sys = {func, NULL, NSPECIES, &dummy_param};
    gsl_odeiv2_driver *d = gsl_odeiv2_driver_alloc_y_new(&sys,
                                                         gsl_odeiv2_step_rk8pd,
                                                         1e-6,
                                                         1e-6,
                                                         0.0);
    int i,j,status;
    double t = NSTART;
    //  double t1 = NSTOP;
    //double foobar; /* trash parameter */
    double ti;
    double y[NSPECIES];

    double freq = FREQ;
    double doser = DOSER;
    double simtime = SIMTIME;
    double resol = RESOL;
    double period = 0;
    double dpulse = 0;
    double tick = 0;
    double ta = 1e-3; /* plotting clocks*/
    int pulse_left = 0;
    int pulse_counter = 0;
    signed long int tick_counter = 0;
    char *fname;
    int c;
    int flagp = 0, flagl = 0, flagc = 1; /* default: print each tick */

    /* parse options */
    opterr = 0;
    while ((c = getopt(argc, argv, "f:d:o:r:t:")) != -1) {
        switch(c) {
        case 'l':
            flagl = 1;
            flagp = 0;
            flagc = 0;
            break;
        case 'p':
            flagl = 0;
            flagp = 1;
            flagc = 0;
            break;
        case 'c':
            flagl = 0;
            flagp = 0;
            flagc = 1;
            break;
        case 'f':
            sscanf(optarg, "%lf", &freq);
            break;
        case 'd':
            sscanf(optarg, "%lf", &doser);
            break;
        case 'o':
            fname = optarg; /* not implemented */
            break;
        case 'r':
            sscanf(optarg, "%lf", &resol);
            break;
        case 't':
            sscanf(optarg, "%lf", &simtime);
            break;
        case '?':
            printf("Options:\n");
            printf("  -p output per pulse\n");
            printf("  -c output per tick\n");
            printf("  -l only end of output\n");
            printf("  -f frequency [Hz]\n");
            printf("  -d dose rate [Gy/min]\n");
            printf("  -o output file\n");
            printf("  -r tick resolution [sec]\n");
            printf("  -t simulation time [sec]\n");
            printf("\n");
            exit(0);
            break;
        default:
            printf ("?? no handle ?? %i\n", c);
            exit(-1);
        }
    }

    period = 1 / freq;
    dpulse = (doser / 60.0) / freq * EVJ; // in eV per pulse per liter

    if (period < resol) {
        printf(" *** Error: ");
        printf("Period must be larger than tick resolution.\n");
        exit(-1);
    }

    /* find ticksize closest to the requested one */
    tick = period/(round(period/ resol));

    /* number of pulses to be simulated */
    pulse_left = (int) ((simtime-RSTART) * freq);

    /* print a header */
    printf("# Frequency    : %.3e Hz \n", freq);
    printf("# Period       : %.3e sec\n", period);
    printf("#\n");
    printf("# Tick size    : %.3e sec\n", tick);
    printf("# Time for sim : %.3e sec\n", simtime);
    printf("#\n");
    printf("# Dose rate    : %.3e Gy/min\n", doser );
    printf("# Pulse size   : %.3e eV/l/pulse\n", dpulse );
    printf("# Pulse count  : %i pulses\n", pulse_left);
    printf("#\n");
    printf("# Total delivered dose for this simulation\n");
    printf("#              : %.6e Gy \n", pulse_left * dpulse / EVJ);
    printf("#\n");

    /* copy start conditions into working array */
    for (i = 0; i < NSPECIES; i++)
        y[i] = ystart[i];

    /* Print first line with starting conditions. */
    nbeprint(t, y, NSPECIES);

    /* loop over all pulses */
    while (pulse_left != 0) {

        /* check if we are at a pulse time step */
        if (t >= ((pulse_counter * period) + RSTART)) {
            if (flagp) /* print per pulse (pre pulse) */
                nbeprint(t, y, NSPECIES);

            printf("# Pulse! %i \n", pulse_left);
            for (j=0; j < NSPECIES; j++)
                y[j] += dpulse * gval[j] * 0.01 / NA;
            pulse_left--;
            pulse_counter++;
        }

        tick_counter++;
        ti = tick_counter * tick;

        status = gsl_odeiv2_driver_apply(d, &t, ti, y);

        if (status != GSL_SUCCESS) {
            printf("error, return value = %d\n", status);
            break;
        }

        if (flagc) /* print per tick */
            nbeprint(t, y, NSPECIES);

    } /* end pulses left iterator */

    if (flagl) /* print last */
        nbeprint(t, y, NSPECIES);

    gsl_odeiv2_driver_free(d);
    return 0;
}

max_stars_repo_name: nbassler/simwater

max_issues_repo_name: nbassler/simwater

max_forks_repo_name: nbassler/simwater


text:
#ifndef libCV_ellipse
#define libCV_ellipse

#include <objc/Object.h>
#include "cvtypes.h"
#include "util.h"
#include "Array.h"
#include "math.h"
#include <gsl/gsl_math.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_eigen.h>

@interface Ellipse : Object {

	f64 x_axis;
	f64 y_axis;
	f64 x_center;
	f64 y_center;
	f64 rotation;

	f64 general[6];

	f64 bounds[4];

}

+(id) newAtX:(f64)xc andY:(f64)yc withXAxis:(f64)xa andYAxis:(f64)ya rotatedBy:(f64)phi;
+(id) newWithA:(f64)a b:(f64)b c:(f64)c d:(f64)d e:(f64)e f:(f64)f;

+(id) newFromPoints:(ArrayP)points;

-(void) drawInArray:(ArrayP)array;

-(void) toGeneralConic;
-(void) toGeneralEllipse;
-(void) findBounds;

// Get and set general ellipse parameters

-(f64) x_axis;
-(f64) y_axis;
-(f64) major_axis;
-(f64) minor_axis;
-(f64) rotation;
-(f64) majorRotation;
-(f64) minorRotation;
-(f64) x_center;
-(f64) y_center;

-(void) setX_axis:(f64)newx;
-(void) setY_axis:(f64)newy;
-(void) setX_center:(f64)newx;
-(void) setY_center:(f64)newy;
-(void) setRotation:(f64)newr;
-(void) setX_axis:(f64)xa y_axis:(f64)ya x_center:(f64)xc y_center:(f64)yc rotation:(f64)r;

// Get and set general conic parameters

-(f64) A;
-(f64) B;
-(f64) C;
-(f64) D;
-(f64) E;
-(f64) F;

-(f64 *) general;

-(void) setA:(f64)A;
-(void) setB:(f64)B;
-(void) setC:(f64)C;
-(void) setD:(f64)D;
-(void) setE:(f64)E;
-(void) setF:(f64)F;
-(void) setGeneralA:(f64)A B:(f64)B C:(f64)C D:(f64)D E:(f64)E F:(f64)F;

-(void) printInfoTo:(FILE *)fp;

-(u08) isValid;
-(id) free;

@end

enum ellipseindices {
	AX,
	BX,
	CX,
	DX,
	EX,
	FX,
};

enum ellipsebounds {
	XHI,
	XLO,
	YHI,
	YLO,
};

typedef Ellipse * EllipseP;

void originCenteredLSQEllipse( f64 x_values[], f64 y_values[], u64 number_of_points, f64 ellipse[] );
void generateBoundEllipses( f64 e[], f64 b1[], f64 b2[], f64 treshold );
f64 ellipseCircumference( f64 e[] );
EllipseP ellipseRANSAC( ArrayP edges, f64 fit_treshold, f64 min_percent_inliers, f64 certainty_probability, u64 max_iterations );

#endif


text:
// Copyright (c) 2021 Stig Rune Sellevag
//
// This file is distributed under the MIT License. See the accompanying file
// LICENSE.txt or http://www.opensource.org/licenses/mit-license.php for terms
// and conditions.

#ifndef SCILIB_LINALG_LAPACK_TYPES_H
#define SCILIB_LINALG_LAPACK_TYPES_H

//------------------------------------------------------------------------------
// Define integer type for use with BLAS, LAPACK and Intel MKL.

#ifdef USE_MKL
#include <mkl.h>
#define BLAS_INT MKL_INT
#else
#include <cblas.h>
#ifdef blasint
#define BLAS_INT blasint
#else
#define BLAS_INT int
#endif
#endif

//------------------------------------------------------------------------------
// Define complex type for use with LAPACK.

#ifndef USE_MKL
#include <complex>
#define lapack_complex_float std::complex<float>
#define lapack_complex_double std::complex<double>
#endif

#endif // SCILIB_LINALG_LAPACK_TYPES_H

max_stars_repo_name: stigrs/scilib

max_issues_repo_name: stigrs/scilib

max_forks_repo_name: stigrs/scilib


text:
/**
 */
#include <time.h>
#include <math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_sort_vector.h>
#include <gsl/gsl_interp.h>
#include <gsl/gsl_sort_double.h>
#include "inout_aper.h"
#include "aXe_grism.h"
#include "spce_sect.h"
#include "spc_back.h"
#include "spce_PET.h"
#include "spc_wl_calib.h"
#include "aXe_errors.h"
#include "fringe_conf.h"
#include "spc_resp.h"
#include "spce_pathlength.h"
#include "aper_conf.h"
#include "specmodel_utils.h"
#include "model_utils.h"
#include "spc_fluxcube.h"
#include "fringe_conf.h"
#include "dirimage_model.h"

#define MAX(x,y) (((x)>(y))?(x):(y))
#define MIN(x,y) (((x)<(y))?(x):(y))
#define SQR(x) ((x)*(x))

int
compute_dirimage_model(char dirim_file[], char conf_file[], char tpass_file[],
		       char specmod_file[], char objmod_file[], char aper_file[],
		       const double model_scale, const double tel_area, const double lambda_psf,
		       observation *obs, char map_file[])
{
  object          **oblist;
  dirobject       **dirlist;
  spectral_models *spec_mod;
  object_models   *obj_mod;

  interpolator    *tpass;

  gsl_matrix *dirimage_matrix=NULL;

  px_point npixels;

  // load the object list
  fprintf (stdout, "aXe_DIRIMAGE: Loading object aperture list...");
  oblist = file_to_object_list_seq (aper_file, obs);
  fprintf (stdout,"%d objects loaded.\n",object_list_size(oblist));

  // check whether highres models
  // are given
  if (strlen(specmod_file) > 0)
    // load the spectral models
    spec_mod = load_spectral_models(specmod_file);
  else
    // or set the struct to NULL
    spec_mod = NULL;

  // check whether direct emission models
  // are given
  if (strlen(objmod_file) > 0)
    obj_mod = load_object_models(objmod_file);
  else
    obj_mod = NULL;

  // get the sensitivity curve of the total passband
  tpass =  get_filter_sensitivity(tpass_file, tel_area);

  // get the image dimensions
  npixels = get_npixel(obs);

  // create the list of direct objects to be simulated
  // the interpolation type is fixed to linear
  dirlist = oblist_to_dirlist2(dirim_file, conf_file, npixels, oblist,
			       spec_mod, obj_mod, model_scale, 1);

  // determine the XOFF and YOFF values
  // for the various beams
  fill_xy_offsets(dirlist, conf_file);

  // create and fill the pixel matrix for the direct images
  dirimage_matrix = make_dirimage(oblist, dirlist, npixels, lambda_psf, tpass);

  // store the contamination image
  gsl_to_FITSimage (dirimage_matrix, map_file, 1, NULL);

  // check if memory must be released;
  // do it if necessary
  if (spec_mod != NULL)
    free_spectral_models(spec_mod);
  if (obj_mod != NULL)
    free_object_models(obj_mod);
  if (oblist !=NULL)
    free_oblist (oblist);

  // always release these
  // objects
  free_dirlist(dirlist);
  gsl_matrix_free(dirimage_matrix);
  free_interp(tpass);

  // return always '1'
  return 1;
}

gsl_matrix *
make_dirimage(object **oblist, dirobject **dirlist, const px_point npixels,
	      const double lambda_psf, interpolator *tpass)
{
  dirobject       *actdir;
  beam             actbeam;

  gsl_matrix      *dirimage_matrix;

  double cps        = 0.0;
  double sval       = 0.0;
  double value      = 0.0;

  int ii=0;
  int nx=0;
  int ny=0;

  d_point dpixel;

   // allocate memory for the image matrix
  dirimage_matrix = gsl_matrix_alloc(npixels.x, npixels.y);
  gsl_matrix_set_all(dirimage_matrix,0.0);

  // go over each beam model
  ii = 0;
  while (oblist[ii] != NULL)
    {

      // get the direct object for the actual model spectrum
      actdir = get_dirobject_from_list(dirlist, oblist[ii]->ID);

      // get the beam for the actual model spectrum
      actbeam = oblist[ii]->beams[0];

      // report onto the screen
      fprintf(stdout, "aXe_DIRIMAGE: modelling object %i ...", oblist[ii]->ID);

      // get the integrated direct image intensity in cps
      cps = get_cps_for_dirobject(tpass, actdir);

      // correct the direct image positions for
      // the offset values introduced by building
      // the direct image objects around the reference
      // position, which is shifted from the true
      // direct image position by XOFF and YOFF as
      // given in the aXe configuration file
      actdir->ix_min -= (int)floor(actdir->xy_off[actbeam.ID].x + 0.5);
      actdir->ix_max -= (int)floor(actdir->xy_off[actbeam.ID].x + 0.5);
      actdir->iy_min -= (int)floor(actdir->xy_off[actbeam.ID].y + 0.5);
      actdir->iy_max -= (int)floor(actdir->xy_off[actbeam.ID].y + 0.5);
      actbeam.refpoint.x -= actdir->xy_off[actbeam.ID].x;
      actbeam.refpoint.y -= actdir->xy_off[actbeam.ID].y;

      // go over each pixel in the direct object area
      for (nx=actdir->ix_min; nx<=actdir->ix_max; nx++)
	{
	  // make sure to be inside the image
	  if (nx < 0 || nx >=  (int)dirimage_matrix->size1)
	    // skip column if not
	    continue;
	  for (ny=actdir->iy_min; ny<=actdir->iy_max; ny++)
	    {
	      // make sure to be inside the image
	      if (ny < 0 || ny >=  (int)dirimage_matrix->size2)
		// skip element if not
		continue;

	      // fill the dpixel structure
	      dpixel.x = (double)nx;
	      dpixel.y = (double)ny;

	      // check whether there is
	      // an image template
	      if (actdir->dirim)
		// get the pixel intensity from the image template
		sval = get_diremission_value(actdir->dirim, dpixel.x - actbeam.refpoint.x, dpixel.y - actbeam.refpoint.y);
	      else
		// get the Gaussian pixel intensity;
		//
		// do a subsampling over the pixel
		// to get a more appropriate value for the
		// emission val
		sval = get_sub_emodel_value(dpixel, actbeam, actdir->drzscale);

	      // compute and set the new pixel values
	      value = gsl_matrix_get(dirimage_matrix, nx, ny) + sval*cps;
	      gsl_matrix_set(dirimage_matrix, nx, ny, value);
	    }
	}

      // release the space for the various structures
      fprintf(stdout, " Done\n");
      ii++;
    }

  // return the image matrix
  return dirimage_matrix;
}

double
get_cps_for_dirobject(interpolator *tpass, dirobject *actdir)
{
  interpolator *combine;

  double cps=0.0;

  // create aa new multiplicator
  // interpolator
  combine = combine_tpass_SED(tpass, actdir);

  // integrate over the interpolator
  cps = integrate_interpolator(combine);

  // release memory
  free_interp(combine);

  // return the cps value
  return cps;
}

double
integrate_interpolator(interpolator *combine)
{
  double integral=0.0;

  int index;

  for (index=1; index < combine->nvals-1; index++)
    // add the increment
    integral +=  combine->yvals[index] * (combine->xvals[index+1] - combine->xvals[index-1]);

  // dont forget the start and end piece
  integral +=  combine->yvals[0]                * (combine->xvals[1]                - combine->xvals[0]);
  integral +=  combine->yvals[combine->nvals-1] * (combine->xvals[combine->nvals-1] - combine->xvals[combine->nvals-2]);

  // then return the half
  return integral/2.0;
}

interpolator *
combine_tpass_SED(interpolator *tpass, dirobject *actdir)
{
  interpolator *combine;

  int index=0;
  int n_additional=0;
  int act_index;

  gsl_vector     *indep_data;

  // go over all SED data points
  for (index = 0; index < actdir->SED->npoints; index++)
    {
      // transform to Angstrom
      actdir->SED->wavelength[index] *= 10.0;

      // check whether the independent data point
      // falls in the area of the sensitivity
      if (actdir->SED->wavelength[index] < tpass->xmax
	  && actdir->SED->wavelength[index] >  tpass->xmin)
	// enhance the number
	// of additional data points
	n_additional += 1;
    }

  // allocate memory for the new independent value list
  indep_data   = gsl_vector_alloc(tpass->nvals + n_additional);

  // go over all sensitivity data points and fill in
  // the independent values
  for (index = 0; index < tpass->nvals; index++)
    gsl_vector_set(indep_data, index, tpass->xvals[index]);

  // set the counter to the next free index
  act_index = tpass->nvals;

  // go over all SED data points
  for (index = 0; index < actdir->SED->npoints; index++)
    // check whether the independent data point
    // falls in the area of the sensitivity
    if (actdir->SED->wavelength[index] < tpass->xmax
	&& actdir->SED->wavelength[index] >  tpass->xmin)
      {
	// fill in the additional independent data point
	// from the SED
	gsl_vector_set(indep_data, act_index, actdir->SED->wavelength[index]);

	// enhance the counter
        act_index += 1;
      }

  // sort the vector
  gsl_sort_vector(indep_data);

  // compose a new interpolator from the independent values
  combine = get_combined_tpass_SED(indep_data, tpass, actdir);

  // release memory in the weight vector
  gsl_vector_free(indep_data);

  // go over all SED data points
  for (index = 0; index < actdir->SED->npoints; index++)
    // transform to nm
    actdir->SED->wavelength[index] /= 10.0;

  // return the interpolator
  return combine;
}

interpolator *
get_combined_tpass_SED(gsl_vector *indep_data, interpolator *tpass, dirobject *actdir)
{
  interpolator *combine;

  gsl_vector_int *indep_weight;

  int n_new=0;
  int index;
  int act_index;

  double *xvals_new;
  double *yvals_new;

  // allocat memory for the weight vector
  indep_weight = gsl_vector_int_alloc(indep_data->size);

  // set all the weights
  gsl_vector_int_set_all(indep_weight, 1);

  // go over the independent data
  for (index = 1; index < (int)indep_data->size; index++)
    // check whether the current entry is equal the previous one
    if (gsl_vector_get(indep_data, index) ==  gsl_vector_get(indep_data, index-1))
      // set the weight of the current entry to zero
      gsl_vector_int_set(indep_weight,index, 0);

  // initialize the
  // final number
  n_new = 0;

  // go over the weight array
  for (index=0; index < (int)indep_weight->size; index++)
    // just count the weights
    n_new += gsl_vector_int_get(indep_weight, index);

  // allocate memory for the interpolator data
  xvals_new = (double *)malloc(n_new * sizeof(double));
  yvals_new = (double *)malloc(n_new * sizeof(double));

  // go over the weight array
  act_index=0;
  for (index=0; index < (int)indep_weight->size; index++)
    // check for weight
    if (gsl_vector_int_get(indep_weight, index))
      {
	// take the independent value from the vector
	xvals_new[act_index] = gsl_vector_get(indep_data, index);

	// compute the new dependent value bye multiplying
	// the brightness with the sensitivity
	yvals_new[act_index] = eval_interp(tpass, xvals_new[act_index])
	  * get_flux_from_SED(actdir->SED, xvals_new[act_index]);

	// enhance the counter
	act_index++;
      }

  // create a new interpolator
  combine = create_interp(n_new, TPASS_INTERP_TYPE, xvals_new, yvals_new);

  // release memory in the weight vector
  gsl_vector_int_free(indep_weight);

  // return the interpolator
  return combine;
}

/*
 * Function: get_filter_sensitivity
 * The function reads in a total passband file from fits format
 * and stroes the data as an interpolation function. The dependent
 * data values are transformed from throughput to sensitivity,
 * and the interpolation function is returned.
 *
 * Parameters:
 * @param tpass_file  - the full name of the total passband
 * @param tel_are     - the collecting area of the telescope
 *
 * Returns:
 * @return tpass    - the filter sensitivity
 */
interpolator *
get_filter_sensitivity(const char tpass_file[], const double tel_area)
{
  interpolator    *tpass;

  double h_erg  = 6.6260693E-27;
  double c_cm   = 2.99792458E+10;

  double factor = 0.0;

  int index=0;

  fprintf(stdout, "Load Total bandpass table :%s\n", tpass_file);
  tpass =  create_interp_ftable(tpass_file, 2,"WAVELENGTH","THROUGHPUT", TPASS_INTERP_TYPE);

  for (index=0; index < tpass->nvals; index++)
    {
      // compute the conversion factor
      factor = tel_area / (h_erg * c_cm / (1.0E-08 * tpass->xvals[index]));

      // apply the conversion factor
      tpass->yvals[index] = tpass->yvals[index] * factor;
    }

  // return the filter sensitivity
  return tpass;
}

max_stars_repo_name: sosey/pyaxe

max_issues_repo_name: sosey/pyaxe

max_forks_repo_name: sosey/pyaxe


text:
#pragma once

#include "datastructures/PointBuffer.h"
#include "io/PointsPersistence.h"
#include "math/AABB.h"
#include "point_source/PointSource.h"
#include "pointcloud/FileStats.h"
#include "pointcloud/PointAttributes.h"
#include "tiling/Sampling.h"
#include "util/Definitions.h"
#include "util/Transformation.h"
#include <reflection/StaticReflection.h>
#include <threading/Semaphore.h>
#include <threading/TaskSystem.h>

#include <atomic>
#include <deque>
#include <gsl/gsl>
#include <string>
#include <thread>

#include <taskflow/taskflow.hpp>

struct ProgressReporter;
struct TilingAlgorithmBase;
struct ThroughputSampler;

/**
 * Which tiling strategy to use?
 */
enum class TilingStrategy
{
  /**
   * Use the accurate strategy that samples from the root node. This will be
   * slower than the 'Fast' strategy, that skips the first couple of levels and
   * reconstructs them afterwards
   */
  Accurate,
  /**
   * Use the fast strategy that skips the first couple of levels and starts
   * tiling deeper in the octree to enable increased parallelism. The skipped
   * levels are reconstructed and thus contain duplicated data
   */
  Fast
};

struct FixedThreadCount
{
  uint32_t num_threads_for_reading;
  uint32_t num_threads_for_indexing;
};

struct AdaptiveThreadCount
{
  uint32_t num_threads;
};

/**
 * Thread number configuration. Either a fixed number of read and index threads, or
 * a total fixed number of threads that are scheduled adaptively
 */
using ThreadConfig = std::variant<FixedThreadCount, AdaptiveThreadCount>;

struct TilerMetaParameters
{
  float spacing_at_root;
  uint32_t max_depth;
  size_t max_points_per_node;
  size_t batch_read_size;
  size_t internal_cache_size;
  bool shift_points_to_origin;
  bool create_journal;
  TilingStrategy tiling_strategy;
  std::variant<FixedThreadCount, AdaptiveThreadCount> thread_count;
};

/**
 * Abstract command for reading points from a file.
 */
struct ReadCommand
{
  const fs::path* file_path;
  size_t to_read_count;
};

struct Tiler
{
  Tiler(DatasetMetadata dataset_metadata,
        TilerMetaParameters meta_parameters,
        SamplingStrategy sampling_strategy,
        ProgressReporter* progress_reporter,
        MultiReaderPointSource point_source,
        PointsPersistence& persistence,
        const PointAttributes& input_attributes,
        fs::path output_directory);
  ~Tiler();

  /**
   * Run the tiler. Returns the total number of points that were processed
   */
  size_t run();

private:
  void swap_point_buffers(size_t produced_points_count);

  bool build_execution_graph_for_reading(tf::Taskflow& tf,
                                         uint32_t num_read_threads,
                                         ThroughputSampler& throughput_sampler);
  void execute_read_commands(const std::vector<ReadCommand>& read_commands,
                             util::Range<PointBuffer::PointIterator> read_destination);

  void build_execution_graph_for_indexing(tf::Taskflow& tf,
                                          uint32_t num_indexing_threads,
                                          ThroughputSampler& throughput_sampler);

  void create_read_commands();
  void adjust_read_thread_count(size_t num_read_threads);
  uint32_t max_read_parallelism() const;

  void estimate_read_throughput(ThroughputSampler& sampler, size_t num_points_in_last_cycle) const;
  void estimate_index_throughput(ThroughputSampler& sampler, size_t num_points_in_last_cycle) const;

  DatasetMetadata _dataset_metadata;
  TilerMetaParameters _meta_parameters;
  SamplingStrategy _sampling_strategy;
  ProgressReporter* _progress_reporter;
  MultiReaderPointSource _point_source;
  PointsPersistence& _persistence;

  AABB _bounds;

  const PointAttributes& _input_attributes;
  fs::path _output_directory;

  PointBuffer _points_cache_for_producers, _points_cache_for_consumers;
  size_t _produced_points_count;

  std::deque<ReadCommand> _remaining_read_commands;
  std::vector<ReadCommand> _next_read_commands_per_thread;

  std::unique_ptr<TilingAlgorithmBase> _tiling_algorithm;

  Semaphore _producers, _consumers;

  std::chrono::high_resolution_clock::time_point _begin_read_cycle_time;
  std::chrono::high_resolution_clock::time_point _begin_index_cycle_time;
};

max_stars_repo_name: igd-geo/schwarzwald

max_issues_repo_name: igd-geo/schwarzwald

max_forks_repo_name: igd-geo/schwarzwald


text:
/*
 * Implement Heap sort -- direct and indirect sorting
 * Based on descriptions in Sedgewick "Algorithms in C"
 *
 * Copyright (C) 1999  Thomas Walter
 *
 * 18 February 2000: Modified for GSL by Brian Gough
 *
 * This is free software; you can redistribute it and/or modify it
 * under the terms of the GNU General Public License as published by the
 * Free Software Foundation; either version 2, or (at your option) any
 * later version.
 *
 * This source 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.
 */

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

static inline void swap (void *base, size_t size, size_t i, size_t j);
static inline void downheap (void *data, const size_t size, const size_t N, size_t k, gsl_comparison_fn_t compare);

/* Inline swap function for moving objects around */

static inline void
swap (void *base, size_t size, size_t i, size_t j)
{
  register char *a = size * i + (char *) base;
  register char *b = size * j + (char *) base;
  register size_t s = size;

  if (i == j)
    return;

  do
    {
      char tmp = *a;
      *a++ = *b;
      *b++ = tmp;
    }
  while (--s > 0);
}

#define CMP(data,size,j,k) (compare((char *)(data) + (size) * (j), (char *)(data) + (size) * (k)))

static inline void
downheap (void *data, const size_t size, const size_t N, size_t k, gsl_comparison_fn_t compare)
{
  while (k <= N / 2)
    {
      size_t j = 2 * k;

      if (j < N && CMP (data, size, j, j + 1) < 0)
        {
          j++;
        }

      if (CMP (data, size, k, j) < 0)
        {
          swap (data, size, j, k);
        }
      else
        {
          break;
        }

      k = j;
    }
}

void
gsl_heapsort (void *data, size_t count, size_t size, gsl_comparison_fn_t compare)
{
  /* Sort the array in ascending order. This is a true inplace
     algorithm with N log N operations. Worst case (an already sorted
     array) is something like 20% slower */

  size_t N;
  size_t k;

  if (count == 0)
    {
      return;                   /* No data to sort */
    }

  /* We have n_data elements, last element is at 'n_data-1', first at
     '0' Set N to the last element number. */

  N = count - 1;

  k = N / 2;
  k++;                          /* Compensate the first use of 'k--' */
  do
    {
      k--;
      downheap (data, size, N, k, compare);
    }
  while (k > 0);

  while (N > 0)
    {
      /* first swap the elements */
      swap (data, size, 0, N);

      /* then process the heap */
      N--;

      downheap (data, size, N, 0, compare);
    }
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


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

/**
 * C BLAS function wrappers.
 */

inline void copy_vector_cblas(
  const double* const x, double* const y,
  const unsigned int rows )
{
  cblas_dcopy( rows, x, 1, y, 1 );
}

inline double inner_product_cblas(
  const double* const x, const double* const y,
  const unsigned int rows )
{
  return cblas_ddot( rows, x, 1, y, 1 );
}

inline void outer_product_cblas(
  const double* const x, const double* const y,
  const double alpha, const unsigned int rows_x, const unsigned int rows_y,
  double* m )
{
  cblas_dger( CblasRowMajor, rows_x, rows_y, alpha, x, 1, y, 1, m, rows_y);
}

inline void vector_vector_sum_cblas(
  const double* const x, double* const y,
  const double alpha, const unsigned int rows )
{
  cblas_daxpy( rows, alpha, x, 1, y, 1);
}

/**
 * Helper functions.
 */

inline void copy_matrix(
  const double* const m0, double* const m1,
  const unsigned int rows, const unsigned int cols )
{
  for ( unsigned int i = 0; i < rows; ++i )
  {
    for ( unsigned int j = 0; j < cols; ++j )
    {
      m1[ i * cols + j ] = m0[ i * cols + j ];
    }
  }
}

inline void matrix_matrix_sum(
  const double* const m0, const double* const m1,
  const double alpha, const unsigned int rows, const unsigned int cols,
  double* m2 )
{
  for ( unsigned int i = 0; i < rows; ++i )
  {
    for ( unsigned int j = 0; j < cols; ++j )
    {
      m2[ i * cols + j ] = m0[ i * cols + j ] + alpha * m1[ i * cols + j ];
    }
  }
}

inline void zero_vector(
  const unsigned int rows,
  double* x )
{
  memset(x, 0, rows * sizeof( double ) );
}

inline void zero_matrix(
  const unsigned int rows, const unsigned int cols,
  double* m )
{
  memset(m, 0, rows * cols * sizeof( double ) );
}

inline void zero_diagonal(
  const unsigned int rows, const unsigned int cols,
  double* m )
{
  for ( unsigned int i = 0; i < rows; ++i )
  {
    m[ i * cols + i ] = 0;
  }
}

inline double rand_double( void )
{
  return ( double ) rand() / ( double )( RAND_MAX );
}

inline unsigned int rand_int( const unsigned int max )
{
  return rand() % max;
}

inline unsigned int logistic( const double h )
{
  return 1. / ( 1. + exp( -h ) ) > rand_double();
}

/**
 * Train a fully visible Boltzmann machine with CD-1.
 */
void ctrain(
  double* W, double* b, unsigned int rows, unsigned int cols,
  double* data, unsigned int episodes,
  double epsilon_w, double epsilon_b,
  unsigned int batchsize, unsigned int seed,
  const double momentum_constant, const double decay_constant )
{
  double dW[ rows * cols ];
  double db[ cols ];
  double dW_minus[ rows * cols ];
  double db_minus[ cols ];
  double u;

  const double effective_epsilon_w = 1 / ( double )( batchsize ) * epsilon_w;
  const double effective_epsilon_b = 1 / ( double )( batchsize ) * epsilon_b;

  srand( seed );

  // initialize momentum terms to zero at start of learning
  zero_matrix( rows, cols, dW_minus );
  zero_vector( rows, db_minus );

  for ( unsigned int i = 0; i < episodes / batchsize; ++i )
  {
    // initialize weights and bias change to zero at start of batch
    zero_matrix( rows, cols, dW );
    zero_vector( rows, db );

    for ( unsigned int j = 0; j < batchsize; ++j )
    {
      // compute \Delta W += <ss>_0
      outer_product_cblas( &data[ i * batchsize * cols + j * cols], &data[ i * batchsize * cols + j * cols], 1., rows, cols, dW );

      // compute \Delta b += <s>_0
      vector_vector_sum_cblas( &data[ i * batchsize * cols + j * cols], db, 1, rows );

      // perform a single Gibbs step (update all units once) directly
      // on data vector to avoid copying
      for ( unsigned int k = 0; k < cols; ++k )
      {
        u = inner_product_cblas( &W[ k * cols ], &data[ i * batchsize * cols + j * cols], cols );
        data[ i * batchsize * cols + j * cols + k ] = logistic( u + b[ k ] );
      }

      // compute \Delta W -= <ss>_1
      outer_product_cblas( &data[ i * batchsize * cols + j * cols ], &data[ i * batchsize * cols + j * cols ], -1., rows, cols, dW );

      // compute \Delta b -= <s>_1
      vector_vector_sum_cblas( &data[ i * batchsize * cols + j * cols ], db, -1, rows );
    }

    // add momentum
    if ( momentum_constant > 0. )
    {
      matrix_matrix_sum( dW, dW_minus, momentum_constant, rows, cols, dW );
      vector_vector_sum_cblas( db_minus, db, momentum_constant, rows );
      copy_matrix( dW, dW_minus, rows, cols );
      copy_vector_cblas( db, db_minus, rows );
    }

    // decay weights
    if ( decay_constant > 0. )
    {
      // need to devide by learning rate, since dW is multiplied by it
      // below
      matrix_matrix_sum( dW, W, -decay_constant / effective_epsilon_w, rows, cols, dW );
    }

    // updates weights and biases
    matrix_matrix_sum( W, dW, effective_epsilon_w, rows, cols, W );
    vector_vector_sum_cblas( db, b, effective_epsilon_b, rows );

    // set diagonal of weight matrix to zero
    zero_diagonal( rows, cols, W );
  }
}

/**
 * Sample from a fully visible Boltzmann machine.
 */
void csample(
  const double* W, const double* b, const unsigned int rows, const unsigned int cols,
  const unsigned int episodes, const unsigned int seed,
  double* s, double* samples )
{
  double u;
  srand( seed );

  unsigned int i = 0;
  while(1)
  {
    // perform a single Gibbs step (update all units once on average)
    for ( unsigned int k = 0; k < cols; ++k )
    {
      const unsigned int j = rand_int( cols );
      u = inner_product_cblas( &W[ j * cols ], s, cols );
      s[ j ] = logistic( u + b[ j ] );
    }
    copy_vector_cblas( s, &samples[ i * rows], rows );
    ++i;
    if( i == episodes)
    {
      return;
    }
  }
}

max_stars_repo_name: jakobj/python-vbm

max_issues_repo_name: jakobj/python-vbm

max_forks_repo_name: jakobj/python-vbm


text:
/* ******************************************************************************
 *
 *
 * This program and the accompanying materials are made available under the
 * terms of the Apache License, Version 2.0 which is available at
 * https://www.apache.org/licenses/LICENSE-2.0.
 *
 *  See the NOTICE file distributed with this work for additional
 *  information regarding copyright ownership.
 * 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.
 *
 * SPDX-License-Identifier: Apache-2.0
 ******************************************************************************/

//
//  @author raver119@gmail.com
//

#ifndef LIBND4J_BLAS_HELPER_H
#define LIBND4J_BLAS_HELPER_H

// work around conflict with OpenBLAS
struct bfloat16;
#define BFLOAT16 BFLOAT16

#include <cblas.h>
#include <helpers/logger.h>
#include <types/float16.h>

#ifdef _WIN32
#define CUBLASWINAPI __stdcall
#define CUSOLVERAPI __stdcall
#else
#define CUBLASWINAPI
#define CUSOLVERAPI
#endif

namespace sd {
typedef enum {
  CUBLAS_STATUS_SUCCESS = 0,
  CUBLAS_STATUS_NOT_INITIALIZED = 1,
  CUBLAS_STATUS_ALLOC_FAILED = 3,
  CUBLAS_STATUS_INVALID_VALUE = 7,
  CUBLAS_STATUS_ARCH_MISMATCH = 8,
  CUBLAS_STATUS_MAPPING_ERROR = 11,
  CUBLAS_STATUS_EXECUTION_FAILED = 13,
  CUBLAS_STATUS_INTERNAL_ERROR = 14,
  CUBLAS_STATUS_NOT_SUPPORTED = 15,
  CUBLAS_STATUS_LICENSE_ERROR = 16
} cublasStatus_t;

typedef enum { CUBLAS_OP_N = 0, CUBLAS_OP_T = 1, CUBLAS_OP_C = 2 } cublasOperation_t;

struct cublasContext;
typedef struct cublasContext *cublasHandle_t;

typedef enum {
  CUDA_R_16F = 2,  /* real as a half */
  CUDA_C_16F = 6,  /* complex as a pair of half numbers */
  CUDA_R_32F = 0,  /* real as a float */
  CUDA_C_32F = 4,  /* complex as a pair of float numbers */
  CUDA_R_64F = 1,  /* real as a double */
  CUDA_C_64F = 5,  /* complex as a pair of double numbers */
  CUDA_R_8I = 3,   /* real as a signed char */
  CUDA_C_8I = 7,   /* complex as a pair of signed char numbers */
  CUDA_R_8U = 8,   /* real as a unsigned char */
  CUDA_C_8U = 9,   /* complex as a pair of unsigned char numbers */
  CUDA_R_32I = 10, /* real as a signed int */
  CUDA_C_32I = 11, /* complex as a pair of signed int numbers */
  CUDA_R_32U = 12, /* real as a unsigned int */
  CUDA_C_32U = 13  /* complex as a pair of unsigned int numbers */
} cublasDataType_t;

typedef void (*CblasSgemv)(CBLAS_ORDER Layout, CBLAS_TRANSPOSE TransA, int M, int N, float alpha, float *A, int lda,
                           float *X, int incX, float beta, float *Y, int incY);

typedef void (*CblasDgemv)(CBLAS_ORDER Layout, CBLAS_TRANSPOSE TransA, int M, int N, double alpha, double *A, int lda,
                           double *X, int incX, double beta, double *Y, int incY);

typedef void (*CblasSgemm)(CBLAS_ORDER Layout, CBLAS_TRANSPOSE TransA, CBLAS_TRANSPOSE TransB, int M, int N, int K,
                           float alpha, float *A, int lda, float *B, int ldb, float beta, float *C, int ldc);

typedef void (*CblasDgemm)(CBLAS_ORDER Layout, CBLAS_TRANSPOSE TransA, CBLAS_TRANSPOSE TransB, int M, int N, int K,
                           double alpha, double *A, int lda, double *B, int ldb, double beta, double *C, int ldc);

typedef void (*CblasSgemmBatch)(CBLAS_ORDER Layout, CBLAS_TRANSPOSE *TransA_Array, CBLAS_TRANSPOSE *TransB_Array,
                                int *M_Array, int *N_Array, int *K_Array, float *alpha_Array, float **A_Array,
                                int *lda_Array, float **B_Array, int *ldb_Array, float *beta_Array, float **C_Array,
                                int *ldc_Array, int group_count, int *group_size);

typedef void (*CblasDgemmBatch)(CBLAS_ORDER Layout, CBLAS_TRANSPOSE *TransA_Array, CBLAS_TRANSPOSE *TransB_Array,
                                int *M_Array, int *N_Array, int *K_Array, double *alpha_Array, double **A_Array,
                                int *lda_Array, double **B_Array, int *ldb_Array, double *beta_Array, double **C_Array,
                                int *ldc_Array, int group_count, int *group_size);

#ifdef LAPACK_ROW_MAJOR
#undef LAPACK_ROW_MAJOR
#endif

#ifdef LAPACK_COL_MAJOR
#undef LAPACK_COL_MAJOR
#endif
enum LAPACK_LAYOUT { LAPACK_ROW_MAJOR = 101, LAPACK_COL_MAJOR = 102 };

typedef int (*LapackeSgesvd)(LAPACK_LAYOUT matrix_layout, char jobu, char jobvt, int m, int n, float *a, int lda,
                             float *s, float *u, int ldu, float *vt, int ldvt, float *superb);

typedef int (*LapackeDgesvd)(LAPACK_LAYOUT matrix_layout, char jobu, char jobvt, int m, int n, double *a, int lda,
                             double *s, double *u, int ldu, double *vt, int ldvt, double *superb);

typedef int (*LapackeSgesdd)(LAPACK_LAYOUT matrix_layout, char jobz, int m, int n, float *a, int lda, float *s,
                             float *u, int ldu, float *vt, int ldvt);
typedef int (*LapackeDgesdd)(LAPACK_LAYOUT matrix_layout, char jobz, int m, int n, double *a, int lda, double *s,
                             double *u, int ldu, double *vt, int ldvt);

typedef cublasStatus_t(CUBLASWINAPI *CublasSgemv)(cublasHandle_t handle, cublasOperation_t trans, int m, int n,
                                                  float *alpha, /* host or device pointer */
                                                  float *A, int lda, float *x, int incx,
                                                  float *beta, /* host or device pointer */
                                                  float *y, int incy);

typedef cublasStatus_t(CUBLASWINAPI *CublasDgemv)(cublasHandle_t handle, cublasOperation_t trans, int m, int n,
                                                  double *alpha, /* host or device pointer */
                                                  double *A, int lda, double *x, int incx,
                                                  double *beta, /* host or device pointer */
                                                  double *y, int incy);

typedef cublasStatus_t(CUBLASWINAPI *CublasHgemm)(cublasHandle_t handle, cublasOperation_t transa,
                                                  cublasOperation_t transb, int m, int n, int k,
                                                  __half *alpha, /* host or device pointer */
                                                  __half *A, int lda, __half *B, int ldb,
                                                  __half *beta, /* host or device pointer */
                                                  __half *C, int ldc);

typedef cublasStatus_t(CUBLASWINAPI *CublasSgemm)(cublasHandle_t handle, cublasOperation_t transa,
                                                  cublasOperation_t transb, int m, int n, int k,
                                                  float *alpha, /* host or device pointer */
                                                  float *A, int lda, float *B, int ldb,
                                                  float *beta, /* host or device pointer */
                                                  float *C, int ldc);

typedef cublasStatus_t(CUBLASWINAPI *CublasDgemm)(cublasHandle_t handle, cublasOperation_t transa,
                                                  cublasOperation_t transb, int m, int n, int k,
                                                  double *alpha, /* host or device pointer */
                                                  double *A, int lda, double *B, int ldb,
                                                  double *beta, /* host or device pointer */
                                                  double *C, int ldc);

typedef cublasStatus_t(CUBLASWINAPI *CublasSgemmEx)(cublasHandle_t handle, cublasOperation_t transa,
                                                    cublasOperation_t transb, int m, int n, int k,
                                                    float *alpha, /* host or device pointer */
                                                    void *A, cublasDataType_t Atype, int lda, void *B,
                                                    cublasDataType_t Btype, int ldb,
                                                    float *beta, /* host or device pointer */
                                                    void *C, cublasDataType_t Ctype, int ldc);

typedef cublasStatus_t(CUBLASWINAPI *CublasHgemmBatched)(cublasHandle_t handle, cublasOperation_t transa,
                                                         cublasOperation_t transb, int m, int n, int k,
                                                         __half *alpha, /* host or device pointer */
                                                         __half *Aarray[], int lda, __half *Barray[], int ldb,
                                                         __half *beta, /* host or device pointer */
                                                         __half *Carray[], int ldc, int batchCount);

typedef cublasStatus_t(CUBLASWINAPI *CublasSgemmBatched)(cublasHandle_t handle, cublasOperation_t transa,
                                                         cublasOperation_t transb, int m, int n, int k,
                                                         float *alpha, /* host or device pointer */
                                                         float *Aarray[], int lda, float *Barray[], int ldb,
                                                         float *beta, /* host or device pointer */
                                                         float *Carray[], int ldc, int batchCount);

typedef cublasStatus_t(CUBLASWINAPI *CublasDgemmBatched)(cublasHandle_t handle, cublasOperation_t transa,
                                                         cublasOperation_t transb, int m, int n, int k,
                                                         double *alpha, /* host or device pointer */
                                                         double *Aarray[], int lda, double *Barray[], int ldb,
                                                         double *beta, /* host or device pointer */
                                                         double *Carray[], int ldc, int batchCount);

typedef enum {
  CUSOLVER_STATUS_SUCCESS = 0,
  CUSOLVER_STATUS_NOT_INITIALIZED = 1,
  CUSOLVER_STATUS_ALLOC_FAILED = 2,
  CUSOLVER_STATUS_INVALID_VALUE = 3,
  CUSOLVER_STATUS_ARCH_MISMATCH = 4,
  CUSOLVER_STATUS_MAPPING_ERROR = 5,
  CUSOLVER_STATUS_EXECUTION_FAILED = 6,
  CUSOLVER_STATUS_INTERNAL_ERROR = 7,
  CUSOLVER_STATUS_MATRIX_TYPE_NOT_SUPPORTED = 8,
  CUSOLVER_STATUS_NOT_SUPPORTED = 9,
  CUSOLVER_STATUS_ZERO_PIVOT = 10,
  CUSOLVER_STATUS_INVALID_LICENSE = 11
} cusolverStatus_t;

typedef enum { CUSOLVER_EIG_TYPE_1 = 1, CUSOLVER_EIG_TYPE_2 = 2, CUSOLVER_EIG_TYPE_3 = 3 } cusolverEigType_t;

typedef enum { CUSOLVER_EIG_MODE_NOVECTOR = 0, CUSOLVER_EIG_MODE_VECTOR = 1 } cusolverEigMode_t;

struct cusolverDnContext;
typedef struct cusolverDnContext *cusolverDnHandle_t;

typedef cusolverStatus_t(CUSOLVERAPI *CusolverDnSgesvdBufferSize)(cusolverDnHandle_t handle, int m, int n, int *lwork);

typedef cusolverStatus_t(CUSOLVERAPI *CusolverDnDgesvdBufferSize)(cusolverDnHandle_t handle, int m, int n, int *lwork);

typedef cusolverStatus_t(CUSOLVERAPI *CusolverDnSgesvd)(cusolverDnHandle_t handle, signed char jobu, signed char jobvt,
                                                        int m, int n, float *A, int lda, float *S, float *U, int ldu,
                                                        float *VT, int ldvt, float *work, int lwork, float *rwork,
                                                        int *info);

typedef cusolverStatus_t(CUSOLVERAPI *CusolverDnDgesvd)(cusolverDnHandle_t handle, signed char jobu, signed char jobvt,
                                                        int m, int n, double *A, int lda, double *S, double *U, int ldu,
                                                        double *VT, int ldvt, double *work, int lwork, double *rwork,
                                                        int *info);

enum BlasFunctions {
  GEMV = 0,
  GEMM = 1,
};

class BlasHelper {
 private:
  bool _hasHgemv = false;
  bool _hasHgemm = false;
  bool _hasHgemmBatch = false;

  bool _hasSgemv = false;
  bool _hasSgemm = false;
  bool _hasSgemmBatch = false;

  bool _hasDgemv = false;
  bool _hasDgemm = false;
  bool _hasDgemmBatch = false;

  CblasSgemv cblasSgemv;
  CblasDgemv cblasDgemv;
  CblasSgemm cblasSgemm;
  CblasDgemm cblasDgemm;
  CblasSgemmBatch cblasSgemmBatch;
  CblasDgemmBatch cblasDgemmBatch;
  LapackeSgesvd lapackeSgesvd;
  LapackeDgesvd lapackeDgesvd;
  LapackeSgesdd lapackeSgesdd;
  LapackeDgesdd lapackeDgesdd;

  CublasSgemv cublasSgemv;
  CublasDgemv cublasDgemv;
  CublasHgemm cublasHgemm;
  CublasSgemm cublasSgemm;
  CublasDgemm cublasDgemm;
  CublasSgemmEx cublasSgemmEx;
  CublasHgemmBatched cublasHgemmBatched;
  CublasSgemmBatched cublasSgemmBatched;
  CublasDgemmBatched cublasDgemmBatched;
  CusolverDnSgesvdBufferSize cusolverDnSgesvdBufferSize;
  CusolverDnDgesvdBufferSize cusolverDnDgesvdBufferSize;
  CusolverDnSgesvd cusolverDnSgesvd;
  CusolverDnDgesvd cusolverDnDgesvd;

 public:
  static BlasHelper &getInstance();

  void initializeFunctions(sd::Pointer *functions);
  void initializeDeviceFunctions(sd::Pointer *functions);

  template <typename T>
  bool hasGEMV();

  template <typename T>
  bool hasGEMM();

  bool hasGEMM(const sd::DataType dtype);
  bool hasGEMV(const sd::DataType dtype);

  template <typename T>
  bool hasBatchedGEMM();

  CblasSgemv sgemv();
  CblasDgemv dgemv();

  CblasSgemm sgemm();
  CblasDgemm dgemm();

  CblasSgemmBatch sgemmBatched();
  CblasDgemmBatch dgemmBatched();

  LapackeSgesvd sgesvd();
  LapackeDgesvd dgesvd();

  LapackeSgesdd sgesdd();
  LapackeDgesdd dgesdd();

  // destructor
  ~BlasHelper() noexcept;
};
}  // namespace sd

#endif


text:
/*
ODE: a program to get optime Runge-Kutta and multi-steps methods.

Copyright 2011-2019, Javier Burguete Tolosa.

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

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

	2. Redistributions in binary form must reproduce the above copyright notice,
		this list of conditions and the following disclaimer in the
		documentation and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY Javier Burguete Tolosa ``AS IS'' AND ANY EXPRESS OR
IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT
SHALL Javier Burguete Tolosa OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/

/**
 * \file ode.c
 * \brief Source file with the main function.
 * \author Javier Burguete Tolosa.
 * \copyright Copyright 2011-2019.
 */
#define _GNU_SOURCE
#include <stdio.h>
#include <stdlib.h>
#include <getopt.h>
#include <math.h>
#include <string.h>
#include <alloca.h>
#include <time.h>
#include <unistd.h>
#include <libxml/parser.h>
#include <glib.h>
#include <gsl/gsl_rng.h>
#include <libintl.h>
#if HAVE_MPI
#include <mpi.h>
#endif
#include "config.h"
#include "utils.h"
#include "optimize.h"
#include "steps.h"
#include "rk.h"

// OPTIMES
//
// RK_2_2 (cfl=1, cfl/steps=1/2)
// t1=1 t2=1
// a10=1 a20=a21=1/2
// b10=1 b20=b21=1/2
// c10=1 c20=0 c21=1
//
// RK_3_2 (cfl=2, cfl/steps=2/3)
// t1=1/2 t2=t3=1
// a10=1 a20=0 a21=1 a30=1/3 a31=0 a32=2/3
// b10=1/2 b20=b21=1/2 b30=b31=b32=1/3
// c10=1/2 c20=0 c21=1/2 c30=0 c31=0 c32=1/2
//
// RK_4_2 (cfl=3, cfl/steps=3/4)
// t1=1/3 t2=2/3 t3=t4=1
// a10=1 a20=0 a21=1 a30=a31=0 a32=1 a40=1/4 a41=a42=0 a43=3/4
// b10=1/3 b20=b21=1/2 b30=b31=b32=1/3
// c10=1/3 c20=0 c21=1/3 c30=c31=0 c32=1/3 c40=c41=c42=0 c43=1/3
//
// RK_5_2 (cfl=4, cfl/steps=4/5)
// t1=1/4 t2=1/2 t3=3/4 t4=t5=1
// a10=1 a20=0 a21=1 a30=a31=0 a32=1 a40=a41=a42=0 a43=1
// a50=1/5 a51=a52=a53=0 a54=4/5
// b10=b20=b21=b30=b31=b32=b40=b41=b42=b43=1/4 b50=b51=b52=b53=b54=1/5
// c10=1/4 c20=0 c21=1/4 c30=c31=0 c32=1/4 c40=c41=c42=0 c43=1/4
// c50=c51=c52=c53=0 c54=1/4
//
// RK_6_2 (cfl=5, cfl/steps=5/6)
// t1=1/5 t2=2/5 t3=3/5 t4=4/5 t5=t6=1
// a10=1 a20=0 a21=1 a30=a31=0 a32=1 a40=a41=a42=0 a43=1 a50=a51=a52=a53=0 a54=1
// a60=1/6 a61=a62=a63=a64=0 a65=5/6
// b10=b20=b21=b30=b31=b32=b40=b41=b42=b43b50=b51=b52=b53=b54=1/5
// b60=b61=b62=b63=b64=b65=1/6
// c10=1/5 c20=0 c21=1/5 c30=c31=0 c32=1/5 c40=c41=c42=0 c43=1/5
// c50=c51=c52=c53=0 c54=1/5 c60=c61=c62=c63=c64=0 c65=1/5
//
// RK_7_2 (cfl=6, cfl/steps=6/7)
// t1=1/6 t2=2/6 t3=3/6 t4=4/6 t5=5/6 t6=t7=1
// a10=1 a20=0 a21=1 a30=a31=0 a32=1 a40=a41=a42=0 a43=1 a50=a51=a52=a53=0 a54=1
// a60=a61=a62=a63=a64=0 a65=1 a70=1/7 a71=a72=a73=a74=a75=0 a76=6/7
// b10=b20=b21=b30=b31=b32=b40=b41=b42=b43b50=b51=b52=b53=b54=b60=b61=b62=b63
// =b64=b65=1/6 b70=b71=b72=b73=b74=b75=b76=1/7
// c10=1/6 c20=0 c21=1/6 c30=c31=0 c32=1/6 c40=c41=c42=0 c43=1/6
// c50=c51=c52=c53=0 c54=1/6 c60=c61=c62=c63=c64=0 c65=1/6
// c70=c71=c73=c74=c75=0 c76=1/6
//
// STEPS_3_2 (cfl=1/2)
// a0=3/4 a1=0 a2=1/4
// b0=3/2 b1=b2=0
// c0=2 c1=c2=0
//
// STEPS_4_2 (cfl=2/3)
// a0=8/9 a1=a2=0 a3=1/9
// b0=4/3 b1=b2=b3=0
// c0=3/2 c1=c2=c3=0
//
// STEPS_5_2 (cfl=3/4)
// a0=15/16 a1=a2=a3=0 a4=1/16
// b0=5/4 b1=b2=b3=b4=0
// c0=4/3 c1=c2=c3=c4=0
//
// STEPS_6_2 (cfl=4/5)
// a0=24/25 a1=a2=a3=a4=0 a5=1/25
// b0=6/5 b1=b2=b3=b4=b5=0
// c0=5/4 c1=c2=c3=c4=c5=0
//
// STEPS_7_2 (cfl=5/6)
// a0=35/36 a1=a2=a3=a4=0 a5=1/36
// b0=7/6 b1=b2=b3=b4=b5=0
// c0=6/5 c1=c2=c3=c4=c5=0
//
// STEPS_8_2 (cfl=6/7)
// a0=48/49 a1=a2=a3=a4=0 a5=1/49
// b0=8/7 b1=b2=b3=b4=b5=0
// c0=7/6 c1=c2=c3=c4=c5=0
//
// RK_3_3 (cfl=1, cfl/steps=1/3)
// t1=t3=1 t2=1/2
// a10=1 a20=3/4 a21=1/4 a30=1/3 a31=0 a32=2/3
// b10=1 b20=b21=1/4 b30=b31=1/6 b32=2/3
// c10=1 c20=0 c21=1 c30=0 c31=0 c32=1
//
// RK_4_3 (cfl=2, cfl/steps=1/2)
// t1=t3=1/2 t2=1
// a10=1 a20=0 a21=1 a30=2/3 a31=0 a32=1/3 a40=a41=a42=0 a43=1
// b10=b20=b21=1/2 b30=b31=b32=b40=b41=b42=1/6 b43=1/2
// c10=c20=c21=1/2 c30=0 c31=c32=c40=c41=c42=c43=1/2
//
// STEPS_4_3 (cfl=1/3)
// a0=16/27 a1=a2=0 a3=11/27
// b0=16/9 b1=b2=0 b3=4/9
// c0=3 c1=c2=0 c3=12/11
//
// STEPS_5_3 (cfl=1/2)
// a0=25/32 a1=a2=a3=0 a4=7/32
// b0=25/16 b1=b2=b3=0 b4=5/16
// c0=2 c1=c2=c3=0 c4=10/7
//
// STEPS_6_3 (cfl=0.58282164314264111327)
// a0=0.85070887167257828565    // (1-a4-a5)
// a1=a2=a3=0
// a4=0.030664864534382320725 // (sqrt(10*a5-1)-10*a5+1)/8
// a5=0.11862626379303939363
//      (63*sqrt(2614)+3221)^(1/3)-5*(63*sqrt(2614)+3221)^(-1/3)-13)/45
// b0=1.4596384360152765194     // a0*c0
// b1=b2=b3=0
// b4=0.052614491749197671863   // a4*c4
// b5=0.20353784933825205974    // a5*c5
// c0=c4=c5=1.715790777102726251        // (1+4*a4+5*a5)
// c1=c2=c3=0
//
// RK_4_4
// b20=b30=b31=0
// t1=t2=b10=b21=1/2
// t3=t4=b32=1
// b40=b43=1/6
// b41=b42=1/3
//
// RK_5_4 (cfl=5.7383119999072648693e-01)
// a10=1
// a20=7.6291847536448934231e-01
// a21=2.3708152463551065768e-01
// a30=2.8277566215397998339e-03
// a31=9.4483076489543642366e-01
// a32=5.2341478483023776481e-02
// a40=6.6309349322657261432e-01
// a41=4.4608078655530437600e-04
// a42=3.3644655897571642926e-01
// a43=1.3867011155652078012e-05
// a50=6.1508837844244168818e-01
// a51=1.6215445203250317099e-01
// a52=1.5031540295088349761e-03
// a53=2.9995002701206228898e-03
// a54=2.1825451522542568299e-01
// b10=4.5307898578512145396e-01
// b20=1.0742850906542086858e-01
// b21=4.1289584245869493973e-01
// b30=4.3724780003845571699e-01
// b31=1.6681424138984873981e+00
// b32=9.0588475806896936944e-02
// b40=3.6664090929447934446e-02
// b41=1.3964776663947759666e-01
// b42=5.8495388577489235762e-01
// b43=1.9355255172522962676e-05
// b50=1.8734964237542205069e-01
// b51=3.0001844590932289148e-01
// b52=1.2903718344972795025e-01
// b53=3.4166232498545846666e-03
// b54=3.8017810501567252292e-01
// c10=4.5307898578512145396e-01
// c20=1.5535519922690283242e-05
// c21=1.7415774725317857109e+00
// c30=1.2525364871787093798e+00
// c31=1.7426727581493662919e+00
// c32=1.7307206145557779180e+00
// c40=4.7046995817421687883e-04
// c41=1.5854730818354572705e+00
// c42=1.7386197420604550316e+00
// c43=1.3957769958693602267e+00
// c50=1.6974106727520458494e-01
// c51=1.6275559956612365566e+00
// c52=7.2955698585694530597e-01
// c53=1.1376557995384390896e+00
// c54=1.7419025884665109550e+00
//
// STEPS_6_4 (cfl=1.6475925238473621578e-01)
// a0=3.4246085571701207596e-01
// a1=a2=0
// a3=1.9179825943473607269e-01
// a4=9.3562124939009442453e-02
// a5=3.721787599092424089e-01
// b0=2.078553105578055306e+00
// b1=b2=b5=0
// b3=1.1641122222796929277e+00
// b4=5.6787174974870979874e-01
// c0=c3=c4=6.0694618695213438363e+00
// c1=c2=c5=0

/**
 * enum to define the method type.
 */
enum
{
  METHOD_TYPE_RUNGE_KUTTA = 1,  ///< Runge-Kutta method.
  METHOD_TYPE_STEPS = 2,        ///< multi-steps method.
} MethodType;

/**
 * enum to define the program error code.
 */
enum
{
  ERROR_CODE_NARGS = 1,         ///< bad command line arguments number.
  ERROR_CODE_BAD_DOC = 2,       ///< bad XML file.
  ERROR_CODE_NO_XML_ROOT = 3,   ///< no XML root node.
  ERROR_CODE_BAD_RK = 4,        ///< bad Runge-Kutta method.
  ERROR_CODE_BAD_STEPS = 5,     ///< bad multi-steps method.
  ERROR_CODE_UNKNOWN_METHOD = 6,        ///< unknown method.
  ERROR_CODE_UNKNOWN_OPTION = 7,        ///< unknown option.
} ErrorCode;

unsigned int nsteps;            ///< steps number.
unsigned int order;             ///< accuracy order.

/**
 * Main function
 *
 * \return 0 on succes, error code on error.
 */
int
main (int argn,                 ///< arguments number.
      char **argc)              ///< argument chains array.
{
  const struct option options[] = {
    {"help", no_argument, NULL, 'h'},
    {"seed", required_argument, NULL, 's'},
    {"threads", required_argument, NULL, 't'},
    {NULL, 0, NULL, 0}
  };
  const char *usage = _("Usage is:\n"
                        "./ode "
                        "[-t --threads threads_number] "
                        "[-s --seed random_seed] "
                        "input_file [variables_file]");
  xmlDoc *doc;
  xmlNode *node;
  gsl_rng *rng0, **rng;
  time_t d0;
  clock_t t0;
  unsigned long int seed = 7l;
  int o, option_index;
  unsigned int i, j, k, h = 0;

#if HAVE_MPI
  // Init MPI
  MPI_Init (&argn, &argc);
#endif

  // Enabling spaces in XML files
  xmlKeepBlanksDefault (0);

  // Setting the threads as the number of processors
  nthreads = sysconf (_SC_NPROCESSORS_CONF);

  // Parsing command line options
  while (1)
    {
      o = getopt_long (argn, argc, "hs:t:", options, &option_index);
      if (o == -1)
        break;
      switch (o)
        {
        case 's':
          seed = atol (optarg);
          break;
        case 't':
          nthreads = atoi (optarg);
          break;
        case 'h':
          printf ("%s\n", usage);
          h = 1;
          break;
        default:
          show_error (_("Unknown option"));
          return ERROR_CODE_UNKNOWN_OPTION;
        }
    }
  argn -= optind;
  if (argn != 1 && argn != 2)
    {
      if (h)
        return 0;
      show_error (usage);
      return ERROR_CODE_NARGS;
    }

  // Init the clock
  t0 = clock ();
  d0 = time (NULL);

#if HAVE_MPI

  // Nodes number
  MPI_Comm_size (MPI_COMM_WORLD, &nnodes);

  // Actual node
  MPI_Comm_rank (MPI_COMM_WORLD, &rank);

#else

  nnodes = 1;
  rank = 0;

#endif

  // Select the numerical model
  printf ("Rank=%d nnodes=%d nthreads=%u\n", rank, nnodes, nthreads);
  printf ("Selecting method optind=%d\n", optind);
  doc = xmlParseFile (argc[optind]);
  if (!doc)
    {
      show_error (_("Unable to parse the input file"));
      return ERROR_CODE_BAD_DOC;
    }
  node = xmlDocGetRootElement (doc);
  if (!node)
    {
      show_error (_("No XML root node"));
      return ERROR_CODE_NO_XML_ROOT;
    }

  // Init a random numbers generator per node and thread
  printf ("Initing random numbers\n");
  rng = (gsl_rng **) g_slice_alloc (nnodes * nthreads * sizeof (gsl_rng *));
  rng0 = gsl_rng_alloc (gsl_rng_ranlxs2);
  gsl_rng_set (rng0, seed);
  for (i = k = 0; (int) i < nnodes; ++i)
    for (j = 0; j < nthreads; ++j, ++k)
      {
        rng[k] = gsl_rng_alloc (gsl_rng_taus2);
        gsl_rng_set (rng[k], gsl_rng_get (rng0));
      }

  if (argn == 2)
    file_variables = fopen (argc[++optind], "w");

  j = rank * nthreads;
  if (!xmlStrcmp (node->name, XML_RUNGE_KUTTA))
    {
      if (!rk_run (node, rng))
        {
          show_error (error_message);
          return ERROR_CODE_BAD_RK;
        }
    }
  else if (!xmlStrcmp (node->name, XML_STEPS))
    {
      if (!steps_run (node, rng))
        {
          show_error (error_message);
          return ERROR_CODE_BAD_STEPS;
        }
    }
  else
    {
      show_error (_("Unknown method type"));
      return ERROR_CODE_UNKNOWN_METHOD;
    }

  printf ("cpu time=%lg real time=%lu\n",
          (clock () - t0) / ((double) CLOCKS_PER_SEC), time (NULL) - d0);

  // Free memory
  xmlFreeDoc (doc);
  if (file_variables)
    fclose (file_variables);
  j = nnodes * nthreads;
  for (i = 0; i < j; ++i)
    gsl_rng_free (rng[i]);
  g_slice_free1 (j * sizeof (gsl_rng *), rng);
  gsl_rng_free (rng0);

#if HAVE_MPI
  // Close MPI
  MPI_Finalize ();
#endif

  return 0;
}

max_stars_repo_name: jburguete/ode

max_issues_repo_name: jburguete/ode

max_forks_repo_name: jburguete/ode


text:
#ifndef __INTERP_2D_H__
#define __INTERP_2D_H__

#include <string.h>
//#include <gsl/gsl_vector.h>
//#include <gsl/gsl_matrix.h>
//#include <gsl/gsl_math.h>
#include <gsl/gsl_spline.h>

#include "cubic_bspline_2d_coeffs.h"
#include "cubic_bspline_2d_interpol.h"

enum
{
    INTERP_2D_LINEAR=0,
    INTERP_2D_CUBIC_BSPLINE=1
};

typedef struct
{
    int size1;
    int size2;
    double * xa;
    double * ya;
    double * za;
    int type;
}interp_2d;

interp_2d * interp_2d_alloc(int size1, int size2);
void interp_2d_free(interp_2d * i2d);

void interp_2d_init(interp_2d * i2d, const double * xa, const double * ya, const double * za, int type);

double interp_2d_eval(interp_2d  * i2d, double x, double y, gsl_interp_accel * accx, gsl_interp_accel * accy);
void interp_2d_eval_grad(interp_2d * i2d, double x, double y, double * grad, gsl_interp_accel * accx, gsl_interp_accel * accy);
double interp_2d_eval_cubic_bspline(interp_2d * i2d, double x, double y, gsl_interp_accel * accx,gsl_interp_accel * accy);

#endif

max_stars_repo_name: fardal/galpy

max_issues_repo_name: fardal/galpy

max_forks_repo_name: fardal/galpy


text:
//
// Created by robin on 2018/9/1.
// Copyright (c) 2018 Robin. All rights reserved.
//

#pragma once

#include <gsl/gsl>

namespace satz {

template <typename P>
const P zero = P::ZERO;

template <typename P>
const P one = P::ONE;

template <typename P>
const P X = P::X;

/**
 * @brief Compute $b^e$.
 * @param b base
 * @param e exponent
 * @param acc
 * @pre $e \ge 0, acc = 1$.
 */
template <typename P, typename N>
P pow (P b, N e, P acc = one<P>) {
  Expects(e >= 0);

  while (e) {
    if (e % 2) acc = acc * b;
    b = b * b;
    e = e / 2;
  }

  return acc;
}

/**
 * @brief Compute $(b^e) mod m$.
 * @param b base
 * @param e exponent
 * @param m modulus
 * @param acc
 * @pre $e >= 0, acc = 1$
 */
template <typename P, typename N>
P mod_pow (P b, N e, const P& m, P acc = one<P>) {
  Expects(e >= 0);

  while (e) {
    if (e % 2) acc = (acc * b) % m;
    b = (b * b) % m;
    e = e / 2;
  }

  return acc;
}

template <typename P>
P gcd (P a, P b) {
  while (b) {
    std::tie(a, b) = std::make_pair(b, a % b);
  }
  return a;
}

/**
 * @brief Compute $(x^{2^p} - x) mod m$.
 * @param p
 * @param m modulus
 */
template <typename P, typename N>
P reduce_exponent (N p, const P& m) {
  auto r = X<P>;
  for (; p > 0; --p) r = (r * r) % m;
  return (r - X<P>) % m;
}

/**
 * @brief Given a polynomial p over GF(2) return whether it is irreducible
 *    with Ben-Or's algorithm.
 * @param p polynomial over GF(2)
 * @return true if p is irreducible.
 * @pre p is non-constant.
 */
template <typename P>
bool ben_or_test (const P& p) {
  Expects(p.degree() > 0);

  auto             d = p.degree();
  for (decltype(d) i = 1; i <= d / 2; ++i) {
    auto b = reduce_exponent(i, p);
    auto g = gcd(p, b);
    if (g != one<P>) return false;
  }
  return true;
}

/**
 * @brief Given a polynomial p over GF(2), return whether it is irreducible.
 * @param p polynomial over GF(2)
 * @return true if p is irreducible.
 */
template <typename P>
bool is_irreducible (const P& p) {
  if (p.degree() > 0) return ben_or_test(p);
  else return false;
}

}

max_stars_repo_name: lie-yan/rabin-fingerprint

max_issues_repo_name: lie-yan/rabin-fingerprint

max_forks_repo_name: lie-yan/rabin-fingerprint