text:
static const char help[] =
"Solves a 2D dam-saturation problem.  No options.\n"
"The exact soluution is not known, but\n"
"a coarse-grid discrete solution can be checked against that source.\n"
"Note Poisson2DFunctionLocal() sets-up this unconstrained problem:\n"
"    - u_xx - u_yy = - 1\n"
"while we want this complementarity problem:\n"
"    F(u) = - u_xx - u_yy + 1 >= 0\n"
"    u >= 0\n"
"    u F(u) = 0.\n"
"As with obstacle.c, this is solved (default) by -snes_type vinewtonrsls.\n"
"Reference:  pages 667-668 of Brandt & Cryer (1983).\n\n";

/*
note PFAS is not implemented in PETSc, but the following runs for X = 1,2,3,4,5,6
quickly solve the same problems as in Brandt & Cryer Table 4.2:
   s ./dam -snes_monitor -pc_type mg -snes_grid_sequence X

on ed-galago I can go up to X = 11 giving 4097 x 6145 grid and serial runtime of 297 seconds
(the next step runs out of memory)

on parallel I am getting error messages with vinewtonrsls + mg:
    mpiexec -n 4 ./dam -snes_converged_reason -snes_grid_sequence 5 -snes_type vinewtonrsls -pc_type mg
but this works with either vinewtonssls or another PC (gamg or bjacobi+ilu or asm+lu or etc.)

run as:
    ./dam -da_refine 1 -snes_view_solution :foo.m:ascii_matlab -snes_converged_reason -snes_rtol 1.0e-12
Nonlinear solve converged due to CONVERGED_FNORM_RELATIVE iterations 4
done on   5 x   7 grid

then in Matlab/Octave you can compare to Brandt & Cryer Table 4.1:
>> foo
>> format long g
>> u = flipud(reshape(Vec_0x84000000_0,5,7)')    % use loaded name here
u =
           0                     0                     0                     0             0
           8      2.53716015237691                     0                     0             0
          32       18.148640609503      6.78414305345868                     0             0
          72      47.2732592321936      24.9879316043211      7.91201621181146             0
         128      89.9564647149903       53.982307919772      22.6601332428759             0
         200      146.570291707977      94.3247021169195      44.7462088399388             0
         288                   218                   148                    78             8

regarding computing seepage face height, runs
    ./dam -snes_converged_reason -pc_type mg -snes_rtol 1.0e-8 -snes_grid_sequence X
gives:
    X    grid       height
    6    129x193    8.8750000
    7    257x385    8.7500000
    8    513x769    8.7500000
    9    1025x1537  8.7343750
    10   2049x3073  8.7187500
    11   4097x6145  8.7109375
but note these numbers depend (at the second digit, even) on the value of
"wetthreshold" in GetSeepageFaceHeight()
*/

#include <petsc.h>
#include "../../ch6/poissonfunctions.h"

typedef struct {
    double  a, y1, y2;
} DamCtx;

double g_fcn(double x, double y, double z, void *ctx) {
    PoissonCtx *user = (PoissonCtx*)ctx;
    DamCtx     *dctx = (DamCtx*)(user->addctx);
    const double a = dctx->a,
                 y1 = dctx->y1,
                 y2 = dctx->y2,
                 tol = a * 1.0e-8;
    // see (4.2) in Brandt & Cryer for following:
    if (x < tol) {
       return (y1 - y) * (y1 - y) / 2.0;                      // AB
    } else if (x > a - tol) {
       if (y < y2) {
           return (y2 - y) * (y2 - y) / 2.0;                  // CD
       } else {
           return 0.0;                                        // DF
       }
    } else if (y < tol) {
       return (y1 * y1 * (a - x) + y2 * y2 * x) / (2.0 * a);  // BC
    } else if (y > y1 - tol) {
       return 0.0;                                            // FA
    } else {
       return NAN;
    }
}

double f_fcn(double x, double y, double z, void *ctx) {
    return -1.0;
}

extern PetscErrorCode FormBounds(SNES, Vec, Vec);
extern PetscErrorCode GetSeepageFaceHeight(DMDALocalInfo*, Vec, double*, DamCtx*);

int main(int argc,char **argv) {
  PetscErrorCode ierr;
  DM             da, da_after;
  SNES           snes;
  Vec            u;
  PoissonCtx     user;
  DamCtx         dctx;
  DMDALocalInfo  info;
  double         height;

  PetscInitialize(&argc,&argv,NULL,help);

  dctx.a  = 16.0;  // a, y1, y2 from Brandt & Cryer
  dctx.y1 = 24.0;
  dctx.y2 = 4.0;
  user.cx = 1.0;
  user.cy = 1.0;
  user.cz = 1.0;
  user.g_bdry = &g_fcn;
  user.f_rhs = &f_fcn;
  user.addctx = &dctx;

  ierr = DMDACreate2d(PETSC_COMM_WORLD,
      DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR,
      3,4,                       // override with -da_refine or -da_grid_x,_y
      PETSC_DECIDE,PETSC_DECIDE, // num of procs in each dim
      1,1,NULL,NULL,             // dof = 1 and stencil width = 1
      &da);CHKERRQ(ierr);
  ierr = DMSetFromOptions(da); CHKERRQ(ierr);
  ierr = DMSetUp(da); CHKERRQ(ierr);
  ierr = DMDASetUniformCoordinates(da,0.0,dctx.a,0.0,dctx.y1,-1.0,-1.0);CHKERRQ(ierr);
  ierr = DMSetApplicationContext(da,&user);CHKERRQ(ierr);

  ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr);
  ierr = SNESSetDM(snes,da);CHKERRQ(ierr);
  ierr = SNESSetApplicationContext(snes,&user);CHKERRQ(ierr);

  ierr = SNESSetType(snes,SNESVINEWTONRSLS);CHKERRQ(ierr);
  ierr = SNESVISetComputeVariableBounds(snes,&FormBounds);CHKERRQ(ierr);

  ierr = DMDASNESSetFunctionLocal(da,INSERT_VALUES,
             (DMDASNESFunction)Poisson2DFunctionLocal,&user); CHKERRQ(ierr);
  ierr = DMDASNESSetJacobianLocal(da,
             (DMDASNESJacobian)Poisson2DJacobianLocal,&user); CHKERRQ(ierr);
  ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);

  ierr = DMCreateGlobalVector(da,&u);CHKERRQ(ierr);
  // initial iterate has u=g on boundary and u=0 in interior
  ierr = InitialState(da, ZEROS, PETSC_TRUE, u, &user); CHKERRQ(ierr);

  /* solve */
  ierr = SNESSolve(snes,NULL,u);CHKERRQ(ierr);
  ierr = VecDestroy(&u); CHKERRQ(ierr);
  ierr = DMDestroy(&da); CHKERRQ(ierr);

  // report seepage face
  ierr = SNESGetSolution(snes,&u); CHKERRQ(ierr); /* do not destroy u */
  ierr = SNESGetDM(snes,&da_after); CHKERRQ(ierr);
  ierr = DMDAGetLocalInfo(da_after,&info); CHKERRQ(ierr);
  ierr = GetSeepageFaceHeight(&info,u,&height,&dctx); CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD,
      "done on %3d x %3d grid; computed seepage face height = %.7f\n",
      info.mx,info.my,height); CHKERRQ(ierr);

  SNESDestroy(&snes);
  return PetscFinalize();
}

// for call-back: tell SNESVI we want  0 <= u < +infinity
PetscErrorCode FormBounds(SNES snes, Vec Xl, Vec Xu) {
    PetscErrorCode ierr;
    ierr = VecSet(Xl,0.0);CHKERRQ(ierr);
    ierr = VecSet(Xu,PETSC_INFINITY);CHKERRQ(ierr);
    return 0;
}

PetscErrorCode GetSeepageFaceHeight(DMDALocalInfo *info, Vec u, double *height, DamCtx *dctx) {
    PetscErrorCode ierr;
    MPI_Comm       comm;
    const double   dy = dctx->y1 / (PetscReal)(info->my-1),
                   wetthreshhold = 1.0e-6;  // what does "u>0" mean?
    int            j;
    double         **au, locwetmax = - PETSC_INFINITY;
    ierr = DMDAVecGetArrayRead(info->da,u,&au); CHKERRQ(ierr);
    if (info->xs+info->xm == info->mx) { // do we even own (part of) the x=a side of the rectangle?
        for (j=info->ys; j<info->ys+info->ym; j++) {
           if (au[j][info->mx-2] > wetthreshhold) {  // is the first inter point wet?
               locwetmax = PetscMax(j*dy,locwetmax);
           }
        }
    }
    ierr = DMDAVecRestoreArrayRead(info->da,u,&au); CHKERRQ(ierr);
    ierr = PetscObjectGetComm((PetscObject)(info->da),&comm); CHKERRQ(ierr);
    ierr = MPI_Allreduce(&locwetmax,height,1,MPI_DOUBLE,MPI_MAX,comm); CHKERRQ(ierr);
    *height -= dctx->y2;   // height is segment ED in figure
    return 0;
}


text:
#ifndef PyGSL_RNG_H
#define PyGSL_RNG_H 1

#include <pygsl/intern.h>
#include <gsl/gsl_rng.h>
typedef struct {
  PyObject_HEAD
  gsl_rng * rng;
} PyGSL_rng;

/*
 * Get a gsl_rng object from a PyGSL rng wrapper.
 */
PyGSL_API_EXTERN gsl_rng *
PyGSL_gsl_rng_from_pyobject(PyObject * object);
#ifndef _PyGSL_API_MODULE
#define PyGSL_gsl_rng_from_pyobject \
(*(gsl_rng *  (*) (PyObject *)) PyGSL_API[PyGSL_gsl_rng_from_pyobject_NUM])
#endif /* _PyGSL_API_MODULE */

#define PyGSL_RNG_Check(op) \
   ((op)->ob_type == (PyTypeObject *)PyGSL_API[PyGSL_RNG_ObjectType_NUM])

#define import_pygsl_rng() \
{ \
   PyObject *pygsl = NULL, *c_api = NULL, *md = NULL; \
   if ( \
      (pygsl = PyImport_ImportModule("pygsl.rng"))         != NULL && \
      (md = PyModule_GetDict(pygsl))                       != NULL && \
      (c_api = PyDict_GetItemString(md, "_PYGSL_RNG_API")) != NULL && \
      (PyCObject_Check(c_api))                                        \
     ) { \
	 PyGSL_API = (void **)PyCObject_AsVoidPtr(c_api); \
   } else { \
        PyGSL_API = NULL; \
   } \
   /* fprintf(stderr, "PyGSL_API points to %p\n", (void *) PyGSL_API); */ \
}
#endif  /* PyGSL_RNG_H */

ext: h

max_stars_repo_name: juhnowski/FishingRod

max_issues_repo_name: juhnowski/FishingRod

max_forks_repo_name: juhnowski/FishingRod


text:
#pragma once

#include <functional>
#include <memory>
#include <vector>

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

#include <QObject>

#include "HypCommands.h" // Hyperion::RECORD_LENGTH

class DeviceLink;
class Recording;

class HypReader : public QObject
{
    Q_OBJECT
public:
    explicit HypReader(std::function<void(const QString&)> msgCBck,
                       std::function<void(size_t)> finishedCBck,
                       QObject* parent = nullptr);

    void DownloadFromDevice();
    void EraseDevice();

    bool LoadFromFile(const QString& filePath);
    bool SaveToFile  (const QString& filePath);

    size_t GetNumRecordings() const; // { return m_recordings.size(); }
    const Recording& GetRecording(size_t idx) const;

    Hyperion::DeviceType GetDeviceType() const;
    int         GetFirmwareVersionX100() const;

signals:
    void SeriesEnded();
    void DownloadFinish(bool success);

private:
    void EndSeries();
    void FinishDownload(bool success);

    // These get called from DeviceLink worker thread.
    // They emit the data to be processed on UI thread.
    void MarkSeriesEnd();
    void DownloadFinished(bool success);
    void ReceiveDataChunk(const gsl::span<uint8_t>& data);

//data:
    std::function<void(const QString&)> m_msgCBck;
    std::function<void(size_t)>         m_finishedCBck;
    std::unique_ptr<DeviceLink>         m_deviceReader;

    std::vector<std::vector<uint8_t>> m_downloadedRawData;
    std::vector<Recording>            m_recordings;

public:
    static const QString BS;
};


text:
/* histogram/gsl_histogram2d.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 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 2 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#ifndef __GSL_HISTOGRAM2D_H__
#define __GSL_HISTOGRAM2D_H__

#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_types.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

typedef struct {
  size_t nx, ny ;
  double * xrange ;
  double * yrange ;
  double * bin ;
} gsl_histogram2d ;

typedef struct {
  size_t nx, ny ;
  double * xrange ;
  double * yrange ;
  double * sum ;
} gsl_histogram2d_pdf ;

GSL_EXPORT gsl_histogram2d * gsl_histogram2d_alloc (const size_t nx, const size_t ny);
GSL_EXPORT gsl_histogram2d * gsl_histogram2d_calloc (const size_t nx, const size_t ny);
GSL_EXPORT gsl_histogram2d * gsl_histogram2d_calloc_uniform (const size_t nx, const size_t ny,
                                                             const double xmin, const double xmax,
                                                             const double ymin, const double ymax);

GSL_EXPORT void gsl_histogram2d_free (gsl_histogram2d * h);

GSL_EXPORT int gsl_histogram2d_increment (gsl_histogram2d * h, double x, double y);
GSL_EXPORT int gsl_histogram2d_accumulate (gsl_histogram2d * h,
                                           double x, double y, double weight);
GSL_EXPORT int gsl_histogram2d_find (const gsl_histogram2d * h,
                                     const double x, const double y, size_t * i, size_t * j);

GSL_EXPORT double gsl_histogram2d_get (const gsl_histogram2d * h, const size_t i, const size_t j);
GSL_EXPORT int gsl_histogram2d_get_xrange (const gsl_histogram2d * h, const size_t i,
                                           double * xlower, double * xupper);
GSL_EXPORT int gsl_histogram2d_get_yrange (const gsl_histogram2d * h, const size_t j,
                                           double * ylower, double * yupper);

GSL_EXPORT double gsl_histogram2d_xmax (const gsl_histogram2d * h);
GSL_EXPORT double gsl_histogram2d_xmin (const gsl_histogram2d * h);
GSL_EXPORT size_t gsl_histogram2d_nx (const gsl_histogram2d * h);

GSL_EXPORT double gsl_histogram2d_ymax (const gsl_histogram2d * h);
GSL_EXPORT double gsl_histogram2d_ymin (const gsl_histogram2d * h);
GSL_EXPORT size_t gsl_histogram2d_ny (const gsl_histogram2d * h);

GSL_EXPORT void gsl_histogram2d_reset (gsl_histogram2d * h);

GSL_EXPORT
gsl_histogram2d *
gsl_histogram2d_calloc_range(size_t nx, size_t ny,
                             double *xrange, double *yrange);

GSL_EXPORT
int
gsl_histogram2d_set_ranges_uniform (gsl_histogram2d * h,
                                    double xmin, double xmax,
                                    double ymin, double ymax);

GSL_EXPORT
int
gsl_histogram2d_set_ranges (gsl_histogram2d * h,
                            const double xrange[], size_t xsize,
                            const double yrange[], size_t ysize);

GSL_EXPORT
int
gsl_histogram2d_memcpy(gsl_histogram2d *dest, const gsl_histogram2d *source);

GSL_EXPORT
gsl_histogram2d *
gsl_histogram2d_clone(const gsl_histogram2d * source);

GSL_EXPORT
double
gsl_histogram2d_max_val(const gsl_histogram2d *h);

GSL_EXPORT
void
gsl_histogram2d_max_bin (const gsl_histogram2d *h, size_t *i, size_t *j);

GSL_EXPORT
double
gsl_histogram2d_min_val(const gsl_histogram2d *h);

GSL_EXPORT
void
gsl_histogram2d_min_bin (const gsl_histogram2d *h, size_t *i, size_t *j);

GSL_EXPORT
double
gsl_histogram2d_xmean (const gsl_histogram2d * h);

GSL_EXPORT
double
gsl_histogram2d_ymean (const gsl_histogram2d * h);

GSL_EXPORT
double
gsl_histogram2d_xsigma (const gsl_histogram2d * h);

GSL_EXPORT
double
gsl_histogram2d_ysigma (const gsl_histogram2d * h);

GSL_EXPORT
double
gsl_histogram2d_cov (const gsl_histogram2d * h);

GSL_EXPORT
double
gsl_histogram2d_sum (const gsl_histogram2d *h);

GSL_EXPORT
int
gsl_histogram2d_equal_bins_p(const gsl_histogram2d *h1,
                             const gsl_histogram2d *h2) ;

GSL_EXPORT
int
gsl_histogram2d_add(gsl_histogram2d *h1, const gsl_histogram2d *h2);

GSL_EXPORT
int
gsl_histogram2d_sub(gsl_histogram2d *h1, const gsl_histogram2d *h2);

GSL_EXPORT
int
gsl_histogram2d_mul(gsl_histogram2d *h1, const gsl_histogram2d *h2);

GSL_EXPORT
int
gsl_histogram2d_div(gsl_histogram2d *h1, const gsl_histogram2d *h2);

GSL_EXPORT
int
gsl_histogram2d_scale(gsl_histogram2d *h, double scale);

GSL_EXPORT
int
gsl_histogram2d_shift(gsl_histogram2d *h, double shift);

GSL_EXPORT int gsl_histogram2d_fwrite (FILE * stream, const gsl_histogram2d * h) ;
GSL_EXPORT int gsl_histogram2d_fread (FILE * stream, gsl_histogram2d * h);
GSL_EXPORT int gsl_histogram2d_fprintf (FILE * stream, const gsl_histogram2d * h,
                                        const char * range_format,
                                        const char * bin_format);
GSL_EXPORT int gsl_histogram2d_fscanf (FILE * stream, gsl_histogram2d * h);

GSL_EXPORT gsl_histogram2d_pdf * gsl_histogram2d_pdf_alloc (const size_t nx, const size_t ny);
GSL_EXPORT int gsl_histogram2d_pdf_init (gsl_histogram2d_pdf * p, const gsl_histogram2d * h);
GSL_EXPORT void gsl_histogram2d_pdf_free (gsl_histogram2d_pdf * p);
GSL_EXPORT int gsl_histogram2d_pdf_sample (const gsl_histogram2d_pdf * p,
                                           double r1, double r2,
                                           double * x, double * y);

__END_DECLS

#endif /* __GSL_HISTOGRAM2D_H__ */

max_stars_repo_name: dynaryu/vaws

max_issues_repo_name: dynaryu/vaws

max_forks_repo_name: dynaryu/vaws


text:
#pragma once

#include <petsc.h>
#include <memory>

#include "partition/rank_subset.h"

namespace TimeSteppingScheme
{

namespace NestedMatVecUtility
{

//! from a Petsc Vec with nested type (nestedVec) create a new Petsc Vec (singleVec) that contains all values at once. If the singleVec already exists, do not create again, only copy the values.
void createVecFromNestedVec(Vec nestedVec, Vec &singleVec, std::shared_ptr<Partition::RankSubset> rankSubset);

//! copy the values from a singleVec back to the nested Petsc Vec (nestedVec)
void fillNestedVec(Vec singleVec, Vec nestedVec);

//! from a Petsc Mat with nested type (nestedMat) create a new Petsc Mat (singleMat) that contains all values at once. If the singleMat already exists, do not create again, only copy the values.
void createMatFromNestedMat(Mat nestedMat, Mat &singleMat, std::shared_ptr<Partition::RankSubset> rankSubset);

}
}  // namespace

max_stars_repo_name: maierbn/opendihu

max_issues_repo_name: maierbn/opendihu

max_forks_repo_name: maierbn/opendihu


text:
/*
 * BIG_DATA.c
 *
 *  Created on: 3 Sep 2020
 *      Author: heine
 */

#include <gsl/gsl_sf_bessel.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_linalg.h>
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <time.h>

#include "kalman_filter.h"
#include "particle_filters.h"

#define N_SETTINGS 10
#define N_SCALES 11

void init_with_zeros(double* array, long length);
void scale_sample_sizes(double scaler, long sample_size_collection[N_SETTINGS][2], int n_settings);
void particle_allocation_time_matching(int bpf_sample_size, double std_sig, double obs_std, int Nobs, int data_length, double x0, double *y, double *Rinv1, double *Rinv0, gsl_matrix *R);

int main(void) {

  gsl_rng * rng = gsl_rng_alloc(gsl_rng_taus);
  gsl_rng_set(rng,clock());

  double scaler;
  int N_top = 50;
  int N_bpf = 250;
  int run_time_matching = 0;
  int data_length = 50;
  int Nobs = 500;
  double sig_std = .1;
  double obs_std = 1;
  double x0 = 0;
  double p0 = sig_std * sig_std;
  double sample_size_coeff = (double)68000 / (double) 250;

  /* -----------------------------------------------------------------------------------
   *
   * Create data
   *
   -----------------------------------------------------------------------------------*/
  double *x = (double*) malloc(data_length * sizeof(double));
  double *y = (double*) malloc(data_length * Nobs * sizeof(double));
  double *noise = (double*) malloc(Nobs * sizeof(double));
  gsl_matrix * A;

  double x_prev = 0;

  /* Create the measurement covariance square root */
  A = gsl_matrix_alloc(Nobs, Nobs);
  for (int i = 0; i < Nobs; i++) {
    for (int j = 0; j < Nobs; j++) {
      /* Uncomment for identity matrix */
      /*
       *       if (i == j) {
       *         gsl_matrix_set(A, i, j, 1);
       *       } else {
       *         gsl_matrix_set(A, i, j, 0);
       } */
//      gsl_matrix_set(A, i, j, gsl_ran_gaussian(rng, 1));
      gsl_matrix_set(A, i, j, gsl_rng_uniform(rng));
    }
  }

  /* Calculate the observation covariance */
  gsl_matrix * R = gsl_matrix_alloc(Nobs, Nobs);
  double tmp;
  double decay = 2; // Decay of correlations
  FILE * R_out = fopen("R_matrix.txt", "w");
  for (int i = 0; i < Nobs; i++) {
    for (int j = 0; j < Nobs; j++) {
      tmp = 0;
      for (int k = 0; k < Nobs; k++) {
        tmp += gsl_matrix_get(A, i, k) * gsl_matrix_get(A, j, k);
      }
      tmp = obs_std * obs_std * exp(-decay * (double) abs(i - j)) * tmp;
      gsl_matrix_set(R, i, j, tmp);
      fprintf(R_out, "%e ", tmp);
    }
    fprintf(R_out, "\n");
  }
  fclose(R_out);

  // Copy R to A
  for (int i = 0; i < Nobs; i++) {
    for (int j = 0; j < Nobs; j++) {
      gsl_matrix_set(A, i, j, gsl_matrix_get(R, i, j));
    }
  }
  gsl_linalg_cholesky_decomp1(A);

  printf("Obsevation covariance written in R_matrix.txt\n");

  /* Main loop for generating the signal and the data */
  FILE *data_file = fopen("data.txt", "w");
  for (int n = 0; n < data_length; n++) {

    x[n] = x_prev + gsl_ran_gaussian(rng, sig_std);
    x_prev = x[n];
    fprintf(data_file, "%i %e", n, x[n]);

    /* Uncorrelated observation noise */
    for (int i = 0; i < Nobs; i++) {
      noise[i] = gsl_ran_gaussian(rng, obs_std);
    }
    /* Make the noise correlated by multiplying by A */
    for (int i = 0; i < Nobs; i++) {
      y[n * Nobs + i] = x[n];
      for (int j = 0; j <= i; j++)
        y[n * Nobs + i] += gsl_matrix_get(A, i, j) * noise[j];
      fprintf(data_file, " %e,", y[n * Nobs + i]);
    }
    fprintf(data_file, "\n");

  }
  fclose(data_file);
  printf("Signal and data of length %i generated and written in data.txt\n", data_length);

  /*
   * Calculate the inverse observation covariances
   */
  printf("Computing the observation covariance inverse...\n");
  fflush(stdout);
  clock_t start = clock();

  /* Observation noise covariance */
  gsl_vector * eval = gsl_vector_alloc(Nobs);
  gsl_matrix * evec = gsl_matrix_alloc(Nobs, Nobs);
  double * Rinv1 = (double*) malloc(Nobs * Nobs * sizeof(double));
  gsl_eigen_symmv_workspace * workspace = gsl_eigen_symmv_alloc(Nobs);

  /* Compute the eigen decomposition of R to get the inverse */
  gsl_eigen_symmv(R, eval, evec, workspace);
  FILE* rinv_out = fopen("R_inv.txt", "w");
  for (int i = 0; i < Nobs; i++) {
    for (int k = 0; k < Nobs; k++) {
      tmp = 0;
      for (int j = 0; j < Nobs; j++) {
        tmp += gsl_matrix_get(evec, i, j) * gsl_matrix_get(evec, k, j)
        / gsl_vector_get(eval, j);

      }
      Rinv1[i * Nobs + k] = tmp;
      fprintf(rinv_out, "%e, ", tmp);
    }
    fprintf(rinv_out, "\n");
  }
  fclose(rinv_out);
  printf("Observation covariance inverse computed and written in R_inv.txt\n");

  gsl_matrix_free(evec);
  gsl_vector_free(eval);
  gsl_eigen_symmv_free(workspace);

  /* Limit the observation correlation to a band and ... */
  int band = 0; // Note we only consider diagonal, i.e. band = 0.
  gsl_matrix * R_low = gsl_matrix_alloc(Nobs, Nobs);
  for (int i = 0; i < Nobs; i++) {
    for (int j = 0; j < Nobs; j++) {
      tmp = abs(i - j) <= band ? gsl_matrix_get(R, i, j) : 0;
      gsl_matrix_set(R_low, i, j, tmp);
    }
  }
  /* ...take inverse */
  gsl_vector * eval_low = gsl_vector_alloc(Nobs);
  gsl_matrix * evec_low = gsl_matrix_alloc(Nobs, Nobs);
  gsl_eigen_symmv_workspace * workspace_low = gsl_eigen_symmv_alloc(Nobs);
  double * Rinv0 = (double*) malloc(Nobs * Nobs * sizeof(double));

  gsl_eigen_symmv(R_low, eval_low, evec_low, workspace_low);

  FILE * band_R_out = fopen("R_band_inv.txt", "w");
  for (int i = 0; i < Nobs; i++) {
    for (int k = 0; k < Nobs; k++) {
      tmp = 0;
      for (int j = 0; j < Nobs; j++) {
        tmp += gsl_matrix_get(evec_low, i, j) * gsl_matrix_get(evec_low, k, j)
        / gsl_vector_get(eval_low, j);

      }
      Rinv0[i * Nobs + k] = tmp;
      fprintf(band_R_out, "%e ", tmp);
    }
    fprintf(band_R_out, "\n");
  }
  fclose(band_R_out);
  printf(
         "Inverse of the diagonal covariance approximation computed and written in R_band_inv.txt\n");

  gsl_matrix_free(evec_low);
  gsl_vector_free(eval_low);
  gsl_eigen_symmv_free(workspace_low);

  printf("elapsed %5.2f sec\n", (double) (clock() - start) / CLOCKS_PER_SEC);

  if (run_time_matching == 1) {
    /*
     * Particle allocation time matching
     *
     * This part of the code should be run when trying to find the correct N0 sample sizes to match
     * the BPF running time
     */

    particle_allocation_time_matching(N_bpf, sig_std, obs_std, Nobs, data_length, x0, y, Rinv1, Rinv0, R);

    printf("Particle allocation time matching done! Exiting...\n");
    return 0; // the code exits
  }

  /* ------------------------------------------------------------
   * Kalman filter
   * --------------------------------------------------------- */
  double *x_hat = (double*) malloc(data_length * sizeof(double));
  double *p_hat = (double*) malloc(data_length * sizeof(double));

  kalman_filter(x, y, Rinv1, sig_std, obs_std, x_hat, p_hat, data_length, x0, p0, Nobs);

  /* Write KF results in a file */
  FILE *kf_out = fopen("kf_out.txt", "w");
  for (int i = 0; i < data_length; i++)
    fprintf(kf_out, "%i %f %f %f\n", i, x[i], x_hat[i], p_hat[i]);
  fclose(kf_out);

  /*
   * Top level loop
   * */
  double scalers[N_SCALES] = {0.5,1,2,3,4,5,6,7,8,9,10};
  short N_levels = 2;
  long sample_sizes[2] = {0,0};
  long sample_size[1] = { N_bpf }; // This is for one level BPF

  double *x_hat_bpf = (double*) malloc(data_length * sizeof(double));
  double *errors = (double *) malloc((N_SETTINGS + 1) * sizeof(double));
  double *times = (double *) malloc((N_SETTINGS + 1) * sizeof(double));
  init_with_zeros(errors, N_SETTINGS + 1);
  init_with_zeros(times, N_SETTINGS + 1);

  double error;
  FILE * error_out = fopen("error.txt", "w");

  double filtering_time[4];
  double worst_case_sign_ratio[1] = { 100 }; // This is basically sustitute for "Inf"
  long unilevel_sample_size = N_bpf;
  long error_matched_unilevel_sample_size = 7 * N_bpf;
  int N1, N0;
  double step;

  /*
   * MAIN LOOP FOR TOP LEVEL MONTE CARLO
   */
  for (int i = 0; i < N_top; i++) {

    printf("Top level iteration %i\n", i);

    for (int k = 0; k < N_SCALES; k++) {

      scaler = scalers[k];

      printf("Scaler %f\n", scaler);

      /* Run the basic BPF (time matched) */
      sample_size[0] = (long) round(unilevel_sample_size * scaler);
      bootstrapfilter(y, sig_std, obs_std, data_length, x0, Nobs, sample_size[0],
                      x_hat_bpf, filtering_time, Rinv1, 1);

      /* Sum of squared errors over the whole signal */
      error = 0;
      for (int n = 0; n < data_length; n++) {
        error += (x_hat[n] - x_hat_bpf[n]) * (x_hat[n] - x_hat_bpf[n]);
      }

      fprintf(error_out, "%e %i %lu %f %i %e %e %f %e %f\n",
              error / (double) data_length, 0, sample_size[0], filtering_time[0], 0,
              filtering_time[1], filtering_time[2], worst_case_sign_ratio[0],
              filtering_time[3], scaler);
      fflush(error_out);

      if(scaler < 8) {
        /* Run the basic BPF (error matched) */
        sample_size[0] = (long) round(error_matched_unilevel_sample_size * scaler);
        bootstrapfilter(y, sig_std, obs_std, data_length, x0, Nobs, sample_size[0],
                        x_hat_bpf, filtering_time, Rinv1, 1);

        /* Sum of squared errors over the whole signal */
        error = 0;
        for (int n = 0; n < data_length; n++) {
          error += (x_hat[n] - x_hat_bpf[n]) * (x_hat[n] - x_hat_bpf[n]);
        }

        fprintf(error_out, "%e %i %lu %f %i %e %e %f %e %f\n",
                error / (double) data_length, 0, sample_size[0], filtering_time[0], -1,
                filtering_time[1], filtering_time[2], worst_case_sign_ratio[0],
                filtering_time[3], scaler);
        fflush(error_out);
      }

      /* Iterate over different configurations */
      N1 = 0;
      for (int j = 0; j < N_SETTINGS; j++) {

        step = (double)(N_bpf - 5) / (double)(N_SETTINGS-1);
        N1 = (int) j * step;
        N0 = (int) floor((N_bpf-N1)*sample_size_coeff);
        worst_case_sign_ratio[0] = 100;

        sample_sizes[0] = N0 * scaler;
        sample_sizes[1] = N1 * scaler;

        printf("N0 = %lu, N1 = %lu, N = %lu\n", sample_sizes[0], sample_sizes[1], sample_sizes[0]+sample_sizes[1]);

        MLbootstrapfilter(y, sig_std, obs_std, data_length, x0, Nobs, N_levels,
                          sample_sizes, x_hat_bpf, filtering_time, Rinv0, Rinv1,
                          worst_case_sign_ratio);

        /* Sum of squared errors over the whole signal */
        error = 0;
        for (int n = 0; n < data_length; n++) {
          error += (x_hat[n] - x_hat_bpf[n]) * (x_hat[n] - x_hat_bpf[n]);
        }

        fprintf(error_out, "%e %lu %lu %f %i %e %e %f %e %f\n",
                error / (double) data_length, sample_sizes[0],
                sample_sizes[1], filtering_time[0], j+1, filtering_time[1],
                filtering_time[2], worst_case_sign_ratio[0], filtering_time[3], scaler);

      }
      fflush(error_out);
    }

  }

  fclose(error_out);

  printf("Results of all BPF runs written in error.txt\n");
  printf("Error summary written in error_summary.txt\n");

  return 0;
}

void init_with_zeros(double *array, long length) {
  for (int i = 0; i < length; i++)
    array[i] = 0;
}

void scale_sample_sizes(double scaler, long sample_size_collection[N_SETTINGS][2], int n_settings) {

  for (int i = 0; i < n_settings; i++)
    for (int j = 0; j < 1; j++)
      sample_size_collection[i][j] = (long) round(scaler * sample_size_collection[i][j]);

}

void particle_allocation_time_matching(int bpf_sample_size, double std_sig, double obs_std, int Nobs, int data_length, double x0, double *y, double *Rinv1, double *Rinv0, gsl_matrix *R) {

  int N1_step = 20;
  int N0_step = 10000;
  int N0_init = 10;
  int N_top = 5;
  double time_per_iteration;
  double time_tolerance = 1.00; //
  double filtering_time[1];
  double worst_case_sign_ratio[1] = { 100 };

  double *x_hat_bpf = (double*) malloc(data_length * sizeof(double));

  /* Find the reference time per iteration for BPF */
  double bpf_time_per_iteration = 0;

  for(int nmc = 0; nmc < N_top; nmc++) {
    bootstrapfilter(y, std_sig, obs_std, data_length, x0, Nobs, bpf_sample_size,
                    x_hat_bpf, filtering_time, Rinv1, 1);
    bpf_time_per_iteration += filtering_time[0] / (double) N_top;
  }

  printf("Time per iteration (BPF): %0.5f\n",bpf_time_per_iteration);
  fflush(stdout);

  long sample_sizes[2] = {0,0};
  int N0;
  short direction;

  FILE* out = fopen("time_matched_sample_sizes.txt","w");

  /* For given L0 mesh, iterate over N1 */
  for(int N1 = 0; N1 <= bpf_sample_size ; N1 = N1 + N1_step) {

    N0 = N0_init;
    N0_step = 10000;

    /* For given L0 mesh and N1 iterate over increasing N0 until MLBPF time per iteration exceeds
     BPF time per iteration */
    time_per_iteration = 0;

    direction = 1; // up
    while(1) {

      sample_sizes[0] = N0;
      sample_sizes[1] = N1;

      /* Estimate the time per iteration by averaging N_top runs */
      time_per_iteration = 0;
      for(int nmc = 0; nmc < N_top; nmc++) {

        MLbootstrapfilter(y, std_sig, obs_std, data_length, x0, Nobs, 2,
                          sample_sizes, x_hat_bpf, filtering_time, Rinv0, Rinv1,
                          worst_case_sign_ratio);

        time_per_iteration += filtering_time[0] / (double) N_top;
      }

      if (direction == 1 && time_per_iteration > bpf_time_per_iteration * time_tolerance && N0_step > 10) {
        direction = -1; // start going down
        N0_step /= 2; // halve the step

      } else if (direction == -1 && time_per_iteration < bpf_time_per_iteration * 1.00 && N0_step > 10) {
        direction = 1;
        N0_step /= 2;
      } else if (N0_step < 10) {
        break;
      }
      N0 = N0 + direction * N0_step;

      if(N0 < 1) { // Make sure N0 remains positive and go up if this happens
        N0 = 1;
        direction = 1;
      }

    }
    fprintf(out,"%i %i %i %i\n",N0-N0_step,N1,0,0);
    fflush(out);
  }

  fclose(out);
  free(x_hat_bpf);
}

ext: c

max_stars_repo_name: heinekmp/MLBPF

max_issues_repo_name: heinekmp/MLBPF

max_forks_repo_name: heinekmp/MLBPF


text:
/* -*- mode: C; c-basic-offset: 4 -*- */
/* ex: set shiftwidth=4 tabstop=4 expandtab: */
/*
 * Copyright (c) 2011, Georgia Tech Research Corporation
 * All rights reserved.
 *
 * Author(s): Neil T. Dantam <ntd@gatech.edu>
 * Georgia Tech Humanoid Robotics Lab
 * Under Direction of Prof. Mike Stilman <mstilman@cc.gatech.edu>
 *
 *
 * This file is provided under the following "BSD-style" License:
 *
 *
 *   Redistribution and use in source and binary forms, with or
 *   without modification, are permitted provided that the following
 *   conditions are met:
 *
 *   * Redistributions of source code must retain the above copyright
 *     notice, this list of conditions and the following disclaimer.
 *
 *   * Redistributions in binary form must reproduce the above
 *     copyright notice, this list of conditions and the following
 *     disclaimer in the documentation and/or other materials provided
 *     with the distribution.
 *
 *   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND
 *   CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,
 *   INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
 *   MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 *   DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
 *   CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
 *   SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
 *   LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
 *   USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED
 *   AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
 *   LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
 *   ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 *   POSSIBILITY OF SUCH DAMAGE.
 *
 */

#include <amino.h>
#include <cblas.h>
#include "reflex.h"

#define SYMM_PART CblasUpper
#define SYMM_PARTC "U"

AA_API void rfx_lqg_kf_predict_cov( size_t n, const double *A, const double *V, double *P )
{
    // P = A * P * A**T + V
    int ni = (int)n;
    double *T = (double*)aa_mem_region_local_tmpalloc( 2 * n*n * sizeof(P[0]) );

    // T := A*P
    cblas_dsymm( CblasColMajor, CblasRight, SYMM_PART,
                 ni, ni,
                 1.0, P, ni,
                 A, ni,
                 0.0, T, ni );

    // P := (A*P) * A**T + V
    dlacpy_( SYMM_PARTC, &ni, &ni,
             V, &ni,
             P, &ni );
    cblas_dgemm( CblasColMajor, CblasNoTrans, CblasTrans,
                 ni, ni, ni,
                 1.0, T, ni,
                 A, ni,
                 1.0, P, ni );
}

int rfx_lqg_kf_correct_gain
( size_t n_x, size_t n_z, const double *C, const double *P, const double *W, double *K )
{
    int nxi = (int)n_x;
    int nzi = (int)n_z;

    double *ptr = (double*) aa_mem_region_local_alloc( sizeof(ptr[0]) *
                                                       (n_z*n_x + n_z*n_z) );

    double *CP = ptr;
    double *Kp = &ptr[n_z*n_x];

    // P is symmetric, so P*C^T == (C*P^T)^T == (C*P)^T

    // Kp := C * P * C**T + W = C*(C*P)^T
    // Kp^T := C*(C * P)^T
    cblas_dsymm( CblasColMajor, CblasRight, SYMM_PART,
                 nzi, nxi,
                 1.0, P, nxi,
                 C, nzi,
                 0.0, CP, nzi );

    dlacpy_( SYMM_PARTC, &nzi, &nzi,
             W, &nzi,
             Kp, &nzi );
    cblas_dgemm( CblasColMajor, CblasNoTrans, CblasTrans,
                 nzi, nzi, nxi,
                 1.0, C, nzi,
                 CP, nzi,
                 1.0, Kp, nzi );

    // K := P * C**T * Kp^-1
    // K Kp := P * C**T
    // CP = Kp^T K^T
    // Kp is symmetric
    // CP = Kp K^T
    int info;

    // Cholesky
    dpotrf_( SYMM_PARTC, &nzi,
             Kp, &nzi,
             &info );

    // Solve it!
    dpotrs_( SYMM_PARTC, &nzi, &nxi,
             Kp, &nzi,
             CP, &nzi,
             &info );

    // Transpose
    for( size_t j = 0; j < n_z; j ++ )
        for( size_t i = 0; i < n_x; i ++ )
            K[j*n_x+i] = CP[i*n_z+j];

    int i = 0;

    return i;
}

void rfx_lqg_kf_correct_cov
( size_t n_x, size_t n_z, const double *C, double *P, double *K )
{
    int nxi = (int)n_x;
    int nzi = (int)n_z;

    double *ptr = (double*)aa_mem_region_local_tmpalloc( 2 * n_x*n_x * sizeof(P[0]) );
    double *KC = ptr;
    double *P1 = &ptr[n_x*n_x];

    // P := (I - K*C) * P
    cblas_dgemm( CblasColMajor, CblasNoTrans, CblasNoTrans,
                 nxi, nxi, nzi,
                 -1.0, K, nxi,
                 C, nzi,
                 0.0, KC, nxi );
    // KC += I
    for( size_t i = 0; i < n_x*n_x; i = i + n_x + 1 )
        KC[i] += 1;

    cblas_dsymm( CblasColMajor, CblasRight, SYMM_PART,
                 nxi, nxi,
                 1.0, P, nxi,
                 KC, nxi,
                 0.0, P1, nxi );
    dlacpy_( SYMM_PARTC, &nxi, &nxi, P1, &nxi, P, &nxi );
}

int rfx_lqg_ekf_predict
( void *cx, size_t n_x, double *x, const double *u, double *P, const double *V,
  rfx_lqg_ekf_process_fun process )
{
    double F[n_x*n_x];
    int i = process( cx, x, u, F );
    rfx_lqg_kf_predict_cov( n_x, F, V, P );
    return i;
}

int rfx_lqg_ekf_correct
( void *cx, size_t n_x, size_t n_z, double *x, const double *z, double *P, const double *W,
  rfx_lqg_ekf_measure_fun measure, rfx_lqg_ekf_innovate_fun innovate, rfx_lqg_ekf_update_fun update )
{
    int i;

    double H[n_z*n_x];
    double K[n_x*n_z];
    {
        double y[n_z];
        i = measure(cx, x, y, H);
        if(i) return i;

        i = rfx_lqg_kf_correct_gain( n_x, n_z, H, P, W, K );
        if(i) return i;

        i = innovate(cx, x, z, y);
        if(i) return i;

        {
            double Ky[n_x];
            cblas_dgemv( CblasColMajor, CblasNoTrans,
                         (int)n_x, (int)n_z,
                         1.0, K, (int)n_x,
                         y, 1,
                         0.0, Ky, 1 );

            i = update(cx, x, Ky);
        }
    }

    rfx_lqg_kf_correct_cov( n_x, n_z, H, P, K );

    return i;
}

AA_API void rfx_lqg_init( rfx_lqg_t *lqg, size_t n_x, size_t n_u, size_t n_z ) {
    memset( lqg, 0, sizeof(*lqg) );
    lqg->n_x = n_x;
    lqg->n_u = n_u;
    lqg->n_z = n_z;

    lqg->x = AA_NEW0_AR( double, lqg->n_x );
    lqg->u = AA_NEW0_AR( double, lqg->n_u );
    lqg->z = AA_NEW0_AR( double, lqg->n_z );

    lqg->A = AA_NEW0_AR( double, lqg->n_x * lqg->n_x );
    lqg->B = AA_NEW0_AR( double, lqg->n_x * lqg->n_u );
    lqg->C = AA_NEW0_AR( double, lqg->n_z * lqg->n_x );

    lqg->P = AA_NEW0_AR( double, lqg->n_x * lqg->n_x );
    lqg->V = AA_NEW0_AR( double, lqg->n_x * lqg->n_x );
    lqg->W = AA_NEW0_AR( double, lqg->n_z * lqg->n_z );

    lqg->Q = AA_NEW0_AR( double, lqg->n_x * lqg->n_x );
    lqg->R = AA_NEW0_AR( double, lqg->n_u * lqg->n_u );

    lqg->K = AA_NEW0_AR( double, lqg->n_z * lqg->n_x );
    lqg->L = AA_NEW0_AR( double, lqg->n_u * lqg->n_x );

}
AA_API void rfx_lqg_destroy( rfx_lqg_t *lqg ) {

    free(lqg->x);
    free(lqg->u);
    free(lqg->z);

    free(lqg->A);
    free(lqg->B);
    free(lqg->C);

    free(lqg->Q);
    free(lqg->R);

    free(lqg->L);
    free(lqg->K);

}

/*
 * X = op(A)*op(B)*op(C) + beta*X
 *
 * m*n = (m*k) * (k*n)
 *
 * X:     m*n
 * op(A): m*k
 * op(B): k*p
 * op(C): p*n
 */

static inline void matmul3( size_t m, size_t n, size_t k, size_t p,
                            enum CBLAS_TRANSPOSE transA,
                            enum CBLAS_TRANSPOSE transB,
                            enum CBLAS_TRANSPOSE transC,
                            const double *A, size_t lda,
                            const double *B, size_t ldb,
                            const double *C, size_t ldc,
                            double beta, double *X, size_t ldx) {
    double T[m*p];
    // T := A*B
    cblas_dgemm( CblasColMajor, transA, transB,
                 (int)m, (int)p, (int)k,
                 1.0, A, (int)lda,
                 B, (int)ldb,
                 0.0, T, (int)m );

    // X := (A*B) * C + beta*X
    cblas_dgemm( CblasColMajor, CblasNoTrans, transC,
                 (int)m, (int)n, (int)p,
                 1.0, T, (int)m,
                 C, (int)ldc,
                 beta, X, (int)ldx );
}

static void kf_innovate( rfx_lqg_t *lqg, double *xh ) {

    // xh = x + K * (z - C*x)

    // T := z - C*x
    double T[lqg->n_z];
    memcpy(T, lqg->z, sizeof(T));
    cblas_dgemv( CblasColMajor, CblasNoTrans,
                 (int)lqg->n_z, (int)lqg->n_x,
                 -1.0, lqg->C, (int)lqg->n_z,
                 lqg->x, 1,
                 1.0, T, 1 );
    // xh = K * (z - C*x) + xh
    cblas_dgemv( CblasColMajor, CblasNoTrans,
                 (int)lqg->n_x, (int)lqg->n_z,
                 1.0, lqg->K, (int)lqg->n_x,
                 T, 1,
                 1.0, xh, 1 );
}

AA_API void rfx_lqg_kf_predict( rfx_lqg_t *lqg ) {
    // x = A*x + B*u
    {
        double x[lqg->n_x];
        aa_lsim_dstep( lqg->n_x, lqg->n_u,
                       lqg->A, lqg->B,
                       lqg->x, lqg->u,
                       x );
        memcpy(lqg->x, x, sizeof(x));
    }
    // P = A * P * A**T + V
    rfx_lqg_kf_predict_cov( lqg->n_x, lqg->A, lqg->V, lqg->P );
}

AA_API void rfx_lqg_kf_correct( rfx_lqg_t *lqg ) {

    rfx_lqg_kf_correct_gain( lqg->n_x, lqg->n_z,
                             lqg->C, lqg->P, lqg->W, lqg->K );

    // x = x + K * (z - C*x)
    kf_innovate( lqg, lqg->x );

    // P = (I - K*C) * P
    rfx_lqg_kf_correct_cov( lqg->n_x, lqg->n_z,
                            lqg->C, lqg->P, lqg->K );
}

// kalman-bucy gain
void rfx_lqg_kbf_gain( rfx_lqg_t *lqg )
{
    rfx_lqg_kbf_gain_work( lqg->n_x, lqg->n_z,
                           lqg->A, lqg->C, lqg->V, lqg->W, lqg->P, lqg->K );
}

// kalman-bucy gain
void rfx_lqg_kbf_gain_work
( size_t n_x, size_t n_z,
  const double *A, const double *C, const double *V, const double *W, double *P, double *K )
{
    // dP = A*P + P*A' - P*C'*W^{-1}*C*P + V
    // K = P * C' * W^{-1}

    double Winv[n_z*n_z];
    memcpy(Winv, W, sizeof(Winv));
    aa_la_inv( n_z, Winv);

    {
        double ric_B[n_x*n_x];

        // ric_B := C' * W^{-1} * C
        matmul3( n_x, n_x, n_z, n_z,
                 CblasTrans, CblasNoTrans, CblasNoTrans,
                 C, n_z,
                 Winv, n_z,
                 C, n_z,
                 0.0, ric_B, n_x );

        double ric_A[n_x*n_x];
        aa_la_transpose2(n_x, n_x, A, ric_A );

        // solve ARE with dP = 0, result is P
        //aa_tick("laub: ");
        aa_la_care_laub( n_x, n_x, n_x,
                         ric_A, ric_B, V, P );
        //aa_tock();
    }

    // K = P * C' * W^{-1}
    matmul3( n_x, n_z, n_x, n_z,
             CblasNoTrans, CblasTrans, CblasNoTrans,
             P, n_x,
             C, n_z,
             Winv, n_z,
             0.0, K, n_x );

}

void rfx_lqg_lqr_gain( rfx_lqg_t *lqg ) {
    // A'*S + S*A - S*B*R^{-1}*B'*S + Q = 0

    double Rinv[lqg->n_u*lqg->n_u];
    memcpy(Rinv, lqg->R, sizeof(Rinv));
    aa_la_inv( lqg->n_u, Rinv);

    double S[lqg->n_x*lqg->n_x];
    {
        double ric_B[lqg->n_x*lqg->n_x];

        // ric_B := B * R^{-1} * B'
        matmul3( lqg->n_x, lqg->n_x, lqg->n_u, lqg->n_u,
                 CblasNoTrans, CblasNoTrans, CblasTrans,
                 lqg->B, lqg->n_x,
                 Rinv, lqg->n_u,
                 lqg->B, lqg->n_x,
                 0.0, ric_B, lqg->n_x );

        // solve ARE with dS = 0, result is S
        aa_la_care_laub( lqg->n_x, lqg->n_x, lqg->n_x,
                         lqg->A, ric_B, lqg->Q, S );
    }

    // L = R^{-1}*B'*S
    matmul3( lqg->n_u, lqg->n_x, lqg->n_u, lqg->n_x,
             CblasNoTrans, CblasTrans, CblasNoTrans,
             Rinv, lqg->n_u,
             lqg->B, lqg->n_x,
             S, lqg->n_x,
             0.0, lqg->L, lqg->n_u );
}

static void kbf_step_fun( const void *cx,
                          double t, const double *restrict x,
                          double *restrict dx ) {
    // dx = A*x + B*u + K(z-Cx)
    (void)t;
    rfx_lqg_t *lqg = (rfx_lqg_t*)cx;

    // dx := A*x + B*u
    aa_lsim_dstep( lqg->n_x, lqg->n_u,
                   lqg->A, lqg->B,
                   x, lqg->u,
                   dx );
    // dx := K * (z - C*x) + dx
    kf_innovate( lqg, dx );
}

void rfx_lqg_kbf_step1( rfx_lqg_t *lqg, double dt ) {
    // dx = A*x + B*u + K(z-Cx)
    double dx[lqg->n_x];
    kbf_step_fun( lqg, 0, lqg->x, dx );

    // x := dt*dx + x
    cblas_daxpy( (int)lqg->n_x, dt, dx, 1, lqg->x, 1 );
}

void rfx_lqg_kbf_step4( rfx_lqg_t *lqg, double dt ) {
    double x1[lqg->n_x];
    // runge-kutta 4
    aa_odestep_rk4( lqg->n_x, kbf_step_fun, lqg,
                    0, dt, lqg->x, x1 );
    memcpy( lqg->x, x1, sizeof(x1) );
}

AA_API void rfx_lqg_lqr_ctrl( rfx_lqg_t *lqg ) {
    cblas_dgemv( CblasColMajor, CblasNoTrans,
                 (int)lqg->n_u, (int)lqg->n_x,
                 -1.0, lqg->L, (int)lqg->n_u,
                 lqg->x, 1,
                 0.0, lqg->u, 1 );
}

void rfx_lqg_sys( const void *cx,
                  double t, const double *restrict x,
                  double *restrict dx ) {
    rfx_lqg_t * lqg = (rfx_lqg_t*)cx;
    (void)t;
    aa_lsim_dstep( lqg->n_x, lqg->n_u,
                   lqg->A, lqg->B,
                   x, lqg->u,
                   dx );
}

//AA_API void rfx_lqg_observe_euler( rfx_lqg_t *lqg, double dt, aa_mem_region_t *reg ) {
    // --- calculate kalman gain ---

    // --- predict from previous input ---
    // dx = A*x + B*u + K * ( z - C*xh )
    //double *dx = (double*)aa_mem_region_alloc(reg, sizeof(double)*lqg->n_x);
    // zz := z
    //double *zz = (double*)aa_mem_region_alloc(reg, sizeof(double)*lqg->n_z);

    // x = x + dx * dt
    //aa_la_axpy( lqg->n_x, dt, dx, lqg->x );
//}

//AA_API void rfx_lqg_ctrl( rfx_lqg_t *lqg, double dt, aa_mem_region_t *reg ) {
    // compute optimal gain
    //  -dS = A'*S + S*A - S*B*R^{-1}*B' + Q
    // solve ARE, result is S
    // L = R^{-1} * B' * S  : (nu*nu) * (nu*nx) * (nx*nx)

    // compute current input
    // u = -Lx
//}

/* int rfx_lqg_duqu_predict */
/* ( double dt, double *S, double *dx, double *P, const double *V ) */
/* { */
/*     return 0; */
/* } */

max_stars_repo_name: golems/reflex

max_issues_repo_name: golems/reflex

max_forks_repo_name: golems/reflex


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

#ifndef SQ_INCLUDE_GUARD_parser_FilterSpec_h_
#define SQ_INCLUDE_GUARD_parser_FilterSpec_h_

#include "core/Primitive.h"
#include "core/typeutil.h"

#include <compare>
#include <gsl/gsl>
#include <iosfwd>
#include <optional>
#include <variant>

namespace sq::parser {

/**
 * Represents the lack of a filter specification for a field access.
 */
struct NoFilterSpec {
  SQ_ND auto operator<=>(const NoFilterSpec &) const = default;
};
std::ostream &operator<<(std::ostream &os, NoFilterSpec nlfs);

/**
 * Represents access of an indexed element in a list of results for a field
 * access.
 */
struct ElementAccessSpec {
  gsl::index index_;
  SQ_ND auto operator<=>(const ElementAccessSpec &) const = default;
};
std::ostream &operator<<(std::ostream &os, ElementAccessSpec leas);

/**
 * Represents a Python-style slice of a list of results for a field access.
 */
struct SliceSpec {
  std::optional<gsl::index> start_;
  std::optional<gsl::index> stop_;
  std::optional<gsl::index> step_;
  SQ_ND auto operator<=>(const SliceSpec &) const = default;
};
std::ostream &operator<<(std::ostream &os, SliceSpec lss);

enum class ComparisonOperator {
  GreaterThanOrEqualTo,
  GreaterThan,
  LessThanOrEqualTo,
  LessThan,
  Equals
};

/**
 * Represents a comparison to determine whether to keep a field or not.
 */
struct ComparisonSpec {
  std::string member_;
  ComparisonOperator op_;
  Primitive value_;
  SQ_ND auto operator<=>(const ComparisonSpec &) const = default;
};

std::ostream &operator<<(std::ostream &os, const ComparisonOperator &op);
std::ostream &operator<<(std::ostream &os, const ComparisonSpec &cs);

using FilterSpec =
    std::variant<NoFilterSpec, ElementAccessSpec, SliceSpec, ComparisonSpec>;

} // namespace sq::parser

#endif // SQ_INCLUDE_GUARD_parser_FilterSpec_h_

max_stars_repo_name: jonathanhaigh/sq

max_issues_repo_name: jonathanhaigh/sq

max_forks_repo_name: jonathanhaigh/sq


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

#include "SM.h"
#include "THDM.h"
#include <complex>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_integration.h>

using namespace std;

/**
* @brief Calculates the decay modes of 2HDM Higgs bosons
*
* Given a THDM object, a DecayTable can be generated. From this table, the
* Higgs boson decay widths and branching ratios are obtained. For the complete
* list of available decay modes, we refer to the complete documentation, or
* the list of member methods below.
*/
class DecayTable {

 public:

  /**
  * @brief Default constructor
  *
  * This default constructor takes a THDM object as argument for which
  * the decays are to be calculated. %SM properties are taken from the SM
  * object in the THDM.
  *
  * @param mod Two-Higgs doublet model for which to calculate decay modes
  */
  DecayTable(THDM mod);

  /**
  * @brief Sets underlying 2HDM
  *
  * This method sets the THDM underlying the DecayTable.
  *
  * @param model Two-Higgs doublet model for which to calculate decay modes
  */
	void set_model(THDM model);

  /**
  * @brief Underlying 2HDM
  *
  * Use to obtain the underlying THDM object
  *
  * @returns The THDM object on which this DecayTable operates
  */
	THDM get_model();

  /**
  * @brief Decay width \f$\Gamma(h\to u_1\overline{u}_2) \f$
  *
  * This method calculates the on-shell decay width for the decay of Higgs
  * boson \a h to a pair of up-type quarks. QCD corrections have been included.
  *
  * @param h  Index of Higgs boson (1,2,3 = h,H,A)
  * @param u1 Index of up-type quark (1,2,3 = \f$ u,c,t \f$)
  * @param u2 Index of up-type antiquark (1,2,3 = \f$ \bar{u},\bar{c},\bar{t} \f$)
  *
  * @returns The decay width in GeV
  */
  double  get_gamma_huu(int h, int u1, int u2);

  /**
  * @brief Decay width \f$ \Gamma(h\to d_1\overline{d}_2) \f$
  *
  * This method calculates the on-shell decay width for the decay of Higgs
  * boson \a h to a pair of down-type quarks. QCD corrections have been included.
  *
  * @param h  Index of Higgs boson (1,2,3 = h,H,A)
  * @param d1 Index of down-type quark (1,2,3 = \f$ d,s,b \f$)
  * @param d2 Index of down-type antiquark (1,2,3 = \f$ \bar{d},\bar{s},\bar{b} \f$)
  *
  * @returns The decay width in GeV
  */
  double  get_gamma_hdd(int h, int d1, int d2);

  /**
  * @brief Decay width \f$ \Gamma(h\to l_1\overline{l}_2) \f$
  *
  * This method calculates the on-shell decay width for the decay of Higgs
  * boson \a h to a pair of charged leptons.
  *
  * @param h  Index of Higgs boson (1,2,3 = h,H,A)
  * @param l1 Index of lepton (1,2,3 = \f$ e,\mu,\tau \f$)
  * @param l2 Index of antilepton (1,2,3 = \f$ \bar{e},\bar{\mu},\bar{\tau} \f$)
  *
  * @returns The decay width in GeV
  */
  double  get_gamma_hll(int h, int l1, int l2);

  /**
  * @brief Decay width \f$ \Gamma(h\to u\overline{d}) \f$
  *
  * This method calculates the on-shell decay width for the decay of Higgs
  * boson \a h to a pair of quarks. QCD corrections have been included.
  *
  * @param h  Index of Higgs boson (4 = H+)
  * @param d  Index of down-type antiquark (1,2,3 = \f$ \bar{d},\bar{s},\bar{b} \f$)
  * @param u  Index of up-type quark (1,2,3 = \f$ u,c,t \f$)
  *
  * @returns The decay width in GeV
  */
  double  get_gamma_hdu(int h, int d, int u);

  /**
  * @brief Decay width \f$ \Gamma(h\to l\overline{\nu}_{l'}) \f$
  *
  * This method calculates the on-shell decay width for the decay of Higgs
  * boson \a h to lepton-neutrino.
  *
  * @param h  Index of Higgs boson (4 = H+)
  * @param l  Index of charged lepton (1,2,3 = \f$ e^+,\mu^+,\tau^+ \f$)
  * @param n  Index of neutrino (1,2,3 = \f$ \nu_e,\nu_\mu,\nu_\tau \f$)
  *
  * @returns The decay width in GeV
  */
  double  get_gamma_hln(int h, int l, int n);

  /**
  * @brief Decay width \f$ \Gamma(h\to H_1 H_2) \f$
  *
  * This method calculates the on-shell decay width for the decay of Higgs
  * boson \a h to a pair of Higgs bosons
  *
  * @param h  Index of decaying Higgs boson (1,2,3,4 = h,H,A,H+)
  * @param h1 Index of first Higgs boson (1,2,3,4 = h,H,A,H+)
  * @param h2 Index of second Higgs boson (1,2,3,4 = h,H,A,H+)
  *
  * @returns The decay width in GeV
  */
  double  get_gamma_hhh(int h, int h1, int h2);

  /**
  * @brief Decay width \f$ \Gamma(h\to VV) \f$
  *
  * This method calculates the decay width for the Higgs boson \a h
  * to a pair of vector bosons. The decay mode with one vector
  * boson off-shell is included.
  *
  * @param h  Index of decaying Higgs boson (1,2,3,4 = h,H,A,H+)
  * @param v  Index of vector bosons (1,2,3 = \f$\gamma \f$,Z,W)
  *
  * @returns The decay width in GeV
  */
  double  get_gamma_hvv(int h, int v);

  /**
  * @brief Decay width \f$ \Gamma(h\to VH) \f$
  *
  * This method calculates the decay width for the Higgs boson \a h
  * to one massive vector and one Higgs boson. Decay with the vector
  * boson off-shell is included.
  *
  * @param h  Index of decaying Higgs boson (1,2,3,4 = h,H,A,H+)
  * @param V  Index of vector boson (2,3 = \f$\gamma \f$,Z,W)
  * @param H  Index of final-state Higgs boson (1,2,3,4 = h,H,A,H+)
	*
  * @returns The decay width in GeV
  */
  double  get_gamma_hvh(int H, int V, int h);

  /**
  * @brief Decay width \f$ \Gamma(h\to gg) \f$
  *
  * This method calculates the decay width for the Higgs boson \a h
  * to a pair of gluons. LO QCD corrections are included.
  *
  * @param h  Index of decaying Higgs boson (1,2,3,4 = h,H,A,H+)
	*
  * @returns The decay width in GeV
  */
  double  get_gamma_hgg(int h);

  /**
  * @brief Decay width \f$ \Gamma(h\to \gamma\gamma) \f$
  *
  * This method calculates the decay width for the neutral
  * Higgs boson \a h to a pair of photons.
  *
  * @param h  Index of decaying Higgs boson (1,2,3 = h,H,A)
	*
  * @returns The decay width in GeV
  */
	double  get_gamma_hgaga(int h);

  /**
  * @brief Decay width \f$ \Gamma(h\to Z\gamma) \f$
  *
  * This method calculates the decay width for the neutral
  * Higgs boson \a h to a Z and a photon.
  *
  * @param h  Index of decaying Higgs boson (1,2,3 = h,H,A)
	*
  * @returns The decay width in GeV
  */
	double  get_gamma_hZga(int h);

  /**
  * @brief Decay width \f$ \Gamma(h\to WZ) \f$
  *
  * This method calculates the decay width for the neutral
  * Higgs boson \a h to a Z and a photon.
  *
  * @param h  Index of decaying Higgs boson (4 = H+-)
    *
  * @returns The decay width in GeV
  */
    double  get_gamma_hWZ(int h);

  /**
  * @brief Decay width \f$ \Gamma(h\to W\gamma) \f$
  *
  * This method calculates the decay width for the neutral
  * Higgs boson \a h to a Z and a photon.
  *
  * @param h  Index of decaying Higgs boson (4 = H+-)
    *
  * @returns The decay width in GeV
  */
    double  get_gamma_hWga(int h);

  /**
  * @brief Total width \f$ \Gamma_h \f$
  *
  * Calculates the total decay width of Higgs boson \a h
  *
  * @param h  Index of Higgs boson (1,2,3,4 = h,H,A,H+)
	*
  * @returns Total width in GeV
  */
  double  get_gammatot_h(int h);

  /**
  * @brief Total width \f$ \Gamma_V \f$
  *
  * Returns the total decay width of vector boson \a v
  *
  * @param v  Index of vector boson (2,3 = Z,W)
	*
  * @returns Total width in GeV
  */
  double  get_gammatot_v(int v);

  /**
  * @brief Total width \f$ \Gamma_t \f$
  *
  * Returns the total decay width of the top quark
  *
  * @returns Total width in GeV
  */
  double  get_gammatot_top();

  /**
  * @brief Decay width for \f$ t \to H^+X \f$
  *
  * Returns the decay width of the top quark in the charged Higgs mode
  *
  * @returns Decay width in GeV
  *
  * @param u Index of decaying quark (1,2,3 = \f$ u,c,t \f$)
  * @param h Index of Higgs boson (4 = H+)
	* @param d Index of down-type quark (1,2,3 = \f$ d,s,b \f$)
  */
  double  get_gamma_uhd(int u, int h, int d);

  double  get_gamma_uhu(int u1, int h, int u2);

  /**
  *	@brief Prints the decay modes of a Higgs boson
  *
  * The decay modes of Higgs boson \a h are printed to stdout
  *
  * @param h Index of Higgs boson (1,2,3,4 = h,H,A,H+)
  */
  void    print_decays(int h);

  /**
  *	@brief Prints the total width of a Higgs boson
  *
  * The total decay width of Higgs boson \a h are printed to stdout
  *
  * @param h Index of Higgs boson (1,2,3,4 = h,H,A,H+)
  */
  void    print_width(int h);

  /**
  *	@brief Prints the decay modes of the top quark
  *
  * The decay modes of the top quark are printed to stdout
  */
  void    print_top_decays();

  /**
  *	@brief Prints decay information for a Higgs boson in LesHouches format
  *
  * The decay information for Higgs boson \a h are printed to a file in
  * LesHouches format
  *
  * @param output The name of the output file to write
  * @param h 			Index of Higgs boson (1,2,3,4 = h,H,A,H+)
  * @param full 	If \a true, all decay modes are printed, if \a false only the total width
  */
  void    print_decays_LesHouches(FILE* output, int h, bool full);

  /**
  *	@brief Prints decay information for the top quark in LesHouches format
  *
  * The decay information for the top quark are printed to a file in
  * LesHouches format
  *
  * @param output The name of the output file to write
  * @param full 	If \a true, all decay modes are printed, if \a false only the total width
  */
  void    print_top_decays_LesHouches(FILE* output, bool full);

  /**
  *     @brief Turns QCD corrections on or off
  *
  * This method is used to turn QCD corrections on or off. If the output is
  * meant to be used with the MadGraph/MadEvent 2HDMC model
  * QCD corrections should be turned off to get a consistent result.
  *
  * @param set  If \a true QCD corrections are turned on,
  *             if \a false QCD corrections are turned off
  */
  void    set_qcd(bool set);

  static double  DHp(double ui, double uj, double xi, double xj, double sqL);
  static double  DHm(double ui, double uj, double xi, double xj, double sqL);
  static double  BHp(double ui, double uj, double xi, double xj, double sqL);

 private:
  THDM model;
  SM sm;

  complex <double> F_sf(double t);
  complex <double> F_pf(double t);
  complex <double> F_0(double t);
  complex <double> F_1(double t);
  complex <double> ftau(double t);
  complex <double> gtau(double t);
  complex <double> I_1(double tau, double lambda);
  complex <double> I_2(double tau, double lambda);
  complex <double> FF_s(double tau, double lambda);
  complex <double> FF_p(double tau, double lambda);
  complex <double> FW(double tau, double lambda);
  complex <double> FHp(double tau, double lambda);

  double gammatot_h[5];
  double gamma_uhd[5][5][5];
  double gamma_uhu[5][5][5];
  double gamma_hdd[5][5][5];
  double gamma_huu[5][5][5];
  double gamma_hdu[5][5][5];
  double gamma_hll[5][5][5];
  double gamma_hln[5][5][5];
  double gamma_hgg[5];
  double gamma_hgaga[5];
  double gamma_hZga[5];
  double gamma_hWga[5]; // For H+- only
  double gamma_hWZ[5]; // For H+- only
  double gamma_hvv[5][5];
  double gamma_hvh[5][5][5];
  double gamma_hhh[5][5][5];

  double  br(double dG, double G);

  double  hvv_onshell(int h, int V, double M);
  double  hvv_offshell(int h, int V, double M);
  double  hvv_all(int h, int V, double M);
  double  hff_onshell(double M, double m1, double m1run, double m2, double m2run, complex<double> cs, complex<double> cp, int Nc, int h, bool tt);
  double  hpff_onshell(double M, double m1, double m1run, double m2, double m2run, complex<double> cs, complex<double> cp, int Nc, int h);
  double  htt_onshell(double M, double m1, double m2, int Nc, int h);
  double  htt_offshell(double M, double m1, double m2, int Nc, int h);
  double  htb_offshell(double M, double m1, double m1run, double m2, double m2run, complex<double> cs, complex<double> cp, int Nc);

  double  hvh_onshell(int H, int V, int h, double M);
  double  hvh_offshell(int H, int V, int h, double M);
  double  hdu_offshell(int h, int d, int u, double M);
  double  hgaga(int h);
  double  hZga(int h);
  double  hWga(int h);
  double  hWZ(int h);
  double  hgg(int h);
  double  PS2(double M, double m1, double m2);
  double  LAM(double m12, double m22, double m32);

  double  interp(double R, double x, double y, double c);

  bool    qcd_on;
  int     err_code;

  void    print_decay_LesHouches(FILE* output, double br, int id1, int id2);
  void    print_decay(const char *h, const char *id1, const char *id2, double g, double br);
  void    print_decays(FILE* output, int h, bool full, bool lh);

};

#endif


text:
/* Easel's interfaces to the GNU Scientific Library
 */
#ifdef HAVE_LIBGSL
#include "esl_config.h"

#include <stdlib.h>
#include "easel/easel.h"
#include "easel/dmatrix.h"
#include <gsl/gsl_math.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_eigen.h>

int
esl_GSL_MatrixInversion(ESL_DMATRIX *A, ESL_DMATRIX **ret_Ai)
{
  ESL_DMATRIX     *Ai;          /* RETURN: A^{-1}             */
  gsl_matrix_view  Av;          /* gsl view of matrix A       */
  gsl_matrix      *LU;          /* LU decomposition of A      */
  gsl_matrix      *Aiv;         /* gsl version of A^{-1}      */
  gsl_permutation *permute;
  int              signum;
  int              i,j;

  Ai = esl_dmx_Alloc(A->n, A->m);

  /* Invert U to get Ui, using LU decomposition.
   */
  Av      = gsl_matrix_view_array(A->mx[0], A->n, A->n);
  LU      = gsl_matrix_alloc(A->n, A->n);
  Aiv     = gsl_matrix_alloc(A->n, A->n); /* U^{-1}: inverse of U    */
  permute = gsl_permutation_alloc(A->n);
  gsl_matrix_memcpy(LU, &Av.matrix);

  if (gsl_linalg_LU_decomp(LU, permute, &signum) != 0) ESL_EXCEPTION(eslEUNKNOWN, "gsl failed");
  if (gsl_linalg_LU_invert(LU, permute, Aiv) != 0)     ESL_EXCEPTION(eslEUNKNOWN, "gsl failed");

  gsl_matrix_free(LU);
  gsl_permutation_free(permute);

  /* recover the matrix from gsl.
   */
  for (i = 0; i < A->n; i++)
    for (j = 0; j < A->n; j++)
      Ai->mx[i][j] = gsl_matrix_get(Aiv, i, j);
  gsl_matrix_free(Aiv);

  ret->Ai = Ai;
  return eslOK;
}

#endif /*HAVE_LIBGSL*/

/*****************************************************************
 * @LICENSE@
 *
 * SVN $Id$
 * SVN $URL$
 *****************************************************************/

max_stars_repo_name: STRIDES-Codes/Exploring-the-Microbiome-

max_issues_repo_name: STRIDES-Codes/Exploring-the-Microbiome-

max_forks_repo_name: STRIDES-Codes/Exploring-the-Microbiome-


text:
/*
 * Copyright (c) 2016-2021 lymastee, All rights reserved.
 * Contact: lymastee@hotmail.com
 *
 * This file is part of the gslib project.
 *
 * 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.
 */

#pragma once

#ifndef scene_048bb735_09b4_43d2_9ff5_5e40441c230b_h
#define scene_048bb735_09b4_43d2_9ff5_5e40441c230b_h

#include <gslib/std.h>
#include <ariel/config.h>
#include <ariel/rendersys.h>
#include <ariel/framesys.h>
#include <ariel/widget.h>

__ariel_begin__

class rose;
class stage;
typedef list<stage*> stages;

class __gs_novtable stage abstract:
    public frame_listener
{
public:
    stage();
    virtual ~stage() {}
    virtual const gchar* get_name() const = 0;
    virtual bool setup() = 0;
    virtual void draw() = 0;

public:
    template<class _sty> inline
    _sty* as() { return static_cast<_sty*>(this); }
    template<class _sty> inline
    const _sty* as_const() const { return static_cast<const _sty*>(this); }

protected:
    stage*              _prev_presentation_stage;
    stage*              _next_presentation_stage;
    stage*              _prev_notification_stage;
    stage*              _next_notification_stage;

protected:
    void on_next_draw();
    void on_next_frame_start();
    void on_next_frame_end();
    bool on_next_event(const frame_event& event);

public:
    stage* prev_presentation_stage() const { return _prev_presentation_stage; }
    stage* next_presentation_stage() const { return _next_presentation_stage; }
    stage* prev_notification_stage() const { return _prev_notification_stage; }
    stage* next_notification_stage() const { return _next_notification_stage; }

public:
    static void connect_presentation_order(stage* stg1, stage* stg2);
    static void connect_notification_order(stage* stg1, stage* stg2);
    static void connect_presentation_order(stage* stg1, stage* stg2, stage* stg3);
    static void connect_notification_order(stage* stg1, stage* stg2, stage* stg3);
};

class ui_stage:
    public stage
{
public:
    ui_stage(rose* r) { _rose = r; }
    virtual const gchar* get_name() const override { return _t("ui"); }
    virtual bool setup() override;
    virtual void draw() override;
    virtual void on_frame_start() override { on_next_frame_start(); }
    virtual void on_frame_end() override { on_next_frame_end(); }
    virtual bool on_frame_event(const frame_event& event) override;

public:
    wsys_manager* get_wsys_manager() { return &_wsys_manager; }

protected:
    rose*               _rose;
    wsys_manager        _wsys_manager;

private:
    void setup_draw();
};

class scene:
    public frame_listener
{
public:
    static scene* get_singleton_ptr()
    {
        static scene inst;
        return &inst;
    }

private:
    scene();

public:
    ~scene();
    stage* get_stage(const gchar* name);
    wsys_manager* get_ui_system() const { return _uisys; }
    fontsys* get_fontsys() const { return _fontsys; }
    rendersys* get_rendersys() const { return _rendersys; }
    void set_rendersys(rendersys* rsys) { _rendersys = rsys; }
    void set_rose(rose* ptr) { _rose = ptr; }
    void set_fontsys(fontsys* fsys);
    void setup();
    void destroy();
    void destroy_all_stages();
    void add_stage(stage* stg) { _stages.push_back(stg); }
    void set_presentation_before(stage* pos, stage* tar);
    void set_presentation_after(stage* pos, stage* tar);
    void set_notification_before(stage* pos, stage* tar);
    void set_notification_after(stage* pos, stage* tar);

protected:
    rendersys*          _rendersys;
    rose*               _rose;
    stages              _stages;
    wsys_manager*       _uisys;
    fontsys*            _fontsys;
    stage*              _notify;
    stage*              _present;

protected:
    void draw();

public:
    virtual void on_frame_start() override;
    virtual void on_frame_end() override;
    virtual bool on_frame_event(const frame_event& event) override;
};

__ariel_end__

#endif

max_stars_repo_name: lymastee/gslib

max_issues_repo_name: lymastee/gslib

max_forks_repo_name: lymastee/gslib


text:
/* histogram/test_trap.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 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 2 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

#include <config.h>
#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_histogram.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_ieee_utils.h>

#define N 397

void my_error_handler (const char *reason, const char *file,
		       int line, int err);
int status = 0;

int
main (void)
{
  gsl_histogram *h;
  double result, lower, upper;
  size_t i;

  gsl_set_error_handler (&my_error_handler);

  gsl_ieee_env_setup ();

  status = 0;
  h = gsl_histogram_calloc (0);
  gsl_test (!status, "gsl_histogram_calloc traps zero-length histogram");
  gsl_test (h != 0,
	    "gsl_histogram_calloc returns NULL for zero-length histogram");

  status = 0;
  h = gsl_histogram_calloc_uniform (0, 0.0, 1.0);
  gsl_test (!status,
	    "gsl_histogram_calloc_uniform traps zero-length histogram");
  gsl_test (h != 0,
     "gsl_histogram_calloc_uniform returns NULL for zero-length histogram");

  status = 0;
  h = gsl_histogram_calloc_uniform (10, 1.0, 1.0);
  gsl_test (!status,
	    "gsl_histogram_calloc_uniform traps equal endpoints");
  gsl_test (h != 0,
	    "gsl_histogram_calloc_uniform returns NULL for equal endpoints");

  status = 0;
  h = gsl_histogram_calloc_uniform (10, 2.0, 1.0);
  gsl_test (!status,
	    "gsl_histogram_calloc_uniform traps invalid range");
  gsl_test (h != 0,
	    "gsl_histogram_calloc_uniform returns NULL for invalid range");

  h = gsl_histogram_calloc_uniform (N, 0.0, 1.0);

  status = gsl_histogram_accumulate (h, 1.0, 10.0);
  gsl_test (status != GSL_EDOM, "gsl_histogram_accumulate traps x at xmax");

  status = gsl_histogram_accumulate (h, 2.0, 100.0);
  gsl_test (status != GSL_EDOM, "gsl_histogram_accumulate traps x above xmax");

  status = gsl_histogram_accumulate (h, -1.0, 1000.0);
  gsl_test (status != GSL_EDOM, "gsl_histogram_accumulate traps x below xmin");

  status = gsl_histogram_increment (h, 1.0);
  gsl_test (status != GSL_EDOM, "gsl_histogram_increment traps x at xmax");

  status = gsl_histogram_increment (h, 2.0);
  gsl_test (status != GSL_EDOM, "gsl_histogram_increment traps x above xmax");

  status = gsl_histogram_increment (h, -1.0);
  gsl_test (status != GSL_EDOM, "gsl_histogram_increment traps x below xmin");

  result = gsl_histogram_get (h, N);
  gsl_test (result != 0, "gsl_histogram_get traps index at n");

  result = gsl_histogram_get (h, N + 1);
  gsl_test (result != 0, "gsl_histogram_get traps index above n");

  status = gsl_histogram_get_range (h, N, &lower, &upper);
  gsl_test (status != GSL_EDOM,
	    "gsl_histogram_get_range traps index at n");

  status = gsl_histogram_get_range (h, N + 1, &lower, &upper);
  gsl_test (status != GSL_EDOM,
	    "gsl_histogram_get_range traps index above n");

  status = 0;
  gsl_histogram_find (h, -0.01, &i);
  gsl_test (status != GSL_EDOM, "gsl_histogram_find traps x below xmin");

  status = 0;
  gsl_histogram_find (h, 1.0, &i);
  gsl_test (status != GSL_EDOM, "gsl_histogram_find traps x at xmax");

  status = 0;
  gsl_histogram_find (h, 1.1, &i);
  gsl_test (status != GSL_EDOM, "gsl_histogram_find traps x above xmax");

  gsl_histogram_free (h);

  exit (gsl_test_summary ());
}

void
my_error_handler (const char *reason, const char *file, int line, int err)
{
  if (0)
    printf ("(caught [%s:%d: %s (%d)])\n", file, line, reason, err);
  status = 1;
}


text:
// Use the geodesic solver code to create a gsl implementation so you can
// compare the speeds.  The template needs to get close to ,ideally better than,
// exactly to the gsl run times

#include <stdio.h>
#include <iostream>
#include <chrono>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv2.h>

const double alpha = 1.32754125e20;
int odefunc (double x, const double y[], double f[], void *params)
{
    f[0] = y[3];
    f[1] = y[4];
    f[2] = y[5];
    double r = sqrt(y[0]*y[0]+y[1]*y[1]+y[2]*y[2]);

    f[3] = -(alpha*y[0]) / (r*r*r);
    f[4] = -(alpha*y[1]) / (r*r*r);
    f[5] = -(alpha*y[2]) / (r*r*r);
    return GSL_SUCCESS;
}

int main ()
{
    double AU = 1.496e11; // astronomical unit
    double V = 30000.0; // initial velocity
    // initial conditions
    double y0[] = {0.0,AU,0.0,-V,0.0,0.0};
    double T0 = 0.0;
    double Tf = 7.2e6;
    double hi = 7.2e3/3;
    int dim = 6;
    gsl_odeiv2_system sys = {odefunc, NULL, dim, NULL};

    gsl_odeiv2_driver * d = gsl_odeiv2_driver_alloc_y_new (&sys, gsl_odeiv2_step_rkf45, hi, 1e-5, 0.0);
    auto Nruns = 5000;
    auto totalDuration = 0;
    auto t1 = std::chrono::high_resolution_clock::now();
    double T;
    for (int i = 1; i <= Nruns; i++)
    {
        T0 = 0.0;
        T = Tf;
        y0[0]=0.0;
        y0[1] = AU;
        y0[2] = 0.0;
        y0[3] = -V;
        y0[4] = 0.0;
        y0[5] = 0.0;
        gsl_odeiv2_driver_apply (d, &T0, T, y0);
    }

    auto t2 = std::chrono::high_resolution_clock::now();
    auto duration = std::chrono::duration_cast<std::chrono::microseconds>( t2 - t1 ).count();
    totalDuration += duration;
    std::cout << "GSL runs" << std::endl;
    std::cout << totalDuration/Nruns << std::endl;
    for(int j = 0 ; j < 6 ; j++)
    {
        std::cout << y0[j] << ",";
    }
    std::cout << std::endl;
    gsl_odeiv2_driver_free (d);
    return 0;
}


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

/*
 * The code in this module is based on IQPACK, specifically the LGPL
 * implementation found in HERMITE_RULE:
 * https://people.sc.fsu.edu/~jburkardt/c_src/hermite_rule/hermite_rule.html
 */

#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_sf_gamma.h>

static int
jacobi_check(const size_t n, const gsl_integration_fixed_params * params)
{
  (void) n;

  if (fabs(params->b - params->a) <= GSL_DBL_EPSILON)
    {
      GSL_ERROR("|b - a| too small", GSL_EDOM);
    }
  else if (params->a >= params->b)
    {
      GSL_ERROR("lower integration limit must be smaller than upper limit", GSL_EDOM);
    }
  else if (params->alpha <= -1.0 || params->beta <= -1.0)
    {
      GSL_ERROR("alpha and beta must be > -1", GSL_EDOM);
    }
  else
    {
      return GSL_SUCCESS;
    }
}

static int
jacobi_init(const size_t n, double * diag, double * subdiag, gsl_integration_fixed_params * params)
{
  const double absum = params->beta + params->alpha;
  const double abdiff = params->beta - params->alpha;
  const double a2b2 = absum * abdiff; /* beta^2 - alpha^2 */
  size_t i;

  /* construct the diagonal and subdiagonal elements of Jacobi matrix */
  diag[0] = abdiff/(absum + 2.0);
  subdiag[0] = 2.0*sqrt((params->alpha + 1.0)*(params->beta + 1.0)/(absum + 3.0))/(absum + 2.0);
  for (i = 1; i < n; i++)
    {
      diag[i] = a2b2 / ( (absum + 2.0*i) * (absum + 2.0*i + 2.0) );
      subdiag[i] = sqrt ( 4.0*(i + 1.0) * (params->alpha + i + 1.0) * (params->beta + i + 1.0) * (absum + i + 1.0) / ( pow((absum + 2.0*i + 2.0), 2.0) - 1.0 ) ) / ( absum + 2.0*i + 2.0 );
    }

  params->zemu = pow(2.0, absum + 1.0) * gsl_sf_gamma(params->alpha + 1.0) * gsl_sf_gamma(params->beta + 1.0) / gsl_sf_gamma(absum + 2.0);
  params->shft = 0.5*(params->b + params->a);
  params->slp = 0.5*(params->b - params->a);
  params->al = params->alpha;
  params->be = params->beta;

  return GSL_SUCCESS;
}

static const gsl_integration_fixed_type jacobi_type =
{
  jacobi_check,
  jacobi_init
};

const gsl_integration_fixed_type *gsl_integration_fixed_jacobi = &jacobi_type;

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
////
//
// Copyright (c) 2010-
//      NAKASATO, Naohito
//      All rights reserved.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions
// are met:
// 1. Redistributions of source code must retain the above copyright
//    notice, this list of conditions and the following disclaimer.
// 2. Redistributions in binary form must reproduce the above copyright
//    notice, this list of conditions and the following disclaimer in the
//    documentation and/or other materials provided with the distribution.
//
// THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
// ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
// FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
// DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
// OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
// HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
// LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
// OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
// SUCH DAMAGE.
//
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <malloc.h>
#include <string.h>
#include <sys/time.h>
#include <sys/resource.h>
#include <cblas.h>
#include "accllib.h"

#include "kernel_TN.h"
#include "kernel_NN.h"
#include "kernel_NT.h"
#include "kernel_TT.h"

int dev = 0;

double e_time(void)
{
  static struct timeval now;
  gettimeofday(&now, NULL);
  return (double)(now.tv_sec  + now.tv_usec/1000000.0);
}

double *tmp_A, *tmp_B, *tmp_C, *tmp_D, *tmp_E;

void get_matrix(int dev, double *src, int isrc, double *dma, int idma, int s)
{
  ACCLmemcpy(dev, isrc, idma);
  memcpy(src, dma, s);
}

void set_matrix(int dev, double *src, int isrc, double *dma, int idma, int s)
{
  memcpy(dma, src, s);
  ACCLmemcpy(dev, idma, isrc);
}

void my_dgemm(int dev, int mod, int n)
{
  set_matrix(dev, tmp_A, 0, tmp_E, 6, sizeof(double)*n*n);
  set_matrix(dev, tmp_B, 1, tmp_E, 6, sizeof(double)*n*n);
  set_matrix(dev, tmp_C, 2, tmp_E, 6, sizeof(double)*n*n);

  ACCLrun_on_domain(dev, mod, n/4, n/4);

  get_matrix(dev, tmp_C, 2, tmp_E, 6, sizeof(double)*n*n);
}

void check(double *p1, double *p2, int n)
{
  double norm = 0.0;
  for(int i = 0; i < n*n; i++) {
    double dx = fabs(p1[i] - p2[i]);
    if (dx > norm) {
      norm = dx;
    }
  }
  printf("N = %5i : norm %e\t", n, norm);
}

void out(int n, double t_noio, double t_io, int niter)
{
  fprintf(stdout, "%g %g %g %g\n",
	  t_noio/niter, niter*(double)n*(double)n*(double)n*2.0/t_noio/1.0e9,
	  t_io/niter,   niter*(double)n*(double)n*(double)n*2.0/t_io/1.0e9
	  );
}

int main(int narg, char **argv)
{
  int n;
  int ss;
  int flag = 0;
  int test_mode = 0;
  int mod;
  double alpha, beta;

  ACCLopen();
  ACCLallocate(dev);

  if ( ACCLloadkernel(dev, 0, d_TN, 2) < 0) exit(-1);
  if ( ACCLloadkernel(dev, 1, d_NN, 2) < 0) exit(-1);
  if ( ACCLloadkernel(dev, 2, d_NT, 2) < 0) exit(-1);
  if ( ACCLloadkernel(dev, 3, d_TT, 2) < 0) exit(-1);

  if (narg == 1) {
    n = 512;
  } else {
    n = atoi(argv[1]);
  }
  if (n % 64 != 0) {
    printf("N is not multiple of 64! Error!\n");
    exit(-1);
  }
  if (narg >= 3) {
    test_mode = atoi(argv[2]);
    switch(test_mode) {
    case 0:
      break;
    case 1:
      break;
    case 2:
      break;
    case 3:
      break;
    default:
      printf("Error!!!!\n");
      exit(-1);
    }
  }
  if (narg == 4) flag = 1;

  alpha = 1.0;
  beta = 1.0;

  ss = n*n;

  mod = test_mode;
  ACCLallocatememory(dev, mod, 0, ACCL_DOUBLE_2_2DIM_GPU, n/2, n);
  ACCLallocatememory(dev, mod, 1, ACCL_DOUBLE_2_2DIM_GPU, n/2, n);
  ACCLallocatememory(dev, mod, 2, ACCL_DOUBLE_2_2DIM_GPU | ACCL_GLOBAL_BUFFER, n/2, n);
  ACCLallocateresource0(dev, 3, ACCL_FLOAT_4_1DIM_WRITE, 1, 1);
  ACCLrebindmemory(dev, mod, 3, "cb0");
  ACCLallocateresource0(dev, 4, ACCL_INT_4_1DIM_WRITE, 1, 1);
  ACCLrebindmemory(dev, mod, 4, "cb1");
  ACCLallocateresource0(dev, 5, ACCL_DOUBLE_2_1DIM_WRITE, 1, 1);
  ACCLrebindmemory(dev, mod, 5, "cb2");

  srand(time(NULL));
  tmp_A = (double *)memalign(256, ss*sizeof(double));
  tmp_B = (double *)memalign(256, ss*sizeof(double));
  tmp_C = (double *)memalign(256, ss*sizeof(double));
  tmp_D = (double *)memalign(256, ss*sizeof(double));
  tmp_E = (double *)memalign(256, ss*sizeof(double));
  ACCLallocateresource_pinned(dev, 6, ACCL_DOUBLE_2_2DIM_CPU, n/2, n, tmp_E, ss*sizeof(double));

  for(int i = 0; i < ss; i++) {
    tmp_A[i] = (double)rand()/(RAND_MAX);
    tmp_B[i] = (double)rand()/(RAND_MAX);
    tmp_C[i] = (double)rand()/(RAND_MAX);
  }

  if (flag == 1) {
    memcpy(tmp_D, tmp_C, ss*sizeof(double));

    switch(test_mode) {
    case 0:
      cblas_dgemm(CblasRowMajor, CblasTrans, CblasNoTrans, n, n, n, alpha, tmp_A, n, tmp_B, n, beta, tmp_D, n);
      break;
    case 1:
      cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, n, n, n, alpha, tmp_A, n, tmp_B, n, beta, tmp_D, n);
      break;
    case 2:
      cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasTrans, n, n, n, alpha, tmp_A, n, tmp_B, n, beta, tmp_D, n);
      break;
    case 3:
      cblas_dgemm(CblasRowMajor, CblasTrans, CblasTrans, n, n, n, alpha, tmp_A, n, tmp_B, n, beta, tmp_D, n);
      break;
    }
  }

  float cb[4] = {(float)(n/2), (float)n, 0.0, 0.0};
  ACCLwritememory_float(dev, 3, 4, cb);
  int cb1[4] = {0, n/2, n, 3*n/2};
  ACCLwritememory_int(dev, 4, 4, cb1);
  double cb2[2] = {beta, alpha};
  ACCLwritememory_double(dev, 5, 2, cb2);

  my_dgemm(dev, mod, n);

  if (flag == 1) {
    check(tmp_D, tmp_C, n);
  }

  {
    const int ccc = 5;
    double dum = e_time();
    for(int i = 0; i < ccc; i++) {
      ACCLrun_on_domain(dev, mod, n/4, n/4);
    }
    dum = e_time() - dum;

    double dum2 = e_time();
    for(int i = 0; i < ccc; i++) my_dgemm(dev, mod, n);
    dum2 = e_time() - dum2;

    out(n, dum, dum2, ccc);
  }

  ACCLcleanup(dev);

  return 0;
}


text:
#ifndef __RNG_H
#define __RNG_H

// #include <Rcpp.h>
#include <RcppGSL.h>
#include <gsl/gsl_rng.h>

// [[Rcpp::depends(RcppGSL)]]

using namespace Rcpp;

// Setup RNG for GSL
extern gsl_rng *r_RNG;

unsigned long int random_seed();

//' Generate one sample from a Dirichlet distribution.
//'
//' ## RNG
//'
//' this function uses GSL's RNG seed, unaffected by R's RNG.
//'
//' @param alpha the Dirichlet hyperparameter
//' @param seed the RNG seed: if 0 (default), generate a time-based seed
//' @return a numeric vector
//' @export
//' @family RNG functions
// [[Rcpp::export]]
Rcpp::NumericVector rdirichlet_cpp(const Rcpp::NumericVector &alpha, const unsigned long int seed = 0);

//' Generate from a Dirichlet distribution using the stick breaking definition (safer).
//'
//' ## RNG
//'
//' this function uses R's RNG seed.
//'
//' @param n how many samples to generate
//' @param alpha the Dirichlet hyperparameter, with p entries
//' @return a numeric matrix, n*p
//' @export
//' @family RNG functions
// [[Rcpp::export]]
Rcpp::NumericMatrix rdirichlet_beta_cpp(const unsigned int n, Rcpp::NumericVector alpha);

//' Generate a Dirichlet-Dirichlet-Gamma population (unsafe).
//'
//' Generate samples from m sources and p parameters, n sample per source.
//' The between-source alpha hyperparameter used to generate the source parameters is mandatory.
//'
//' ## RNG
//'
//' this function uses GSL's RNG seed, unaffected by R's RNG.
//'
//' @param n number of samples per source
//' @param m number of sources
//' @param alpha_0 between-source Gamma hyperparameter, a scalar
//' @param beta_0 between-source Gamma hyperparameter, a scalar
//' @param nu_0 between-source Dirichlet hyperparameter, a numeric vector
//' @export
//' @return a matrix with n*m rows
//' @inheritParams rdirichlet_cpp
//' @family RNG functions
// [[Rcpp::export]]
RcppGSL::Matrix rdirdirgamma_cpp(
      const unsigned int &n, const unsigned int &m,
      const double &alpha_0, const double &beta_0,
      const Rcpp::NumericVector &nu_0,
      const unsigned int seed = 0
);

//' Generate a Dirichlet-Dirichlet-Gamma population (safer).
//'
//' Generate samples from m sources and p parameters, n sample per source.
//' The between-source alpha hyperparameter used to generate the source parameters is mandatory.
//'
//' ## RNG
//'
//' this function uses R's RNG seed.
//'
//' @export
//' @return a matrix with n*m rows
//' @inheritParams rdirdirgamma_cpp
//' @family RNG functions
// [[Rcpp::export]]
Rcpp::NumericMatrix rdirdirgamma_beta_cpp(
      const unsigned int &n, const unsigned int &m,
      const double &alpha_0, const double &beta_0,
      const Rcpp::NumericVector nu_0
);

#endif

max_stars_repo_name: lgaborini/rdirdirgamma

max_issues_repo_name: lgaborini/rdirdirgamma

max_forks_repo_name: lgaborini/rdirdirgamma


text:
#include <math.h>
#include <lapacke.h>
#include <complex.h>

//#include "WignerD_fftw.h"

#include <stdio.h>

#define max(a,b) \
   ({ __typeof__ (a) _a = (a); \
       __typeof__ (b) _b = (b); \
     _a > _b ? _a : _b; })

#define min(x,y) (((x) < (y)) ? (x) : (y))

#define MAXFILENAME 100

lapack_complex_double *Iy_Matrix(lapack_int Ndim)
{
double J = (Ndim-1)/2.;
lapack_complex_double *AB = calloc(Ndim*2,sizeof(lapack_complex_double));

lapack_complex_double zero = lapack_make_complex_double(0.0,0.0);
lapack_complex_double myimag = lapack_make_complex_double(0.0,1.0);
//prepare the upper triangle of Iy
//the array is computed in FORTRAN style (transposed)
   double m;
   for (int i=0; i<Ndim; i++){
     m = i-J;
     AB[2*i] =  sqrt((J+m)*(J-m+1.)) / 2.*(-myimag);
	 AB[2*i+1] =  zero;
	 //printf("Hal %f \n" , cimag(AB[i]));
    }

  return AB;
}

lapack_complex_double *Iy_Eigsystem_Lapacke(lapack_int Ndim)
{

char JOBZ, UPLO;
   lapack_int KD, LDAB, LDZ, INFO;

   JOBZ='V'; //= 'V':  Compute eigenvalues and eigenvectors.
   UPLO='U'; // = 'U':  Upper triangle of A is stored;
   KD=1;  //number of superdiagonals
   LDAB=KD+1;  //The leading dimension of the array AB = 2
   LDZ=max(1,Ndim); //The leading dimension of the array eigenvectors

   //int rworksize=max(1,3*Ndim-2);

   //Initialize the array that holds the eigenvalues
    double *Eigenvalue;
	lapack_complex_double *Eigenvec, *WORK; //AB is transposed because of FORTRAN array struct.

	Eigenvec = calloc(Ndim*Ndim,sizeof(lapack_complex_double));
	WORK = calloc(Ndim,sizeof(lapack_complex_double));
	Eigenvalue = calloc(Ndim,sizeof(double));

    lapack_complex_double *AB = Iy_Matrix(Ndim);  //Create the upper triangle matrix

    INFO = LAPACKE_zhbev(LAPACK_COL_MAJOR, JOBZ, UPLO, Ndim,
                          KD, AB,
                          LDAB, Eigenvalue, Eigenvec,
                          LDZ );

	free(AB);
	free(Eigenvalue);
	free(WORK);
	return Eigenvec;

}

void PrecalculateEigenvec(int Ndim)
{
    lapack_complex_double *Eigenvec = Iy_Eigsystem_Lapacke(Ndim);

	char filename[MAXFILENAME+1];

	snprintf(filename, MAXFILENAME, "../../Calculated/Eigenvectors/EigenvectorsD%d.dat", Ndim);

	FILE* file = fopen(filename, "wb");

	if (file == NULL){
		printf("CANNOT OPEN FILE FOR WRITING: %s\n",filename);
		abort();
	}
	//else printf("WRITING TO FILE: %s\n",filename);

    const int size = sizeof(lapack_complex_double);
	const int fftdim = 2*(Ndim-1)+1;

    for (int m=0; m<Ndim; m++){
	for (int mu=0; mu<Ndim; mu++){
        if (fwrite(&Eigenvec[m + Ndim*mu], 1, size, file) != size){
			printf("ERROR WHILE WRITING TO FILE: %s\n",filename);
			abort();
		}
	}}
	free(Eigenvec);
	fclose(file);
}

int main(){
    int MAXDIM = 120;

	for (int Ndim=2; Ndim<=MAXDIM; Ndim++){
		PrecalculateEigenvec(Ndim);
		printf("Eigenvectors calculated for dimension %d\n", Ndim);
	}
	return 0;
}

max_stars_repo_name: BalintKoczor/fast-spherical-phase-space

max_issues_repo_name: BalintKoczor/fast-spherical-phase-space

max_forks_repo_name: BalintKoczor/fast-spherical-phase-space


text:
/*****************************************************************\
           __
          / /
		 / /                     __  __
		/ /______    _______    / / / / ________   __       __
	   / ______  \  /_____  \  / / / / / _____  | / /      / /
	  / /      | / _______| / / / / / / /____/ / / /      / /
	 / /      / / / _____  / / / / / / _______/ / /      / /
	/ /      / / / /____/ / / / / / / |______  / |______/ /
   /_/      /_/ |________/ / / / /  \_______/  \_______  /
                          /_/ /_/                     / /
			                                         / /
		                                            /_/

  ---------------------------------------------------------------

  Copyright (c) 2007-2011 - Rodrigo Braz Monteiro.
  This file is subject to the terms of halley_license.txt.

\*****************************************************************/

#pragma once

#include <gsl/gsl>
#include "halley/maths/vector2.h"
#include "halley/text/halleystring.h"
#include "halley/resources/resource.h"
#include "halley/file/path.h"
#include "halley/maths/rect.h"
#include "halley/maths/colour.h"

namespace Halley {
	class ResourceDataStatic;
	class ResourceLoader;

	class Image final : public Resource {
	public:
		enum class Format {
			Undefined,
			Indexed,
			RGB,
			RGBA,
			RGBAPremultiplied,
			SingleChannel
		};

		Image(Format format = Format::RGBA, Vector2i size = {});
		Image(gsl::span<const gsl::byte> bytes, Format format = Format::Undefined);
		explicit Image(const ResourceDataStatic& data);
		Image(const ResourceDataStatic& data, const Metadata& meta);
		Image(Image&& other);
		~Image();

		std::unique_ptr<Image> clone();

		void setSize(Vector2i size);

		void load(gsl::span<const gsl::byte> bytes, Format format = Format::Undefined);
		Bytes savePNGToBytes(bool allowDepthReduce = true) const;
		Bytes saveQOIToBytes() const;
		static Vector2i getImageSize(gsl::span<const gsl::byte> bytes);
		static Format getImageFormat(gsl::span<const gsl::byte> bytes);

		static bool isQOI(gsl::span<const gsl::byte> bytes);
		static bool isPNG(gsl::span<const gsl::byte> bytes);

		gsl::span<unsigned char> getPixelBytes();
		gsl::span<const unsigned char> getPixelBytes() const;
		gsl::span<const unsigned char> getPixelBytesRow(int x0, int x1, int y) const;

		int getPixel4BPP(Vector2i pos) const;
		int getPixelAlpha(Vector2i pos) const;
		gsl::span<int> getPixels4BPP();
		gsl::span<const int> getPixels4BPP() const;
		gsl::span<const int> getPixelRow4BPP(int x0, int x1, int y) const;
		size_t getByteSize() const;

		static unsigned int convertRGBAToInt(unsigned int r, unsigned int g, unsigned int b, unsigned int a=255);
		static void convertIntToRGBA(unsigned int col, unsigned int& r, unsigned int& g, unsigned int& b, unsigned int& a);
		static Colour4c convertIntToColour(unsigned int col);

		unsigned int getWidth() const { return w; }
		unsigned int getHeight() const { return h; }
		Vector2i getSize() const { return Vector2i(int(w), int(h)); }

		int getBytesPerPixel() const;
		Format getFormat() const;

		Rect4i getTrimRect() const;
		Rect4i getRect() const;

		void clear(int colour);
		void blitFrom(Vector2i pos, gsl::span<const unsigned char> buffer, size_t width, size_t height, size_t pitch, size_t srcBpp);
		void blitFromRotated(Vector2i pos, gsl::span<const unsigned char> buffer, size_t width, size_t height, size_t pitch, size_t bpp);
		void blitFrom(Vector2i pos, const Image& srcImg, bool rotated = false);
		void blitFrom(Vector2i pos, const Image& srcImg, Rect4i srcArea, bool rotated = false);
		void blitDownsampled(Image& src, int scale);

		void drawImageAlpha(const Image& src, Vector2i pos, uint8_t opacity = 255);
		void drawImageAdd(const Image& src, Vector2i pos, uint8_t opacity = 255);
		void drawImageLighten(const Image& src, Vector2i pos, uint8_t opacity = 255);

		static std::unique_ptr<Image> loadResource(ResourceLoader& loader);
		constexpr static AssetType getAssetType() { return AssetType::Image; }
		void reload(Resource&& resource) override;

		Image& operator=(const Image& o) = delete;
		Image& operator=(Image&& o) = default;

		void serialize(Serializer& s) const;
		void deserialize(Deserializer& s);

		void preMultiply();

		ResourceMemoryUsage getMemoryUsage() const override;

	private:
		std::unique_ptr<unsigned char, void(*)(unsigned char*)> px;
		size_t dataLen = 0;
		unsigned int w = 0;
		unsigned int h = 0;
		Format format = Format::Undefined;
	};

	template <>
	struct EnumNames<Image::Format> {
		constexpr std::array<const char*, 6> operator()() const {
			return{{
				"undefined",
				"indexed",
				"rgb",
				"rgba",
				"rgba_premultiplied",
				"single_channel"
			}};
		}
	};
}

max_stars_repo_name: amrezzd/halley

max_issues_repo_name: amrezzd/halley

max_forks_repo_name: amrezzd/halley


text:
#ifndef SmtkMatcher_h
#define SmtkMatcher_h

/**
 * @file
 * $Revision: 1.1 $
 * $Date: 2009/09/09 23:42:41 $
 *
 *   Unless noted otherwise, the portions of Isis written by the USGS are
 *   public domain. See individual third-party library and package descriptions
 *   for intellectual property information, user agreements, and related
 *   information.
 *
 *   Although Isis has been used by the USGS, no warranty, expressed or
 *   implied, is made by the USGS as to the accuracy and functioning of such
 *   software and related material nor shall the fact of distribution
 *   constitute any such warranty, and no responsibility is assumed by the
 *   USGS in connection therewith.
 *
 *   For additional information, launch
 *   $ISISROOT/doc//documents/Disclaimers/Disclaimers.html
 *   in a browser or see the Privacy & Disclaimers page on the Isis website,
 *   http://isis.astrogeology.usgs.gov, and the USGS privacy and disclaimers on
 *   http://www.usgs.gov/privacy.html.
 */

#include <memory>

#include <QSharedPointer>
#include "SmtkStack.h"
#include "Gruen.h"
#include "SmtkPoint.h"

#include "GSLUtility.h"
#include <gsl/gsl_rng.h>

namespace Isis {

/**
 * @brief Workhorse of stereo matcher
 *
 * This class provides stereo matching functionality to the SMTK toolkit.  It
 * registers points, clones them by adjusting parameters to nearby point
 * locations and manages point selection processes.
 *
 * The Gruen algorithm is initialized here and maintained for use in the
 * stereo matching process.
 *
 * @author 2011-05-28 Kris Becker
 *
 * @internal
 *   @history 2012-12-20 Debbie A. Cook - Removed unused Projection.h
 *                           References #775.
 *   @history 2017-08-18 Summer Stapleton, Ian Humphrey, Tyler Wilson -
 *                           Changed auto_ptr reference to QSharedPointer
 *                           so this class compiles under C++14.
 *                           References #4809.
 */
class SmtkMatcher {
  public:
    SmtkMatcher();
    SmtkMatcher(const QString &regdef);
    SmtkMatcher(const QString &regdef, Cube *lhImage, Cube *rhImage);
    SmtkMatcher(Cube *lhImage, Cube *rhImage);
    ~SmtkMatcher();

    void setImages(Cube *lhImage, Cube *rhImage);
    void setGruenDef(const QString &regdef);

    bool isValid(const Coordinate &pnt);
    bool isValid(const SmtkPoint &spnt);

    /** Return pattern chip */
    Chip *PatternChip() const {
      validate();
      return (m_gruen->PatternChip());
    }

    /** Return search chip */
    Chip *SearchChip() const {
      validate();
      return (m_gruen->SearchChip());
    }

    /** Returns the fit chip */
    Chip *FitChip() const {
      validate();
      return (m_gruen->FitChip());
    }

    void setWriteSubsearchChipPattern(const QString &fileptrn = "SmtkMatcher");

    SmtkQStackIter FindSmallestEV(SmtkQStack &stack);
    SmtkQStackIter FindExpDistEV(SmtkQStack &stack, const double &seedsample,
                                 const double &minEV, const double &maxEV);

    SmtkPoint Register(const Coordinate &lpnt,
                       const AffineRadio &affrad = AffineRadio());
    SmtkPoint Register(const PointPair &pnts,
                       const AffineRadio &affrad = AffineRadio());
    SmtkPoint Register(const SmtkPoint &spnt,
                       const AffineRadio &affrad = AffineRadio());
    SmtkPoint Register(const PointGeometry &lpg, const PointGeometry &rpg,
                       const AffineRadio &affrad  = AffineRadio());

    SmtkPoint Create(const Coordinate &left, const Coordinate &right);
    SmtkPoint Clone(const SmtkPoint &point, const Coordinate &left);

    inline BigInt OffImageErrorCount() const { return (m_offImage);  }
    inline BigInt SpiceErrorCount() const { return (m_spiceErr);  }

    /** Return Gruen template parameters */
    PvlGroup RegTemplate() { return (m_gruen->RegTemplate());  }
    /** Return Gruen registration statistics */
    Pvl RegistrationStatistics() { return (m_gruen->RegistrationStatistics()); }

  private:
    SmtkMatcher &operator=(const SmtkMatcher &matcher); // Assignment disabled
    SmtkMatcher(const SmtkMatcher &matcher);            // Copy const disabled

    Cube     *m_lhCube;                // Left image cube (not owned)
    Cube     *m_rhCube;                // Right image cube (not owned)
    QSharedPointer<Gruen> m_gruen;     // Gruen matcher
    BigInt  m_offImage;                // Offimage counter
    BigInt  m_spiceErr;                // SPICE distance error
    bool    m_useAutoReg;              // Select AutoReg features
    const gsl_rng_type * T;            // GSL random number type
    gsl_rng * r;                       // GSL random number generator

    void randomNumberSetup();
    bool validate(const bool &throwError = true) const;

    inline Camera &lhCamera() { return (*m_lhCube->camera());   }
    inline Camera &rhCamera() { return (*m_rhCube->camera());   }

    Coordinate getLineSample(Camera &camera, const Coordinate &geom);
    Coordinate getLatLon(Camera &camera, const Coordinate &pnt);

    bool inCube(const Camera &camera, const Coordinate &point) const;

    SmtkPoint makeRegisteredPoint(const PointGeometry &left,
                                  const PointGeometry &right, Gruen *gruen);
};
} // namespace Isis

#endif