text:
/**
 * Subroutines to calculate the
 * various contamination models
 *
 */

#include <math.h>

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

#include "fitsio.h"
#include "inout_aper.h"
#include "aXe_grism.h"
#include "spce_PET.h"
#include "spc_wl_calib.h"
#include "aXe_errors.h"
#include "fringe_conf.h"
#include "spc_resp.h"
#include "spce_pathlength.h"
#include "aper_conf.h"
#include "trace_conf.h"
#include "specmodel_utils.h"
#include "model_utils.h"
#include "spc_fluxcube.h"

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

/**
 * The function creates and loads a new flux_cube structure
 * The data is taken from the file specified in the input.
 *
 * @param  fcube_file - the file name of the fluxcube
 *
 * @return ret        - the new flux_cube structure
 */
flux_cube *
load_fluxcube(const char fcube_file[])
{
  flux_cube *fcube;

  //gsl_matrix *test;

  int nflux, n_ext=0;
  int i;

  // get the number of extensions
  n_ext = FITSextnum(fcube_file);
  if (n_ext <3)
    aXe_message(aXe_M_FATAL, __FILE__, __LINE__,
                "aXe_PETCONT: " "the fluxcube file %s has only %i extensions!\n",
                fcube_file, n_ext);

  // determine the number of flux images and allocate
  // the space for the fluxcube structure
  nflux = n_ext-2;
  fcube = alloc_fluxcube(nflux);

  // load XOFFS and YOFFS
  load_offsets(fcube_file, fcube);

  // load the segmentation image into the structure
  fcube->segmentation = load_segmentation(fcube_file);
  //  gsl_to_FITSimage (fcube->segmentation, "gogo.fits", 1, "MOO");

  // load the fluximages
  for (i=0; i < nflux; i++)
    {
      fcube->fluxims[i] = load_fluximage(fcube_file, i+3);
    }

  // fill the number of fluximages
  fcube->n_fimage = nflux;

  // order the fluximages and store the
  // vector with the ordered indices
  fcube->fimage_order = order_fluxims(fcube);

  return fcube;
}

/**
 * The function allocates space for a new
 * fluxcube structure.
 *
 * @param  nflux - the number of fluximages
 *
 * @return ret   - the fluxcube structure
 */
flux_cube *
alloc_fluxcube(const int nflux)
{
   flux_cube *fcube;

   // allocate space for the flux_cube structure
   fcube = (flux_cube *)malloc(sizeof(flux_cube));

   // allocate space for the array of fluximages
   fcube->fluxims      = (flux_image **)malloc(nflux * sizeof(flux_image *));

   // return the result
   return fcube;
}

/**
 * The function fills the keywords XOFFS and YOFFS
 * from the primary header of the fitscube file
 * into the according items of the fitscube structure.
 * In case that those keywords are not present,
 * they are set to 0.
 *
 * @param  fcube_file - the file name of the fluxcube
 * @param  fcube      - the fluxcube structure
 *
 * @return 1          - always true
 */
int
load_offsets(const char fcube_file[], flux_cube *fcube)
{
  char tmp[MAXCHAR];

  double dtmp=0.0;

  // avoid problems with the 'const' declaration
  strcpy(tmp, fcube_file);

  // read XOFFS from the primary header
  dtmp = get_float_from_keyword(tmp, 1, "XOFFS");

  // store the keyword value or 0 (if not present) in the fcube structure
  if (isnan(dtmp))
    fcube->xoffs = 0;
  else
    fcube->xoffs = (int)dtmp;

  // read XOFFS from the primary header
  dtmp = get_float_from_keyword(tmp, 1, "YOFFS");

  // store the keyword value or 0 (if not present) in the fcube structure
  if (isnan(dtmp))
    fcube->yoffs = 0;
  else
    fcube->yoffs = (int)dtmp;

  return 1;
}

/**
 * The function creates and loads the segmentation image
 * of a flux_cube file into a gsl_matrix_int structure.
 *
 * @param  fcube_file - the file name of the fluxcube
 *
 * @return ret        - pointer to the gsl_matrix_int structure
 */
gsl_matrix_int *
load_segmentation(const char fcube_file[])
{

  gsl_matrix_int *segmentation;
  gsl_matrix     *dvals;

  int i,j;

  double diff;

  // load the segmentation image into a gsl_matrix_(double)
  dvals = FITSimage_to_gsl(fcube_file, 2, 1);

  // add 0.5 to get proper rounding with (int)
  gsl_matrix_add_constant (dvals, 0.5);

  // allocate space for the integer matrix
  segmentation = gsl_matrix_int_alloc(dvals->size1,dvals->size2);

  // fill the integer matrix with values from the double matrix
  for (i=0; i < (int)dvals->size1; i++)
    {
      for (j=0; j < (int)dvals->size2; j++)
        {
          gsl_matrix_int_set(segmentation, i, j, (int)gsl_matrix_get(dvals, i, j));

          diff = gsl_matrix_get(dvals, i, j) - (double)gsl_matrix_int_get(segmentation, i, j);
          if (diff > 0.6 || diff < 0.4)
            fprintf(stdout, "diff %f; %f --> %i, (%i, %i)\n", diff, gsl_matrix_get(dvals, i, j), gsl_matrix_int_get(segmentation, i, j), i, j);
        }
    }

  // free the space for the double matrix
  gsl_matrix_free(dvals);

  return segmentation;
}

/**
 * The function creates and loads a fluximage structure
 * from an extention of a fluxcube file.
 *
 * @param  fcube_file - the file name of the fluxcube
 * @param  hdunum     - the extension number to load in a fluximage
 *
 * @return fimage     - pointer to the fluximage structure loaded
 */
flux_image *
load_fluximage(const char fcube_file[], int hdunum)
{
  flux_image *fimage;

  char tmp[MAXCHAR];

  double wavelength;

  strcpy(tmp, fcube_file);

  // allocate space for the fluximage structure
  fimage = (flux_image *)malloc(sizeof(flux_image));

  // load the array from the specified extension number
  fimage->flux = FITSimage_to_gsl(fcube_file, hdunum, 1);

  // get the keyword "WAVELENG' from the specified extension number
  wavelength = (double)get_float_from_keyword(tmp, hdunum, "WAVELENG");

  // report an error in case that the extension does not have that keyword.
  // Otherwise store its content in the fluximage structure
  if (isnan(wavelength))
    aXe_message(aXe_M_FATAL, __FILE__, __LINE__,
                "aXe_PETCONT: " "Missing keyword WAVELENGTH in fluxcube %s, extension %i!\n",
                fcube_file, hdunum);
  else
    fimage->wavelength = wavelength;

  return fimage;
}

/**
 * The function analyzes the fluximages in a
 * fluxcube structure and sorts their indices in
 * ascending wavelength. The vector with the
 * ordered indices is given back.
 *
 * @param  fcube - the pointer to the fluxcube structure
 *
 */
gsl_vector_int *
order_fluxims(flux_cube *fcube)
{
  gsl_vector     *fimage_wavelength;

  gsl_vector_int *fimage_order;

  int i, j, k;

  // allocate space for the order vector
  fimage_order = gsl_vector_int_alloc(fcube->n_fimage);

  // put a default order in the vector
  for (i=0; i<fcube->n_fimage; i++)
    gsl_vector_int_set(fimage_order, i, i);

  // order the wavelength
  if (fcube->n_fimage > 1)
    {

      // create a vector for the wavelengths
      fimage_wavelength = gsl_vector_alloc(fcube->n_fimage);

      // set the first entry for the wavelength
      gsl_vector_set(fimage_wavelength, 0, fcube->fluxims[0]->wavelength);

      // go over all fluximages
      for (i=1; i<fcube->n_fimage; i++)
        {

          // find the correct position of the fluximage
          j=0;
          while (j < i && gsl_vector_get(fimage_wavelength, j) <= fcube->fluxims[i]->wavelength)
            j++;

          // check whether the correct poition of the fluximage is at the end
          if (j == i)
            {

              // append the fluximage at the end
              gsl_vector_set(fimage_wavelength, j, fcube->fluxims[i]->wavelength);
              gsl_vector_int_set(fimage_order, j, i);
            }
          else
            {

              // move all fluximages one index up
              for (k=i; k > j; k--)
                {
                  gsl_vector_set(fimage_wavelength, k, gsl_vector_get(fimage_wavelength, k-1));
                  gsl_vector_int_set(fimage_order, k, gsl_vector_int_get(fimage_order, k-1));
                }

              // place the flux image at the correct position
              gsl_vector_set(fimage_wavelength, j, fcube->fluxims[i]->wavelength);
              gsl_vector_int_set(fimage_order, j, i);
            }
        }

      // release the space for the wavelength vector
      gsl_vector_free(fimage_wavelength);
    }

  // return the order vector
  return fimage_order;
}

/**
 * The function releases the space allocated
 * for a fluxcube structure.
 *
 * @param  fcube - the pointer to the fluxcube structure
 *
 */
void
free_fluxcube(flux_cube *fcube)
{
  int i;

  // free the segmentation matrix
  gsl_matrix_int_free(fcube->segmentation);

  // free the wavelength order vector
  gsl_vector_int_free(fcube->fimage_order);

  // go over the fluximages, and free them
  for (i=0; i < fcube->n_fimage; i++)
    free_fluximage(fcube->fluxims[i]);

  // free the fluximage vector
  free(fcube->fluxims);

  // free the fluximage order vector
  //  gsl_vector_int_free(fcube->fimage_order);

  // free the fluxcube
  free(fcube);

  // set the fluxcube to NULL
  fcube=NULL;
}

/**
 * The function releases the space allocated
 * for a fluximage structure.
 *
 * @param  fimage - the pointer to the fluximage structure
 *
 */
void
free_fluximage(flux_image *fimage)
{
  // free the flux matrix
  gsl_matrix_free(fimage->flux);

  // free the fluximage
  free(fimage);

  // set the fluximage to NULL
  fimage=NULL;
}

/**
 * The function transforms a pixel coordinate point in a
 * flt image into an pixel coordinate point in the associated
 * fluxcube image
 *
 * @param  fcube       - the fluxcube structure
 * @param  flt_point   - pixel coordinate point in the flt-image
 *
 * @return fcube_point - the pixel coordinate in the fluxcube
 */
d_point
flt_to_fcube_trans(const flux_cube *fcube, const d_point flt_point)
{
  d_point fcube_point;

  fcube_point.x = flt_point.x - (double)fcube->xoffs;
  fcube_point.y = flt_point.y - (double)fcube->yoffs;

  return fcube_point;
}

/**
 * The function transforms an image coordinate point in a
 * fluxcube into an image coo point in the associated
 * flt image
 *
 * @param  fcube       - the fluxcube structure
 * @param  fcube_point - image-coo in the fcube
 *
 * @return flt_point   - the image coo in the flt
 */
d_point
fcube_to_flt_trans(const flux_cube *fcube, const d_point fcube_point)
{
  d_point flt_point;

  flt_point.x = fcube_point.x + (double)fcube->xoffs;
  flt_point.y = fcube_point.y + (double)fcube->yoffs;

  return flt_point;
}

/**
 * The function converts the information stored in the segmenation
 * matrix of a fluxcube and the object list to a list of direct objects
 *
 * @param  fcube   - the fluxcube structure
 * @param  oblist  - the object list
 *
 * @return dirlist - the dirobject list created
 */
dirobject **
fluxcube_to_dirlist(const flux_cube *fcube, object  **oblist)
{
  dirobject **dirlist;
  dirobject  *actdir;

  //px_point fcube_point;
  px_point flt_point;

  d_point tmp1, tmp2;

  int iact=0;
  int nobjects=0;
  int ndir=0;
  int i, j;
  int objindex;

  //int itmp=0;

  // determine the number of objects in the object list
  nobjects = object_list_size(oblist);

  // allocate space for the dirobject list
  dirlist = (dirobject **) malloc((nobjects+1) * sizeof(dirobject *));
  if (dirlist == NULL)
    aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
                 "fluxcube_to_dirlist:" " Could not allocate"
                 " memory for pointers to %i dirobject objects", nobjects+1);

  // initialize the number of direct object
  // and set the first item to NULL
  ndir=0;
  dirlist[ndir]=NULL;

  //  if (nobjects < 1)
  //    return dirlist;

  // go over each column in the segmentation image
  for (i=0; i < (int)fcube->segmentation->size1; i++)
    {
      // go over each row in the segmentation image
      for (j=0; j < (int)fcube->segmentation->size2; j++)
        {

          // get the pixel value in the segmentation image
          iact = gsl_matrix_int_get(fcube->segmentation, i, j);

          // do something if the pixel is part of an object
          if (iact > 0)
            {

              // transform the cooridnates from the fcube-system
              // into the flt system
              tmp1.x = (double)i;
              tmp1.y = (double)j;
              tmp2 = fcube_to_flt_trans(fcube, tmp1);
              flt_point.x = (int)tmp2.x;
              flt_point.y = (int)tmp2.y;

              // try to get a dirobject for that object (=pixel value)
              actdir = get_dirobject_from_list(dirlist, iact);

              // check whether a dirobject just exists
              if (actdir != NULL)
                {
                  // udate an existing object
                  update_dirobject(actdir, flt_point);
                }
              else
                {

                  // find the the object index for a (hypotetical) new dirobject
                  objindex = find_object_in_object_list(oblist, iact);
                  // check whether object exists in oblist
                  if (objindex >-1)
                    {
                      // create a new dirobject and enhance the counter
                      dirlist[ndir++] = dirobject_from_segpoint(flt_point, oblist[objindex]);

                      // set the actually last dirobject to NULL
                      // to get a properly closed list at all time
                      dirlist[ndir]=NULL;
                    }
                }
            }
        }
    }

  // return the direct object list
  return dirlist;
}

/**
 * The function creates a a dirobject on the basis
 * of an object and the pixel coos of a point
 * which is part of that object.
 *
 * @param  flt_point - the fluxcube structure
 * @param  actobject - the object list
 *
 * @return dirobject - the dirobject created
 */
dirobject *
dirobject_from_segpoint(const px_point flt_point, const object  *actobject)
{

  dirobject  *actdir;

  // allocate space for the dirobject
  actdir = (dirobject *) malloc (sizeof (dirobject));
  if (actdir == NULL)
    aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
                 "fill_dirobject:" " Could not allocate"
                 " memory for a dirobject object");

  // transfer refpoint and ID
  actdir->ID = actobject->ID;
  actdir->refpoint.x = actobject->beams[0].refpoint.x;
  actdir->refpoint.y = actobject->beams[0].refpoint.y;

  // derive and store min/max in x/y for the corners
  actdir->ix_min = flt_point.x;
  actdir->ix_max = flt_point.x;
  actdir->iy_min = flt_point.y;
  actdir->iy_max = flt_point.y;

  // derive the distortion scales along the major an minor axis
  actdir->drzscale.x = 1.0;
  actdir->drzscale.y = 1.0;

  // set the variable
  actdir->bb_sed = 0;

  // make it a dummy
  actdir->SED = NULL;

  // return the new dirobject
  return actdir;
}

/**
 * The function updates a dirobject with a point
 * belonging to the same object. The new point
 * may extend the object size.
 *
 * @param actdir    - the direct object
 * @param flt_point - the new coordinate point
 *
 */
void
update_dirobject(dirobject *actdir, const px_point flt_point)
{
  // check whether the new points extends the
  // old min/max values in x or y.
  // correct the values if necessary.
  actdir->ix_min = MIN(actdir->ix_min, flt_point.x);
  actdir->ix_max = MAX(actdir->ix_max, flt_point.x);
  actdir->iy_min = MIN(actdir->iy_min, flt_point.y);
  actdir->iy_max = MAX(actdir->iy_max, flt_point.y);
}

/**
 * The function computes and stores at the according position
 * the xy offsets of a direct image. This is only important
 * in the fluxcube emission model, since the image coordinates
 * in fluxcubes are 'filter' coordinates and do not take into
 * account the offsets introduce by the grism.
 * For each object and beam those offsets are evaluated
 * and stored in the direct object structure.
 *
 * @param dirlist   - the direct object
 * @param CONF_file - the new coordinate point
 *
 */
void fill_xy_offsets(dirobject **dirlist, char CONF_file[])
{
  aperture_conf   *conf;

  gsl_vector *x_coeffs;
  gsl_vector *y_coeffs;

  d_point m_point;

  //double xoffs, yoffs;

  int beamID=0;
  int i;

  // load the configuration file
  conf = get_aperture_descriptor (CONF_file);

  // go over each beam defined in the configuration
  // file
  for (beamID=0; beamID < conf->nbeams; beamID++)
    {

      // determine the coeeficients for the offsets
      // in x and in y
      x_coeffs = get_beam_trace_xoff (CONF_file, beamID);
      y_coeffs = get_beam_trace_yoff (CONF_file, beamID);

      // go over each direct object
      i=0;
      while (dirlist[i] !=NULL)
        {

          // get the mean direct object position
          m_point = get_dirobject_meanpos(dirlist[i]);

          // evaluate the coefficients for the 2D variable offsets
          // at the mean position
          dirlist[i]->xy_off[beamID].x = eval_trace_off_at_pos (x_coeffs, m_point, beamID);
          dirlist[i]->xy_off[beamID].y = eval_trace_off_at_pos (y_coeffs, m_point, beamID);

          // iterate the counter
          i++;
        }

      // free the vectors for the coefficients
      gsl_vector_free(x_coeffs);
      gsl_vector_free(y_coeffs);
    }
        // release memory
        free_aperture_conf(conf);
}

/**
 * The function fills the flux information from a point in
 * a fluxcube structure into the flux part of a direct object.
 *
 * @param fcube  - the fluxcube structure
 * @param point  - the position to take the flux from
 * @param actdir - the direct object to update
 *
 */
void fill_fluxvalues(const flux_cube *fcube, const px_point point,
                     dirobject *actdir, const int inter_type)
{

  energy_distrib *sed;
  double *sed_wavs;
  double *sed_flux;

  int i;
  int iact;

  sed = (energy_distrib*) malloc(sizeof(energy_distrib));
  sed_wavs = (double*) malloc(fcube->n_fimage*sizeof(double));
  sed_flux = (double*) malloc(fcube->n_fimage*sizeof(double));

  if (actdir->SED)
    free_enerdist (actdir->SED);

  // go over each fluximage in the fluxcube structure
  for (i=0; i < fcube->n_fimage; i++)
    {

      // determine the index of the next fluximage
      iact = gsl_vector_int_get(fcube->fimage_order, i);

      // transfer the wavelength and flux from the fluxcube
      sed_wavs[i] = fcube->fluxims[iact]->wavelength;
      sed_flux[i] = gsl_matrix_get(fcube->fluxims[iact]->flux, point.x, point.y);
    }

  sed->npoints    = fcube->n_fimage;
  sed->wavelength = sed_wavs ;
  sed->flux       = sed_flux;

  if (fcube->n_fimage > 1)
    {
      if (inter_type == 2)
        sed->interp     = gsl_interp_alloc (gsl_interp_polynomial, (size_t)fcube->n_fimage);
      else if (inter_type == 3)
        sed->interp     = gsl_interp_alloc (gsl_interp_cspline, (size_t)fcube->n_fimage);
      else
        sed->interp     = gsl_interp_alloc (gsl_interp_linear, (size_t)fcube->n_fimage);

      sed->accel      = gsl_interp_accel_alloc ();
      gsl_interp_init (sed->interp, sed->wavelength, sed->flux, (size_t)sed->npoints);
    }
  else
    {
      sed->interp = NULL;
      sed->accel= NULL;
    }

  // transfer the SED;
  // announce the SED
  actdir->SED    = sed;
  actdir->bb_sed = 1;
}

/**
 * The function prints the main properties of a
 * fluxcube structure onto the screen.
 * This is mostly for debuggin/development
 * purposes, and not for production
 *
 * @param fcube      - the fluxcube structure
 *
 */
void
print_fluxcube(const flux_cube *fcube)
{
  int i;

  fprintf(stdout, "XOFFS,YOFFS: %i,%i\n", fcube->xoffs, fcube->yoffs);
  fprintf(stdout, "Number of flux-extension: %i\n", fcube->n_fimage);

  for (i=0; i<fcube->n_fimage; i++)
    fprintf(stdout, "Wavelength for flux extension %i: %f\n", i,
            fcube->fluxims[gsl_vector_int_get(fcube->fimage_order,i)]->wavelength);
}

max_stars_repo_name: sosey/pyaxe

max_issues_repo_name: sosey/pyaxe

max_forks_repo_name: sosey/pyaxe


text:
/*
 *  gslrandomgen.h
 *
 *  This file is part of NEST.
 *
 *  Copyright (C) 2004 The NEST Initiative
 *
 *  NEST 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.
 *
 *  NEST 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 NEST.  If not, see <http://www.gnu.org/licenses/>.
 *
 */

/*
 *  Interface to GSL Random Number Generators
 *
 */

#ifndef GSLRANDOMGEN_H
#define GSLRANDOMGEN_H

// C++ includes:
#include <cassert>
#include <list>
#include <string>

// Generated includes:
#include "config.h"

// Includes from librandom:
#include "random_datums.h"
#include "randomgen.h"

// Includes from sli:
#include "dictdatum.h"

// essential GSL includes or replacements
// GSL Versions < 1.2 have weak MT seeding
#ifdef HAVE_GSL

// "Real" version in presence of GSL

// External includes:
#include <gsl/gsl_rng.h>

namespace librandom
{

/**
 * class GslRandomGen
 * C++ wrapper for GSL/GSL-style generators.
 * @note
 * This class should only be used within librandom.
 *
 * @ingroup RandomNumberGenerators
 */

class GslRandomGen : public RandomGen
{
  friend class GSL_BinomialRandomDev;

public:
  explicit GslRandomGen( const gsl_rng_type*, //!< given RNG, given seed
    unsigned long );

  ~GslRandomGen();

  //! Add all GSL RNGs to rngdict
  static void add_gsl_rngs( Dictionary& );

  RngPtr
  clone( unsigned long s )
  {
    return RngPtr( new GslRandomGen( rng_type_, s ) );
  }

private:
  void seed_( unsigned long );
  double drand_( void );

private:
  gsl_rng_type const* rng_type_;
  gsl_rng* rng_;
};

inline void
GslRandomGen::seed_( unsigned long s )
{
  gsl_rng_set( rng_, s );
}

inline double
GslRandomGen::drand_( void )
{
  return gsl_rng_uniform( rng_ );
}

//! Factory class for GSL-based random generators
class GslRNGFactory : public GenericRNGFactory
{
public:
  GslRNGFactory( gsl_rng_type const* const );
  RngPtr create( unsigned long ) const;

private:
  //! GSL generator type information
  gsl_rng_type const* const gsl_rng_;
};
}

#else

// NO GSL Available---Implement class as empty shell
namespace librandom
{

class GslRandomGen : public RandomGen
{
public:
  //! Add all GSL RNGs to rngdict
  //! Do nothing if GSL not available
  static void
  add_gsl_rngs( Dictionary& )
  {
  }

private:
  GslRandomGen()
  {
    assert( false );
  }
  ~GslRandomGen()
  {
    assert( false );
  }
};
}

#endif

#endif

ext: h

max_stars_repo_name: OpenHEC/SNN-simulator-on-PYNQcluster

max_issues_repo_name: zlchai/SNN-simulator-on-PYNQcluster

max_forks_repo_name: OpenHEC/SNN-simulator-on-PYNQcluster


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

#ifndef SCILIB_LINALG_BLAS2_MATRIX_VECTOR_PRODUCT_H
#define SCILIB_LINALG_BLAS2_MATRIX_VECTOR_PRODUCT_H

#ifdef USE_MKL
#include <mkl.h>
#else
#include <cblas.h>
#endif

#include <scilib/mdarray.h>
#include <scilib/linalg_impl/lapack_types.h>
#include <experimental/mdspan>
#include <cassert>
#include <complex>
#include <type_traits>
#include <iostream>

namespace Sci {
namespace Linalg {

namespace stdex = std::experimental;

template <class T_a,
          stdex::extents<>::size_type nrows_a,
          stdex::extents<>::size_type ncols_a,
          class Layout_a,
          class Accessor_a,
          class T_x,
          stdex::extents<>::size_type ext_x,
          class Layout_x,
          class Accessor_x,
          class T_y,
          stdex::extents<>::size_type ext_y,
          class Layout_y,
          class Accessor_y>
    requires(!std::is_const_v<T_y>)
inline void matrix_vector_product(
    stdex::mdspan<T_a, stdex::extents<nrows_a, ncols_a>, Layout_a, Accessor_a>
        a,
    stdex::mdspan<T_x, stdex::extents<ext_x>, Layout_x, Accessor_x> x,
    stdex::mdspan<T_y, stdex::extents<ext_y>, Layout_y, Accessor_y> y)
{
    static_assert(x.static_extent(0) == a.static_extent(1));
    using size_type = stdex::extents<>::size_type;

    for (size_type i = 0; i < a.extent(0); ++i) {
        y(i) = T_y{0};
        for (size_type j = 0; j < a.extent(1); ++j) {
            y(i) += a(i, j) * x(j);
        }
    }
}

template <class Layout>
inline void matrix_vector_product(Sci::Matrix_view<double, Layout> a,
                                  Sci::Vector_view<double, Layout> x,
                                  Sci::Vector_view<double, Layout> y)
{
    static_assert(x.static_extent(0) == a.static_extent(1));

    constexpr double alpha = 1.0;
    constexpr double beta = 0.0;

    const BLAS_INT m = static_cast<BLAS_INT>(a.extent(0));
    const BLAS_INT n = static_cast<BLAS_INT>(a.extent(1));
    const BLAS_INT incx = static_cast<BLAS_INT>(x.stride(0));
    const BLAS_INT incy = static_cast<BLAS_INT>(y.stride(0));

    auto matrix_layout = CblasRowMajor;
    BLAS_INT lda = n;

    if constexpr (std::is_same_v<Layout, stdex::layout_left>) {
        matrix_layout = CblasColMajor;
        lda = m;
    }
    cblas_dgemv(matrix_layout, CblasNoTrans, m, n, alpha, a.data(), lda,
                x.data(), incx, beta, y.data(), incy);
}

template <class Layout>
inline void matrix_vector_product(Sci::Matrix_view<const double, Layout> a,
                                  Sci::Vector_view<const double, Layout> x,
                                  Sci::Vector_view<double, Layout> y)
{
    static_assert(x.static_extent(0) == a.static_extent(1));

    constexpr double alpha = 1.0;
    constexpr double beta = 0.0;

    const BLAS_INT m = static_cast<BLAS_INT>(a.extent(0));
    const BLAS_INT n = static_cast<BLAS_INT>(a.extent(1));
    const BLAS_INT incx = static_cast<BLAS_INT>(x.stride(0));
    const BLAS_INT incy = static_cast<BLAS_INT>(y.stride(0));

    auto matrix_layout = CblasRowMajor;
    BLAS_INT lda = n;

    if constexpr (std::is_same_v<Layout, stdex::layout_left>) {
        matrix_layout = CblasColMajor;
        lda = m;
    }
    cblas_dgemv(matrix_layout, CblasNoTrans, m, n, alpha, a.data(), lda,
                x.data(), incx, beta, y.data(), incy);
}

#ifdef USE_MKL
// Does not work with OpenBLAS version 0.2.14.1
template <class Layout>
inline void
matrix_vector_product(Sci::Matrix_view<std::complex<double>, Layout> a,
                      Sci::Vector_view<std::complex<double>, Layout> x,
                      Sci::Vector_view<std::complex<double>, Layout> y)
{
    static_assert(x.static_extent(0) == a.static_extent(1));

    constexpr std::complex<double> alpha = {1.0, 0.0};
    constexpr std::complex<double> beta = {0.0, 0.0};

    const BLAS_INT m = static_cast<BLAS_INT>(a.extent(0));
    const BLAS_INT n = static_cast<BLAS_INT>(a.extent(1));
    const BLAS_INT incx = static_cast<BLAS_INT>(x.stride(0));
    const BLAS_INT incy = static_cast<BLAS_INT>(y.stride(0));

    auto matrix_layout = CblasRowMajor;
    BLAS_INT lda = n;

    if constexpr (std::is_same_v<Layout, stdex::layout_left>) {
        matrix_layout = CblasColMajor;
        lda = m;
    }
    cblas_zgemv(matrix_layout, CblasNoTrans, m, n, &alpha, a.data(), lda,
                x.data(), incx, &beta, y.data(), incy);
}

template <class Layout>
inline void
matrix_vector_product(Sci::Matrix_view<const std::complex<double>, Layout> a,
                      Sci::Vector_view<const std::complex<double>, Layout> x,
                      Sci::Vector_view<std::complex<double>, Layout> y)
{
    static_assert(x.static_extent(0) == a.static_extent(1));

    constexpr std::complex<double> alpha = {1.0, 0.0};
    constexpr std::complex<double> beta = {0.0, 0.0};

    const BLAS_INT m = static_cast<BLAS_INT>(a.extent(0));
    const BLAS_INT n = static_cast<BLAS_INT>(a.extent(1));
    const BLAS_INT incx = static_cast<BLAS_INT>(x.stride(0));
    const BLAS_INT incy = static_cast<BLAS_INT>(y.stride(0));

    auto matrix_layout = CblasRowMajor;
    BLAS_INT lda = n;

    if constexpr (std::is_same_v<Layout, stdex::layout_left>) {
        matrix_layout = CblasColMajor;
        lda = m;
    }
    cblas_zgemv(matrix_layout, CblasNoTrans, m, n, &alpha, a.data(), lda,
                x.data(), incx, &beta, y.data(), incy);
}
#endif

template <class T, class Layout, class Allocator>
inline Sci::Vector<T, Layout, Allocator>
matrix_vector_product(const Sci::Matrix<T, Layout, Allocator>& a,
                      const Sci::Vector<T, Layout, Allocator>& x)
{
    Sci::Vector<T, Layout, Allocator> res(a.extent(0));
    matrix_vector_product(a.view(), x.view(), res.view());
    return res;
}

template <class T, class Layout, class Allocator>
inline void matrix_vector_product(const Sci::Matrix<T, Layout, Allocator>& a,
                                  const Sci::Vector<T, Layout, Allocator>& x,
                                  Sci::Vector<T, Layout, Allocator>& res)
{
    assert(res.size() == a.extent(0));
    matrix_vector_product(a.view(), x.view(), res.view());
}

} // namespace Linalg
} // namespace Sci

#endif // SCILIB_LINALG_BLAS2_MATRIX_VECTOR_PRODUCT_H

max_stars_repo_name: stigrs/scilib

max_issues_repo_name: stigrs/scilib

max_forks_repo_name: stigrs/scilib


text:
/****************************************************************************
 *                                                                          *
 * Copyright (C) 2001 ~ 2016 Neutrino International Inc.                    *
 *                                                                          *
 * Author : Brian Lin <lin.foxman@gmail.com>, Skype: wolfram_lin            *
 *                                                                          *
 * QtFFT acts as an interface between Qt and FFT libraries.                 *
 * Please keep QtFFT as simple as possible.                                 *
 *                                                                          *
 * Qt Version   : 5.4.1                                                     *
 * CIOS Version : 1.6.0                                                     *
 *                                                                          *
 ****************************************************************************/

#ifndef QT_FFT_H
#define QT_FFT_H

#include <QtCore>
#include <QtNetwork>
#include <QtSql>
#ifndef QT_STATIC
#include <QtScript>
#endif

#include <Essentials>

#ifdef QT_QTGMP_LIB
#include <QtGMP>
#endif

#ifdef QT_QTGSL_LIB
#include <QtGSL>
#endif

#ifdef QT_QTCUDA_LIB
#include <QtCUDA>
#endif

extern "C" {

#include <gsl/gsl_wavelet.h>
#include <gsl/gsl_wavelet2d.h>
#include <gsl/gsl_fft.h>
#include <gsl/gsl_fft_complex.h>
#include <gsl/gsl_fft_complex_float.h>
#include <gsl/gsl_fft_halfcomplex.h>
#include <gsl/gsl_fft_halfcomplex_float.h>
#include <gsl/gsl_fft_real.h>
#include <gsl/gsl_fft_real_float.h>
#include <QtFFT/fftw3.h>

}

QT_BEGIN_NAMESPACE

#ifndef QT_STATIC
#    if defined(QT_BUILD_QTFFT_LIB)
#      define Q_FFT_EXPORT Q_DECL_EXPORT
#    else
#      define Q_FFT_EXPORT Q_DECL_IMPORT
#    endif
#else
#    define Q_FFT_EXPORT
#endif

namespace N
{

class Q_FFT_EXPORT Wavelet        ;
class Q_FFT_EXPORT Hankel         ;
class Q_FFT_EXPORT AbstractFFT    ;
class Q_FFT_EXPORT CooleyFFT      ;
class Q_FFT_EXPORT KissFFT        ;
class Q_FFT_EXPORT GSLFFT         ;
class Q_FFT_EXPORT FFTW           ;
class Q_FFT_EXPORT AccelerateFFT  ;

class Q_FFT_EXPORT Wavelet
{
  public:

    gsl_wavelet           * wavelet   ;
    gsl_wavelet_workspace * work      ;
    gsl_wavelet_direction   direction ;

    explicit Wavelet             (void) ;
    virtual ~Wavelet             (void) ;

    bool    allocate             (const gsl_wavelet_type * T,size_t k) ;
    bool    workspace            (size_t n) ;
    QString Name                 (void) ;

    int Transform                (double * data,size_t stride,size_t n) ;
    int TransformForward         (double * data,size_t stride,size_t n) ;
    int TransformInverse         (double * data,size_t stride,size_t n) ;
    int Transform2d              (double * data,size_t tda,size_t size1,size_t size2) ;
    int Transform2dForward       (double * data,size_t tda,size_t size1,size_t size2) ;
    int Transform2dInverse       (double * data,size_t tda,size_t size1,size_t size2) ;
    int TransformMatrix          (gsl_matrix * m) ;
    int TransformMatrixForward   (gsl_matrix * m) ;
    int TransformMatrixInverse   (gsl_matrix * m) ;
    int nsTransform              (double * data,size_t tda,size_t size1,size_t size2) ;
    int nsTransformForward       (double * data,size_t tda,size_t size1,size_t size2) ;
    int nsTransformInverse       (double * data,size_t tda,size_t size1,size_t size2) ;
    int nsTransformMatrix        (gsl_matrix * m) ;
    int nsTransformMatrixForward (gsl_matrix * m) ;
    int nsTransformMatrixInverse (gsl_matrix * m) ;

  protected:

  private:

};

class Q_FFT_EXPORT Hankel
{
  public:

    gsl_dht * hankel ;

    explicit Hankel   (void) ;
    explicit Hankel   (int size) ;
    explicit Hankel   (int size,double nu, double xmax) ;
             Hankel   (const Hankel & hankel) ;
    virtual ~Hankel   (void) ;

    bool   Allocate   (int size) ;
    bool   New        (int size,double nu, double xmax) ;
    int    Initialize (double nu,double xmax) ;
    int    Apply      (double * In,double * Out) ;
    double x          (int n) ;
    double k          (int n) ;

  protected:

  private:

};

class Q_FFT_EXPORT AbstractFFT
{
  public:

    typedef enum         {
      fftForward  = 0    ,
      fftBackward = 1    }
      TransformDirection ;

    TransformDirection direction  ;
    Cpp::ValueTypes    vType      ;
    int                sourceType ; // 0 - Undecided , 1 - Real , 2 - Complex
    int                targetType ; // 0 - Undecided , 1 - Real , 2 - Complex
    int                dimension  ;

    explicit AbstractFFT        (void);
    virtual ~AbstractFFT        (void);

    virtual int  type           (void) const = 0 ;
    virtual bool isSupported    (Cpp::ValueTypes type,int dimension) const = 0 ;
    virtual int  AllocateSource (CUIDs elements,QByteArray & Source) = 0 ;
    virtual int  AllocateTarget (CUIDs elements,QByteArray & Target) = 0 ;
    virtual bool Prepare        (CUIDs        elements        ,
                                 VarArgs    & Paraments       ,
                                 QByteArray & Source          ,
                                 QByteArray & Target    ) = 0 ;
    virtual void Execution      (void) = 0 ;

  protected:

  private:

};

class Q_FFT_EXPORT CooleyFFT : public AbstractFFT
{
  public:

    explicit CooleyFFT          (void) ;
    virtual ~CooleyFFT          (void) ;

    virtual int  type           (void) const ;
    virtual bool isSupported    (Cpp::ValueTypes T,int dimension) const ;
    virtual int  AllocateSource (CUIDs elements,QByteArray & Source) ;
    virtual int  AllocateTarget (CUIDs elements,QByteArray & Target) ;
    virtual bool Prepare        (CUIDs        elements  ,
                                 VarArgs    & Paraments ,
                                 QByteArray & Source    ,
                                 QByteArray & Target  ) ;
    virtual void Execution      (void) ;

  protected:

  private:

};

class Q_FFT_EXPORT KissFFT : public AbstractFFT
{
  public:

    explicit KissFFT            (void) ;
    virtual ~KissFFT            (void) ;

    virtual int  type           (void) const ;
    virtual bool isSupported    (Cpp::ValueTypes T,int dimension) const ;
    virtual int  AllocateSource (CUIDs elements,QByteArray & Source) ;
    virtual int  AllocateTarget (CUIDs elements,QByteArray & Target) ;
    virtual bool Prepare        (CUIDs        elements  ,
                                 VarArgs    & Paraments ,
                                 QByteArray & Source    ,
                                 QByteArray & Target  ) ;
    virtual void Execution      (void) ;

  protected:

  private:

};

class Q_FFT_EXPORT GSLFFT : public AbstractFFT
{
  public:

    explicit GSLFFT             (void) ;
    virtual ~GSLFFT             (void) ;

    virtual int  type           (void) const ;
    virtual bool isSupported    (Cpp::ValueTypes T,int dimension) const ;
    virtual int  AllocateSource (CUIDs elements,QByteArray & Source) ;
    virtual int  AllocateTarget (CUIDs elements,QByteArray & Target) ;
    virtual bool Prepare        (CUIDs        elements  ,
                                 VarArgs    & Paraments ,
                                 QByteArray & Source    ,
                                 QByteArray & Target  ) ;
    virtual void Execution      (void) ;

  protected:

  private:

};

class Q_FFT_EXPORT FFTW : public AbstractFFT
{
  public:

    explicit       FFTW           (void) ;
    virtual       ~FFTW           (void) ;

    virtual int    type           (void) const ;
    virtual bool   isSupported    (Cpp::ValueTypes T,int dimension) const ;
    virtual int    AllocateSource (CUIDs elements,QByteArray & Source) ;
    virtual int    AllocateTarget (CUIDs elements,QByteArray & Target) ;
    virtual bool   Prepare        (CUIDs        elements  ,
                                   VarArgs    & Paraments ,
                                   QByteArray & Source    ,
                                   QByteArray & Target  ) ;
    virtual void   Execution      (void) ;

    void           Cleanup        (void);

    fftw_complex * array          (int size);
    fftw_complex * fromByteArray  (int & size,QByteArray & Array);
    QByteArray     toByteArray    (fftw_complex * complex,int size);
    void           Execute        (const fftw_plan plan);
    void           Destroy        (fftw_plan plan);

    fftw_plan      Frequency      (int N,double * data,fftw_complex * dout,unsigned int flags = FFTW_MEASURE) ;
    void           Retrieve       (int component,int index,int length,double * data,fftw_complex * dout) ;

    void           Execute        (const fftwf_plan plan);
    void           Destroy        (fftwf_plan plan);

    #ifndef DISABLE_FFTW_LDOUBLE
    void           Execute        (const fftwl_plan plan);
    void           Destroy        (fftwl_plan plan);
    #endif

  protected:

    fftw_plan  fftplan  ;
    fftwf_plan fftfplan ;
    #ifndef DISABLE_FFTW_LDOUBLE
    fftwl_plan fftlplan ;
    #endif

  private:

    int fftSize                   (int fft,int length) ;
    int fftSize                   (int fft,CUIDs elements) ;

};

class Q_FFT_EXPORT AccelerateFFT : public AbstractFFT
{ // Apple iPhone OS Accelerate Framework
  public:

    explicit     AccelerateFFT  (void) ;
    virtual     ~AccelerateFFT  (void) ;

    virtual int  type           (void) const ;
    virtual bool isSupported    (Cpp::ValueTypes T,int dimension) const ;
    virtual int  AllocateSource (CUIDs elements,QByteArray & Source) ;
    virtual int  AllocateTarget (CUIDs elements,QByteArray & Target) ;
    virtual bool Prepare        (CUIDs        elements  ,
                                 VarArgs    & Paraments ,
                                 QByteArray & Source    ,
                                 QByteArray & Target  ) ;
    virtual void Execution      (void) ;

  protected:

  private:

};

}

QT_END_NAMESPACE

Q_DECLARE_METATYPE(N::Wavelet)
Q_DECLARE_METATYPE(N::Hankel)
Q_DECLARE_METATYPE(N::CooleyFFT)
Q_DECLARE_METATYPE(N::KissFFT)
Q_DECLARE_METATYPE(N::GSLFFT)
Q_DECLARE_METATYPE(N::FFTW)
Q_DECLARE_METATYPE(N::AccelerateFFT)

Q_DECLARE_INTERFACE(N::AbstractFFT  , "com.neutrino.math.fft" )

#endif

max_stars_repo_name: Vladimir-Lin/QtFFT

max_issues_repo_name: Vladimir-Lin/QtFFT

max_forks_repo_name: Vladimir-Lin/QtFFT


text:
/**
 * Original Author:
 * * author: Jochen K"upper
 * * created: April 2002
 *
 * author: Pierre Schnizer
 * created: December 2003
 * file: pygsl/src/simanmodule.c
 * $Id: simanmodule.c,v 1.9 2008/10/25 20:45:04 schnizer Exp $
 *
 * Jochen K"upper wrote the original version of this module. In December 2003 I
 * rewrote it. Now I only support the variable type as it is more pythonic.
 *
 */

#include <Python.h>
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <setjmp.h>
#include <gsl/gsl_siman.h>
#include <gsl/gsl_nan.h>
#include <pygsl/error_helpers.h>
#include <pygsl/general_helpers.h>
#include <pygsl/rng.h>
#include <pygsl/rng_helpers.h>

/*
 * Common to all objects of one problem
 * Currently the individual method pointers are not used. Instead the method
 * is resolved every time. Do you know which method is faster? Is it advicable
 * to resolve it once and store it here? Pierre
 */
typedef struct{
	/*
	  PyObject * efunc;
	  PyObject * step;
	  PyObject * metric;
	  PyObject * print;
	*/
	PyObject * rng;
	jmp_buf  buffer;
}pygsl_siman_func_t;

/*
 * The Linked list keeping the reference to the individual objects.
 *
 * Necessary as I have to use longjmp as the current implementation does
 * not allow error propagation yet!
 */
struct _pygsl_siman_t{
	pygsl_siman_func_t *func;
	PyObject * x;
	struct _pygsl_siman_t * prev;
	struct _pygsl_siman_t * next;
};
typedef struct _pygsl_siman_t pygsl_siman_t;

static const char module_doc[] = "C Implementation needed for the siman module.";
static PyObject *module = NULL;
static const char filename[] = __FILE__;

/*
 * Naming of the various methods the object has provide as callbacks.
 */
static const char EFunc_name[]  = "EFunc";
static const char Metric_name[] = "Metric";
static const char Step_name[]   = "Step";
static const char Clone_name[]  = "Clone";
static const char Print_name[]  = "Print";

/*
 * Get a callable method from a object.
 *
 * Basically it is a flat wrapper around PyObject_GetAttrString but
 * checks if the method is callable and adds a Traceback frame
 *
 * Flag usage:
 *           flag == 1 Must exist and must be callable
 *           flag == 2 If it exists it must be callable
 *
 * For the traceback frame all arguments following the module are needed.
 * If you do not know the module, you can pass a NULL pointer.
 */
static PyObject *
PyGSL_get_callable_method(PyObject *o, const char * attr, int flag, PyObject *module,
			  const char * filename, const char * func_name,
			  int lineno)
{
	PyObject * method = NULL;

	FUNC_MESS_BEGIN();
	method = PyObject_GetAttrString(o, (char *) attr);
	if(method == NULL){
		if(flag == 1){
			PyGSL_add_traceback(module, filename, func_name, lineno);
		}else if(flag == 2){
			/* Clear the error otherwise it will show up later on! */
			PyErr_Clear();
		}
		return NULL;
	}
	if(!(PyCallable_Check(method))) {
		/* I must add the method name! I must change it to an more descriptive exception! */
		PyGSL_add_traceback(module, (const char *) filename, func_name, lineno);
		PyErr_SetString(PyExc_TypeError, "Found a attribute which was not callable!"
				"XXX must add the method name!");
		return NULL;
	}
	DEBUG_MESS(2, "Found a method at %p", (void *) method);
	FUNC_MESS_END();
	return method;
}

/* This function type should return the energy of a configuration XP. */
static double
PyGSL_siman_efunc(void *xp)
{

	PyObject *result = NULL, *callback = NULL, *arglist = NULL;
	PyGSL_error_info info;
	pygsl_siman_t *x;
	int flag=GSL_EFAILED;
	double value;
	/* static char *functionname  = __FUNCTION__; */

	FUNC_MESS_BEGIN();

	assert(xp);
	x = (pygsl_siman_t *) xp;

	DEBUG_MESS(2, "Found a pygsl_siman_t at %p and a pygsl_siman_func_t at %p and x at %p",
		   (void *)x, (void *) x->func, (void *) x->x);

	assert(x);
	assert(x->func);

	callback = PyGSL_get_callable_method(x->x, EFunc_name, 1, module, filename, __FUNCTION__, __LINE__);
	if(callback == NULL)
		goto fail;

	info.callback = callback;
	info.message  = __FUNCTION__;
	info.error_description = "and the description ???";
	info.argnum = 1;

	arglist = PyTuple_New(0);
	result = PyEval_CallObject(callback, arglist);
	Py_DECREF(arglist);
	if((flag = PyGSL_CHECK_PYTHON_RETURN(result, 1, &info)) != GSL_SUCCESS){
		PyGSL_add_traceback(module, filename, __FUNCTION__, __LINE__);
		goto fail;
	}
	if((flag = PyGSL_PYFLOAT_TO_DOUBLE(result, &value, &info)) != GSL_SUCCESS){
		PyGSL_add_traceback(module, filename, __FUNCTION__, __LINE__);
		goto fail;
	}
	Py_DECREF(result);
	FUNC_MESS_END();
	return value;

  fail:
	FUNC_MESS("In Fail");
	Py_XDECREF(result);
	longjmp(x->func->buffer, flag);
	return GSL_NAN;
}

/*
 *  This function type should modify the configuration XP using a random step
 *  taken from the generator R, up to a maximum distance of STEP_SIZE.
 */
static void
PyGSL_siman_step(const gsl_rng *r, void *xp, double step_size)
{

	PyObject *result = NULL, *arglist = NULL, *callback = NULL;
	PyGSL_error_info info;
	pygsl_siman_t *x;
	int flag=GSL_EFAILED;

	/* static const char * functionname  = __FUNCTION__; */

	FUNC_MESS_BEGIN();
	x = (pygsl_siman_t *) xp;
	DEBUG_MESS(2, "Found x at %p", xp);

	callback = PyGSL_get_callable_method(x->x, Step_name, 1, module, filename, __FUNCTION__, __LINE__);
	if(callback == NULL)
		goto fail;

	info.callback = callback;
	info.message  = __FUNCTION__;
	info.error_description = "???";
	info.argnum = 1;

	assert(PyGSL_RNG_Check(x->func->rng));
	assert(((PyGSL_rng *) x->func->rng)->rng == r);

	/* create argument list */
	arglist = PyTuple_New(2);
	PyTuple_SET_ITEM(arglist, 0, x->func->rng);
	Py_INCREF(x->func->rng); /* Don't forget tuple is owner! */
	PyTuple_SET_ITEM(arglist, 1, PyFloat_FromDouble(step_size));

	result = PyEval_CallObject(callback, arglist);
	Py_DECREF(arglist);
	if((flag = PyGSL_CHECK_PYTHON_RETURN(result, 0, &info)) != GSL_SUCCESS){
		PyGSL_add_traceback(module, filename, __FUNCTION__, __LINE__);
		goto fail;
	}
	Py_DECREF(result);
	FUNC_MESS_END();
	return;

  fail:
	FUNC_MESS("In Fail");
	Py_XDECREF(result);
	longjmp(x->func->buffer, flag);
	return;
}

/* This function type should return the distance between two configurations XP
   and YP. */
static double
PyGSL_siman_metric(void *xp, void *yp)
{

	PyObject *result = NULL, *arglist = NULL, *callback = NULL;
	PyGSL_error_info info;
	pygsl_siman_t *x, *y;
	int flag=GSL_EFAILED;
	double value;
	/* static const char * functionname = __FUNCTION__; */

	FUNC_MESS_BEGIN();
	x = (pygsl_siman_t *) xp;
	y = (pygsl_siman_t *) yp;

	DEBUG_MESS(2, "Found x at (%p,%p) and y at (%p %p)",
		   (void *) x, (void *)x->x, (void *) y, (void *) y->x);

	assert(x);
	assert(y);
	assert(x->x);
	assert(y->x);

	callback = PyGSL_get_callable_method(x->x, Metric_name, 1, module, filename, __FUNCTION__, __LINE__);
	if(callback == NULL)
		goto fail;

	info.callback = callback;
	info.message  = __FUNCTION__;
	info.error_description = "???";
	info.argnum = 1;

	arglist = PyTuple_New(1);
	PyTuple_SET_ITEM(arglist, 0, y->x);
	Py_INCREF(y->x); /* Tuple is owner! */

	result = PyEval_CallObject(callback, arglist);
	Py_XDECREF(arglist);

	if((flag = PyGSL_CHECK_PYTHON_RETURN(result, 0, &info)) != GSL_SUCCESS){
		PyGSL_add_traceback(module, filename, __FUNCTION__, __LINE__);
		goto fail;
	}
	if((flag = PyGSL_PYFLOAT_TO_DOUBLE(result, &value, &info)) != GSL_SUCCESS){
		PyGSL_add_traceback(module, filename, __FUNCTION__, __LINE__);
		goto fail;
	}
	Py_DECREF(result);
	FUNC_MESS_END();
	return value;

  fail:
	FUNC_MESS("In Fail");
	Py_XDECREF(result);
	longjmp(x->func->buffer, flag);
	return GSL_NAN;
}

/* This function type should print the contents of the configuration XP. */
static void
PyGSL_siman_print(void *xp)
{
	PyObject *result = NULL, *callback=NULL, *arglist=NULL;
	PyGSL_error_info info;
	pygsl_siman_t *x;
	int flag=GSL_EFAILED;
	/* static const char * functionname  = __FUNCTION__; */

	FUNC_MESS_BEGIN();
	x = (pygsl_siman_t *) xp;

	callback = PyGSL_get_callable_method(x->x, Print_name, 1, module, filename, __FUNCTION__, __LINE__);
	if(callback == NULL)
		goto fail;

	info.callback = callback;
	info.message  = __FUNCTION__;
	info.error_description = "what goes here ???";
	info.argnum = 1;

	arglist = PyTuple_New(0);
	result = PyEval_CallObject(callback, arglist);
	Py_DECREF(arglist);

	if((flag = PyGSL_CHECK_PYTHON_RETURN(result, 0, &info)) != GSL_SUCCESS){
		PyGSL_add_traceback(module, (const char*)filename, __FUNCTION__, __LINE__);
		goto fail;
	}
	Py_DECREF(result);
	FUNC_MESS_END();
	return;

  fail:
	FUNC_MESS("In Fail");
	Py_XDECREF(result);
	longjmp(x->func->buffer, flag);
	return;

}

/*
 * As Python is generating the new objects I can not skip the reference
 * counting, but must add a reference for each object in construct and
 * depose it here and in the siman_destroy function.
 *
 * siman_solve will call siman_release_x anyway at the end to clear up
 * all objects. Necessary to non existing error propagation in the siman
 * module.
 */
static void
PyGSL_siman_copy(void *src, void *dst)
{
	PyObject *callback=NULL, *new = NULL, *arglist = NULL;
	PyGSL_error_info info;
	pygsl_siman_t *x, *y;
	int flag=GSL_EFAILED;
	/* static const char * functionname = __FUNCTION__; */

	FUNC_MESS_BEGIN();
	x = (pygsl_siman_t *) src;
	y = (pygsl_siman_t *) dst;

	DEBUG_MESS(2, "Got source at %p, Destination at %p", (void *) x, (void *) y);
	assert(x->x);
	callback = PyGSL_get_callable_method(x->x, Clone_name, 1, module, filename, __FUNCTION__, __LINE__);
	if(callback == NULL)
		goto fail;

	arglist = PyTuple_New(0);
	new = PyEval_CallObject(callback, arglist);
	Py_DECREF(arglist);

	info.callback = callback;
	info.message  = __FUNCTION__;
	info.error_description = "???";
	info.argnum = 1;

	if((flag = PyGSL_CHECK_PYTHON_RETURN(new, 1, &info)) != GSL_SUCCESS){
		PyGSL_add_traceback(module, (const char*)filename, __FUNCTION__, __LINE__);
		goto fail;
	}
	Py_XDECREF(y->x);
	/* Py_INCREF(new);  Necessary? I don't think so. */
	y->x = new;
	FUNC_MESS_END();
	return;

  fail:
	FUNC_MESS("Fail");
	Py_XDECREF(new);
	longjmp(x->func->buffer, flag);

}

static void *
PyGSL_siman_copy_construct(void *new)
{
	pygsl_siman_t * ret, * n, *p;
	int flag = GSL_ENOMEM;

	FUNC_MESS_BEGIN();
	n   = (pygsl_siman_t *) new;

	/* The pointers next and prev are checked against NULL */
	ret = (pygsl_siman_t *) calloc(1, sizeof(pygsl_siman_t));
	DEBUG_MESS(2, "New was %p, Constructed a new object at %p", new, (void *) ret);
	if(ret == NULL){
		pygsl_error("Could not allocate the object for the linked list",
			  filename, __LINE__ - 3, GSL_ENOMEM);
		goto fail;
	}
	/* Put the Link to the old object so that I can clone when I need to copy */
	ret->x    = n->x;
	/* Eventually I will dispose the object */
	Py_INCREF(ret->x);
	/* Pointer to the func struct. This information is the same for all objects */
	ret->func = n->func;

	/* Find the first open object in the linked list ... */
	p = n;
	while(p->next != NULL){
		p = p->next;
	}
	DEBUG_MESS(2, "I found a open object at %p", (void *) p);
	/* and connect the links */
	p->next = ret;
	ret->prev = p;

	FUNC_MESS_END();
	return ret;

  fail:
	FUNC_MESS("Fail");
	longjmp(n->func->buffer, flag);
	return NULL;

}

static void
PyGSL_siman_destroy(void * old)
{
	pygsl_siman_t * o;

	FUNC_MESS_BEGIN();
	o = (pygsl_siman_t *) old;
	assert(o);

	/* fprintf(stderr, "Destroying: Previous object = %p, Next Object = %p\n",
	   (void *)o->prev, (void *)o->next); */

	/* Reconnect the linked list */
	if (o->prev && o->next){
		/* Connect the both */
		o->prev->next = o->next;
		o->next->prev = o->prev;
	}
	else if (o->prev && o->next == NULL){
		/* Prev last element. Terminate the list */
		o->prev->next = NULL;
	}
	else if (o->prev == NULL && o->next == NULL){
		/* Last Element, better to leave it */
		DEBUG_MESS(2, "I do not dispose the last element %p!", (void *) o);
		return;
	}

	/* Dispose the object */
	Py_XDECREF(o->x);
	free(o);
	FUNC_MESS_END();
}

/* Clean up of the linked list of objects */
int
PyGSL_siman_release_x(pygsl_siman_t * myargs, pygsl_siman_t * x)
{
	pygsl_siman_t *p=NULL;
	FUNC_MESS_BEGIN();

	p = myargs;
	/* fprintf(stderr, "Releasing list!\n"); */
	while(1){
		/* fprintf(stderr, "Previous object = %p, Next Object = %p\n",
		   (void *)p->prev, (void *)p->next); */
		/* Don't delete the object containing the result! */
		if(p != x){
			/* fprintf(stderr, "Deleting object at %p\n", (void *) p); */
			PyGSL_siman_destroy((void *) p);
		}
		if(p->next == NULL)
			break;
		p = p->next;
	}
	FUNC_MESS_END();
	return GSL_SUCCESS;
}

static const char pygsl_siman_solve_doc[] =
"Simulated annealing driver.\n\
\n\
Usage:\n\
result = solve(r, x0, ...)\n\
\n\
Input:\n\
    r   ... a random generator from pygsl.rng\n\
    x0  ... a configuration. It must be an object providing the following\n\
            methods:\n\
            EFunc()\n\
            Metric()\n\
            Step()\n\
            Clone()\n\
               If you want to use the print functionality you must provide\n\
               the following method:\n\
            Print()\n\
\n\
Output:\n\
    result ... a object of type x0 with the final value.\n\
\n\
Keywords:\n\
    n_tries       = 200  ... how many points to try for each step\n\
    iters_fixed_T = 10   ... how many iterations at each temperature?\n\
    step_size     = 10   ... max step size in the random walk\n\
\n\
                      parameters for the Boltzmann distribution\n\
    k             = 1.0    ... Boltzmann constant\n\
    t_initial     = 0.002  ... initial temperature\n\
    mu_t          = 1.005  ... damping factor for the temperature\n\
    t_min         = 2.0e-6\n\
\n\
    do_print      = 0      ... print the status of the annealing process\n\
                               (== 0: do not print)\n\
                               ( > 0: print)\n\
";

/* wrapper functions */
static PyObject *
PyGSL_siman_solve(PyObject *self, PyObject *args, PyObject *kw)
{
	PyObject *result = NULL;
	PyObject *efunc = NULL, *step = NULL, *metric = NULL, *clone = NULL,
		 *print = NULL, *r_o = NULL, *x_o=NULL;

	gsl_rng *rng = NULL;

	gsl_siman_print_t    a_print= PyGSL_siman_print;
	gsl_siman_params_t   params = {200, 10, 10.0, 1.0, 0.002, 1.005, 2.0e-6};

	pygsl_siman_func_t   myargs_func = {NULL};
	pygsl_siman_t        myargs = {NULL, NULL, NULL, NULL};

	/* static const char  * functionname = __FUNCTION__; */

	int flag=GSL_EFAILED, do_print=0;
	void * x0 = NULL;

	static const char * kwlist[] = {"rng", "x0", "n_tries", "iters_fixed_T", "step_size", "k",
					"t_initial", "mu_t",  "t_min", "do_print", NULL};

	FUNC_MESS_BEGIN();
	/* python arguments are (rng, x0, settings) */
	if(! PyArg_ParseTupleAndKeywords(args, kw, "OO|iidddddi",  (char **) kwlist, &r_o, &x_o,
			      &params.n_tries, &params.iters_fixed_T, &params.step_size,
			      &params.k, &params.t_initial, &params.mu_t, &params.t_min, &do_print))
		return NULL;

	/* The following methods must exist */
	efunc  = PyGSL_get_callable_method(x_o, EFunc_name,  1, module, filename, __FUNCTION__, __LINE__);
	step   = PyGSL_get_callable_method(x_o, Step_name,   1, module, filename, __FUNCTION__, __LINE__);
	metric = PyGSL_get_callable_method(x_o, Metric_name, 1, module, filename, __FUNCTION__, __LINE__);
	clone  = PyGSL_get_callable_method(x_o, Clone_name,  1, module, filename, __FUNCTION__, __LINE__);
	if( efunc == NULL || step == NULL || metric == NULL || clone == NULL){
		return NULL;
	}

	/* optional print */
	if(do_print == 0){
		a_print = NULL;
	} else {
		print  = PyGSL_get_callable_method(x_o, Print_name,  1, module, filename, __FUNCTION__, __LINE__);
		if(print == NULL){
			DEBUG_MESS(2, "Did not get a print method! print = %p", print);
			a_print = NULL;
			return NULL;
		}
	}

	rng = PyGSL_gsl_rng_from_pyobject(r_o);
	if(rng == NULL){
		return NULL;
	}

	/* initialize/assign functions */
	Py_INCREF(x_o);
	/*
	  myargs_func.efunc  = efunc;
	  myargs_func.step   = step;
	  myargs_func.metric = metric;
	  myargs_func.print  = print;
	*/
	myargs_func.rng    = r_o;

	myargs.func = &myargs_func;
	myargs.x   = x_o;
	myargs.prev = NULL;
	myargs.next = NULL;

	x0 = (void *) &myargs;
	DEBUG_MESS(2, "x0 @ %p; myargs at %p; myargs_func at %p", x0, (void *) &myargs, (void *) &myargs_func);
	DEBUG_MESS(2, "Found a pygsl_siman_t at %p and a pygsl_siman_func_t at %p",
		   (void *) x0,
		   (void *) (((pygsl_siman_t *) x0)->func));

	if((flag = setjmp(myargs_func.buffer)) == 0){
		FUNC_MESS("Starting siman");
		gsl_siman_solve(rng, x0, PyGSL_siman_efunc, PyGSL_siman_step,
				PyGSL_siman_metric, a_print, PyGSL_siman_copy,
				PyGSL_siman_copy_construct, PyGSL_siman_destroy,
				0, /* Only variable mode supported by this wrapper. */
				params);
		FUNC_MESS("End siman");
	}else{
		PyGSL_add_traceback(module, filename, __FUNCTION__, __LINE__);
		goto fail;
	}
	Py_DECREF(x_o);
	DEBUG_MESS(2, "I found x0 at %p", x0);
	result = ((pygsl_siman_t *) x0)->x;

	PyGSL_siman_release_x(&myargs, x0);
	FUNC_MESS_END();
	return result;

  fail:
	FUNC_MESS("In Fail");
	PyGSL_siman_release_x(&myargs, x0);
	Py_XDECREF(x_o);
	PyGSL_error_flag(flag);
	return NULL;
}

/* module initialization */
static PyMethodDef simanMethods[] = {
	{"solve", (PyCFunction) PyGSL_siman_solve,
	 METH_VARARGS | METH_KEYWORDS, (char *) pygsl_siman_solve_doc},
	{NULL, NULL} /* Sentinel */
};

DL_EXPORT(void) init_siman(void)
{
	PyObject *m = NULL;
	FUNC_MESS_BEGIN();
	m = Py_InitModule("_siman", simanMethods);
	module = m;
	init_pygsl();
	import_pygsl_rng();
	FUNC_MESS_END();
	return;
}

/*
 * Local Variables:
 * mode: C
 * c-file-style: "python"
 * End:
 */

ext: c

max_stars_repo_name: juhnowski/FishingRod

max_issues_repo_name: juhnowski/FishingRod

max_forks_repo_name: juhnowski/FishingRod


text:
/////////////////
//example25.5.c
/////////////////
#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv.h>
#include <math.h>

typedef struct{
  double *a; //y=a[0]x[0]+a[1]x[1]+...a[ka]y[ka]+b[0]u[0]+b[1]u[1]+...+b[kb]u[kb]
  //  double *a; //y[t+1]=a[0]y[t]+a[1]y[t-1]+...a[ka]y[t-ka]+b[0]u[t]+b[1]u[t-1]+...+b[kb]u[t-kb]
  double *b;
  double *x;
  double *u;
  int k,kx,ku;
  double xmax;
  double umax;
  double y;
  int _dim;
  double *_y;
  double *_y_err;
  //  double *_dydt_in;
  //  double *_dydt_out;
  double h;
  double _t;
} LINEAR;

//#define dim_linear 4
//int linearfunc (double t, const double y[], double f[], void *params)
//{
//  LINEAR *c = (LINEAR *)params;
//  f[0] = y[1];
//  f[1] = (-(2*c->dr+c->C)*y[1] -c->ddX*cos(y[0]) -c->g*sin(y[0]))/c->r;
//  f[2] = y[3];
//  f[3] = (c->F+c->T*sin(y[0]))/c->M;
//  return GSL_SUCCESS;
//}
////////////////
#define square(x) ((x)*(x))
LINEAR linear;
int PLANT_NO=1;
int initialize()
{
  if(PLANT_NO==1){
    linear.kx=1;
    linear.ku=1;
    linear.k=linear.kx+linear.ku;
    linear.a=(double*)malloc(sizeof(double)*linear.k);
    linear.b=&linear.a[linear.kx];
    linear.x=(double*)malloc(sizeof(double)*linear.k+1);
    linear.u=&linear.x[linear.kx];

    linear.a[0]=0.9;
    linear.b[0]=1;
  }
  else if(PLANT_NO==2){
    linear.kx=3;
    linear.ku=1;
    linear.k=linear.kx+linear.ku;
    linear.a=(double*)malloc(sizeof(double)*linear.k);
    linear.b=&linear.a[linear.kx];
    linear.x=(double*)malloc(sizeof(double)*linear.k+1);
    linear.u=&linear.x[linear.kx];

    linear.a[0]=1;
    linear.a[1]=1;
    linear.a[2]=1;
    linear.b[0]=1;
  }
  linear.umax=6.0;//check//  linear.Fmax=30;//check
  int i;for(i=0;i<linear.kx;i++) linear.x[i]=0;
#ifdef LINEARSUB
  //  linear.Fmax=_linear_Fmax;//check//  linear.Fmax=30;//check
  //  linear.xmax=20.0;//check//  linear.Fmax=30;//check
  if(_AP_umax>0) AP_u_max=linear.umax=_AP_umax;
  else AP_u_max=linear.umax;

  AP_u_min=-AP_u_max;
  rr=AP_r=_AP_r;//10
  rr_kyoyou=_rr_kyoyou;
  C_MODE=11;
#else
  linear.h=0.01;
#endif
  linear._t=0;
  //  int i;for(i=0;i<linear._dim;i++) linear._dydt_in[i]=0;

  return(0);
}
#ifndef LINEARSUB
char *fn="linear2io.dat";
FILE *fp;
#endif//#ifndef LINEARSUB

double plant(double uu)
{
  //store previous x and u
  int i;
  for(i=0;i<linear.kx-1;i++) linear.x[i+1]=linear.x[i];
  for(i=0;i<linear.ku-1;i++) linear.u[i+1]=linear.u[i];
  //store former output
  linear.x[0]=linear.y;

  //current input
  if(uu>linear.umax) linear.u[0]=linear.umax;
  else if(uu<-linear.umax) linear.u[0]=-linear.umax;
  else linear.u[0]=uu;
  // current output
  linear.y=0;
  for(i=0;i<linear.k;i++) linear.y+=linear.a[i]*linear.x[i];
  linear._t++;

  //  fprintf(stdout,"y%g=(%g)*(%g)+(%g)*(%g) %g %g %g\n",linear.y,linear.a[0],linear.x[0],linear.a[1],linear.x[1],uu,linear.u[0],linear.x[1]);

#ifndef LINEARSUB
  fprintf(fp,"%g %g %g %g", linear._t,linear.y,linear.u[0],linear.x[0]);//linear.a,linear.da
  //  fprintf(fp,"%.7e %.7e %.7e %.7e", linear._t,linear.y,linear.x[0],linear.u[0]);//linear.a,linear.da
  fprintf(fp,"\n");
#endif
  return(linear.y);
  //  return(linear.X);
}

#ifndef LINEARSUB
int main()
{
  ////////////////////////////////////////////////////
  /// method 1                                     ///
  /// input ddX                                    ///
  /// output a,x,y,F                               ///
  ////////////////////////////////////////////////////
  //  double h = 0.001;//,hh=0.1; h = 0.0001;//,hh=0.1;
  double t0=5;     //stationary
  double t1=5.0+t0;//accelerate(speed up)
  double t2=2.0+t1;//free run
  double t3=5.0+t2;//decelerate(slow down)
  double t4=20+t3;   //free run

  initialize();
  linear.h=1;
  int M=(t4/linear.h)+1;
  double *uu=(double*)malloc(sizeof(double)*M);
  double uumax=20;//

  int n;
  double t;
  ////////////////////////////////////////////////////
  ///   set ddX                                    ///
  ////////////////////////////////////////////////////
  int n0 =t0/linear.h+0.5, n4=t4/linear.h+0.5;
  for(t=0;t<t4;t+=linear.h){
    n=t/linear.h;
    if(t<=t0) uu[n]=0; //zero
    else if(t<=t1) uu[n]=uumax/(t1-t0); //accelerate
    else if(t<=t2) uu[n]=0;            //free move
    else if(t<=t3) uu[n]=-uumax/(t3-t2);//decelerate
    else uu[n]=0;                      //free move
    //    fprintf(stdout,"uu[%d]=%e\n",n,uu[n]);
  }
  ////////////////////////////////////////////////////
  ///   Solve by the Runge Kutta Method of GSL     ///
  ////////////////////////////////////////////////////
  fp=fopen(fn,"w");
  fprintf(fp,"#%.7e %d %d %.7e %.7e %.7e %.7e #h,n0,n4,a[0]-a[3]\n",linear.h,n0,n4,linear.a[0],linear.a[1],linear.a[2],linear.a[3]);
  fprintf(fp,"#t,y,x,u\n");
  for(t=0;t<t4;t+=linear.h){
    n=t/linear.h;
    plant(uu[n]);
  }
  fclose(fp);
  fprintf(stdout,"Results are stored in '%s'.\n",fn);

  return 0;
}

#endif //#ifndef SUB


text:
#include <math.h>
#include <stdio.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_cblas.h>

#include "motionControl.h"

void printMatrix(gsl_matrix* m, int row, int col) {
    for(int i = 0; i < row; i++) {
        for(int j = 0; j < col; j++) {
            printf("%g ", gsl_matrix_get(m, i, j));
        }
        printf("\n");
    }
}

gsl_matrix* measurementMatrix() {
    gsl_matrix *m = gsl_matrix_alloc (3, 3);

    gsl_matrix_set (m, 0, 0, -.5);
    gsl_matrix_set (m, 0, 1, .866);
    gsl_matrix_set (m, 0, 2, (.866*6.9282 - -.5*4));
    gsl_matrix_set (m, 1, 0, -.5);
    gsl_matrix_set (m, 1, 1, -.866);
    gsl_matrix_set (m, 1, 2, (.866*6.9282 - -.5*4));
    gsl_matrix_set (m, 2, 0, 1);
    gsl_matrix_set (m, 2, 1, 0);
    gsl_matrix_set (m, 2, 2, 8);

    return m;
}

gsl_matrix* rotationMatrix(float omegaWorld, float omegaBody) {
    gsl_matrix *m = gsl_matrix_alloc (3, 3);

    float theta = omegaWorld - omegaBody; // TODO check if radians or degrees

    gsl_matrix_set (m, 0, 0, cos(theta));
    gsl_matrix_set (m, 0, 1, sin(theta));
    gsl_matrix_set (m, 0, 2, 0);
    gsl_matrix_set (m, 1, 0, -sin(theta));
    gsl_matrix_set (m, 1, 1, cos(theta));
    gsl_matrix_set (m, 1, 2, 0);
    gsl_matrix_set (m, 2, 0, 0);
    gsl_matrix_set (m, 2, 1, 0);
    gsl_matrix_set (m, 2, 2, 1);

    return m;
}

gsl_matrix* pointToVelocity(float cX, float cY, float cOmega, float dX, float dY, float dOmega, float time) { //d stands for desired
    gsl_matrix *m = gsl_matrix_alloc (3, 1);

    float vX = (dX - cX) / time;
    float vY = (dY - cY) / time;
    float omega = (dOmega - cOmega) / time;

    gsl_matrix_set (m, 0, 0, vX);
    gsl_matrix_set (m, 1, 0, vY);
    gsl_matrix_set (m, 2, 0, omega);

    return m;
}

gsl_matrix* linearToAngular(gsl_matrix* R, gsl_matrix* M, gsl_matrix* V) {
    gsl_matrix_transpose(R);

    int signum;

    gsl_matrix * inverseR = gsl_matrix_alloc (3, 3);
    gsl_permutation * perm = gsl_permutation_alloc (3);

    gsl_linalg_LU_decomp (M, perm, &signum);
    gsl_linalg_LU_invert (M, perm, inverseR);

    gsl_matrix_mul_elements (inverseR, M);
    gsl_matrix* result = inverseR;

    float omega1 = rowSum(result, V, 0);
    float omega2 = rowSum(result, V, 1);
    float omega3 = rowSum(result, V, 2);

    gsl_matrix * final = gsl_matrix_alloc (3, 1);
    gsl_matrix_set (final, 0, 0, omega1);
    gsl_matrix_set (final, 1, 0, omega2);
    gsl_matrix_set (final, 2, 0, omega3);

    gsl_matrix_free(inverseR);
    gsl_permutation_free(perm);
    return final;
}

float rowSum(gsl_matrix* result, gsl_matrix* V, int row) {
    float sum = 0;
    for(int i = 0; i < 3; i++) {
        sum += gsl_matrix_get(result, row, i) * gsl_matrix_get(V, i, 0);
    }
    return sum;
}

int main (void) {
    gsl_matrix* M = measurementMatrix();
    gsl_matrix* R = rotationMatrix(1, 2);
    gsl_matrix* V = pointToVelocity(0, 0, 0, 1, 1, 0, 1);

    gsl_matrix* final = linearToAngular(R, M, V);

    printMatrix(final, 3, 1);
    gsl_matrix_free (M);
    gsl_matrix_free (R);
    gsl_matrix_free (V);
    gsl_matrix_free (final);

    return 0;
}


text:
/*
 * Copyright (c) 2018-2021 Aleksas Mazeliauskas, Stefan Floerchinger,
 *                    Eduardo Grossi, and Derek Teaney
 * All rights reserved.
 *
 * FastReso is distributed under MIT license;
 * see the LICENSE file that should be present in the root
 * of the source distribution, or alternately available at:
 * https://github.com/amazeliauskas/FastReso/
 */
#ifndef FASTFO_TParticle_h
#define FASTFO_TParticle_h

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

#include <vector>
#include <string>
#include <iostream>
#include <grid_params.h>
#include <gsl/gsl_interp.h>
#include <gsl/gsl_spline.h>
//! Enum class for deciding on Bose-Einstein, Fermi-Dirac or Botlzmann particle dsitribution
enum class EParticleType : int { kBoson = 1, kFermion = -1, kBoltzman = 0 };
enum class EFjIndex : int {
  kFeq1 = 0,
  kFeq2 = 1,
  kFshear1 = 2,
  kFshear2 = 3,
  kFshear3 = 4,
  kFbulk1  = 5,
  kFbulk2  = 6,
  kFtemp1  = 7,
  kFtemp2  = 8,
  kFvel1  = 9,
  kFvel2  = 10,
  kFvel3  = 11 };

//! Particle class with basic properties and universal decay spectra components
class TParticle {
  protected:
    //! Particle properties
    std::string fParticleName;  // particle name
    std::string fDescription;  // particle name
    std::string fDescriptionHeader;  // particle name
    std::string fHeader = "1:pbar [GeV]\t2:m [GeV]";
    std::vector<std::string> fComponentNames = {"Feq 1", "Feq 2","Fshear 1","Fshear 2","Fshear 3","Fbulk 1","Fbulk 2","Ftemp 1","Ftemp 2","Fvel 1","Fvel 2","Fvel 3"};
    //std::string fHeader = "1:pbar [GeV]\t2:m [GeV]\t3:feq 1\t4:feq 2\t5:fshear 1\t6:fshear 2\t7:fshear 3\t8:fbulk 1\t9:fbulk 2\t10:ftemperature 1\t11:ftemperature 2\t12:fvelocity 1\t13:fvelocity 2\t14:fvelocity 3";

        // particle name
    double fMass;   // [GeV] resonance mass
    double fGamma;  // [GeV] resonance width (currently not used)
    double fIsospin; // isospin J |J,m>
    int fPDGCode;   // number code according to PDG
    int fNu;        // degeneracy of a particle
    double fQB;     // baryon charge
    double fQS;     // strange charge
    double fQC;     // charme charge
    EParticleType fParticleType; // particle statistic (bose-einstein, fermi-dirac, boltzmann)

    double fTfo;
    double fNtotal; //integrated particle number
    double fNtotal_buffer;

    //! Grid and interpolator for decay spectra components
    const gsl_interp_type   *fSpline_type = gsl_interp_cspline;
    double fPbar_arr[grid_params::fNpbar]; // array of fluid restframe momentum
    double fFj_arr[grid_params::fNf][grid_params::fNpbar]; // array of scalar functions Fj
    double fFj_arr_buffer[grid_params::fNf][grid_params::fNpbar]; // buffer array of scalar functions Fj
    gsl_spline *fSpline_Fj[grid_params::fNf] ;
    gsl_interp_accel *fPbar_acc;
    //! Bool variables to control data read in/out-put
    bool fIsLocked = false; //if true, prevent modifying fFj_arr data
    bool fIsModified[grid_params::fNf]; //if true, require initilization of interpolator
    //! After any update of the grid, don't forget to reinitialized the interpolator!
    //! initiliazes the interpolator
    void init_grid(std::string comps="Feq Fshear Fbulk Ftemp Fvel");
    void init(int j);
  public:

    std::vector <EFjIndex> fComponents;
    //! Contructor with particle properties
    ~TParticle() ;
    //! return parameters about the particle
    std::string getName() {return fParticleName;};
    double getM() {return fMass;};
    double getGamma() {return fGamma;};
    double getQB() {return fQB;};
    double getQS() {return fQS;};
    double getQC() {return fQC;};
    double getIsospin() {return fIsospin;};
    double getN() {return fNtotal;};
    int getNu() {return fNu;};
    int getPDG() {return fPDGCode;};
    EParticleType getType() {return fParticleType;};

    //! increment integrated particle number
    void addN(double n) {
            if (fIsLocked){
        std::cerr << "\033[1mTParticle.h\033[0m : \033[1;31merror\033[0m : attemped to modify locked "<< fParticleName <<" data" << std::endl;
        exit(EXIT_FAILURE);}
      else { fNtotal+=n; fNtotal_buffer+=n; }};
    //! set initial thermal temperature
    void setTfo(double Tfo) {fTfo=Tfo;};
    //! set lock on data modification
    void lock() {fIsLocked=true;};
    bool hasdecayed() {return fIsLocked;};
    //! clean the buffer
    void clean_buffer() {
      for(int i = 0; i <grid_params::fNpbar; i++){
        for (int jj=0; jj<grid_params::fNf; jj++){
          fFj_arr_buffer[jj][i]=0.0;
        }
        fNtotal_buffer=0.0; }
    };

   //! return the energy and momentum of grid point i
   double getEbar(int i) {return sqrt(fMass*fMass+getPbar(i)*getPbar(i));};
   double getPbar(int i) {return fPbar_arr[i];};
   //! Return the number of momentum grid points
    int getNpbar() {return grid_params::fNpbar;};
    //! Return of the initialization/freeze-out temperature
    double getTfo() {return fTfo;};

    //! increment Fj array at cite i by F
    void addFj(int j, int i, double F) {
      if (fIsLocked){
        std::cerr << "\033[1mTParticle.h\033[0m : \033[1;31merror\033[0m : attemped to modify locked "<< fParticleName <<" data" << std::endl;
        exit(EXIT_FAILURE);}
      else {
        fFj_arr[j][i]+=F;
        fFj_arr_buffer[j][i]+=F;
        fIsModified[j]=true;
      }
    };

    //! Returns the universal scalar functions Fj at Ebar
    double get_Fj(int j, double Ebar);

    //Routines to print the data to screen or file
    std::string getDescriptionHeader(){ return fDescriptionHeader;};
    std::string getDescription(){return fDescription;};
    void print();
    void print(std::string tag);
    void print_buffer();
    void print_buffer(std::string tag);
};

#endif

max_stars_repo_name: amazeliauskas/FastReso

max_issues_repo_name: amazeliauskas/FastReso

max_forks_repo_name: amazeliauskas/FastReso


text:
#pragma once

#include <glm/glm.hpp>
#include <gsl/gsl>
#include <optional>

namespace rev {

struct Segment {
    glm::vec3 start;
    glm::vec3 end;
};

struct Ray {
    glm::vec3 origin;
    glm::vec3 direction;
};

struct Sphere {
    glm::vec3 center;
    float radius;
};

template <typename VertexRange, typename SegmentVisitor>
void iteratePolygonVertices(const VertexRange& vertices, SegmentVisitor&& visitor)
{
    Expects(vertices.size() >= 3);
    auto begin = std::begin(vertices);
    auto end = std::end(vertices);

    auto segmentStart = begin;
    auto segmentEnd = begin + 1;

    while (segmentStart != end) {
        if (segmentEnd == end) {
            segmentEnd = begin;
        }

        visitor(Segment{ *segmentStart, *segmentEnd });
        segmentStart++;
        segmentEnd++;
    }
}

class NullPolygonBuilder {
public:
    void addVertex(const glm::vec3& vertex) {}
};

template <typename IteratorType>
class OutputIteratorPolygonBuilder {
public:
    OutputIteratorPolygonBuilder(IteratorType iterator)
        : _iter(iterator)
    {
    }

    void addVertex(const glm::vec3& vertex)
    {
        *_iter = vertex;
        _iter++;
    }

    IteratorType getIterator() const { return _iter; }

private:
    IteratorType _iter;
};

struct AxisAlignedPlane {
    struct Hit {
        glm::vec3 intersectionPoint;
        float t;
    };

    Hit intersect(const Ray& ray) const
    {
        float boundaryDelta = boundary - ray.origin[dimensionIndex];
        float t = boundaryDelta / ray.direction[dimensionIndex];

        glm::vec3 intersectPoint = ray.origin + (t * ray.direction);

        // Avoid floating point error along the boundary dimension.
        intersectPoint[dimensionIndex] = boundary;
        return { intersectPoint, t };
    }

    glm::vec3 intersect(const Segment& segment) const
    {
        glm::vec3 direction = segment.end - segment.start;
        auto [intersectionPoint, t] = intersect(Ray{ segment.start, direction });
        if (!(t > 0.0f)) {
            return segment.start;
        } else if (!(t < 1.0f)) {
            return segment.end;
        }

        for (size_t k = 0; k < 3; k++) {
            if (segment.start[k] == segment.end[k]) {
                intersectionPoint[k] = segment.start[k];
            }
        }
        return intersectionPoint;
    }

    template <typename InputVertexRange, typename LeftBuilder, typename RightBuilder>
    void splitConvexPolygon(const InputVertexRange& inputVertices, LeftBuilder&& leftBuilder,
        RightBuilder&& rightBuilder) const
    {
        iteratePolygonVertices(
            inputVertices, [this, &leftBuilder, &rightBuilder](const Segment& segment) {
                float startPosition = segment.start[dimensionIndex];
                float endPosition = segment.end[dimensionIndex];

                bool emitIntersection = false;
                if (startPosition < boundary) {
                    leftBuilder.addVertex(segment.start);
                    if (endPosition > boundary) {
                        emitIntersection = true;
                    }
                } else if (startPosition > boundary) {
                    rightBuilder.addVertex(segment.start);
                    if (endPosition < boundary) {
                        emitIntersection = true;
                    }
                } else {
                    // Right on the boundary.
                    leftBuilder.addVertex(segment.start);
                    rightBuilder.addVertex(segment.start);
                }

                if (emitIntersection) {
                    auto intersection = intersect(segment);
                    leftBuilder.addVertex(intersection);
                    rightBuilder.addVertex(intersection);
                }
            });
    }

    uint8_t dimensionIndex;
    float boundary;
};

struct AxisAlignedBoundingBox {
    void expandToVertex(const glm::vec3& vertex)
    {
        for (int k = 0; k < 3; k++) {
            if (vertex[k] < minimum[k]) {
                minimum[k] = vertex[k];
            }
            if (vertex[k] > maximum[k]) {
                maximum[k] = vertex[k];
            }
        }
    }

    void expandToBox(const AxisAlignedBoundingBox& other)
    {
        for (int k = 0; k < 3; k++) {
            if (other.minimum[k] < minimum[k]) {
                minimum[k] = other.minimum[k];
            }
            if (other.maximum[k] > maximum[k]) {
                maximum[k] = other.maximum[k];
            }
        }
    }

    void reduceToBox(const AxisAlignedBoundingBox& other)
    {
        for (int k = 0; k < 3; k++) {
            if (other.minimum[k] > minimum[k]) {
                minimum[k] = other.minimum[k];
            }
            if (other.maximum[k] < maximum[k]) {
                maximum[k] = other.maximum[k];
            }
            Expects(!(minimum[k] > maximum[k]));
        }
    }

    float getSurfaceArea() const
    {
        glm::vec3 diagonal = maximum - minimum;
        float surfaceArea = 0.0f;
        surfaceArea += diagonal.x * diagonal.y;
        surfaceArea += diagonal.y * diagonal.z;
        surfaceArea += diagonal.z * diagonal.z;

        return 2.0f * surfaceArea;
    }

    float getVolume() const
    {
        glm::vec3 diagonal = maximum - minimum;
        return diagonal.x * diagonal.y * diagonal.z;
    }

    std::array<AxisAlignedBoundingBox, 2> split(const AxisAlignedPlane& plane) const
    {
        Expects(!(plane.boundary < minimum[plane.dimensionIndex]));
        Expects(!(plane.boundary > maximum[plane.dimensionIndex]));

        std::array<AxisAlignedBoundingBox, 2> splitBoxes{ *this, *this };
        splitBoxes[0].maximum[plane.dimensionIndex] = plane.boundary;
        splitBoxes[1].minimum[plane.dimensionIndex] = plane.boundary;

        return splitBoxes;
    }

    struct Hit {
        glm::vec3 intersectionPoint;
        float t;
        uint8_t planeDimension;
    };

    bool containsPoint(const glm::vec3& point) const
    {
        for (int k = 0; k < 3; k++) {
            if (point[k] < minimum[k]) {
                return false;
            }
            if (point[k] > maximum[k]) {
                return false;
            }
        }

        return true;
    }

    bool intersectsSphere(const Sphere& sphere) const
    {
        float r2 = sphere.radius * sphere.radius;
        for (size_t k = 0; k < 3; k++) {
            if (sphere.center[k] < minimum[k]) {
                float diff = minimum[k] - sphere.center[k];
                r2 -= diff * diff;
            } else if (sphere.center[k] > maximum[k]) {
                float diff = sphere.center[k] - maximum[k];
                r2 -= diff * diff;
            }
        }
        return r2 > 0.0f;
    }

    std::optional<Hit> castExternalRay(const Ray& ray) const
    {
        for (uint8_t k = 0; k < 3; k++) {
            float d = ray.direction[k];
            if (d < 0.0f) {
                float boundary = maximum[k];
                if (ray.origin[k] > boundary) {
                    auto intersection = AxisAlignedPlane{ k, boundary }.intersect(ray);
                    if (containsPoint(intersection.intersectionPoint)) {
                        return Hit{ intersection.intersectionPoint, intersection.t, k };
                    }
                }
            } else if (d > 0.0f) {
                float boundary = minimum[k];
                if (ray.origin[k] < boundary) {
                    auto intersection = AxisAlignedPlane{ k, boundary }.intersect(ray);
                    if (containsPoint(intersection.intersectionPoint)) {
                        return Hit{ intersection.intersectionPoint, intersection.t, k };
                    }
                }
            }
        }
        return std::nullopt;
    }

    Hit castInternalRay(const Ray& ray) const
    {
        for (uint8_t k = 0; k < 3; k++) {
            float d = ray.direction[k];

            float boundary;
            if (d < 0.0f) {
                boundary = minimum[k];
            } else if (d > 0.0f) {
                boundary = maximum[k];
            } else {
                continue;
            }

            auto intersection = AxisAlignedPlane{ k, boundary }.intersect(ray);
            if (containsPoint(intersection.intersectionPoint)) {
                return Hit{ intersection.intersectionPoint, intersection.t, k };
            }
        }
        throw std::runtime_error("Internal ray intersection not found.");
    }

    glm::vec3 minimum{ std::numeric_limits<float>::infinity() };
    glm::vec3 maximum{ -std::numeric_limits<float>::infinity() };
};

template <typename VertexRange>
AxisAlignedBoundingBox smallestBoxContainingVertices(const VertexRange& range)
{
    AxisAlignedBoundingBox box;
    for (const auto& vertex : range) {
        box.expandToVertex(vertex);
    }
    return box;
}

}

max_stars_repo_name: eyebrowsoffire/rev

max_issues_repo_name: eyebrowsoffire/rev

max_forks_repo_name: eyebrowsoffire/rev


text:
#pragma once

#include <unordered_map>
#include <vector>

#include <gsl/gsl>

#include "huffman.h"
#include "interface.h"
#include "minqueue.h"

class statistics {
    public:
    size_t original;
    size_t modified;
};

class huffman_encoder {
    // Statistics of the current encoding.
    statistics stats;

    std::unordered_map<char, std::vector<int>> encoding;
    std::unordered_map<char, int> count_chars(std::vector<char> buf);

    void evaluate(minqueue_node* root, std::vector<int> encoding);

    public:
    huffman_encoder();

    std::vector<int> encode(std::vector<char> buf);

    void display_encoding();
    void display_stats();
};

max_stars_repo_name: joshleeb/Huffman

max_issues_repo_name: joshleeb/Huffman

max_forks_repo_name: joshleeb/Huffman


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

#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_roots.h>

#include "demof.h"
#include "demof.c"

int
main ()
{
  int status;
  int iterations = 0, max_iterations = 100;
  gsl_root_fsolver *s;
  double r = 0, r_expected = sqrt (5.0);
  double x_lower x = 0.0, x_upper = 5.0;
  gsl_function F;
  struct quadratic_params params =
  {1.0, 0.0, -5.0};

  F.function = &quadratic;
  F.params = &params;

  s = gsl_root_fsolver_alloc (gsl_root_fsolver_bisection);
  gsl_root_fsolver_set (s, &F, x);

  printf ("using %s method\n", gsl_root_fsolver_name (s));

  printf ("%5s [%9s, %9s] %9s %9s %10s %9s\n",
          "iter", "lower", "upper", "root", "actual", "err", "err(est)");

  do
    {
      iterations++;
      status = gsl_root_fsolver_iterate (s);
      r = gsl_root_fsolver_root (s);
      x_lower = gsl_root_fsolver_x_lower (s);
      x_upper = gsl_root_fsolver_x_upper (s);
      status = gsl_root_test_interval (x, 0, 0.001);

      if (status == GSL_SUCCESS)
        printf ("Converged:\n");

      printf ("%5d [%.7f, %.7f] %.7f %.7f %+.7f %.7f\n",
              iterations, x_lower, x_upper,
              r, r_expected, r - r_expected, x_upper - x_lower);
    }
  while (status == GSL_CONTINUE && iterations < max_iterations);

}

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim


text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_rng.h>

#define N 65536 //the final number of nodes
#define I 2 // the initial number of nodes

#define MI 0 //the initial trial number
#define MF 33 //(the final trial number) - 1

#define RS 35382 //seed of the random number generator

int main()
{
	FILE *fp;

	double delta=sqrt(1.5);
	int tout[]={512,1024,2048,4096,8192,16384,32768,65536,131072};

	int i,j,k,l,m,t,tt;
	int tmp,flag;

	int **nb; //neighbors of each node

	int x; //new node
	int vdx; //virtual degree of the node x
	int *tmparr;

	int z;
	int nzsize;
	int *nz;

	int *deg; //degree of each node

	double degav; //average degree
	double *knn; //average degree of the nearest neighbors of each node
	double *cc; //local clustering coefficient of each node
	double ccav; //average local clustering coefficient

	char filename[100];

	const gsl_rng_type * TYPE;
	gsl_rng * ran;

	gsl_rng_env_setup();
	TYPE=gsl_rng_default;
	ran=gsl_rng_alloc(TYPE);
	gsl_rng_set(ran,RS);

	tmparr=malloc(sizeof(int)*N);
	nz=malloc(sizeof(int)*N);

	deg=malloc(sizeof(int)*N);
	knn=malloc(sizeof(double)*N);
	cc=malloc(sizeof(double)*N);

	nb=malloc(sizeof(int*)*N);
	for (i=0;i<N;i++){
		nb[i]=malloc(sizeof(int)*N);
	}

	for (m=MI;m<MF;m++){

		//generation of the initial network
		for (i=0;i<I;i++){
			deg[i]=0;
			for (j=0;j<I;j++){
				if (j!=i){
					nb[i][deg[i]]=j;
					deg[i]++;
				}
			}
		}

		tt=0;
		for (t=I;t<N;t++){
			//adding (t+1)-th node: copying degree step
			x=(int)(gsl_rng_uniform(ran)*t);
			vdx=1+(int)(gsl_rng_uniform(ran)*ceil(delta*deg[x]));
			deg[t]=0;

			/*//new OPA pre-adjunction rule
			j=0;
			for (i=0;i<t;i++){
				if (vdx<=ceil(delta*deg[i])){
					tmparr[j]=i;
					j++;
				}
			}
			z=tmparr[(int)(gsl_rng_uniform(ran)*j)];
			nzsize=deg[z]+1;
			for (i=0;i<nzsize-1;i++){
				nz[i]=nb[z][i];
			}
			nz[nzsize-1]=z;
			*/
			/*//OPA pre-adjunction rule
			nzsize=0;
			for (i=0;i<t;i++){
				if (vdx<=ceil(delta*deg[i])){
					nz[nzsize]=i;
					nzsize++;
				}
			}
			*/

			//OPA adjunction rule
			tmp=t;
			for (i=0;i<t;i++){
				if (vdx<=ceil(delta*deg[i]) && deg[i]<tmp)tmp=deg[i];
			}
			nzsize=0;
			for (i=0;i<t;i++){
				if (deg[i]==tmp){
					nz[nzsize]=i;
					nzsize++;
				}
			}

			//choosing the targets of the new node x
			if (vdx>=nzsize){
				for (i=0;i<nzsize;i++){
					nb[nz[i]][deg[nz[i]]]=t; deg[nz[i]]++;
					nb[t][deg[t]]=nz[i]; deg[t]++;
				}
			} else{
				for (j=0;j<vdx;j++){
					i=j+(int)(gsl_rng_uniform(ran)*(nzsize-j));

					nb[nz[i]][deg[nz[i]]]=t; deg[nz[i]]++;
					nb[t][deg[t]]=nz[i]; deg[t]++;

					tmp=nz[j];
					nz[j]=nz[i];
					nz[i]=tmp;
				}
			}

			//output
			if (t+1==tout[tt]){
				//degav
				degav=0.0;
				for (i=0;i<t+1;i++){
					degav+=deg[i];
				}
				degav/=(t+1);

				//knn
				for (i=0;i<t+1;i++){
					knn[i]=0.0;
					for (j=0;j<deg[i];j++){
						knn[i]+=deg[nb[i][j]];
					}
					if (deg[i]>0)knn[i]/=deg[i];
				}

				//cc
				ccav=0.0;
				for (i=0;i<t+1;i++){
					tmp=0;
					for (j=0;j<deg[i];j++){
						for (k=j+1;k<deg[i];k++){
							flag=0;
							for (l=0;l<deg[nb[i][j]];l++){
								if (nb[nb[i][j]][l]==nb[i][k]){
									flag=1;
									break;
								}
							}
							if (flag==1)tmp++;
						}
					}

					if (deg[i]>1)cc[i]=2.0*tmp/(deg[i]*(deg[i]-1));
					else cc[i]=0.0;

					ccav+=cc[i];
				}
				ccav/=(t+1);

				sprintf(filename,"id-deg-knn-cc_t%d_trial%d.txt",t+1,m);
				fp=fopen(filename,"w");
				for (i=0;i<t+1;i++){
					fprintf(fp,"%d %d %.15f %.15f\n",i,deg[i],knn[i],cc[i]);
				}
				fclose(fp);

				printf("m=%d t=%d degav=%f ccav=%f\n",m,t+1,degav,ccav);

				tt++;
			}//end of output

		}//end of t-loop

	}//end of m-loop

	free(tmparr);
	free(nz);
	free(deg);
	free(knn);
	free(cc);
	for (i=0;i<N;i++){
		free(nb[i]);
	}
	free(nb);

	return 0;
}


text:
#ifndef CLASS_H
#define CLASS_H

#include <iostream>
#include <ctime>
#include <fstream>
#include <string>
#include <algorithm>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_cdf.h>
//multivariate normal libraries
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>

using namespace std;
int poisson(double lambda);
double expTrunc(double, double);

class Parameters { //store all parameters in one class to stop them being everywhere. Will also make it easier to change parameters when needed.
  gsl_rng * rando;
  const gsl_rng_type * T;
  double lbda_original,v_to_h_original,sig_original; //original parameters stored so when changed from intervention still have in memory.
  public:
    Parameters(); //imports parameters from text file.
    Parameters(int); //imports parameters from text file and use index to seed random number generator.
    Parameters(double,double); //bite exposure gamma a, bite exposure gamma b
    Parameters(double,double,double); //bite exposure gamma a, bite exposure gamma b, v_to_h ratio.
    ~Parameters(void);
    //functions
    void setParameters(void);
    void importParameters(void);
    void setBR(bool intervention);
    void setVH(bool intervention);
    void setMu(bool intervention);
    double L3Uptake(double);
    double readParamLine(ifstream& myfile);
    double gamma_dist(double, double);
    double normal_dist(double, double);
    double uniform_dist(void);
    double expTrunc(double, double);
    int poisson_dist(double);
    int rmvnorm(const gsl_rng *, const int, const gsl_vector *, const gsl_matrix *, gsl_vector *);
    int setUB(double *, double *, double *);
    double immuneRK4Step(double,double);
    //parameters
    double a;
    double b;
    double shapeRisk; //shape of probability distribution for risk.
    double riskMu1; //mean of risk for age group less than 16
    double riskMu2; //mean of risk for age group  17-29
    double riskMu3; //mean of risk for age group  30+
    double mu; //death rate of worms
    double theta; //immune system response parameter. 0.112
    double gamma; //mf death rate
    double alpha; //mf birth rate
    double v_to_h; //vector to human ratio.
    double xi; //constant bite rate (5760.0 / 30.0)
    double kappas1; //vector uptake and development
    double r1; //vector uptake and development
    double tau;
    double z; //waning immunity
    double nu; //poly-monogamy parameter
    //shared state parameters
    double L3; //larvae density.
    double g; //Proportion of mosquitoes which pick up infection when biting an infected host
    double sig; //death rate of mosquitos
    double psi1; //Proportion of L3 leaving mosquito per bite
    double psi2; //Proportion of L3 leaving mosquito that enter host
    double lbda; //bite rate per mosquito per month
    double dt;
    double lbdaR, v_to_hR, sigR; //reduction in bite-rate and v_to_h ratio.
    int nMDA; //number of rounds of MDA.
    int mdaFreq;
    double covMDA;
    double s2;
    double mfPropMDA; //proportion of mf killed due to MDA
    double wPropMDA; //proportion of worms killed due to MDA
    double sysComp,u0Comp, sigComp; //parameters for systematic non-compliance
    double sysCompN, u0CompN, sigCompN; //parameters for systematic non-compliance of bed-nets.
    double covN; //coverage of bednets.
    double rhoCN; //correlation of non-compliance between bed-nets and MDA (Chemotherapy).
    int mosquitoSpecies; //species of mosquito 0 - Anopheles, 1 - Culex
    double rhoBU; //correlation between non-compliance and bite risk.
    int aWol; //using doxycycline as intervention 0- not used. 1- is used.
    double sN,dN; //probability of mosquito successfully biting given use of bed nets. prob of death due to bed net.
    int IDAControl; //switch to IDA after five rounds of normal MDA treatment. 0- not used, 1 - is used.
} ;

class Host {

    //double dt; //time-step for tau.
    void initialise(void); //initialise state variables.

  public:

    Host(double);
    Host(void);
    void react(void);
    void evolve(double);
    double mfConc(void);
    double biteRate(void);
    double repRate(void);
    void updateRisk(bool);
    //state variables
    int W,WM,WF; //number of worms, number of male worms, number of female worms
    double I; //immune response (assumed to be deterministic)
    double M; //mf produced. //int
    double t; //time of host
    double a; //age of host.
    int aWol; //host is treated with Doxycycline.
    int bedNet; //host used bednet.
    //host parameters
    Parameters *param; //pointer to params.
    double b; //risk of infection
    double pRisk; //p-value for risk of infection (updates b)
    double uComp; //random parameter used to calculate compliance with MDA
    double uCompN; //random parameter used to calculate compliance with bed nets.
    //const double mu = 0.1/(30.0); //death rate of worms
    //double beta = 0.112; //immune system response parameter.
    //double gamma = 0.1/(30.0);
    //double alpha = 0.2/(30.0);
    int births;
    int deaths;
    //double xi = (5760.0 / 30.0) * 0.414 * 0.32 * 0.2;
    //double tau = 0.04 / (365.0); //0.02
    //double L3 = 1.0; //keeps track of L3 in total population.
} ;

class Model {

    Host *host_pop;
    int n;
    bool aWolInt; //undergoing aWol intervention.
    int bedNetInt; //undergoing bednet intervention.
  public:
    Model(int,double,double);
    Model(int, Parameters*);
    Parameters *param;
    void react(void);
    void evolve(double);
    void save(string);
    void outputs(void);
    void saveOngoing(string); //appends WM,WF,M to the bottom of filename
    void evolveAndSaves(double tot_t, string filename); //functionally same as evolve, but saves outputs at regular intervals.
    double L3(void);
    void MDAEvent(void);
    void aWolEvent(bool);
    void bedNetEvent(void);

} ;

#endif

max_stars_repo_name: sempwn/LF-model

max_issues_repo_name: sempwn/LF-model

max_forks_repo_name: sempwn/LF-model


text:
#include "openblas_utest.h"
#include <stdio.h>
#include <stdlib.h>
#include <cblas.h>

#define LAPACK_ROW_MAJOR               101
blasint LAPACKE_dgesvd( blasint matrix_layout, char jobu, char jobvt,
                           blasint m, blasint n, double* a,
                           blasint lda, double* s, double* u, blasint ldu,
                           double* vt, blasint ldvt, double* superb );

#define DATASIZE 100

double s[DATASIZE];
double u[DATASIZE*DATASIZE];
double vt[DATASIZE*DATASIZE];
double X[DATASIZE*DATASIZE];
double superb[DATASIZE];
double tmp[DATASIZE*DATASIZE];
double m[DATASIZE*DATASIZE];

CTEST(kernel_regress,skx_avx)
{
#ifdef BUILD_DOUBLE
    double norm;
    int i, j, info;
    srand(0);
	for (i = 0; i < DATASIZE*DATASIZE; i++) {
        m[i] = (rand()+0.0)/RAND_MAX * 10;
        tmp[i] = m[i];
    }

    info = LAPACKE_dgesvd( LAPACK_ROW_MAJOR, 'A', 'A', DATASIZE, DATASIZE, m, DATASIZE,
                        s, u, DATASIZE, vt, DATASIZE, superb);

	for (i = 0; i < DATASIZE; i++) {
	    for (j = 0; j < DATASIZE; j++) {
            u[i*DATASIZE+j] = u[i*DATASIZE+j]*s[j];
        }
    }
    cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
                DATASIZE, DATASIZE, DATASIZE, 1, u, DATASIZE, vt, DATASIZE, 0, X, DATASIZE);

	for (i = 0; i < DATASIZE*DATASIZE; i++) {
        X[i] = X[i] - tmp[i];
    }

    norm = cblas_dnrm2(DATASIZE*DATASIZE, X, 1);
    ASSERT_DBL_NEAR_TOL(0.0, norm, 1e-10);
#endif
}

max_stars_repo_name: dnoan/OpenBLAS

max_issues_repo_name: dnoan/OpenBLAS

max_forks_repo_name: dnoan/OpenBLAS


text:
/* randist/gauss.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, 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 <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

/* Of the two methods provided below, I think the Polar method is more
 * efficient, but only when you are actually producing two random
 * deviates.  We don't produce two, because then we'd have to save one
 * in a static variable for the next call, and that would screws up
 * re-entrant or threaded code, so we only produce one.  This makes
 * the Ratio method suddenly more appealing.  There are further tests
 * one can make if the log() is slow.  See Knuth for details */

/* Both methods pass the statistical tests; but the polar method
 * seems to be a touch faster on my home Pentium, EVEN though we
 * are only using half of the available random deviates!
 */

/* Polar (Box-Mueller) method; See Knuth v2, 3rd ed, p122 */

double
gsl_ran_gaussian (const gsl_rng * r, const double sigma)
{
  double x, y, r2;

  do
    {
      /* choose x,y in uniform square (-1,-1) to (+1,+1) */

      x = -1 + 2 * gsl_rng_uniform (r);
      y = -1 + 2 * gsl_rng_uniform (r);

      /* see if it is in the unit circle */
      r2 = x * x + y * y;
    }
  while (r2 > 1.0 || r2 == 0);

  /* Box-Muller transform */
  return sigma * y * sqrt (-2.0 * log (r2) / r2);
}

/* Ratio method (Kinderman-Monahan); see Knuth v2, 3rd ed, p130 */
/* K+M, ACM Trans Math Software 3 (1977) 257-260. */

double
gsl_ran_gaussian_ratio_method (const gsl_rng * r, const double sigma)
{
  double u, v, x;

  do
    {
      v = gsl_rng_uniform (r);
      do
	{
	  u = gsl_rng_uniform (r);
	}
      while (u == 0);
      /* Const 1.715... = sqrt(8/e) */
      x = 1.71552776992141359295 * (v - 0.5) / u;
    }
  while (x * x > -4.0 * log (u));

  return sigma * x;
}

double
gsl_ran_gaussian_pdf (const double x, const double sigma)
{
  double u = x / fabs (sigma);
  double p = (1 / (sqrt (2 * M_PI) * fabs (sigma))) * exp (-u * u / 2);
  return p;
}

double
gsl_ran_ugaussian (const gsl_rng * r)
{
  return gsl_ran_gaussian (r, 1.0);
}

double
gsl_ran_ugaussian_ratio_method (const gsl_rng * r)
{
  return gsl_ran_gaussian_ratio_method (r, 1.0);
}

double
gsl_ran_ugaussian_pdf (const double x)
{
  return gsl_ran_gaussian_pdf (x, 1.0);
}