text:
///////////////////////////////////////////////////////////////////////////////
//
// File eigen_mm.h
//
// MIT License

// Copyright (c) 2020 EigenMM

// 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.
//
// Description: EigenMM library definitions
//
///////////////////////////////////////////////////////////////////////////////

#ifndef eigen_mm_H
#define eigen_mm_H

#include <petsc.h>
#include <slepceps.h>

#include <algorithm>
#include <vector>
#include <map>
#include <string>
#include <fstream>
#include <iterator>
#include <iostream>
#include "mpi.h"
#include "tinyxml.h"

#define PI 3.141592653589793238462

class SolverOptions
{
private:
    // communicator options
    PetscInt _nodesperevaluator = 1;
    PetscInt _subproblemsperevaluator = 1;
    PetscInt _totalsubproblems = 1;
    PetscInt _nevaluators = 1;
    PetscInt _nevals = -1;
    PetscInt _taskspernode = 28;

    // eigenvalue partitioning options
    PetscInt _nk = 7;
    PetscInt _nb = 4;
    PetscInt _p = 0;
    PetscInt _nv = 10;
    PetscInt _splitmaxiters = 10;
    PetscInt _raditers = 10;
    PetscReal _splittol = 0.9;
    PetscReal _radtol = 1e-3;
    PetscReal _L = 0.01;
    PetscReal _R = -1.0;

    // save and print options
    bool _terse = false;
    bool _details = false;
    bool _debug = false;
    bool _save_operators = false;
    bool _save_correctness = false;
    bool _save_eigenvalues = false;
    bool _save_eigenbasis = false;
    char _operators_filename[2048];
    char _correctness_filename[2048];
    char _eigenvalues_filename[2048];
    char _eigenbasis_filename[2048];

public:
    SolverOptions() {}

    // set communicator options
    void set_nodesperevaluator(PetscInt v) { _nodesperevaluator = v; }
    void set_subproblemsperevaluator(PetscInt v) { _subproblemsperevaluator = v; }
    void set_totalsubproblems(PetscInt v) { _totalsubproblems = v; }
    void set_nevaluators(PetscInt v) { _nevaluators = v; }
    void set_nevals(PetscInt v) { _nevals = v; }
    void set_taskspernode(PetscInt v) { _taskspernode = v; }

    // set eigenvalue partitioning options
    void set_nk(PetscInt v) { _nk = v; }
    void set_nb(PetscInt v) { _nb = v; }
    void set_p(PetscInt v) { _p = v; }
    void set_nv(PetscInt v) { _nv = v; }
    void set_splitmaxiters(PetscInt v) { _splitmaxiters = v; }
    void set_raditers(PetscInt v) { _raditers = v; }
    void set_splittol(PetscReal v) { _splittol = v; }
    void set_radtol(PetscReal v) { _radtol = v; }
    void set_L(PetscReal v) { _L = v; }
    void set_R(PetscReal v) { _R = v; }

    // set save and print options
    void set_terse(bool v) { _terse = v; }
    void set_details(bool v) { _details = v; }
    void set_debug(bool v) { _debug = v; }
    void set_save_operators(bool v, const char* filename)
    {
        _save_operators = v;
        sprintf(_operators_filename, "%s", filename);
    }
    void set_save_correctness(bool v, const char* filename)
    {
        _save_correctness = v;
        sprintf(_correctness_filename, "%s", filename);
    }
    void set_save_eigenvalues(bool v, const char* filename)
    {
        _save_eigenvalues = v;
        sprintf(_eigenvalues_filename, "%s", filename);
    }
    void set_save_eigenbasis(bool v, const char* filename)
    {
        _save_eigenbasis = v;
         sprintf(_eigenbasis_filename, "%s", filename);
    }

    // get communicator options
    PetscInt nodesperevaluator() { return _nodesperevaluator; }
    PetscInt subproblemsperevaluator() { return _subproblemsperevaluator; }
    PetscInt totalsubproblems() { return _totalsubproblems; }
    PetscInt nevaluators() { return _nevaluators; }
    PetscInt nevals() { return _nevals; }
    PetscInt taskspernode() { return _taskspernode; }

    // get eigenvalue partitioning options
    PetscInt nk() { return _nk; }
    PetscInt nb() { return _nb; }
    PetscInt p() { return _p; }
    PetscInt nv() { return _nv; }
    PetscInt splitmaxiters() { return _splitmaxiters; }
    PetscInt raditers() { return _raditers; }
    PetscReal splittol() { return _splittol; }
    PetscReal radtol() { return _radtol; }
    PetscReal L() { return _L; }
    PetscReal R() { return _R; }

    // get save and print options
    bool terse() { return _terse; }
    bool details() { return _details; }
    bool debug() { return _debug; }
    bool save_operators() { return _save_operators; }
    bool save_correctness() { return _save_correctness; }
    bool save_eigenvalues() { return _save_eigenvalues; }
    bool save_eigenbasis() { return _save_eigenbasis; }
    char* operators_filename() { return _operators_filename; }
    char* correctness_filename() { return _correctness_filename; }
    char* eigenvalues_filename() { return _eigenvalues_filename; }
    char* eigenbasis_filename() { return _eigenbasis_filename; }
};

struct NodeInfo
{
    // world
    PetscInt worldrank, worldsize;

    // evaluator
    PetscInt rank, size, id;
    MPI_Comm comm;

    // additional properties
    PetscInt nevaluators;

    const char* procname;
    PetscInt procid;

    // Results
    PetscInt neval;
    PetscInt neval0;
};

class eigen_mm
{
private:

    EPS eps;
    SolverOptions opts;
    NodeInfo node;
    Mat K, M, V;
    Mat K_global, M_global;
    Vec lambda;
    std::vector<PetscReal> intervals;
    std::vector<PetscReal> residuals;

    void findUpperBound();
    void rescaleInterval();
    void formSubproblems();
    PetscReal count_subintervals(PetscInt n, PetscReal Nhat, PetscInt Nbar,
        std::vector<PetscReal> x, std::vector<PetscInt> &RA, std::vector<PetscInt> &A);
    void merge_approximations(PetscReal Nhat, PetscInt Nbar,
        std::vector<PetscReal> &x, std::vector<PetscInt> &RA,
        std::vector<PetscReal> y, std::vector<PetscInt> RB,
        std::vector<PetscInt> &A);
    PetscInt solveSubproblems();
    void formEigenbasis(PetscInt neval);

    PetscInt solveSubProblem(PetscReal *intervals, int job);
    void splitSubProblem(PetscReal a, PetscReal b, PetscReal *c,
        PetscInt *out_ec_left, PetscInt *out_ec_right);
    void global_refine(PetscInt n,
                       std::vector<PetscReal>  x,
                       std::vector<PetscReal> &y,
                       std::vector<PetscInt>   C,
                       PetscReal Nhat);
    void balance_intervals(PetscReal   a, PetscReal   b, PetscReal   *c,
                           PetscInt reva, PetscInt revb, PetscInt *revc,
                           PetscInt  *Cl, PetscInt  *Cr);
    PetscInt computeDev_approximate(PetscReal a, PetscReal U, PetscBool rl);
    PetscInt computeDev_exact(PetscReal a, PetscBool rl);
    PetscReal computeRadius(Mat &A);
    void countInterval(PetscReal a, PetscReal b, PetscInt *count);

public:

    eigen_mm();
    ~eigen_mm();

    int init(Mat &K_in, Mat &M_in, SolverOptions &opt);
    int solve(Mat *V_out, Vec *lambda_out);

    void exactEigenvalues_square_neumann(int Ne,
        std::vector<PetscReal> &lambda,
        std::vector<PetscReal> &eta1,
        std::vector<PetscReal> &eta2);
    void exactEigenvalues_cube_neumann(int Ne,
        std::vector<PetscReal> &lambda,
        std::vector<PetscReal> &eta1,
        std::vector<PetscReal> &eta2,
        std::vector<PetscReal> &eta3);

};

#endif //eigen_mm_H

max_stars_repo_name: paralab/EigenMM

max_issues_repo_name: paralab/EigenMM

max_forks_repo_name: paralab/EigenMM


text:
/** Fit a set of points to a model that is subject to rotation and translation.
 *
 * The model is defined by a set of primitive faces such as plane, cylinder, sphere, etc.
 * The measured points are first translated by (xt, yt, zt),
 * then rotated about x-y-z axes by g, b, a angles in sequence,
 *
 * These 6 parameters are the output of the fit.
 *
 * Author: Yuan Mei
 */
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_multifit_nlinear.h>

#ifdef LINE_MAX
#undef LINE_MAX
#define LINE_MAX 4096
#endif

#ifndef LINE_MAX
#define LINE_MAX 4096
#endif

/** Check if the character is a field separator. */
#ifndef sepq
#define sepq(a) ((a)==' ' || (a)=='\t')
#endif

#define NDIM   3    /**< number of dimensions */

enum {
    FACE_PLANE=1,
    FACE_CYLINDER=2,
    FACE_SPHERE=3,
    FACE_LINE=12,
    FACE_POINT=13
};

struct face {
    int ftype; /**< type of face: plane, cylinder etc... */
    double x0; /**< origin */
    double y0;
    double z0;
    double nx; /**< normal vector */
    double ny;
    double nz;
    double r;  /**< radius */
};

struct data {
    double pr;   /**< probe radius */
    const struct face *faces;
    size_t n;    /**< number of measured data points */
    size_t *fid; /**< idx of face in the registered faces */
    double *x;   /**< measured point */
    double *y;
    double *z;
    double *w;   /**< weights */
    double *x1;  /**< after transformation */
    double *y1;
    double *z1;
};

/** Distance from measured point to face. */
int dist_f(const gsl_vector *p, void *data, gsl_vector *f)
{
    double pr = ((struct data *)data)->pr;
    const struct face *faces = ((struct data *)data)->faces;
    size_t n = ((struct data *)data)->n;
    size_t *fid = ((struct data *)data)->fid;
    double *x = ((struct data *)data)->x;
    double *y = ((struct data *)data)->y;
    double *z = ((struct data *)data)->z;

    double a  = gsl_vector_get(p, 0), Ca = cos(a), Sa = sin(a);
    double b  = gsl_vector_get(p, 1), Cb = cos(b), Sb = sin(b);
    double g  = gsl_vector_get(p, 2), Cg = cos(g), Sg = sin(g);
    double xt = gsl_vector_get(p, 3);
    double yt = gsl_vector_get(p, 4);
    double zt = gsl_vector_get(p, 5);

    size_t i;
    for (i = 0; i < n; i++) {
        int ftype = faces[fid[i]].ftype;
        double x0 = faces[fid[i]].x0;
        double y0 = faces[fid[i]].y0;
        double z0 = faces[fid[i]].z0;
        double nx = faces[fid[i]].nx;
        double ny = faces[fid[i]].ny;
        double nz = faces[fid[i]].nz;
        double r  = faces[fid[i]].r;

        double x1, y1, z1;
        x1 = Ca*Cb * (x[i]+xt) + (Ca*Sb*Sg-Sa*Cg) * (y[i]+yt) + (Ca*Sb*Cg+Sa*Sg) * (z[i]+zt);
        y1 = Sa*Cb * (x[i]+xt) + (Sa*Sb*Sg+Ca*Cg) * (y[i]+yt) + (Sa*Sb*Cg-Ca*Sg) * (z[i]+zt);
        z1 =  - Sb * (x[i]+xt) +            Cb*Sg * (y[i]+yt) +            Cb*Cg * (z[i]+zt);

        ((struct data *)data)->x1[i] = x1;
        ((struct data *)data)->y1[i] = y1;
        ((struct data *)data)->z1[i] = z1;

        double d, dist=0.0;
        switch (ftype) {
        case FACE_PLANE:
            d = -(nx*x0 + ny*y0 + nz*z0);
            dist = fabs(nx*x1 + ny*y1 + nz*z1 + d)
                / sqrt(nx*nx + ny*ny + nz*nz);
            dist -= pr;
            break;
        case FACE_CYLINDER:
        case FACE_LINE:
        {
            double dx = x1 - x0, dy = y1 - y0, dz = z1 - z0;
            double s1 = dy * nz - dz * ny, s2 = dz*nx - dx*nz, s3 = dx*ny - dy*nx;
            dist = fabs(sqrt((s1*s1 + s2*s2 + s3*s3)/(nx*nx + ny*ny + nz*nz)) - r);
            if (ftype == FACE_CYLINDER)
                dist -= pr;
        }
            break;
        case FACE_SPHERE:
        case FACE_POINT:
            dist = fabs(sqrt((x1-x0)*(x1-x0) + (y1-y0)*(y1-y0) + (z1-z0)*(z1-z0)) - r);
            if (ftype == FACE_SPHERE)
                dist -= pr;
            break;
        default:
            dist = 0.0;
            break;
        }
        gsl_vector_set(f, i, dist);
    }

    return GSL_SUCCESS;
}

void callback(const size_t iter, void *params,
              const gsl_multifit_nlinear_workspace *w)
{
    gsl_vector *f = gsl_multifit_nlinear_residual(w);
    gsl_vector *x = gsl_multifit_nlinear_position(w);
    double rcond=0;

    /* compute reciprocal condition number of J(x) */
    /* cond(J) = %8.4f, 1.0 / rcond */
    gsl_multifit_nlinear_rcond(&rcond, w);

    fprintf(stderr, "iter %2zu: a = %7.4f, b = %7.4f, g = %7.4f, xt = %7.4f, yt = %7.4f, zt = %7.4f, |f(x)| = %.4f\n",
            iter,
            gsl_vector_get(x, 0),
            gsl_vector_get(x, 1),
            gsl_vector_get(x, 2),
            gsl_vector_get(x, 3),
            gsl_vector_get(x, 4),
            gsl_vector_get(x, 5),
            gsl_blas_dnrm2(f));
}

/** Read a long line from file.
 * @param[inout] s string of the line, is allocated when s==NULL and n==0 and grown as needed.
 * @param[inout] n current size of s
 * @return s
 */
static char *file_read_long_line(char **s, size_t *n, FILE *fp)
{
    const int bufsz = LINE_MAX;
    char *p;
    size_t cnt, sz;

    if ( *s == NULL && *n == 0 ) {
        *n = bufsz;
        if ( (*s = calloc(*n, sizeof(char))) == NULL ) exit(-1);
    }
    p = *s;
    sz = *n;
    while ( 1 ) {
        if ( fgets(p, sz, fp) == NULL ) return NULL;
        cnt = strlen(*s);
        if ( (*s)[cnt-1] == '\n' ) {
            break;
        } else { /* line too long, expand the buffer */
            *n += bufsz;
            if ( (*s = realloc(*s, (*n)*sizeof(char))) == NULL ) exit(-1);
            p = *s + cnt;
            sz = bufsz;
        }
    }

    return *s;
}

/** Read face definition file.
 *
 * @param[inout] n number of faces read from file.  If the given value
 * *n > 0, *n is interpreted as the number of elements in the faces
 * array and only up to *n faces will be read from file.  However, if
 * the file contains less than *n faces, *n will be updated to reflect
 * the available number of elements.
 *
 * @param[inout] faces array of faces.  If *faces != NULL, *faces will
 * be used rather than allocated.
 */
int read_faces(const char *fname, size_t *n, struct face **faces)
{
    char *linebuf = NULL;
    size_t linen = 0;
    FILE *fp;
    if ((fp = fopen(fname, "r"))==NULL) {
        perror(fname);
        return -1;
    }
    struct face fc;
    int fid;
    size_t lid = 0;
    ssize_t nelem = -1;
    if (*faces == NULL) {
        /* get number of elements in the file */
        while (file_read_long_line(&linebuf, &linen, fp)) {
            lid++;
            if (linebuf[0] == '#' || linebuf[0] == '\n') continue;
            int ret = sscanf(linebuf, "%d %d %lf %lf %lf %lf %lf %lf %lf", &fid, &fc.ftype,
                             &fc.x0, &fc.y0, &fc.z0, &fc.nx, &fc.ny, &fc.nz, &fc.r);
            if (ret < 9 || fid < 0 || fc.ftype <= 0) {
                fprintf(stderr, "Malformatted face at line %zd\n", lid);
            } else {
                if (fid > nelem) nelem = fid;
            }
        }
        nelem++;
        if (nelem == 0) {
            fprintf(stderr, "No valid face in file.\n");
            return -1;
        } else {
            fprintf(stderr, "%zd faces available in file.\n", nelem);
        }
        *n = nelem;
        if ((*faces = calloc(nelem, sizeof(struct face))) == NULL) {
            perror("calloc *faces");
            return -1;
        }
    }
    rewind(fp);
    nelem = 0;
    lid = 0;
    while (file_read_long_line(&linebuf, &linen, fp) && (nelem < *n)) {
        lid++;
        if (linebuf[0] == '#' || linebuf[0] == '\n') continue;
        int ret = sscanf(linebuf, "%d %d %lf %lf %lf %lf %lf %lf %lf", &fid, &fc.ftype,
                         &fc.x0, &fc.y0, &fc.z0, &fc.nx, &fc.ny, &fc.nz, &fc.r);
        if (ret < 9 || fid < 0 || fc.ftype <= 0) {
            fprintf(stderr, "Malformatted face at line %zd, skipped.\n", lid);
        } else {
            if (fid >= *n) break;
            struct face *fc1 = *faces;
            memcpy(&fc1[fid], &fc, sizeof(fc));
            nelem++;
        }
    }

    fprintf(stderr, "%zd faces constructed.\n", nelem);
    free(linebuf);
    fclose(fp);
    return 0;
}

int read_points(const char *fname, size_t *n, struct data *data)
{
    char *linebuf = NULL;
    size_t linen = 0;
    FILE *fp;
    if ((fp = fopen(fname, "r"))==NULL) {
        perror(fname);
        return -1;
    }

    int fidmax = -1;
    *n = 0;
    size_t lid = 0;
    while (file_read_long_line(&linebuf, &linen, fp)) {
        lid++;
        if (linebuf[0] == '#' || linebuf[0] == '\n') continue;
        int ret = 0;
        for (int i=0; i<linen; i++) { /* count number of ';' */
            if (linebuf[i] == ';') ret++;
        }
        int fid;
        sscanf(linebuf, "%*d;%d", &fid);
        if (ret < 8 || fid < 0) {
            fprintf(stderr, "Malformatted point at line %zd\n", lid);
        } else {
            (*n)++;
            if (fid > fidmax) fidmax = fid;
        }
    }
    fprintf(stderr, "%zd points in file.  fidmax = %d\n", *n, fidmax);
    data->n = *n;
    data->fid = calloc(data->n, sizeof(size_t));
    data->x   = calloc(data->n, sizeof(double));
    data->y   = calloc(data->n, sizeof(double));
    data->z   = calloc(data->n, sizeof(double));
    data->w   = calloc(data->n, sizeof(double));
    data->x1  = calloc(data->n, sizeof(double));
    data->y1  = calloc(data->n, sizeof(double));
    data->z1  = calloc(data->n, sizeof(double));

    rewind(fp);

    lid = 0;
    size_t idx = 0;
    while (file_read_long_line(&linebuf, &linen, fp)) {
        lid++;
        if (linebuf[0] == '#' || linebuf[0] == '\n') continue;
        int fid = 0, ret = 0;
        double x, y, z, s;
        // ret = sscanf(linebuf, "%*d;%d;%*s ;%*d;%lf;%lf;%lf;;%lf", &fid, &x, &y, &z, &s);
        char *buf = linebuf, *endptr;
        for (; *buf != ';' ; buf++){;} buf++; /* find the next character past a ';' */
        fid = strtol(buf, &endptr, 0); if (endptr > buf) ret++; buf = endptr++;
        for (; *buf != ';' ; buf++){;} buf++;
        for (; *buf != ';' ; buf++){;} buf++;
        for (; *buf != ';' ; buf++){;} buf++;
        x = strtod(buf, &endptr); if (endptr > buf) ret++; buf = endptr++;
        for (; *buf != ';' ; buf++){;} buf++;
        y = strtod(buf, &endptr); if (endptr > buf) ret++; buf = endptr++;
        for (; *buf != ';' ; buf++){;} buf++;
        z = strtod(buf, &endptr); if (endptr > buf) ret++; buf = endptr++;
        for (; *buf != ';' ; buf++){;} buf++;
        for (; *buf != ';' ; buf++){;} buf++;
        s = strtod(buf, &endptr); if (endptr > buf) ret++; buf = endptr++;

        if (ret < 5 || fid < 0) {
            fprintf(stderr, "Malformatted point at line %zd, skipped.\n", lid);
        } else {
            data->fid[idx] = fid;
            data->x[idx] = x;
            data->y[idx] = y;
            data->z[idx] = z;
            data->w[idx] = (fabs(s)>1e-8) ? (1.0/(s*s)) : 1.0;
            idx++;
        }
    }

    fprintf(stderr, "%zd points read.\n", idx);
    free(linebuf);
    fclose(fp);
    return 0;
}

int main(int argc, char **argv)
{
    if (argc != 10) {
        fprintf(stderr,
                "Usage: %s faces_file points_file pr a b g xt yt zt\n\n"
                "    pr is the probe radius.\n"
                "    The last 6 parameters are initial guesses.\n"
                "    Points are first translated, then rotated about x(g), y(b), z(a).\n"
                "    The 2nd column in points_file shall be faceid,\n"
                "        which is the first column of faces_file.\n"
                "    The last column in points_file shall be measurement sigma.\n", argv[0]);
        return EXIT_FAILURE;
    }
    size_t nfaces=0;
    struct face *faces=NULL;
    read_faces(argv[1], &nfaces, &faces);

    size_t nd=0;
    struct data data;
    read_points(argv[2], &nd, &data);
    data.faces = faces;

    data.pr = atof(argv[3]);

    #define np 6
    double p_init[np] = {0}; /* starting values */
    for (int i=0; i<np; i++) {
        p_init[i] = atof(argv[4+i]);
    }

    const gsl_multifit_nlinear_type *T = gsl_multifit_nlinear_trust;
    gsl_multifit_nlinear_workspace *w;
    gsl_multifit_nlinear_fdf fdf;
    gsl_multifit_nlinear_parameters fdf_params =
        gsl_multifit_nlinear_default_parameters();

    gsl_vector *f;
    gsl_matrix *J;
    gsl_matrix *covar = gsl_matrix_alloc(np, np);

    gsl_vector_view params = gsl_vector_view_array(p_init, np);
    gsl_vector_view wts = gsl_vector_view_array(data.w, nd);
    gsl_rng *r;
    double chisq, chisq0;
    int status, info;

    const double xtol = 1e-8;
    const double gtol = 1e-8;
    const double ftol = 1e-8;

    gsl_rng_env_setup();
    r = gsl_rng_alloc(gsl_rng_default);

    /* define the function to be minimized */
    fdf.f = dist_f;
    fdf.df = NULL;   /* set to NULL for finite-difference Jacobian */
    fdf.fvv = NULL;  /* not using geodesic acceleration */
    fdf.n = nd;
    fdf.p = np;
    fdf.params = &data; /* `data' supplied to the penalty function */

    /* allocate workspace with default parameters */
    w = gsl_multifit_nlinear_alloc(T, &fdf_params, nd, np);

    /* initialize solver with starting point and weights */
    gsl_multifit_nlinear_winit(&params.vector, &wts.vector, &fdf, w);

    /* compute initial cost function */
    f = gsl_multifit_nlinear_residual(w);
    gsl_blas_ddot(f, f, &chisq0);

    /* solve the system with a maximum of 100 iterations */
    status = gsl_multifit_nlinear_driver(100, xtol, gtol, ftol,
                                         callback, NULL, &info, w);

    /* compute covariance of best fit parameters */
    J = gsl_multifit_nlinear_jac(w);
    gsl_multifit_nlinear_covar(J, 0.0, covar);

    /* compute final cost */
    gsl_blas_ddot(f, f, &chisq);

#define FIT(i) gsl_vector_get(w->x, i)
#define ERR(i) sqrt(gsl_matrix_get(covar,i,i))

    fprintf(stderr, "summary from method '%s/%s'\n",
            gsl_multifit_nlinear_name(w),
            gsl_multifit_nlinear_trs_name(w));
    fprintf(stderr, "number of iterations: %zu\n",
            gsl_multifit_nlinear_niter(w));
    fprintf(stderr, "function evaluations: %zu\n", fdf.nevalf);
    fprintf(stderr, "Jacobian evaluations: %zu\n", fdf.nevaldf);
    fprintf(stderr, "reason for stopping: %s\n",
            (info == 1) ? "small step size" : "small gradient");
    fprintf(stderr, "initial |f(x)| = %f\n", sqrt(chisq0));
    fprintf(stderr, "final   |f(x)| = %f\n", sqrt(chisq));

    {
        double dof = nd - np;
        double c = GSL_MAX_DBL(1, sqrt(chisq / dof));

        fprintf(stderr, "chisq/dof = %g\n", chisq / dof);

        fprintf(stderr, "a  = %16g +/- %g\n", FIT(0), c*ERR(0));
        fprintf(stderr, "b  = %16g +/- %g\n", FIT(1), c*ERR(1));
        fprintf(stderr, "g  = %16g +/- %g\n", FIT(2), c*ERR(2));
        fprintf(stderr, "xt = %16g +/- %g\n", FIT(3), c*ERR(3));
        fprintf(stderr, "yt = %16g +/- %g\n", FIT(4), c*ERR(4));
        fprintf(stderr, "zt = %16g +/- %g\n", FIT(5), c*ERR(5));
    }

    fprintf(stderr, "status = %s\n", gsl_strerror(status));

    /* compute residual of every point */
    dist_f(w->x, &data, f); // pure distance, no weights.
    printf("#    fid   ftype  distance          x1          y1          z1\n");
    for (int i=0; i<nd; i++) {
        printf("%8zd %7d %9.6f %11.5f %11.5f %11.5f\n", data.fid[i], faces[data.fid[i]].ftype,
               gsl_vector_get(f, i), data.x1[i], data.y1[i], data.z1[i]);
    }

    gsl_multifit_nlinear_free(w);
    gsl_matrix_free(covar);
    gsl_rng_free(r);

    free(faces);
    free(data.fid);
    free(data.x);
    free(data.y);
    free(data.z);
    free(data.w);
    free(data.x1);
    free(data.y1);
    free(data.z1);
    return EXIT_SUCCESS;
}

ext: c

max_stars_repo_name: ymei/BodyRotTransFit

max_issues_repo_name: ymei/BodyRotTransFit

max_forks_repo_name: ymei/BodyRotTransFit


text:
/* -----------------------------------------------------------------------------
 * Copyright 2021 Jonathan Haigh
 * SPDX-License-Identifier: MIT
 * ---------------------------------------------------------------------------*/

#ifndef SQ_INCLUDE_GUARD_system_linux_SqParamSchemaImpl_h_
#define SQ_INCLUDE_GUARD_system_linux_SqParamSchemaImpl_h_

#include "core/typeutil.h"
#include "system/SqParamSchema.gen.h"
#include "system/schema.h"

#include <gsl/gsl>

namespace sq::system::linux {

class SqParamSchemaImpl : public SqParamSchema<SqParamSchemaImpl> {
public:
  explicit SqParamSchemaImpl(const ParamSchema &param_schema);

  SQ_ND Result get_name() const;
  SQ_ND Result get_doc() const;
  SQ_ND Result get_index() const;
  SQ_ND Result get_type() const;
  SQ_ND Result get_required() const;
  SQ_ND Result get_default_value() const;
  SQ_ND Result get_default_value_doc() const;

  SQ_ND Primitive to_primitive() const override;

private:
  gsl::not_null<const ParamSchema *> param_schema_;
};

} // namespace sq::system::linux

#endif // SQ_INCLUDE_GUARD_system_linux_SqParamSchemaImpl_h_

ext: h

max_stars_repo_name: jonathanhaigh/sq

max_issues_repo_name: jonathanhaigh/sq

max_forks_repo_name: jonathanhaigh/sq


text:
/*
 * Copyright (c) 2010-2018 The University of Tennessee and The University
 *                         of Tennessee Research Foundation.  All rights
 *                         reserved.
 * Copyright (c) 2013      Inria. All rights reserved.
 *
 * @precisions normal z -> c d s
 *
 */

#include <lapacke.h>
#include "dplasma.h"
#include "dplasmatypes.h"

#include "map2.h"

static int
dplasma_zlacpy_operator( parsec_execution_stream_t *es,
                         const parsec_tiled_matrix_dc_t *descA,
                         const parsec_tiled_matrix_dc_t *descB,
                         const void *_A, void *_B,
                         PLASMA_enum uplo, int m, int n,
                         void *args )
{
    int tempmm, tempnn, ldam, ldbm;
    const parsec_complex64_t *A = (const parsec_complex64_t*)_A;
    parsec_complex64_t       *B = (parsec_complex64_t*)_B;
    (void)es;
    (void)args;

    tempmm = ((m)==((descA->mt)-1)) ? ((descA->m)-(m*(descA->mb))) : (descA->mb);
    tempnn = ((n)==((descA->nt)-1)) ? ((descA->n)-(n*(descA->nb))) : (descA->nb);
    ldam = BLKLDD( descA, m );
    ldbm = BLKLDD( descB, m );

    LAPACKE_zlacpy_work(
        LAPACK_COL_MAJOR, lapack_const( uplo ), tempmm, tempnn, A, ldam, B, ldbm);

    return 0;
}

/**
 *******************************************************************************
 *
 * @ingroup dplasma_complex64
 *
 * dplasma_zlacpy_New - Generates an taskpool that performs a copy of the matrix A
 * into the matrix B.
 *
 * See dplasma_map2_New() for further information.
 *
 * WARNING: The computations are not done by this call.
 *
 *******************************************************************************
 *
 * @param[in] uplo
 *          Specifies which part of matrix A is copied:
 *          = PlasmaUpperLower: All matrix is referenced.
 *          = PlasmaUpper:      Only upper part is refrenced.
 *          = PlasmaLower:      Only lower part is referenced.
 *
 * @param[in] A
 *          Descriptor of the distributed original matrix A. Any tiled matrix
 *          descriptor can be used. However, if the data is stored in column
 *          major, the tile distribution must match the one of the matrix B.
 *
 * @param[in,out] B
 *          Descriptor of the distributed destination matrix B. Any tiled matrix
 *          descriptor can be used, with no specific storage.
 *
 *******************************************************************************
 *
 * @return
 *          \retval NULL if incorrect parameters are given.
 *          \retval The parsec taskpool describing the operation that can be
 *          enqueued in the runtime with parsec_context_add_taskpool(). It, then, needs to be
 *          destroy with dplasma_zlacpy_Destruct();
 *
 *******************************************************************************
 *
 * @sa dplasma_zlacpy
 * @sa dplasma_zlacpy_Destruct
 * @sa dplasma_clacpy_New
 * @sa dplasma_dlacpy_New
 * @sa dplasma_slacpy_New
 *
 ******************************************************************************/
parsec_taskpool_t*
dplasma_zlacpy_New( PLASMA_enum uplo,
                    const parsec_tiled_matrix_dc_t *A,
                    parsec_tiled_matrix_dc_t *B)
{
    parsec_taskpool_t* tp;

    tp = dplasma_map2_New(uplo, PlasmaNoTrans, A, B,
                          dplasma_zlacpy_operator, NULL );

    return tp;
}

/**
 *******************************************************************************
 *
 * @ingroup dplasma_complex64
 *
 *  dplasma_zlacpy_Destruct - Free the data structure associated to an taskpool
 *  created with dplasma_zlacpy_New().
 *
 *******************************************************************************
 *
 * @param[in,out] taskpool
 *          On entry, the taskpool to destroy.
 *          On exit, the taskpool cannot be used anymore.
 *
 *******************************************************************************
 *
 * @sa dplasma_zlacpy_New
 * @sa dplasma_zlacpy
 *
 ******************************************************************************/
void
dplasma_zlacpy_Destruct( parsec_taskpool_t *tp )
{
    dplasma_map2_Destruct(tp);
}

/**
 *******************************************************************************
 *
 * @ingroup dplasma_complex64
 *
 * dplasma_zlacpy - Generates an taskpool that performs a copy of the matrix A
 * into the matrix B.
 *
 * See dplasma_map2() for further information.
 *
 *******************************************************************************
 *
 * @param[in,out] parsec
 *          The parsec context of the application that will run the operation.
 *
 * @param[in] uplo
 *          Specifies which part of matrix A is copied:
 *          = PlasmaUpperLower: All matrix is referenced.
 *          = PlasmaUpper:      Only upper part is refrenced.
 *          = PlasmaLower:      Only lower part is referenced.
 *
 * @param[in] A
 *          Descriptor of the distributed original matrix A. Any tiled matrix
 *          descriptor can be used. However, if the data is stored in column
 *          major, the tile distribution must match the one of the matrix B.
 *
 * @param[in,out] B
 *          Descriptor of the distributed destination matrix B. Any tiled matrix
 *          descriptor can be used, with no specific storage.
 *
 *******************************************************************************
 *
 * @return
 *          \retval -i if the ith parameters is incorrect.
 *          \retval 0 on success.
 *
 *******************************************************************************
 *
 * @sa dplasma_zlacpy_New
 * @sa dplasma_zlacpy_Destruct
 * @sa dplasma_clacpy
 * @sa dplasma_dlacpy
 * @sa dplasma_slacpy
 *
 ******************************************************************************/
int
dplasma_zlacpy( parsec_context_t *parsec,
                PLASMA_enum uplo,
                const parsec_tiled_matrix_dc_t *A,
                parsec_tiled_matrix_dc_t *B)
{
    parsec_taskpool_t *parsec_zlacpy = NULL;

    if ((uplo != PlasmaUpperLower) &&
        (uplo != PlasmaUpper)      &&
        (uplo != PlasmaLower))
    {
        dplasma_error("dplasma_zlacpy", "illegal value of uplo");
        return -2;
    }

    if ( (A->m > B->m) || (A->n > B->n) ) {
        dplasma_error("dplasma_zlacpy", "illegal matrix A (B is smaller than A)");
        return -3;
    }

    parsec_zlacpy = dplasma_zlacpy_New(uplo, A, B);

    if ( parsec_zlacpy != NULL )
    {
        parsec_context_add_taskpool(parsec, (parsec_taskpool_t*)parsec_zlacpy);
        dplasma_wait_until_completion(parsec);
        dplasma_zlacpy_Destruct( parsec_zlacpy );
    }
    return 0;
}

max_stars_repo_name: NLAFET/ABFT

max_issues_repo_name: NLAFET/ABFT

max_forks_repo_name: NLAFET/ABFT


text:
/*
 * C version of Diffusive Nested Sampling (DNest4) by Brendon J. Brewer
 *
 * Yan-Rong Li, liyanrong@mail.ihep.ac.cn
 * Jun 30, 2016
 *
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <stdbool.h>
#include <gsl/gsl_rng.h>

#include "dnestvars.h"

/* output files */
FILE *fsample, *fsample_info;

/* random number generator */
const gsl_rng_type * dnest_gsl_T;
gsl_rng * dnest_gsl_r;

DNestOptions options;
char options_file[STR_MAX_LENGTH];

// sampler
bool save_to_disk;
double compression;
unsigned int regularisation;

void *particles;
int dnest_size_of_modeltype;
int particle_offset_size, particle_offset_double;

LikelihoodType *log_likelihoods;
unsigned int *level_assignments;

// number account of unaccepted times
unsigned int *account_unaccepts;

int size_levels, size_levels_combine;
Level *levels;
Level *copies_of_levels, *levels_combine;
LikelihoodType *all_above;
unsigned int count_saves, num_saves, num_saves_restart;
int dnest_which_particle_update; // which particle to be updated
int dnest_which_level_update;    // which level to be updated;
unsigned long long int count_mcmc_steps;
LikelihoodType *above;
unsigned int size_above, size_all_above;

double post_logz;
int dnest_num_params;
char dnest_sample_postfix[STR_MAX_LENGTH], dnest_sample_tag[STR_MAX_LENGTH], dnest_sample_dir[STR_MAX_LENGTH];
double *dnest_param_range, *dnest_prior_info;
int *dnest_prior_type;
void *dnest_args;

double *limits, *copies_of_limits;

int dnest_thistask, dnest_totaltask;
int *dnest_perturb_accept;
int dnest_root;

int dnest_flag_restart=0, dnest_flag_postprc=0, dnest_flag_sample_info=0, dnest_flag_limits=0;
double dnest_post_temp=1.0;
char file_restart[STR_MAX_LENGTH], file_save_restart[STR_MAX_LENGTH];

//***********************************************
/*                  functions                  */
double mod(double y, double x);
void wrap_limit(double *x, double min, double max);
int mod_int(int y, int x);
int dnest_cmp(const void *pa, const void *pb);

void options_load(char *optfile, DNestOptions *opts);
void setup(int argc, char** argv, DNestFptrSet *fptrset, int num_params,
           double *param_range, int *prior_type, double *prior_info,
           char *sample_dir, char *optfile, DNestOptions *opts, void *args);
void finalise();

void dnest_run();
void dnest_mcmc_run();
void update_particle(unsigned int which);
void update_level_assignment(unsigned int which);
double log_push(unsigned int which_level);
bool enough_levels(Level *l, int size_l);
void do_bookkeeping();
void save_levels();
void save_particle();
void save_limits();
void kill_lagging_particles();
void renormalise_visits();
void recalculate_log_X();
void dnest_postprocess(double temperature, char *optfile, DNestOptions *opts);
void postprocess(double temperature);
void initialize_output_file();
void close_output_file();
void dnest_save_restart();
void dnest_restart();
void dnest_restart_action(int iflag);
void dnest_accept_action();
void dnest_kill_action(int i, int i_copy);
void dnest_from_prior(void *model);
double dnest_perturb(void *model);
void dnest_print_particle(FILE *fp, const void *model);
void dnest_read_particle(FILE *fp, void *model);
/*=====================================================*/
// users responsible for following functions
void (*print_particle)(FILE *fp, const void *model);
void (*read_particle)(FILE *fp, void *model);
void (*from_prior)(void *model);
double (*log_likelihoods_cal)(const void *model);
double (*log_likelihoods_cal_initial)(const void *model);
double (*log_likelihoods_cal_restart)(const void *model);
double (*perturb)(void *model);
void (*restart_action)(int iflag);
void (*accept_action)();
void (*kill_action)(int i, int i_copy);
/*=====================================================*/

max_stars_repo_name: LiyrAstroph/CDNest


text:
// get some values from a GNU scienctific library random number generator.
// CK 2014
#include <stdio.h>
#include <gsl/gsl_rng.h>

void rantest(long seed,int N,double* out){
    gsl_rng* rng = gsl_rng_alloc (gsl_rng_taus); // initialize a tausworthe rng
    gsl_rng_set(rng, seed); // seed it
    double number = -1;

    for(int i=0; i<N; i++){
        number = gsl_rng_uniform(rng); // sample a number
        printf("%ld  %f\n",seed,number);
        out[i] = number;
    }
    puts("---");

    gsl_rng_free(rng);
}

max_stars_repo_name: cknd/synchrony

max_issues_repo_name: cknd/synchrony

max_forks_repo_name: cknd/synchrony


text:
/* roots/test_funcs.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Reid Priedhorsky, 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 <config.h>
#include <math.h>
#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_roots.h>

#include "test.h"

gsl_function create_function (double (*f)(double, void *))
{
  gsl_function F ;
  F.function = f;
  F.params = 0;
  return F ;
}

gsl_function_fdf create_fdf (double (*f)(double, void *),
                             double (*df)(double, void *),
                             void (*fdf)(double, void *, double *, double *))
{
  gsl_function_fdf FDF ;
  FDF.f = f ;
  FDF.df = df ;
  FDF.fdf = fdf ;
  FDF.params = 0 ;
  return FDF ;
}

/* f(x) = x^{20} - 1 */
/* f'(x) = 20x^{19} */
/* zero at x = 1 or -1 */

double
func1 (double x, void *p)
{
  return pow (x, 20.0) - 1;
}

double
func1_df (double x, void * p)
{
  return 20.0 * pow (x, 19.0);
}

void
func1_fdf (double x, void * p, double *y, double *yprime)
{
  *y = func1 (x, p);
  *yprime = 20.0 * pow (x, 19.0);
}

/* f(x) = sqrt(abs(x))*sgn(x) */
/* f'(x) = 1 / sqrt(abs(x) */
/* zero at x = 0 */
double
func2 (double x, void * p)
{
  double delta;

  if (x > 0)
    delta = 1.0;
  else if (x < 0)
    delta = -1.0;
  else
    delta = 0.0;

  return sqrt (fabs (x)) * delta;
}

double
func2_df (double x, void * p)
{
  return 1 / sqrt (fabs (x));
}

void
func2_fdf (double x, void * p, double *y, double *yprime)
{
  *y = func2 (x, p);
  *yprime = 1 / sqrt (fabs (x));
}

/* f(x) = x^2 - 1e-8 */
/* f'(x) = 2x */
/* zero at x = sqrt(1e-8) or -sqrt(1e-8) */
double
func3 (double x, void * p)
{
  return pow (x, 2.0) - 1e-8;
}

double
func3_df (double x, void * p)
{
  return 2 * x;
}

void
func3_fdf (double x, void * p, double *y, double *yprime)
{
  *y = func3 (x, p);
  *yprime = 2 * x;
}

/* f(x) = x exp(-x) */
/* f'(x) = exp(-x) - x exp(-x) */
/* zero at x = 0 */
double
func4 (double x, void * p)
{
  return x * exp (-x);
}

double
func4_df (double x, void * p)
{
  return exp (-x) - x * exp (-x);
}

void
func4_fdf (double x, void * p, double *y, double *yprime)
{
  *y = func4 (x, p);
  *yprime = exp (-x) - x * exp (-x);
}

/* f(x) = 1/(1+exp(x)) */
/* f'(x) = -exp(x) / (1 + exp(x))^2 */
/* no roots! */
double
func5 (double x, void * p)
{
  return 1 / (1 + exp (x));
}

double
func5_df (double x, void * p)
{
  return -exp (x) / pow (1 + exp (x), 2.0);
}

void
func5_fdf (double x, void * p, double *y, double *yprime)
{
  *y = func5 (x, p);
  *yprime = -exp (x) / pow (1 + exp (x), 2.0);
}

/* f(x) = (x - 1)^7 */
/* f'(x) = 7 * (x - 1)^6 */
/* zero at x = 1 */
double
func6 (double x, void * p)
{
  return pow (x - 1, 7.0);
}

double
func6_df (double x, void * p)
{
  return 7.0 * pow (x - 1, 6.0);
}

void
func6_fdf (double x, void * p, double *y, double *yprime)
{
  *y = func6 (x, p);
  *yprime = 7.0 * pow (x - 1, 6.0);
}

/* sin(x) packaged up nicely. */
double
sin_f (double x, void * p)
{
  return sin (x);
}

double
sin_df (double x, void * p)
{
  return cos (x);
}

void
sin_fdf (double x, void * p, double *y, double *yprime)
{
  *y = sin (x);
  *yprime = cos (x);
}

/* cos(x) packaged up nicely. */
double
cos_f (double x, void * p)
{
  return cos (x);
}

double
cos_df (double x, void * p)
{
  return -sin (x);
}

void
cos_fdf (double x, void * p, double *y, double *yprime)
{
  *y = cos (x);
  *yprime = -sin (x);
}

max_stars_repo_name: ruslankuzmin/julia


text:
/**
 *
 * @file core_zlantr.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 Mathieu Faverge
 * @date 2010-11-15
 * @precisions normal z -> c d s
 *
 **/
#include <lapacke.h>
#include <math.h>
#include "common.h"

#define LAPACKE_CORRECT_DLANTR

/***************************************************************************//**
 *
 * @ingroup dplasma_cores_complex64
 *
 *  CORE_zlantr returns the value
 *
 *     zlantr = ( max(abs(A(i,j))), NORM = PlasmaMaxNorm
 *              (
 *              ( norm1(A),         NORM = PlasmaOneNorm
 *              (
 *              ( normI(A),         NORM = PlasmaInfNorm
 *              (
 *              ( normF(A),         NORM = PlasmaFrobeniusNorm
 *
 *  where norm1 denotes the one norm of a matrix (maximum column sum),
 *  normI denotes the infinity norm of a matrix (maximum row sum) and
 *  normF denotes the Frobenius norm of a matrix (square root of sum
 *  of squares). Note that max(abs(A(i,j))) is not a consistent matrix
 *  norm.
 *
 *******************************************************************************
 *
 * @param[in] norm
 *          = PlasmaMaxNorm: Max norm
 *          = PlasmaOneNorm: One norm
 *          = PlasmaInfNorm: Infinity norm
 *          = PlasmaFrobeniusNorm: Frobenius norm
 *
 * @param[in] uplo
 *          Specifies whether the matrix A is upper triangular or lower triangular:
 *          = PlasmaUpper: Upper triangle of A is stored;
 *          = PlasmaLower: Lower triangle of A is stored.
 *
 * @param[in] diag
 *          Specifies whether or not A is unit triangular:
 *          = PlasmaNonUnit: A is non unit;
 *          = PlasmaUnit:    A us unit.
 *
 * @param[in] M
 *          The number of rows of the matrix A. M >= 0.
 *          If uplo == PlasmaUpper, M <= N. When M = 0, CORE_zlantr returns 0.
 *
 * @param[in] N
 *          The number of columns of the matrix A. N >= 0.
 *          If uplo == PlasmaLower, N <= M. When N = 0, CORE_zlantr returns 0.
 *
 * @param[in] A
 *          The LDA-by-N matrix A.
 *
 * @param[in] LDA
 *          The leading dimension of the array A. LDA >= max(1,M).
 *
 * @param[in,out] work
 *          Array of dimension (MAX(1,LWORK)), where LWORK >= M when norm =
 *          PlasmaInfNorm, or LWORK >= N when norm = PlasmaOneNorm; otherwise,
 *          work is not referenced.
 *
 * @param[out] normA
 *          On exit, normA is the norm of matrix A.
 *
 ******************************************************************************/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_zlantr = PCORE_zlantr
#define CORE_zlantr PCORE_zlantr
#endif
void CORE_zlantr(PLASMA_enum norm, PLASMA_enum uplo, PLASMA_enum diag,
                 int M, int N,
                 const PLASMA_Complex64_t *A, int LDA,
                 double *work, double *normA)
{
#if defined(LAPACKE_CORRECT_DLANTR)
    *normA = LAPACKE_zlantr_work(
        LAPACK_COL_MAJOR,
        lapack_const(norm),
        lapack_const(uplo),
        lapack_const(diag),
        M, N, A, LDA, work);
#else
    const PLASMA_Complex64_t *tmpA;
    double value;
    int i, j, imax;
    int idiag = (diag == PlasmaUnit) ? 1 : 0;

    if ( min(M, N) == 0 ) {
        *normA = 0;
        return;
    }

    switch ( norm ) {
    case PlasmaMaxNorm:
        if ( diag == PlasmaUnit ) {
            *normA = 1.;
        } else {
            *normA = 0.;
        }

        if ( uplo == PlasmaUpper ) {
            M = min(M, N);
            for (j = 0; j < N; j++) {
                tmpA = A+(j*LDA);
                imax = min(j+1-idiag, M);

                for (i = 0; i < imax; i++) {
                    value = cabs( *tmpA );
                    *normA = ( value > *normA ) ? value : *normA;
                    tmpA++;
                }
            }
        } else {
            N = min(M, N);
            for (j = 0; j < N; j++) {
                tmpA = A + j * (LDA+1) + idiag;

                for (i = j+idiag; i < M; i++) {
                    value = cabs( *tmpA );
                    *normA = ( value > *normA ) ? value : *normA;
                    tmpA++;
                }
            }
        }
        break;

    case PlasmaOneNorm:
        CORE_ztrasm( PlasmaColumnwise, uplo, diag, M, N,
                     A, LDA, work );
        if ( uplo == PlasmaLower )
            N = min(M,N);

        *normA = 0;
        for (i = 0; i < N; i++) {
            *normA = ( work[i] > *normA ) ? work[i] : *normA;
        }
        break;

    case PlasmaInfNorm:
        CORE_ztrasm( PlasmaRowwise, uplo, diag, M, N,
                     A, LDA, work );
        if ( uplo == PlasmaUpper )
            M = min(M,N);

        *normA = 0;
        for (i = 0; i < M; i++) {
            *normA = ( work[i] > *normA ) ? work[i] : *normA;
        }
        break;

    case PlasmaFrobeniusNorm:
    {
        double scale = 0.;
        double sumsq = 1.;
        CORE_ztrssq( uplo, diag, M, N,
                     A, LDA, &scale, &sumsq );

        *normA = scale * sqrt( sumsq );
    }
    break;
    default:
        coreblas_error(1, "Illegal value of norm");
        return;
    }
#endif
}

max_stars_repo_name: therault/dplasma

max_issues_repo_name: therault/dplasma

max_forks_repo_name: therault/dplasma


text:
/*
 * Copyright (c) 1997-1999 Massachusetts Institute of Technology
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 */

/* This file was automatically generated --- DO NOT EDIT */
/* Generated on Mon Mar  8 17:44:40 EST 1999 */

#include <fftw-int.h>
#include <fftw.h>

/* Generated by: ./genfft -magic-alignment-check -magic-twiddle-load-all -magic-variables 4 -magic-loopi -notwiddle 9 */

/*
 * This function contains 80 FP additions, 40 FP multiplications,
 * (or, 60 additions, 20 multiplications, 20 fused multiply/add),
 * 30 stack variables, and 36 memory accesses
 */
static const fftw_real K939692620 = FFTW_KONST(+0.939692620785908384054109277324731469936208134);
static const fftw_real K342020143 = FFTW_KONST(+0.342020143325668733044099614682259580763083368);
static const fftw_real K984807753 = FFTW_KONST(+0.984807753012208059366743024589523013670643252);
static const fftw_real K173648177 = FFTW_KONST(+0.173648177666930348851716626769314796000375677);
static const fftw_real K642787609 = FFTW_KONST(+0.642787609686539326322643409907263432907559884);
static const fftw_real K766044443 = FFTW_KONST(+0.766044443118978035202392650555416673935832457);
static const fftw_real K866025403 = FFTW_KONST(+0.866025403784438646763723170752936183471402627);
static const fftw_real K500000000 = FFTW_KONST(+0.500000000000000000000000000000000000000000000);

/*
 * Generator Id's :
 * $Id: exprdag.ml,v 1.36 1999/02/19 17:22:11 athena Exp $
 * $Id: fft.ml,v 1.41 1999/02/19 17:22:13 athena Exp $
 * $Id: to_c.ml,v 1.24 1999/02/19 17:22:17 athena Exp $
 */

void fftw_no_twiddle_9(const fftw_complex *input, fftw_complex *output, int istride, int ostride)
{
     fftw_real tmp5;
     fftw_real tmp17;
     fftw_real tmp50;
     fftw_real tmp20;
     fftw_real tmp78;
     fftw_real tmp53;
     fftw_real tmp10;
     fftw_real tmp26;
     fftw_real tmp58;
     fftw_real tmp74;
     fftw_real tmp31;
     fftw_real tmp59;
     fftw_real tmp15;
     fftw_real tmp37;
     fftw_real tmp62;
     fftw_real tmp75;
     fftw_real tmp42;
     fftw_real tmp61;
     ASSERT_ALIGNED_DOUBLE();
     {
	  fftw_real tmp1;
	  fftw_real tmp2;
	  fftw_real tmp3;
	  fftw_real tmp4;
	  ASSERT_ALIGNED_DOUBLE();
	  tmp1 = c_re(input[0]);
	  tmp2 = c_re(input[3 * istride]);
	  tmp3 = c_re(input[6 * istride]);
	  tmp4 = tmp2 + tmp3;
	  tmp5 = tmp1 + tmp4;
	  tmp17 = tmp1 - (K500000000 * tmp4);
	  tmp50 = K866025403 * (tmp3 - tmp2);
     }
     {
	  fftw_real tmp51;
	  fftw_real tmp18;
	  fftw_real tmp19;
	  fftw_real tmp52;
	  ASSERT_ALIGNED_DOUBLE();
	  tmp51 = c_im(input[0]);
	  tmp18 = c_im(input[3 * istride]);
	  tmp19 = c_im(input[6 * istride]);
	  tmp52 = tmp18 + tmp19;
	  tmp20 = K866025403 * (tmp18 - tmp19);
	  tmp78 = tmp51 + tmp52;
	  tmp53 = tmp51 - (K500000000 * tmp52);
     }
     {
	  fftw_real tmp6;
	  fftw_real tmp28;
	  fftw_real tmp9;
	  fftw_real tmp27;
	  fftw_real tmp25;
	  fftw_real tmp29;
	  fftw_real tmp22;
	  fftw_real tmp30;
	  ASSERT_ALIGNED_DOUBLE();
	  tmp6 = c_re(input[istride]);
	  tmp28 = c_im(input[istride]);
	  {
	       fftw_real tmp7;
	       fftw_real tmp8;
	       fftw_real tmp23;
	       fftw_real tmp24;
	       ASSERT_ALIGNED_DOUBLE();
	       tmp7 = c_re(input[4 * istride]);
	       tmp8 = c_re(input[7 * istride]);
	       tmp9 = tmp7 + tmp8;
	       tmp27 = K866025403 * (tmp8 - tmp7);
	       tmp23 = c_im(input[4 * istride]);
	       tmp24 = c_im(input[7 * istride]);
	       tmp25 = K866025403 * (tmp23 - tmp24);
	       tmp29 = tmp23 + tmp24;
	  }
	  tmp10 = tmp6 + tmp9;
	  tmp22 = tmp6 - (K500000000 * tmp9);
	  tmp26 = tmp22 + tmp25;
	  tmp58 = tmp22 - tmp25;
	  tmp74 = tmp28 + tmp29;
	  tmp30 = tmp28 - (K500000000 * tmp29);
	  tmp31 = tmp27 + tmp30;
	  tmp59 = tmp30 - tmp27;
     }
     {
	  fftw_real tmp11;
	  fftw_real tmp39;
	  fftw_real tmp14;
	  fftw_real tmp38;
	  fftw_real tmp36;
	  fftw_real tmp40;
	  fftw_real tmp33;
	  fftw_real tmp41;
	  ASSERT_ALIGNED_DOUBLE();
	  tmp11 = c_re(input[2 * istride]);
	  tmp39 = c_im(input[2 * istride]);
	  {
	       fftw_real tmp12;
	       fftw_real tmp13;
	       fftw_real tmp34;
	       fftw_real tmp35;
	       ASSERT_ALIGNED_DOUBLE();
	       tmp12 = c_re(input[5 * istride]);
	       tmp13 = c_re(input[8 * istride]);
	       tmp14 = tmp12 + tmp13;
	       tmp38 = K866025403 * (tmp13 - tmp12);
	       tmp34 = c_im(input[5 * istride]);
	       tmp35 = c_im(input[8 * istride]);
	       tmp36 = K866025403 * (tmp34 - tmp35);
	       tmp40 = tmp34 + tmp35;
	  }
	  tmp15 = tmp11 + tmp14;
	  tmp33 = tmp11 - (K500000000 * tmp14);
	  tmp37 = tmp33 + tmp36;
	  tmp62 = tmp33 - tmp36;
	  tmp75 = tmp39 + tmp40;
	  tmp41 = tmp39 - (K500000000 * tmp40);
	  tmp42 = tmp38 + tmp41;
	  tmp61 = tmp41 - tmp38;
     }
     {
	  fftw_real tmp76;
	  fftw_real tmp16;
	  fftw_real tmp73;
	  fftw_real tmp77;
	  fftw_real tmp79;
	  fftw_real tmp80;
	  ASSERT_ALIGNED_DOUBLE();
	  tmp76 = K866025403 * (tmp74 - tmp75);
	  tmp16 = tmp10 + tmp15;
	  tmp73 = tmp5 - (K500000000 * tmp16);
	  c_re(output[0]) = tmp5 + tmp16;
	  c_re(output[3 * ostride]) = tmp73 + tmp76;
	  c_re(output[6 * ostride]) = tmp73 - tmp76;
	  tmp77 = K866025403 * (tmp15 - tmp10);
	  tmp79 = tmp74 + tmp75;
	  tmp80 = tmp78 - (K500000000 * tmp79);
	  c_im(output[3 * ostride]) = tmp77 + tmp80;
	  c_im(output[6 * ostride]) = tmp80 - tmp77;
	  c_im(output[0]) = tmp78 + tmp79;
     }
     {
	  fftw_real tmp21;
	  fftw_real tmp54;
	  fftw_real tmp44;
	  fftw_real tmp49;
	  fftw_real tmp48;
	  fftw_real tmp55;
	  fftw_real tmp45;
	  fftw_real tmp56;
	  ASSERT_ALIGNED_DOUBLE();
	  tmp21 = tmp17 + tmp20;
	  tmp54 = tmp50 + tmp53;
	  {
	       fftw_real tmp32;
	       fftw_real tmp43;
	       fftw_real tmp46;
	       fftw_real tmp47;
	       ASSERT_ALIGNED_DOUBLE();
	       tmp32 = (K766044443 * tmp26) + (K642787609 * tmp31);
	       tmp43 = (K173648177 * tmp37) + (K984807753 * tmp42);
	       tmp44 = tmp32 + tmp43;
	       tmp49 = K866025403 * (tmp43 - tmp32);
	       tmp46 = (K766044443 * tmp31) - (K642787609 * tmp26);
	       tmp47 = (K173648177 * tmp42) - (K984807753 * tmp37);
	       tmp48 = K866025403 * (tmp46 - tmp47);
	       tmp55 = tmp46 + tmp47;
	  }
	  c_re(output[ostride]) = tmp21 + tmp44;
	  tmp45 = tmp21 - (K500000000 * tmp44);
	  c_re(output[7 * ostride]) = tmp45 - tmp48;
	  c_re(output[4 * ostride]) = tmp45 + tmp48;
	  c_im(output[ostride]) = tmp54 + tmp55;
	  tmp56 = tmp54 - (K500000000 * tmp55);
	  c_im(output[4 * ostride]) = tmp49 + tmp56;
	  c_im(output[7 * ostride]) = tmp56 - tmp49;
     }
     {
	  fftw_real tmp57;
	  fftw_real tmp69;
	  fftw_real tmp72;
	  fftw_real tmp66;
	  fftw_real tmp64;
	  fftw_real tmp65;
	  fftw_real tmp71;
	  fftw_real tmp70;
	  ASSERT_ALIGNED_DOUBLE();
	  {
	       fftw_real tmp67;
	       fftw_real tmp68;
	       fftw_real tmp60;
	       fftw_real tmp63;
	       ASSERT_ALIGNED_DOUBLE();
	       tmp57 = tmp17 - tmp20;
	       tmp67 = (K173648177 * tmp59) - (K984807753 * tmp58);
	       tmp68 = (K342020143 * tmp62) + (K939692620 * tmp61);
	       tmp69 = tmp67 - tmp68;
	       tmp72 = K866025403 * (tmp67 + tmp68);
	       tmp66 = tmp53 - tmp50;
	       tmp60 = (K173648177 * tmp58) + (K984807753 * tmp59);
	       tmp63 = (K342020143 * tmp61) - (K939692620 * tmp62);
	       tmp64 = tmp60 + tmp63;
	       tmp65 = K866025403 * (tmp63 - tmp60);
	  }
	  c_re(output[2 * ostride]) = tmp57 + tmp64;
	  tmp71 = tmp57 - (K500000000 * tmp64);
	  c_re(output[8 * ostride]) = tmp71 - tmp72;
	  c_re(output[5 * ostride]) = tmp71 + tmp72;
	  c_im(output[2 * ostride]) = tmp66 + tmp69;
	  tmp70 = tmp66 - (K500000000 * tmp69);
	  c_im(output[5 * ostride]) = tmp65 + tmp70;
	  c_im(output[8 * ostride]) = tmp70 - tmp65;
     }
}

fftw_codelet_desc fftw_no_twiddle_9_desc =
{
     "fftw_no_twiddle_9",
     (void (*)()) fftw_no_twiddle_9,
     9,
     FFTW_FORWARD,
     FFTW_NOTW,
     145,
     0,
     (const int *) 0,
};

max_stars_repo_name: jlprieur/jlplib

max_issues_repo_name: jlprieur/jlplib

max_forks_repo_name: jlprieur/jlplib


text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "compearth.h"
#ifdef COMPEARTH_USE_MKL
#ifdef __clang__
#pragma clang diagnostic push
#pragma clang diagnostic ignored "-Wreserved-id-macro"
#pragma clang diagnostic ignored "-Wstrict-prototypes"
#endif
#include <mkl_cblas.h>
#ifdef __clang__
#pragma clang diagnostic pop
#endif
#else
#include <cblas.h>
#endif

static double compearth_norm64f(const int n, const double *__restrict__ x,
                                const enum ceNormType_enum norm,
                                const double p, int *ierr);
/*!
 * @brief Computes the matrix (Frobenius) norm for a matrix
 *
 * @param[in] n       number of matrices
 * @param[in] M       3 x 3 input matrix as a length 9 array [9*n]
 * @param[in] Lnorm   matrix norm
 *                      TWO_NORM (2)
 *                      ONE_NORM (1)
 *                      P_NORM (in this case must set p)
 *                      INFINITY_NORM
 *                      NEGATIVE_INFINITY_NORM
 * @param[in] p       if using a p-norm this is the value for p (> 0)
 *
 * @param[out] mnorm  matrix norms for each matrix
 *
 * @result 0 indicates success
 *
 */
int compearth_normMat(const int n,
                      const double *__restrict__ M,
                      const enum ceNormType_enum Lnorm,
                      const double p,
                      double *__restrict__ mnorm)
{
    int i, ierr;
    ierr = 0;
    for (i=0; i<n; i++)
    {
        mnorm[i] = compearth_norm64f(9, &M[9*i], Lnorm, p, &ierr);
        if (ierr != 0)
        {
            fprintf(stderr, "%s: Error computing matrix norm!\n", __func__);
            mnorm[i] = 0.0;
            break;
        }
    }
    return ierr;
}

/*!
 * @brief Computes the norm of a vector.  This is from ISTI's ISCL.
 *
 * @param[in] n      Length of array x.
 * @param[in] x      Array of dimension [n] of which to compute norm.
 * @param[in] norm   Type of norm to compute.
 * @param[in] p      If performing a P norm then this must be defined to
 *                   a real number greater than or equal to 1.  Otherwise,
 *                   it will not be accessed.
 *
 * @param[out] ierr  0 indicates success
 *
 * @result P-norm of array x.
 *
 * @author Ben Baker
 *
 * @date August 2017 - redefined variables so that this may exist in
 *       compearth without collisions with ISCL.
 *
 */
static double compearth_norm64f(const int n, const double *__restrict__ x,
                                const enum ceNormType_enum norm,
                                const double p, int *ierr)
{
    double xnorm;
    int i;
    *ierr = 0;
    xnorm = 0.0;
    if (n < 1 || x == NULL)
    {
        if (n < 1){fprintf(stderr, "%s: Error no data\n", __func__);}
        if (x == NULL){fprintf(stderr, "%s: Error x is NULL\n", __func__);}
        *ierr = 1;
        return xnorm;
    }
    // 2 norm
    if (norm == CE_TWO_NORM)
    {
        xnorm = cblas_dnrm2(n, x, 1);
    // 1 norm
    }
    else if (norm == CE_ONE_NORM)
    {
        xnorm = cblas_dasum(n, x, 1);
    }
    // p norm
    else if (norm == CE_P_NORM)
    {
        if (p <= 0.0)
        {
            fprintf(stderr, "%s: Invalid p value %f\n", __func__, p);
            *ierr = 1;
            return xnorm;
        }
        #pragma omp simd reduction(+:xnorm)
        for (i=0; i<n; i++)
        {
            xnorm = xnorm + pow(fabs(x[i]), p);
        }
        xnorm = pow(xnorm, 1./p);
    }
    // infinity norm
    else if (norm == CE_INFINITY_NORM)
    {
        xnorm = fabs(x[cblas_idamax(n, x, 1)]);
    }
    // negative infinity norm
    else if (norm == CE_NEGATIVE_INFINITY_NORM)
    {
        xnorm = fabs(x[0]);
        #pragma omp simd reduction(min:xnorm)
        for (i=1; i<n; i++)
        {
            xnorm = fmin(fabs(x[i]), xnorm);
        }
    }
    // not sure - default to 2-norm
    else
    {
        fprintf(stderr, "%s: Defaulting to 2-norm\n", __func__);
        xnorm = cblas_dnrm2(n, x, 1);
    }
    return xnorm;
}

max_stars_repo_name: OUCyf/mtbeach

max_issues_repo_name: carltape/mtbeach

max_forks_repo_name: carltape/mtbeach


text:
#ifndef H_EASY_GSL
#define H_EASY_GSL

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

#ifdef __cplusplus
extern "C" {
#endif

struct egsl_val {
	gsl_matrix * gslm;
	int cid : 16;
	int index : 16;
};

typedef struct egsl_val val;

/* Core functions */

/* Push a new context. */
void egsl_push();
void egsl_push_named(const char*name);
/* Pops a context */
void egsl_pop();
void egsl_pop_named(const char*name);
void egsl_free(void);

double* egsl_atmp(val v, size_t i, size_t j);
val egsl_alloc(size_t rows, size_t columns);
val egsl_alloc_in_context(int cid, size_t rows, size_t cols);
gsl_matrix * egsl_gslm(val v);
/** Creates a copy of v in the previous context.*/
val egsl_promote(val v);

/** Operations among values */
val egsl_scale(double, val);
val egsl_sum(val, val);
val egsl_sum3(val, val, val);
val egsl_mult(val, val);
val egsl_transpose(val);
val egsl_inverse(val);
val egsl_sub(val,val);
val egsl_sum(val v1,val v2);
val egsl_compose_col(val v1, val v2);
val egsl_compose_row(val v1, val v2);
void egsl_add_to(val v1, val v2);
void egsl_add_to_col(val v1, size_t j, val v2);

double egsl_norm(val);

void egsl_symm_eig(val v, double* eigenvalues, val* eigenvectors);

double egsl_atv(val, size_t i);
double egsl_atm(val, size_t i, size_t j);

/* File: egsl_conversions.c
  Conversions */

val egsl_vFa(size_t rows, const double*);
val egsl_vFda(size_t rows, size_t columns, const double*);

/** Copies a VECTOR value into array */
void egsl_v2a(val, double*);
/** Copies a MATRIX value into array (row1 .. rown) */
void egsl_v2da(val, double*);
/** Copies a vector value into a gsl_vector */
void egsl_v2vec(val, gsl_vector*);

val egsl_vFgslv(const gsl_vector*);
val egsl_vFgslm(const gsl_matrix*);

gsl_matrix* egsl_v2gslm(val);

/*/ File: egsl_misc.c
    Miscellaneous useful matrixes. */
val egsl_zeros(size_t rows, size_t columns);
val egsl_ones(size_t rows, size_t columns);
val egsl_vers(double theta);
val egsl_rot(double theta);

/* Misc */
void egsl_print(const char*str, val);
/** Prints eigenvalues and eigenvectors of a symmetric matrix */
void egsl_print_spectrum(const char*s, val v);
void egsl_print_stats(void);

/** Private implementations things */
void egsl_expect_size(val v, size_t rows, size_t cols);
void egsl_error(void);

#ifdef __cplusplus
}
#endif

#endif


text:
#pragma once

#include <gsl.h>
#include <ImmutableRanges.h>
#include <qobject.h>

#include <mitkBaseData.h>
#include <FaceIdentifier.h>

class QWidget;

namespace crimson
{
class QtPropertyStorage;

/*! \brief   An interface for data associated with model faces.
 *
 * This includes boundary conditions and materials.
 */
class FaceData : public mitk::BaseData
{
public:
    mitkClassMacro(FaceData, BaseData);

    /*!
     * \brief   Gets the FaceData object name.
     */
    virtual gsl::cstring_span<> getName() const = 0;

    /*!
     * \brief   Gets the list of face types (as defined by FaceIdentifer::FaceType) that the face
     *  data can be applied to.
     */
    virtual ImmutableRefRange<FaceIdentifier::FaceType> applicableFaceTypes() = 0;

    ///@{
    /*!
     * \brief   Sets the identifiers of the faces the FaceData is applied to.
     */
    virtual void setFaces(ImmutableRefRange<FaceIdentifier> faces) = 0;

    /*!
    * \brief   Gets the identifiers of the faces the FaceData is applied to.
    */
    virtual ImmutableRefRange<FaceIdentifier> getFaces() const = 0;
    ///@}

    /*!
     * \brief   Creates custom editor widget which should be added to the UI.
     */
    virtual QWidget* createCustomEditorWidget() { return nullptr; }

    /*!
     * \brief   Gets property storage for the FaceData properties.
     */
    virtual gsl::not_null<QtPropertyStorage*> getPropertyStorage() = 0;

	/*!
	* \brief  stores a pointer to C++ owned data object. For use in python-based boundary conditions that have the logic
	implemented in C++
	*/
	virtual void setDataObject(QVariantList data){};

	virtual void setDataUID(std::string dataUID){};

	virtual std::string getDataUID(){ return std::string{}; };

protected:
    FaceData() { mitk::BaseData::InitializeTimeGeometry(1); }
    virtual ~FaceData() {}

    FaceData(const Self&) = default;

    void SetRequestedRegion(const itk::DataObject*) override {}
    void SetRequestedRegionToLargestPossibleRegion() override {}
    bool RequestedRegionIsOutsideOfTheBufferedRegion() override { return true; }
    bool VerifyRequestedRegion() override { return true; }
};

} // end namespace crimson

max_stars_repo_name: carthurs/CRIMSONGUI

max_issues_repo_name: carthurs/CRIMSONGUI

max_forks_repo_name: carthurs/CRIMSONGUI


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

#pragma once

#include "arcana/utils/serialization/base_stream.h"
#include "arcana/utils/serialization/serializable.h"
#include "arcana/analysis/determinator.h"
#include "Proxies/Proxy.h"
#include "arcana/analysis/binary_iterator.h"
#include "Mapping/MapPointKeyframeAssociations.h"
#include "BundleAdjustment/BundleAdjust.h"
#include "Proxies/KeyframeFields.h"

#include "Tracking/PoseEstimator.h"
#include "Utils/cv.h"

#include <type_traits>
#include <memory>
#include <gsl/gsl>

namespace mage
{
    class MapPoint;
    class Keyframe;

    template<typename IterT>
    inline void binary_iterate(const Id<MapPoint>& id, mira::binary_iterator<IterT>& iterator)
    {
        iterator.iterate(reinterpret_cast<const IdT&>(id));
    }

    template<typename IterT>
    inline void binary_iterate(const Id<Keyframe>& id, mira::binary_iterator<IterT>& iterator)
    {
        iterator.iterate(reinterpret_cast<const IdT&>(id));
    }

    /*
    Serialize only the proxy properties that support serialization
    */
    template<typename IterT, typename T, typename ...Props>
    inline void binary_iterate(const Proxy<T, Props...>& proxy, mira::binary_iterator<IterT>& iterator)
    {
        iterator.iterate(static_cast<const proxy::Id<T>&>(proxy));

        int unused[] = {
            0, // in case it's just a Proxy<MapPoint>, we can't declare an empty c-style array
            (iterator.iterate(static_cast<const Props&>(proxy)) , 0)...
        };
        (void)unused;
    };

    template<typename IterT>
    inline void binary_iterate(const KeyframeReprojection& data, mira::binary_iterator<IterT>& itr)
    {
        itr.iterate(reinterpret_cast<const IdT&>(data.KeyframeId));
        itr.iterate(data.Pose);
        for (auto& point : data.Points)
        {
            itr.iterate(point.As<Proxy<MapPoint, proxy::Position, proxy::UpdateStatistics, proxy::ViewingData>>());
            itr.iterate(*point.GetRepresentativeDescriptor());
        }
    }

    template<typename IterT>
    inline void binary_iterate(const MapPointAssociation& data, mira::binary_iterator<IterT>& itr)
    {
        itr.iterate(&data, sizeof(MapPointAssociation));
    }

    template<typename IterT>
    inline void binary_iterate(const HistoricalFrame& data, mira::binary_iterator<IterT>& itr)
    {
        itr.iterate(data.UpdatedPose);
        itr.iterate(data.Keyframe);
    }

    template<typename IterT>
    inline void binary_iterate(const AdjustableData& data, mira::binary_iterator<IterT>& itr)
    {
        itr.iterate(data.Keyframes);
        itr.iterate(data.ExternallyTetheredKeyframes);
        itr.iterate(data.MapPointAssociations);
        itr.iterate(data.MapPoints);
    }

    template<typename IterT>
    inline void binary_iterate(const typename MapPointAssociations<MapPointTrackingProxy>::Association& assoc, mira::binary_iterator<IterT>& itr)
    {
        itr.iterate_all(assoc.MapPoint, assoc.Index);
    }

    template<typename IterT>
    inline void binary_iterate(const MapPointKeyframeAssociations& data, mira::binary_iterator<IterT>& itr)
    {
        itr.iterate(data.MapPoint);
        itr.iterate(data.Keyframes);
    }

    template<typename IterT>
    inline void binary_iterate(const KeyframeAssociation& data, mira::binary_iterator<IterT>& itr)
    {
        itr.iterate_all(data.KeyframeId, data.KeypointDescriptorIndex);
    }

    template<typename IterT>
    inline void binary_iterate(const FrameId& data, mira::binary_iterator<IterT>& itr)
    {
        itr.iterate(data.CorrelationId);
        itr.iterate(data.Camera);
    }

    template<typename IterT>
    inline void binary_iterate(const AnalyzedImage& data, mira::binary_iterator<IterT>& itr)
    {
        itr.iterate(data.GetFrameId());
        itr.iterate(data.GetTimeStamp().time_since_epoch());
        itr.iterate(data.GetImageData());
        itr.iterate(data.GetDistortedCalibration());
        itr.iterate(data.GetUndistortedCalibration());
        itr.iterate(data.GetKeyPoints());
    }

    template<typename IterT, typename T>
    inline void binary_iterate(const Tether<T>& data, mira::binary_iterator<IterT>& itr)
    {
        itr.iterate(data.Type());
        itr.iterate(data.Weight());
        itr.iterate(data.OriginId());

        if (data.Type() == TetherType::DISTANCE)
        {
            itr.iterate(data.Distance());
        }

        if(data.Type() == TetherType::SIX_DOF || data.Type() == TetherType::EXTRINSIC)
        {
            itr.iterate(data.Position());
        }

        if (data.Type() == TetherType::THREE_DOF || data.Type() == TetherType::SIX_DOF || data.Type() == TetherType::EXTRINSIC)
        {
            itr.iterate(data.Rotation());
        }
    }

    namespace proxy
    {
        // no-ops that we don't check for determinism
        inline void binary_iterate(const Image&, mira::binary_iterator<mira::determinator::iterator_callback>&) {}
        inline void binary_iterate(const Associations<MapPointTrackingProxy>&, mira::binary_iterator<mira::determinator::iterator_callback>&) {}

        inline void binary_iterate(const UnAssociatedMask& data, mira::binary_iterator<mira::determinator::iterator_callback>& itr)
        {
            itr.iterate(data.GetUnassociatedKeypointCount());
            itr.iterate(data.GetUnassociatedKeypointMask());
        }

        template<typename IterT>
        inline void binary_iterate(const PoseConstraints& data, mira::binary_iterator<IterT>& itr)
        {
            itr.iterate(data.IsImmortal());
            itr.iterate(data.IsFixed());
            itr.iterate(data.GetTethers());
        }

        template<typename IterT>
        inline void binary_iterate(const Descriptor& data, mira::binary_iterator<IterT>& itr)
        {
            // proxy::Descriptor is just a pointer to the descriptor
            itr.iterate(data.GetRepresentativeDescriptor()->Data(), mage::ORBDescriptor::DESCRIPTOR_SIZE_BYTES);
        }

        template<typename IterT>
        inline void binary_iterate(const Associations<mage::MapPointTrackingProxy>& data, mira::binary_iterator<IterT>& itr)
        {
            // not the real size as we don't account for the bimap overhead
            std::vector<mage::MapPointAssociations<mage::MapPointTrackingProxy>::Association> assocs;
            data.GetAssociations(assocs);
            itr.iterate(assocs);
        }
    }
}

namespace cv
{
    template<typename IterT>
    inline void binary_iterate(const cv::KeyPoint& data, mira::binary_iterator<IterT>& itr)
    {
        itr.iterate(&data, sizeof(data));
    }

    template<typename IterT>
    inline void binary_iterate(const cv::Vec3f& data, mira::binary_iterator<IterT>& itr)
    {
        itr.iterate(&data, sizeof(data));
    }
}

namespace Eigen
{
    template<typename IterT, typename T>
    inline void binary_iterate(const Quaternion<T>& data, mira::binary_iterator<IterT>& itr)
    {
        itr.iterate(data.coeffs().data(), data.coeffs().size()*sizeof(T));
    }
}

max_stars_repo_name: syntheticmagus/mageslam

max_issues_repo_name: syntheticmagus/mageslam

max_forks_repo_name: syntheticmagus/mageslam


text:
#ifndef MSHADOW_TENSOR_BASE_H
#define MSHADOW_TENSOR_BASE_H
/*!
 * \file tensor_base.h
 * \brief definitions of base types, macros functions
 *
 * \author Bing Xu, Tianqi Chen
 */
#include <cmath>
#include <cstdio>
#include <cfloat>
#include <climits>
#include <algorithm>
// macro defintiions

/*!\brief if this macro is define to be 1, mshadow should compile without any of other libs */
#ifndef MSHADOW_STAND_ALONE
    #define MSHADOW_STAND_ALONE 0
#endif

/*! \brief whether do padding during allocation */
#ifndef MSHADOW_ALLOC_PAD
    #define MSHADOW_ALLOC_PAD true
#endif

/*!
 * \brief x dimension of data must be bigger pad_size * ratio to be alloced padded memory, otherwise use tide allocation
 *        for example, if pad_ratio=2, GPU memory alignement size is 32, then we will only allocate padded memory if x dimension > 64
 *        set it to 0 then we will always allocate padded memory
 */
#ifndef MSHADOW_MIN_PAD_RATIO
    #define MSHADOW_MIN_PAD_RATIO 2
#endif

#if MSHADOW_STAND_ALONE
   #define MSHADOW_USE_CBLAS 0
   #define MSHADOW_USE_MKL   0
   #define MSHADOW_USE_CUDA  0
#endif

/*! \brief use CBLAS for CBLAS */
#ifndef MSHADOW_USE_CBLAS
   #define MSHADOW_USE_CBLAS 0
#endif
/*! \brief use MKL for BLAS */
#ifndef MSHADOW_USE_MKL
   #define MSHADOW_USE_MKL   1
#endif
/*! \brief use CUDA support, must ensure that the cuda include path is correct, or directly compile using nvcc */
#ifndef MSHADOW_USE_CUDA
  #define MSHADOW_USE_CUDA   1
#endif
/*! \brief use single precition float */
#ifndef MSHADOW_SINGLE_PRECISION
  #define MSHADOW_SINGLE_PRECISION 1
#endif
/*! \brief whether use SSE */
#ifndef MSHADOW_USE_SSE
  #define MSHADOW_USE_SSE 1
#endif
/*! \brief whether use NVML to get dynamic info */
#ifndef MSHADOW_USE_NVML
  #define MSHADOW_USE_NVML 0
#endif
// SSE is conflict with cudacc
#ifdef __CUDACC__
  #undef MSHADOW_USE_SSE
  #define MSHADOW_USE_SSE 0
#endif

#if MSHADOW_USE_CBLAS
extern "C"{
    #include <cblas.h>
}
#elif MSHADOW_USE_MKL
  #include <mkl.h>
  #include <mkl_cblas.h>
  #include <mkl_vsl.h>
  #include <mkl_vsl_functions.h>
#endif

#if MSHADOW_USE_CUDA
  #include <cublas.h>
  #include <curand.h>
#endif

#if MSHADOW_USE_NVML
  #include <nvml.h>
#endif
// --------------------------------
// MSHADOW_XINLINE is used for inlining template code for both CUDA and CPU code.
#ifdef MSHADOW_XINLINE
  #error "MSHADOW_XINLINE must not be defined"
#endif
#ifdef __CUDACC__
  #define MSHADOW_XINLINE inline __attribute__((always_inline)) __device__ __host__
#else
  #define MSHADOW_XINLINE inline __attribute__((always_inline))
#endif
/*! \brief cpu force inline */
#define MSHADOW_CINLINE inline __attribute__((always_inline))

#if defined(__GXX_EXPERIMENTAL_CXX0X) || defined(__GXX_EXPERIMENTAL_CXX0X__) || __cplusplus >= 201103L
  #define MSHADOW_CONSTEXPR constexpr
#else
  #define MSHADOW_CONSTEXPR const
#endif

/*! \brief namespace for mshadow */
namespace mshadow {
    /*! \brief buffer size for each random number generator */
    const unsigned kRandBufferSize = 1000000;
    /*! \brief pi  */
    const float kPi = 3.1415926f;

#if MSHADOW_SINGLE_PRECISION
    /*! \brief type that will be used for content */
    typedef float real_t;
#else
    typedef double real_t;
#endif
    /*! \brief type that will be used for index */
    typedef unsigned index_t;
}; // namespace mshadow

namespace mshadow {
    /*! \brief namespace for operators */
    namespace op {
        // binary operator
        /*! \brief mul operator */
        struct mul{
            /*! \brief map a, b to result using defined operation */
            MSHADOW_XINLINE static real_t Map(real_t a, real_t b) {
                return a * b;
            }
        };
        /*! \brief plus operator */
        struct plus {
            /*! \brief map a, b to result using defined operation */
            MSHADOW_XINLINE static real_t Map(real_t a, real_t b) {
                return a + b;
            }
        };
        /*! \brief minus operator */
        struct minus {
            /*! \brief map a, b to result using defined operation */
            MSHADOW_XINLINE static real_t Map(real_t a, real_t b) {
                return a - b;
            }
        };
        /*! \brief divide operator */
        struct div {
            /*! \brief map a, b to result using defined operation */
            MSHADOW_XINLINE static real_t Map(real_t a, real_t b) {
                return a / b;
            }
        };
        /*! \brief get rhs */
        struct right {
            /*! \brief map a, b to result using defined operation */
            MSHADOW_XINLINE static real_t Map(real_t a, real_t b) {
                return b;
            }
        };
    }; // namespace op

    /*! \brief namespace for savers */
    namespace sv {
        /*! \brief save to saver: = */
        struct saveto {
            /*! \brief save b to a using save method */
            MSHADOW_XINLINE static void Save(real_t& a, real_t b) {
                a  = b;
            }
            /*! \brief helper constant to use BLAS, alpha */
            MSHADOW_CONSTEXPR static real_t kAlphaBLAS = 1.0f;
            /*! \brief helper constant to use BLAS, beta */
            MSHADOW_CONSTEXPR static real_t kBetaBLAS  = 0.0f;
            /*! \brief corresponding binary operator type */
            typedef op::right OPType;
        };
        /*! \brief save to saver: += */
        struct plusto {
            /*! \brief save b to a using save method */
            MSHADOW_XINLINE static void Save(real_t& a, real_t b) {
                a += b;
            }
            /*! \brief helper constant to use BLAS, alpha */
            MSHADOW_CONSTEXPR static real_t kAlphaBLAS = 1.0f;
            /*! \brief helper constant to use BLAS, beta */
            MSHADOW_CONSTEXPR static real_t kBetaBLAS  = 1.0f;
            /*! \brief corresponding binary operator type */
            typedef op::plus OPType;
        };
        /*! \brief minus to saver: -= */
        struct minusto {
            /*! \brief save b to a using save method */
            MSHADOW_XINLINE static void Save(real_t& a, real_t b) {
                a -= b;
            }
            /*! \brief helper constant to use BLAS, alpha */
            MSHADOW_CONSTEXPR static real_t kAlphaBLAS = -1.0f;
            /*! \brief helper constant to use BLAS, beta */
            MSHADOW_CONSTEXPR static real_t kBetaBLAS  = 1.0f;
            /*! \brief corresponding binary operator type */
            typedef op::minus OPType;
        };
        /*! \brief multiply to saver: *= */
        struct multo {
            /*! \brief save b to a using save method */
            MSHADOW_XINLINE static void Save(real_t& a, real_t b) {
                a *= b;
            }
            /*! \brief corresponding binary operator type */
            typedef op::mul OPType;
        };
        /*! \brief divide to saver: /= */
        struct divto {
            /*! \brief save b to a using save method */
            MSHADOW_XINLINE static void Save(real_t& a, real_t b) {
                a /= b;
            }
            /*! \brief corresponding binary operator type */
            typedef op::div OPType;
        };
    }; // namespace sv

    namespace op {
        // unary operator/ function: example
        // these operators can be defined by user, in the same style as binary and unary operator
        // to use, simply write F<op::identity>( src )
        /*! \brief identity function that maps a real number to it self */
        struct identity{
            /*! \brief map a to result using defined operation */
            MSHADOW_XINLINE static real_t Map(real_t a) {
                return a;
            }
        };
    }; // namespace op

    /*! \brief namespace for potential reducer operations */
    namespace red {
        /*! \brief sum reducer */
        struct sum {
            /*! \brief do reduction into dst */
            MSHADOW_XINLINE static void Reduce( volatile real_t& dst,  volatile real_t src ) {
                dst += src;
            }
            /*! \brief calculate gradient of redres with respect to redsrc,  redres: reduced result, redsrc: one of reduction element */
            MSHADOW_XINLINE static real_t PartialGrad( real_t redres, real_t redsrc ) {
                return 1.0f;
            }
            /*! \brief an intial value of reducer */
            MSHADOW_CONSTEXPR static real_t kInitV = 0.0f;
        };
        /*! \brief maximum reducer */
        struct maximum {
            /*! \brief do reduction into dst */
            MSHADOW_XINLINE static void Reduce( volatile real_t& dst,  volatile real_t src ) {
                using namespace std;
                dst = max( dst, src );
            }
            /*! \brief calculate gradient of redres with respect to redsrc,  redres: reduced result, redsrc: one of reduction element */
            MSHADOW_XINLINE static real_t PartialGrad( real_t redres, real_t redsrc ) {
                return redres == redsrc ? 1.0f: 0.0f;
            }
            /*! \brief an intial value of reducer */
#if MSHADOW_SINGLE_PRECISION
            MSHADOW_CONSTEXPR static real_t kInitV = -FLT_MAX;
#else
            MSHADOW_CONSTEXPR static real_t kInitV = -DBL_MAX;
#endif
        };
    };

    /*! \brief namespace for helper utils of the project */
    namespace utils{
        /*! \brief send error message then exit */
        inline void Error( const char *msg ){
            fprintf( stderr, "Error:%s\n",msg );
            exit( -1 );
        }
        /*! \brief assert a expression is true */
        inline void Assert( bool exp ){
            if( !exp ) Error( "AssertError" );
        }
        /*! \brief assert a expression is true */
        inline void Assert( bool exp, const char *msg ){
            if( !exp ) Error( msg );
        }
        /*! \brief warning */
        inline void Warning( const char *msg ){
            fprintf( stderr, "warning:%s\n",msg );
        }
    }; // namespace utils
}; // namespace mshadow
#endif // TENSOR_BASE_H


text:
/*
 *  read_parameters.h
 *  Spike
 *	parse simple name/value pairs
 *
 *  Created by Ben Evans on 28/11/2008.
 *  Copyright 2008 University of Oxford. All rights reserved.
 *
 */

#ifndef _READ_PARAMETERS_H
#define _READ_PARAMETERS_H

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include <ctype.h>
#include <assert.h>

#include "globals.h"
#include "parameters.h"
#include "utils.h"
#include <gsl/gsl_sf_gamma.h> // N choose M funcion

#define MAXLEN 256
#define BUFFER 512 // Make larger for comments or break up comments
#define VECBUFF 4

char * trim(char * string);
int parse_string(PARAMS * params, char * string);
int read_parameters(PARAMS * params, char * paramfile);
int parseIntVector(char * string, int ** array);
int parseFloatVector(char * string, float ** array);
int printParameters(PARAMS * mp, char * paramfile);
void printIntArray(FILE * fp, char * name, int * array, int len);
void printFloatArray(FILE * fp, char * name, float * array, int len);

#endif


text:
// Copyright (c) Microsoft Corporation. All rights reserved.
// Licensed under the MIT License. See LICENSE in the project root for license information.

#pragma once

#define NOMINMAX

#pragma warning(disable : 6221) // Disable implicit cast warning for C++/WinRT headers (tracked by Bug 17528784: C++/WinRT headers trigger C6221 comparing e.code() to int-typed things)

// Disable factory caching in CppWinRT as the global COM pointers that are released during dll/process
// unload are not safe. Setting this makes CppWinRT just call get_activation_factory directly every time.
#define WINRT_DISABLE_FACTORY_CACHE 1

#include "targetver.h"

#include "BuildMacros.h"

#ifndef WIN32_LEAN_AND_MEAN
#define WIN32_LEAN_AND_MEAN
#endif

#include <windows.h>

#include <wrl\implements.h>
#include <wrl\module.h>
#include <wrl\event.h>

#define MUX_ASSERT(X) _ASSERT(X)
#define MUX_ASSERT_MSG(X, MSG) _ASSERT_EXPR(X, MSG)
#define MUX_ASSERT_NOASSUME(X) _ASSERT(X)

#define MUX_FAIL_FAST() RaiseFailFastException(nullptr, nullptr, 0);
#define MUX_FAIL_FAST_MSG(MSG) RaiseFailFastException(nullptr, nullptr, 0);

#ifdef BUILD_WINDOWS
#include <gsl/gsl_util>
#else
#include <gsl_util.h>
#endif

// windows.ui.xaml.h accesses LoadLibrary in its inline declaration of CreateXamlUiPresenter
// Accessing LoadLibrary is not always allowed (e.g. on phone), so we need to suppress that.
// We can do so by making this define prior to including the header.
#define CREATE_XAML_UI_PRESENTER_API
#include <Windows.UI.Xaml.Hosting.ReferenceTracker.h>

#include <WindowsNumerics.h>

#include <strsafe.h>
#include <robuffer.h>

// STL
#include <vector>
#include <map>
#include <functional>

#define _USE_MATH_DEFINES
#include <math.h>

#ifdef BUILD_WINDOWS
#include <staging.h>
#include <featurestaging-xaml.h>

#undef CATCH_RETURN // See our implementation in ErrorHandling.h

#else
#define WI_IS_FEATURE_PRESENT(FeatureName) 1
#endif

#undef GetCurrentTime

#include "CppWinRTIncludes.h"

max_stars_repo_name: DotNetUX/microsoft-ui-xaml

max_issues_repo_name: agneszitte-nventive/microsoft-ui-xaml

max_forks_repo_name: DotNetUX/microsoft-ui-xaml


text:
/* specfunc/bessel_sequence.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 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.
 */

/* Author:  G. Jungman */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_bessel.h>

#define DYDX_p(p,u,x) (-(p)/(x) + (((nu)*(nu))/((x)*(x))-1.0)*(u))
#define DYDX_u(p,u,x) (p)

static int
rk_step(double nu, double x, double dx, double * Jp, double * J)
{
  double p_0 = *Jp;
  double u_0 = *J;

  double p_1 = dx * DYDX_p(p_0, u_0, x);
  double u_1 = dx * DYDX_u(p_0, u_0, x);

  double p_2 = dx * DYDX_p(p_0 + 0.5*p_1, u_0 + 0.5*u_1, x + 0.5*dx);
  double u_2 = dx * DYDX_u(p_0 + 0.5*p_1, u_0 + 0.5*u_1, x + 0.5*dx);

  double p_3 = dx * DYDX_p(p_0 + 0.5*p_2, u_0 + 0.5*u_2, x + 0.5*dx);
  double u_3 = dx * DYDX_u(p_0 + 0.5*p_2, u_0 + 0.5*u_2, x + 0.5*dx);

  double p_4 = dx * DYDX_p(p_0 + p_3, u_0 + u_3, x + dx);
  double u_4 = dx * DYDX_u(p_0 + p_3, u_0 + u_3, x + dx);

  *Jp = p_0 + p_1/6.0 + p_2/3.0 + p_3/3.0 + p_4/6.0;
  *J  = u_0 + u_1/6.0 + u_2/3.0 + u_3/3.0 + u_4/6.0;

  return GSL_SUCCESS;
}

int
gsl_sf_bessel_sequence_Jnu_e(double nu, gsl_mode_t mode, size_t size, double * v)
{
  /* CHECK_POINTER(v) */

  if(nu < 0.0) {
    GSL_ERROR ("domain error", GSL_EDOM);
  }
  else if(size == 0) {
    GSL_ERROR ("error", GSL_EINVAL);
  }
  else {
    const gsl_prec_t goal   = GSL_MODE_PREC(mode);
    const double dx_array[] = { 0.001, 0.03, 0.1 }; /* double, single, approx */
    const double dx_nominal = dx_array[goal];

    const int cnu = (int) ceil(nu);
    const double nu13 = pow(nu,1.0/3.0);
    const double smalls[] = { 0.01, 0.02, 0.4, 0.7, 1.3, 2.0, 2.5, 3.2, 3.5, 4.5, 6.0 };
    const double x_small = ( nu >= 10.0 ? nu - nu13 : smalls[cnu] );

    gsl_sf_result J0, J1;
    double Jp, J;
    double x;
    size_t i = 0;

    /* Calculate the first point. */
    x = v[0];
    gsl_sf_bessel_Jnu_e(nu, x, &J0);
    v[0] = J0.val;
    ++i;

    /* Step over the idiot case where the
     * first point was actually zero.
     */
    if(x == 0.0) {
      if(v[1] <= x) {
        /* Strict ordering failure. */
        GSL_ERROR ("error", GSL_EFAILED);
      }
      x = v[1];
      gsl_sf_bessel_Jnu_e(nu, x, &J0);
      v[1] = J0.val;
      ++i;
    }

    /* Calculate directly as long as the argument
     * is small. This is necessary because the
     * integration is not very good there.
     */
    while(v[i] < x_small && i < size) {
      if(v[i] <= x) {
        /* Strict ordering failure. */
        GSL_ERROR ("error", GSL_EFAILED);
      }
      x = v[i];
      gsl_sf_bessel_Jnu_e(nu, x, &J0);
      v[i] = J0.val;
      ++i;
    }

    /* At this point we are ready to integrate.
     * The value of x is the last calculated
     * point, which has the value J0; v[i] is
     * the next point we need to calculate. We
     * calculate nu+1 at x as well to get
     * the derivative, then we go forward.
     */
    gsl_sf_bessel_Jnu_e(nu+1.0, x, &J1);
    J  = J0.val;
    Jp = -J1.val + nu/x * J0.val;

    while(i < size) {
      const double dv = v[i] - x;
      const int Nd    = (int) ceil(dv/dx_nominal);
      const double dx = dv / Nd;
      double xj;
      int j;

      if(v[i] <= x) {
        /* Strict ordering failure. */
        GSL_ERROR ("error", GSL_EFAILED);
      }

      /* Integrate over interval up to next sample point.
       */
      for(j=0, xj=x; j<Nd; j++, xj += dx) {
        rk_step(nu, xj, dx, &Jp, &J);
      }

      /* Go to next interval. */
      x = v[i];
      v[i] = J;
      ++i;
    }

    return GSL_SUCCESS;
  }
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
#include <math.h>

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

#include "qdm.h"

#define TAU_RESET_L_SIZE 11
#define TAU_RESET_M_SIZE  8
#define TAU_RESET_H_SIZE 11
#define TAU_RESET_SIZE (TAU_RESET_L_SIZE + TAU_RESET_M_SIZE + TAU_RESET_H_SIZE)

#define TAU_RESET_STEP 0.0001
#define TAU_RESET_GAP  0.001

#define TAU_TABLE_SIZE 1048576

#define TAU_INITIAL_GUESS 0.5
#define TAU_EPS 0.001
#define TAU_ITERATION_MAX 1000

static
void
qdm_tau_reset_setup(
    qdm_tau *t
)
{
  t->reset = gsl_vector_alloc(TAU_RESET_SIZE);

  size_t offset = 0;

  gsl_vector_view view;

  view = gsl_vector_subvector(t->reset, offset, TAU_RESET_L_SIZE);
  qdm_vector_set_seq(&view.vector, t->low, t->low + (TAU_RESET_L_SIZE - 1) * TAU_RESET_STEP);
  offset += TAU_RESET_L_SIZE;

  view = gsl_vector_subvector(t->reset, offset, TAU_RESET_M_SIZE);
  qdm_vector_set_seq(&view.vector, t->low + TAU_RESET_GAP, t->high - TAU_RESET_GAP);
  offset += TAU_RESET_M_SIZE;

  view = gsl_vector_subvector(t->reset, offset, TAU_RESET_H_SIZE);
  qdm_vector_set_seq(&view.vector, t->high - (TAU_RESET_H_SIZE - 1) * TAU_RESET_STEP, t->high);
  offset += TAU_RESET_H_SIZE;
}

static
void
qdm_tau_table_setup(
    qdm_tau *t
)
{
  if (!t->use_table) {
    t->ispline_table = NULL;
    t->mspline_table = NULL;

    return;
  }

  t->ispline_table = gsl_matrix_alloc(TAU_TABLE_SIZE, (t->knots->size - t->spline_df) + 1);
  t->mspline_table = gsl_matrix_alloc(TAU_TABLE_SIZE, (t->knots->size - t->spline_df) + 1);

  double tau = 0;
  gsl_vector_view row;
  for (size_t i = 0; i < TAU_TABLE_SIZE; i++) {
    tau = (double)(i) / TAU_TABLE_SIZE;

    row = gsl_matrix_row(t->ispline_table, i);
    qdm_ispline_vector(
        &row.vector,
        tau,
        t->spline_df,
        t->knots
    );

    row = gsl_matrix_row(t->mspline_table, i);
    qdm_mspline_vector(
        &row.vector,
        tau,
        t->spline_df,
        t->knots
    );
  }
}

qdm_tau *
qdm_tau_alloc(
    int use_table,

    double low,
    double high,

    size_t spline_df,
    const gsl_vector *knots
)
{
  qdm_tau *t = malloc(sizeof(qdm_tau));

  if (use_table == 1) {
    t->use_table = true;
  } else {
    t->use_table = false;
  }

  t->low = low;
  t->high = high;

  t->spline_df = spline_df;
  t->knots = qdm_vector_copy(knots);

  qdm_tau_reset_setup(t);
  qdm_tau_table_setup(t);

  return t;
}

void
qdm_tau_free(
    qdm_tau *t
)
{
  if (t == NULL) {
    return;
  }

  t->use_table = 0;

  t->low = 0;
  t->high = 0;

  t->spline_df = 0;
  t->knots = NULL;

  gsl_vector_free(t->reset);
  t->reset = NULL;

  gsl_matrix_free(t->ispline_table);
  t->ispline_table = NULL;

  gsl_matrix_free(t->mspline_table);
  t->mspline_table = NULL;

  free(t);
}

double
qdm_tau_ispline_mmm(
    const qdm_tau *t,
    const double value,
    const gsl_vector *mmm
)
{
  double result = 0;

  if (t->use_table) {
    size_t i = (size_t)fmin(floor(value * TAU_TABLE_SIZE), TAU_TABLE_SIZE - 1);
    gsl_vector_view ispline = gsl_matrix_row(t->ispline_table, i);

    gsl_blas_ddot(&ispline.vector, mmm, &result);
  } else {
    size_t m = (t->knots->size - t->spline_df) + 1;
    double ispline_data[m];
    gsl_vector_view ispline = gsl_vector_view_array(ispline_data, m);

    qdm_ispline_vector(&ispline.vector, value, t->spline_df, t->knots);
    gsl_blas_ddot(&ispline.vector, mmm, &result);
  }

  return result;
}

double
qdm_tau_mspline_mmm(
    const qdm_tau *t,
    const double value,
    const gsl_vector *mmm
)
{
  double result = 0;

  if (t->use_table) {
    size_t i = (size_t)fmin(floor(value * TAU_TABLE_SIZE), TAU_TABLE_SIZE - 1);
    gsl_vector_view mspline = gsl_matrix_row(t->mspline_table, i);

    gsl_blas_ddot(&mspline.vector, mmm, &result);
  } else {
    size_t m = (t->knots->size - t->spline_df) + 1;
    double mspline_data[m];
    gsl_vector_view mspline = gsl_vector_view_array(mspline_data, m);

    qdm_mspline_vector(&mspline.vector, value, t->spline_df, t->knots);
    gsl_blas_ddot(&mspline.vector, mmm, &result);
  }

  return result;
}

double
qdm_tau_find(
    const qdm_tau *t,
    const double v,
    const gsl_vector *mmm
)
{
  double tau = TAU_INITIAL_GUESS;
  double qi_u = 0;
  double qm_u = 0;

  size_t reset_i = 0;
  for (size_t i = 0; i < TAU_ITERATION_MAX; i++) {
    if (reset_i >= t->reset->size) {
      tau = TAU_INITIAL_GUESS;

      break;
    }

    /* Calculate position... */
    qi_u = qdm_tau_ispline_mmm(t, tau, mmm);

    /* Calculate slope... */
    qm_u = qdm_tau_mspline_mmm(t, tau, mmm);

    /* Check if the update is within our desired interval [0, 1]. */
    double update = tau - (qi_u - v) / qm_u;
    if (update < 0 || update > 1) {
      tau = gsl_vector_get(t->reset, reset_i);
      reset_i++;

      continue;
    } else {
      tau = update;
    }

    if (fabs(qi_u - v) <= TAU_EPS) {
      break;
    }
  }

  return tau;
}

max_stars_repo_name: calebcase/qdm

max_issues_repo_name: calebcase/qdm

max_forks_repo_name: calebcase/qdm


text:
#include <gsl/gsl_linalg.h>

int
main (void)
{
  double a_data[] = { 0.18, 0.60, 0.57, 0.96,
                      0.41, 0.24, 0.99, 0.58,
                      0.14, 0.30, 0.97, 0.66,
                      0.51, 0.13, 0.19, 0.85 } ;

  double b_data[] = { 1.0, 2.0, 3.0, 4.0 };

  gsl_matrix_view m
    = gsl_matrix_view_array(a_data, 4, 4);

  gsl_vector_view b
    = gsl_vector_view_array(b_data, 4);

  gsl_vector *x = gsl_vector_alloc (4);

  int s;

  gsl_permutation * p = gsl_permutation_alloc (4);

  gsl_linalg_LU_decomp (&m.matrix, p, &s);

  gsl_linalg_LU_solve (&m.matrix, p, &b.vector, x);

  gsl_vector_fprintf(stdout, x, "%g");

  gsl_permutation_free (p);
}

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 2021 atframework
// Created by owent

#pragma once

#include <config/atframe_utils_build_feature.h>
#include <config/compiler_features.h>

#include <design_pattern/nomovable.h>
#include <design_pattern/noncopyable.h>

#include <time/time_utility.h>

#include <gsl/select-gsl.h>

#include <list>
#include <vector>

#if defined(LIBATFRAME_UTILS_ENABLE_UNORDERED_MAP_SET) && LIBATFRAME_UTILS_ENABLE_UNORDERED_MAP_SET
#  include <unordered_map>
#  include <unordered_set>
#else
#  include <map>
#  include <set>
#endif

#include "atframe/atapp_conf.h"
#include "atframe/etcdcli/etcd_discovery.h"

namespace atapp {
class app;
class atapp_connection_handle;
class atapp_endpoint;
class atapp_connector_impl;

struct atapp_endpoint_bind_helper {
  // This API is used by inner system and will not be exported, do not call it directly
  static LIBATAPP_MACRO_API_SYMBOL_HIDDEN void unbind(atapp_connection_handle &handle, atapp_endpoint &connect);
  // This API is used by inner system and will not be exported, do not call it directly
  static LIBATAPP_MACRO_API_SYMBOL_HIDDEN void bind(atapp_connection_handle &handle, atapp_endpoint &connect);
};

class atapp_endpoint {
 public:
  using handle_set_t = std::unordered_set<atapp_connection_handle *>;
  using handle_set_iterator = handle_set_t::iterator;
  using handle_set_const_iterator = handle_set_t::const_iterator;
  using ptr_t = std::shared_ptr<atapp_endpoint>;
  using weak_ptr_t = std::weak_ptr<atapp_endpoint>;

  struct pending_message_t {
    util::time::time_utility::raw_time_t expired_timepoint;
    int32_t type;
    uint64_t message_sequence;
    std::vector<unsigned char> data;
    std::unique_ptr<atapp::protocol::atapp_metadata> metadata;
  };

  UTIL_DESIGN_PATTERN_NOCOPYABLE(atapp_endpoint)
  UTIL_DESIGN_PATTERN_NOMOVABLE(atapp_endpoint)

 private:
  struct construct_helper_t {};

 public:
  LIBATAPP_MACRO_API atapp_endpoint(app &owner, construct_helper_t &helper);
  static LIBATAPP_MACRO_API ptr_t create(app &owner);
  LIBATAPP_MACRO_API ~atapp_endpoint();

  LIBATAPP_MACRO_API uint64_t get_id() const noexcept;
  LIBATAPP_MACRO_API gsl::string_view get_name() const noexcept;

  UTIL_FORCEINLINE bool has_connection_handle() const noexcept { return !refer_connections_.empty(); }
  LIBATAPP_MACRO_API const etcd_discovery_node::ptr_t &get_discovery() const noexcept;
  LIBATAPP_MACRO_API void update_discovery(const etcd_discovery_node::ptr_t &discovery) noexcept;

  LIBATAPP_MACRO_API void add_connection_handle(atapp_connection_handle &handle);
  LIBATAPP_MACRO_API void remove_connection_handle(atapp_connection_handle &handle);
  LIBATAPP_MACRO_API atapp_connection_handle *get_ready_connection_handle() const noexcept;

  LIBATAPP_MACRO_API int32_t push_forward_message(int32_t type, uint64_t &msg_sequence, const void *data,
                                                  size_t data_size, const atapp::protocol::atapp_metadata *metadata);

  LIBATAPP_MACRO_API int32_t retry_pending_messages(const util::time::time_utility::raw_time_t &tick_time,
                                                    int32_t max_count = 0);
  LIBATAPP_MACRO_API void add_waker(util::time::time_utility::raw_time_t wakeup_time);

  UTIL_FORCEINLINE app *get_owner() const noexcept { return owner_; }

  LIBATAPP_MACRO_API size_t get_pending_message_count() const noexcept;
  LIBATAPP_MACRO_API size_t get_pending_message_size() const noexcept;

 private:
  void reset();
  void cancel_pending_messages();

  void trigger_on_receive_forward_response(atapp_connector_impl *connector, atapp_connection_handle *handle,
                                           int32_t type, uint64_t sequence, int32_t error_code, const void *data,
                                           size_t data_size, const atapp::protocol::atapp_metadata *metadata);

 private:
  bool closing_;
  app *owner_;
  util::time::time_utility::raw_time_t nearest_waker_;
  weak_ptr_t watcher_;
  handle_set_t refer_connections_;
  etcd_discovery_node::ptr_t discovery_;

  std::list<pending_message_t> pending_message_;
  size_t pending_message_size_;
#if defined(LIBATAPP_ENABLE_CUSTOM_COUNT_FOR_STD_LIST) && LIBATAPP_ENABLE_CUSTOM_COUNT_FOR_STD_LIST
  size_t pending_message_count_;
#endif

  friend struct atapp_endpoint_bind_helper;
};
}  // namespace atapp

max_stars_repo_name: atframework/libatapp

max_issues_repo_name: atframework/libatapp

max_forks_repo_name: atframework/libatapp