text:
/**
 * @file mcc_localizer.h
 * @brief grid search based localization
 * @author Kenichi Kumatani
 */

#ifndef MCC_LOCALIZER_H
#define MCC_LOCALIZER_H
#include <stdio.h>
#include <assert.h>
#include <float.h>

#include <gsl/gsl_block.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_fft_real.h>
#include <gsl/gsl_fft_complex.h>
#include <common/refcount.h>
#include "common/jexception.h"

#include "stream/stream.h"
#include "feature/feature.h"
//#include "modulated/modulated.h"

/**
   @brief  construct a search grid for source localization and return the time delay corresponding to each position on the grid.
   @note In order to do efficient grid-search for the source position which provides the maximum/minimum objective function, each cell size of a search grid has to be "reasonable".
   @usage
   After you construct an instance, you have to
   1. set the geometry of the microphone array with setRadius(), setDistanceBtwMicrophones() or setPositionsOfMicrophones().
   2. obtain the time delays at the source position with getTimeDelays(), and
   3. go to the next search candidate on the grid with nextSearchGrid().
 */
class SearchGridBuilder {
 public:
  SearchGridBuilder( int nChan, bool isFarField, unsigned int samplingFreq=16000);
  virtual ~SearchGridBuilder();

  const gsl_vector *getSearchPosition(){return(const gsl_vector *)_hypopos;}
  float maxTimeDelay(){return _maxTimeDelay;}
  size_t chanN(){return(_mpos->size1);}
  unsigned int samplingFrequency(){return(_samplingFreq);}
  void reset();

  virtual const gsl_vector *getTimeDelays(){return NULL;}
  virtual bool nextSearchGrid(){return false;}

protected:
  virtual const gsl_vector *nextSearchGridFF(){return NULL;}
  const gsl_vector *nextSearchGridNF();

protected:
  const bool _isFarField; /* if _isFarField = true, the far-filed is assumed. Otherwise, the near-field is assumed. */
  unsigned int _samplingFreq;
  float _maxTimeDelay;
  gsl_matrix* _mpos;
  gsl_vector *_hypopos;
  gsl_vector *_delays;
  float _constV;
};

typedef refcount_ptr<SearchGridBuilder> SearchGridBuilderPtr;

class SGB4LinearArray : public SearchGridBuilder{
public:
  SGB4LinearArray( int nChan, bool isFarField, unsigned int samplingFreq=16000);
  void setDistanceBtwMicrophones( float distance );
  void setPositionsOfMicrophones( const gsl_matrix* mpos );
  virtual const gsl_vector *getTimeDelays();
  virtual bool nextSearchGrid();

protected:
  const gsl_vector *nextSearchGridFF();
};

typedef Inherit<SGB4LinearArray, SearchGridBuilderPtr> SGB4LinearArrayPtr;

class SGB4CircularArray : public SearchGridBuilder {
public:
  SGB4CircularArray( int nChan, bool isFarField, unsigned int samplingFreq=16000 );
  void setRadius( float radius, float height=0.0 );
  virtual const gsl_vector *getTimeDelays();
  virtual bool nextSearchGrid();

protected:
  const gsl_vector *nextSearchGridFF();
};

typedef Inherit<SGB4CircularArray, SearchGridBuilderPtr> SGB4CircularArrayPtr;

/**
   @brief keeps samples in the previous block.
          this object is used in order to fill the gap between previous and current block data.
*/
class SampleHolder {
public:
  SampleHolder( size_t chanN, size_t maxSampleDelay ):
    _filledNum(0),_buffer(gsl_matrix_float_calloc(chanN,maxSampleDelay))
  {}

  ~SampleHolder()
  {
    gsl_matrix_float_free(_buffer);
  }

  void setSamples( size_t chanX, gsl_vector_float *samples )
  {
    size_t maxSampleDelay = _buffer->size2;

    if( samples->size >= maxSampleDelay ){
      size_t frameY = samples->size - maxSampleDelay;
      for(size_t frameX=0;frameX<maxSampleDelay;frameX++,frameY++){
	gsl_matrix_float_set( _buffer, chanX, frameX, gsl_vector_float_get( samples, frameY ) );
      }
      _filledNum = maxSampleDelay;
    }
    else{
      // shift the elements of the buffer
      size_t frameN = maxSampleDelay - samples->size;
      for(size_t frameX=0;frameX<frameN;frameX++){
	gsl_matrix_float_set( _buffer, chanX, frameX, gsl_matrix_float_get( _buffer, chanX, frameX + samples->size ) );
      }
      // hold the samples of the current block
      size_t frameY = 0;
      for(size_t frameX=frameN;frameX<maxSampleDelay;frameX++,frameY++){
	gsl_matrix_float_set( _buffer, chanX, frameX, gsl_vector_float_get( samples, frameY ) );
      }
      _filledNum += samples->size;
      if(  _filledNum > maxSampleDelay )
	_filledNum = maxSampleDelay;
    }
  }

  float getSample( int chanX, int minusFrameX )
  {
    return gsl_matrix_float_get( _buffer, chanX, _buffer->size2 + minusFrameX );
  }

  int nFilled(){ return _filledNum;}

private:
  int  _filledNum;
  gsl_matrix_float *_buffer; /* the number of channels X samples */
};

typedef refcount_ptr<SampleHolder> SampleHolderPtr;

/**
   @brief hold information for a position estimate.
 */
class SourceCandidate {
public:
  SourceCandidate(size_t chanN, double costV=100000):
    _costV(costV), _chanN(chanN){
    _sampledelay = new int[chanN];
    _position = gsl_vector_calloc( 3 );
    _eigenvalues = gsl_vector_alloc( chanN );
  }
  ~SourceCandidate(){
    delete [] _sampledelay;
    gsl_vector_free(_position);
    gsl_vector_free( _eigenvalues );
  }
  void setSourceInfo( int *tau, const gsl_vector *position, double costV, gsl_vector *eigenvalues ){
    for(size_t i=0;i<_chanN;i++){_sampledelay[i]=tau[i];}
    gsl_vector_memcpy(_position, (gsl_vector *)position);
    _costV = costV;
    gsl_vector_memcpy(_eigenvalues, eigenvalues);
  }
  void setConstV( double costV=100000 ){
    _costV = costV;
  }

  double _costV; /* value of the cost function */
  int *_sampledelay;
  gsl_vector *_position;
  gsl_vector *_eigenvalues;
private:
  size_t _chanN;
};

/**
   @class estimate the source positions which provide the larger multi-channel cross correlation values.
   @usage
   1. setChannel()
   2. next()
   3.

   @note The algorithm implemented here is described in:
          J. Chen, J. Benesty and Y. Huang, "Robust Time Delay Estimation Exploiting Redundancy Among Multiple Microphones", IEEE Trans. SAP, vol.11, Sep. 2003.
*/

class MCCLocalizer : public VectorFeatureStream {
public:
  MCCLocalizer( SearchGridBuilderPtr &sgbPtr, size_t maxSource=1, const String& nm= "MCCSourceLocalizer" );
  ~MCCLocalizer();
  virtual const gsl_vector* next(int frameX = -5);
  virtual void  reset();
  //void setChannel(SampleFeaturePtr& chan);
  void setChannel(VectorFloatFeatureStreamPtr& chan);
  /*@brief obtain a relative delay at a microphone which corresponds to the best candidate */
  int getDelayedSample( int chanX ){return _sourceCandidates[0]->_sampledelay[chanX];}
  /*@brief obtain the maximum MCCC value */
  double getMaxMCCC(){/* mistake: it was getMinMCCC() */
#ifdef _NO_LOG_
    return( 1 -_sourceCandidates[0]->_costV );
#else
    return( 1 -exp(_sourceCandidates[0]->_costV) );
#endif
  }
  /*@brief obtain the best candidate of the position estimates */
  const gsl_vector* getPosition(){return((const gsl_vector* )_sourceCandidates[0]->_position);}

  /*@brief obtain a relative delay at a microphone which corresponds to the N-th best candidate */
  int getNthBestDelayedSample( int nth, int chanX ){return _sourceCandidates[nth]->_sampledelay[chanX];}
  /*@brief obtain the N-th best MCCC value */
  double getNthBestMCCC( int nth ){
#ifdef _NO_LOG_
    return( 1 -_sourceCandidates[nth]->_costV );
#else
    return( 1 - exp(_sourceCandidates[nth]->_costV) );
#endif
  }
  /*@brief obtain theN-th  best candidate of the position estimates */
  const gsl_vector* getNthBestPosition( int nth ){return((const gsl_vector* )_sourceCandidates[nth]->_position);}

  const gsl_vector* getEigenValues(){return ((const gsl_vector*)_sourceCandidates[0]->_eigenvalues);}
  gsl_matrix *getR(){return _R;}

protected:
  bool setIncomingData( int frameX );
  void doEigenValueDecomposition();
  double calcObjectiveFunction( bool normalizeVariance = true );
private:
  gsl_matrix *calcCovarianceMatrix();

private:
  gsl_vector* search( int frameX = -5 );

protected:
  typedef list<VectorFloatFeatureStreamPtr>   _ChannelList; //typedef list<SampleFeaturePtr>   _ChannelList;
  typedef _ChannelList::iterator   _ChannelIterator;
  _ChannelList                     _channelList;

  SearchGridBuilderPtr _sgbPtr;
  vector<const gsl_vector_float *> _blockList;
  int *_tau; /* sample delays */
  gsl_vector *_x;
  gsl_matrix *_R; /* covariance matrix */
  gsl_matrix *_Rcopy;
  double _detR;
  void *_workspace;
  gsl_vector *_eigenvalues;
  SampleHolderPtr _shPtr; /* keep samples of the block processed at a prevous frame */
  vector<SourceCandidate *> _sourceCandidates; /* _sourceCandidates[n] is the N-th best candidate */
};

typedef Inherit<MCCLocalizer, VectorFeatureStreamPtr> MCCLocalizerPtr;

/**
   @class calculate multichannel cross correlation coefficient givin a position or time delays.
*/
class MCCCalculator : public MCCLocalizer {
public:
  MCCCalculator( SearchGridBuilderPtr &sgbPtr, bool normalizeVariance=true, const String& nm= "MCCCalculator" );
  ~MCCCalculator();
  void setTimeDelays( gsl_vector *delays );
  double getMCCC();
  double getCostV();
  virtual const gsl_vector* next(int frameX = -5);
  virtual void  reset();
private:
  gsl_matrix *calcCovarianceMatrix( gsl_vector *delays );

private:
  gsl_vector *_delays;
  bool        _normalizeVariance;
};

typedef Inherit<MCCCalculator, MCCLocalizerPtr> MCCCalculatorPtr;

/**
   @class recursively estimate the source position which provides the multi-channel cross correlation.
   @usage

   @note The algorithm implemented here is described in:
         J. Benesty, J. Chen and Y. Huang, "Time-Delay Estimation via Linear Interpolation and Cross Correlation", IEEE Trans. SAP, vol.12, Sep. 2004.
*/
class RMCCLocalizer : public MCCLocalizer {
public:
  RMCCLocalizer( SearchGridBuilderPtr &sgbPtr, float lambda, size_t maxSource=1, const String& nm= "RMCCSourceLocalizer" );
  ~RMCCLocalizer();
  virtual const gsl_vector* next(int frameX = -5);
  virtual void  reset();

private:
  void calcInverseMatrix();
  void updateParameters();

  gsl_matrix *_invR; /* an inverse matrix of the covariance matrix */
  gsl_vector *_kd; /* a priori Kalman gain vector */

  //gsl_vector *_old_x;
  gsl_matrix *_old_R; /* covariance matrix at a previous frame */
  gsl_matrix *_old_invR;
  gsl_vector *_old_kd; /* a priori Kalman gain vector at a previous frame */
};

typedef Inherit<RMCCLocalizer, MCCLocalizerPtr> RMCCLocalizerPtr;

#endif


text:
/* roots/gsl_roots.h
 *
 * 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.
 */

#ifndef __GSL_ROOTS_H__
#define __GSL_ROOTS_H__

#if !defined( GSL_FUN )
#  if !defined( GSL_DLL )
#    define GSL_FUN extern
#  elif defined( BUILD_GSL_DLL )
#    define GSL_FUN extern __declspec(dllexport)
#  else
#    define GSL_FUN extern __declspec(dllimport)
#  endif
#endif

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

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

__BEGIN_DECLS

typedef struct
  {
    const char *name;
    size_t size;
    int (*set) (void *state, gsl_function * f, double * root, double x_lower, double x_upper);
    int (*iterate) (void *state, gsl_function * f, double * root, double * x_lower, double * x_upper);
  }
gsl_root_fsolver_type;

typedef struct
  {
    const gsl_root_fsolver_type * type;
    gsl_function * function ;
    double root ;
    double x_lower;
    double x_upper;
    void *state;
  }
gsl_root_fsolver;

typedef struct
  {
    const char *name;
    size_t size;
    int (*set) (void *state, gsl_function_fdf * f, double * root);
    int (*iterate) (void *state, gsl_function_fdf * f, double * root);
  }
gsl_root_fdfsolver_type;

typedef struct
  {
    const gsl_root_fdfsolver_type * type;
    gsl_function_fdf * fdf ;
    double root ;
    void *state;
  }
gsl_root_fdfsolver;

GSL_FUN gsl_root_fsolver *
gsl_root_fsolver_alloc (const gsl_root_fsolver_type * T);
GSL_FUN void gsl_root_fsolver_free (gsl_root_fsolver * s);

GSL_FUN int gsl_root_fsolver_set (gsl_root_fsolver * s,
                          gsl_function * f,
                          double x_lower, double x_upper);

GSL_FUN int gsl_root_fsolver_iterate (gsl_root_fsolver * s);

GSL_FUN const char * gsl_root_fsolver_name (const gsl_root_fsolver * s);
GSL_FUN double gsl_root_fsolver_root (const gsl_root_fsolver * s);
GSL_FUN double gsl_root_fsolver_x_lower (const gsl_root_fsolver * s);
GSL_FUN double gsl_root_fsolver_x_upper (const gsl_root_fsolver * s);

GSL_FUN gsl_root_fdfsolver *
gsl_root_fdfsolver_alloc (const gsl_root_fdfsolver_type * T);

GSL_FUN int
gsl_root_fdfsolver_set (gsl_root_fdfsolver * s,
                         gsl_function_fdf * fdf, double root);

GSL_FUN int
gsl_root_fdfsolver_iterate (gsl_root_fdfsolver * s);

GSL_FUN void
gsl_root_fdfsolver_free (gsl_root_fdfsolver * s);

GSL_FUN const char * gsl_root_fdfsolver_name (const gsl_root_fdfsolver * s);
GSL_FUN double gsl_root_fdfsolver_root (const gsl_root_fdfsolver * s);

GSL_FUN int
gsl_root_test_interval (double x_lower, double x_upper, double epsabs, double epsrel);

GSL_FUN int
gsl_root_test_residual (double f, double epsabs);

GSL_FUN int
gsl_root_test_delta (double x1, double x0, double epsabs, double epsrel);

GSL_VAR const gsl_root_fsolver_type  * gsl_root_fsolver_bisection;
GSL_VAR const gsl_root_fsolver_type  * gsl_root_fsolver_brent;
GSL_VAR const gsl_root_fsolver_type  * gsl_root_fsolver_falsepos;
GSL_VAR const gsl_root_fdfsolver_type  * gsl_root_fdfsolver_newton;
GSL_VAR const gsl_root_fdfsolver_type  * gsl_root_fdfsolver_secant;
GSL_VAR const gsl_root_fdfsolver_type  * gsl_root_fdfsolver_steffenson;

__END_DECLS

#endif /* __GSL_ROOTS_H__ */

ext: h

max_stars_repo_name: berkus/music-cs

max_issues_repo_name: berkus/music-cs

max_forks_repo_name: berkus/music-cs


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

void integral_recur ( int nmin, int nmax, double vals[]);

void integral_recur (int nmin, int nmax, double vals[])
{
    double IN = 0;
    int i = nmax - 1;
    vals [nmax] = IN;
    while ( i >= nmin)
    {
        vals[i] = (IN + exp (-1.))/((double) i);
        i = i - 1;
    }
}

ext: c


text:
/* -*- Mode: C++; tab-width: 4; c-basic-offset: 4; indent-tabs-mode: nil -*- */
/*
 *     Copyright 2017-Present Couchbase, Inc.
 *
 *   Use of this software is governed by the Business Source License included
 *   in the file licenses/BSL-Couchbase.txt.  As of the Change Date specified
 *   in that file, in accordance with the Business Source License, use of this
 *   software will be governed by the Apache License, Version 2.0, included in
 *   the file licenses/APL2.txt.
 */
#pragma once

#include <mcbp/protocol/opcode.h>
#include <mcbp/protocol/request.h>
#include <mcbp/protocol/response.h>
#include <mcbp/protocol/status.h>

#include <gsl/gsl-lite.hpp>

namespace cb::mcbp {

/**
 * FrameBuilder allows you to build up a request / response frame in
 * a provided memory area. It provides helper methods which makes sure
 * that the fields is formatted in the correct byte order etc.
 *
 * The frame looks like:
 *
 *
 *     Byte/     0       |       1       |       2       |       3       |
 *        /              |               |               |               |
 *       |0 1 2 3 4 5 6 7|0 1 2 3 4 5 6 7|0 1 2 3 4 5 6 7|0 1 2 3 4 5 6 7|
 *       +---------------+---------------+---------------+---------------+
 *      0| Magic         |               | Framing extras| Key Length    |
 *       +---------------+---------------+---------------+---------------+
 *      4| Extras length |               |                               |
 *       +---------------+---------------+---------------+---------------+
 *      8| Total body length                                             |
 *       +---------------+---------------+---------------+---------------+
 *     12|                                                               |
 *       +---------------+---------------+---------------+---------------+
 *     16|                                                               |
 *       |                                                               |
 *       +---------------+---------------+---------------+---------------+
 *     24| Optional section containing framing extras Max 255 bytes      |
 *       +---------------------------------------------------------------+
 *       | Optional section containing command extras                    |
 *       +---------------------------------------------------------------+
 *       | Optional section containing the key                           |
 *       +---------------------------------------------------------------+
 *       | Optional section containing the value.                        |
 *       | The size of the value is the total body length minus the      |
 *       | size of the optional framing extras, extras and key.          |
 *       +---------------------------------------------------------------+
 *
 */
template <typename T>
class FrameBuilder {
public:
    explicit FrameBuilder(cb::char_buffer backing, bool initialized = false)
        : FrameBuilder(
                  {reinterpret_cast<uint8_t*>(backing.data()), backing.size()},
                  initialized) {
    }

    explicit FrameBuilder(cb::byte_buffer backing, bool initialized = false)
        : buffer(backing) {
        checkSize(sizeof(T));
        if (!initialized) {
            std::fill(backing.begin(), backing.begin() + sizeof(T), 0);
        }
    }

    T* getFrame() {
        return reinterpret_cast<T*>(buffer.data());
    }

    void setMagic(Magic magic) {
        getFrame()->setMagic(magic);
    }

    void setOpcode(ClientOpcode opcode) {
        getFrame()->setOpcode(opcode);
    }

    void setOpcode(ServerOpcode opcode) {
        getFrame()->setOpcode(opcode);
    }

    void setDatatype(Datatype datatype) {
        getFrame()->setDatatype(datatype);
    }

    void setVBucket(Vbid value) {
        getFrame()->setVBucket(value);
    }

    void setOpaque(uint32_t opaque) {
        getFrame()->setOpaque(opaque);
    }

    void setCas(uint64_t val) {
        getFrame()->setCas(val);
    }

    /**
     * Insert/replace the Framing extras section and move the rest of the
     * sections to their new locations.
     */
    void setFramingExtras(cb::const_byte_buffer val) {
        if (!is_alternative_encoding(getFrame()->getMagic())) {
            throw std::logic_error(
                    R"(setFramingExtras: Magic needs to be one of the alternative packet encodings)");
        }
        auto* req = getFrame();
        auto existing = req->getFramingExtras();
        // can we fit the data?
        checkSize(existing.size(), val.size());
        moveAndInsert(existing, val, req->getBodylen() - existing.size());
        req->setFramingExtraslen(gsl::narrow<uint8_t>(val.size()));
    }

    /**
     * Insert/replace the Extras section and move the rest of the sections
     * to their new locations.
     */
    void setExtras(cb::const_byte_buffer val) {
        auto* req = getFrame();
        auto existing = req->getExtdata();
        checkSize(existing.size(), val.size());
        auto move_size = req->getKey().size() + req->getValue().size();
        moveAndInsert(existing, val, move_size);
        req->setExtlen(gsl::narrow<uint8_t>(val.size()));
    }
    void setExtras(std::string_view val) {
        setExtras({reinterpret_cast<const uint8_t*>(val.data()), val.size()});
    }

    /**
     * Insert/replace the Key section and move the value section to the new
     * location.
     */
    void setKey(cb::const_byte_buffer val) {
        auto* req = getFrame();
        auto existing = req->getKey();
        checkSize(existing.size(), val.size());
        moveAndInsert(existing, val, req->getValue().size());
        req->setKeylen(gsl::narrow<uint16_t>(val.size()));
    }

    void setKey(std::string_view val) {
        setKey({reinterpret_cast<const uint8_t*>(val.data()), val.size()});
    }

    /**
     * Insert/replace the value section.
     */
    void setValue(cb::const_byte_buffer val) {
        auto* req = getFrame();
        auto existing = req->getValue();
        checkSize(existing.size(), val.size());
        moveAndInsert(existing, val, 0);
    }

    void setValue(std::string_view val) {
        setValue({reinterpret_cast<const uint8_t*>(val.data()), val.size()});
    }

    /**
     * Try to validate the underlying packet
     */
    void validate() {
        getFrame()->isValid();
    }

protected:
    /**
     * check that the requested size fits into the buffer
     *
     * @param size The total requested size (including header)
     */
    void checkSize(size_t size) {
        if (size > buffer.size()) {
            throw std::logic_error(
                    "FrameBuilder::checkSize: too big to fit in buffer");
        }
    }

    /**
     * check if there is room in the buffer to replace the field with
     * one with a different size
     *
     * @param oldfield
     * @param newfield
     */
    void checkSize(size_t oldfield, size_t newfield) {
        checkSize(sizeof(T) + getFrame()->getBodylen() - oldfield + newfield);
    }

    /**
     * Move the rest of the packet to the new location and insert the
     * new value
     */
    void moveAndInsert(cb::byte_buffer existing,
                       cb::const_byte_buffer value,
                       size_t size) {
        auto* location = existing.data();
        if (size > 0 && existing.size() != value.size()) {
            // we need to move the data following this entry to a different
            // location.
            memmove(location + value.size(), location + existing.size(), size);
        }
        std::copy(value.begin(), value.end(), location);
        auto* req = getFrame();
        // Update the total length
        req->setBodylen(gsl::narrow<uint32_t>(req->getBodylen() -
                                              existing.size() + value.size()));
    }

    cb::byte_buffer buffer;
};

/**
 * Specialized class to build a Request
 */
class RequestBuilder : public FrameBuilder<Request> {
public:
    explicit RequestBuilder(cb::char_buffer backing, bool initialized = false)
        : RequestBuilder(
                  {reinterpret_cast<uint8_t*>(backing.data()), backing.size()},
                  initialized) {
    }
    explicit RequestBuilder(cb::byte_buffer backing, bool initialized = false)
        : FrameBuilder<Request>(backing, initialized) {
    }
    void setVbucket(Vbid vbucket) {
        getFrame()->setVBucket(vbucket);
    }
};

/**
 * Specialized class to build a response
 */
class ResponseBuilder : public FrameBuilder<Response> {
public:
    explicit ResponseBuilder(cb::char_buffer backing, bool initialized = false)
        : ResponseBuilder(
                  {reinterpret_cast<uint8_t*>(backing.data()), backing.size()},
                  initialized) {
    }
    explicit ResponseBuilder(cb::byte_buffer backing, bool initialized = false)
        : FrameBuilder<Response>(backing, initialized) {
    }
    void setStatus(Status status) {
        getFrame()->setStatus(status);
    }
};

} // namespace cb::mcbp


text:
/**
 *
 * @file
 *
 *  PLASMA core_blas kernel
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.7.1
 * @author Mathieu Faverge
 * @date 2010-11-15
 **/
/*
 * @precisions normal z -> c d s
 */
#include <cblas.h>
#include "dplasma.h"
#include "dplasma_cores.h"
#include "dplasma_zcores.h"

/**
 ******************************************************************************
 *
 * @ingroup parsec_complex64
 *
 *  dplasma_core_ztradd adds to matrices together as in PBLAS pztradd.
 *
 *       B <- alpha * op(A)  + beta * B
 *
 *******************************************************************************
 *
 * @param[in] uplo
 *          Specifies the shape of A and B matrices:
 *          = PlasmaUpperLower: A and B are general matrices.
 *          = PlasmaUpper: op(A) and B are upper trapezoidal matrices.
 *          = PlasmaLower: op(A) and B are lower trapezoidal matrices.
 *
 * @param[in] trans
 *          Specifies whether the matrix A is non-transposed, transposed, or
 *          conjugate transposed
 *          = PlasmaNoTrans:   op(A) = A
 *          = PlasmaTrans:     op(A) = A'
 *          = PlasmaConjTrans: op(A) = conj(A')
 *
 * @param[in] M
 *          Number of rows of the matrices A and B.
 *
 * @param[in] N
 *          Number of columns of the matrices A and B.
 *
 * @param[in] alpha
 *          Scalar factor of A.
 *
 * @param[in] A
 *          Matrix of size LDA-by-N.
 *
 * @param[in] LDA
 *          Leading dimension of the array A. LDA >= max(1,M)
 *
 * @param[in] beta
 *          Scalar factor of B.
 *
 * @param[in,out] B
 *          Matrix of size LDB-by-N.
 *          On exit, B = alpha * op(A) + beta * B
 *
 * @param[in] LDB
 *          Leading dimension of the array B. LDB >= max(1,M)
 *
 *******************************************************************************
 *
 * @return
 *          \retval PLASMA_SUCCESS successful exit
 *          \retval <0 if -i, the i-th argument had an illegal value
 *
 ******************************************************************************/
int dplasma_core_ztradd(PLASMA_enum uplo, PLASMA_enum trans, int M, int N,
                              parsec_complex64_t  alpha,
                        const parsec_complex64_t *A, int LDA,
                              parsec_complex64_t  beta,
                              parsec_complex64_t *B, int LDB)
{
    static parsec_complex64_t zone = (parsec_complex64_t)1.;
    int j;

    if (uplo == PlasmaUpperLower){
        return dplasma_core_zgeadd( trans, M, N, alpha, A, LDA, beta, B, LDB );
    }

    if ((uplo != PlasmaUpper) &&
        (uplo != PlasmaLower))
    {
        dplasma_error("dplasma_core_ztradd", "illegal value of trans");
        return -1;
    }

    if ((trans != PlasmaNoTrans) &&
        (trans != PlasmaTrans)   &&
        (trans != PlasmaConjTrans))
    {
        dplasma_error("dplasma_core_ztradd", "illegal value of trans");
        return -2;
    }

    if (M < 0) {
        dplasma_error("dplasma_core_ztradd", "Illegal value of M");
        return -3;
    }
    if (N < 0) {
        dplasma_error("dplasma_core_ztradd", "Illegal value of N");
        return -4;
    }
    if ( ((trans == PlasmaNoTrans) && (LDA < dplasma_imax(1,M)) && (M > 0)) ||
         ((trans != PlasmaNoTrans) && (LDA < dplasma_imax(1,N)) && (N > 0)) )
    {
        dplasma_error("dplasma_core_ztradd", "Illegal value of LDA");
        return -7;
    }
    if ( (LDB < dplasma_imax(1,M)) && (M > 0) ) {
        dplasma_error("dplasma_core_ztradd", "Illegal value of LDB");
        return -9;
    }

    if (uplo == PlasmaLower) {
        switch( trans ) {
#if defined(PRECISION_z) || defined(PRECISION_c)
        case PlasmaConjTrans:
            for (j=0; j<N; j++, M--, A+=LDA+1, B+=LDB+1) {
                for(int i=0; i<M; i++) {
                    B[i] = beta * B[i] + alpha * conj(A[LDA*i]);
                }
            }
            break;
#endif /* defined(PRECISION_z) || defined(PRECISION_c) */

        case PlasmaTrans:
            for (j=0; j<N; j++, M--, A+=LDA+1, B+=LDB+1) {
                if (beta != zone) {
                    cblas_zscal(M, CBLAS_SADDR(beta), B, 1);
                }
                cblas_zaxpy(M, CBLAS_SADDR(alpha), A, LDA, B, 1);
            }
            break;

        case PlasmaNoTrans:
        default:
            for (j=0; j<N; j++, M--, A+=LDA+1, B+=LDB+1) {
                if (beta != zone) {
                    cblas_zscal(M, CBLAS_SADDR(beta), B, 1);
                }
                cblas_zaxpy(M, CBLAS_SADDR(alpha), A, 1, B, 1);
            }
        }
    }
    else {
        switch( trans ) {
#if defined(PRECISION_z) || defined(PRECISION_c)
        case PlasmaConjTrans:
            for (j=0; j<N; j++, A++, B+=LDB) {
                for(int i=0; i<=j; i++) {
                    B[i] = beta * B[i] + alpha * conj(A[LDA*i]);
                }
            }
            break;
#endif /* defined(PRECISION_z) || defined(PRECISION_c) */

        case PlasmaTrans:
            for (j=0; j<N; j++, A++, B+=LDB) {
                if (beta != zone) {
                    cblas_zscal(j+1, CBLAS_SADDR(beta), B, 1);
                }
                cblas_zaxpy(j+1, CBLAS_SADDR(alpha), A, LDA, B, 1);
            }
            break;

        case PlasmaNoTrans:
        default:
            for (j=0; j<N; j++, A+=LDA, B+=LDB) {
                if (beta != zone) {
                    cblas_zscal(j+1, CBLAS_SADDR(beta), B, 1);
                }
                cblas_zaxpy(j+1, CBLAS_SADDR(alpha), A, 1, B, 1);
            }
        }
    }

    return 0;
}

max_stars_repo_name: NLAFET/ABFT

max_issues_repo_name: NLAFET/ABFT

max_forks_repo_name: NLAFET/ABFT


text:
/* 生成sample for Regression */
/* theta12, m21 fixed */
/* theta23, theta13, m31 flat distribution in 3 sigma range */
/* deltacp flat distribution in 0~360 */

/* 使用方式(產100K組) :  ./sample_regression <filename.h5> <num> */
/* Caution: `num` should not greater than 1000. */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <sys/time.h>
#include <time.h>

#include <globes/globes.h>   /* GLoBES library */

#include <gsl/gsl_math.h>    /* GNU Scientific library (required for root finding) */
#include <gsl/gsl_roots.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_deriv.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_spline.h>
#include "hdf5.h"

#define degree M_PI/180

/***************************************************************************
 *                            M A I N   P R O G R A M                      *
 ***************************************************************************/

float keithRandom() {
	// Random number function based on the GNU Scientific Library
	// Returns a random float between 0 and 1, exclusive; e.g., (0,1)
	const 	gsl_rng_type * T;
	gsl_rng * r;
	gsl_rng_env_setup();
	struct 	timeval tv; // Seed generation based on time
	gettimeofday(&tv,0);
	unsigned long mySeed = tv.tv_sec + tv.tv_usec;
	T = gsl_rng_default; // Generator setup
	r = gsl_rng_alloc (T);
	gsl_rng_set(r, mySeed);
	double u = gsl_rng_uniform(r); // Generate it!
	gsl_rng_free (r);
	return (float)u;
}

double randn (double mu, double sigma) {
	/*mu is the central value, sigma is the width of gaussian distribution. */
	double U1, U2, W, mult;
	static double X1, X2;
	static int call = 0;

	if (call == 1) {
		call = !call;
		return (mu + sigma * (double) X2);
	}

	do{
		U1 = -1 + ((double) rand () / RAND_MAX) * 2;
		U2 = -1 + ((double) rand () / RAND_MAX) * 2;
		W = pow (U1, 2) + pow (U2, 2);
	} while (W >= 1 || W == 0);

	mult = sqrt ((-2 * log (W)) / W);
	X1 = U1 * mult;
	X2 = U2 * mult;

	call = !call;

	return (mu + sigma * (double) X1);
}

double TCRandom(double mu, double sigma) {
	// Random number function based on the GNU Scientific Library
	// Returns a random float between 0 and 1, exclusive; e.g., (0,1)
	const 	gsl_rng_type * T;
	gsl_rng * r;
	gsl_rng_env_setup();
	struct 	timeval tv; // Seed generation based on time
	gettimeofday(&tv,0);
	unsigned long mySeed = tv.tv_sec + tv.tv_usec;
	T = gsl_rng_default; // Generator setup
	r = gsl_rng_alloc (T);
	gsl_rng_set(r, mySeed);
	double u = mu + gsl_ran_gaussian(r, sigma); // Generate it!
	gsl_rng_free (r);

	return u;
}

int case_judger (int mode, double value) {
	/*
	 * mode 1: Octant;
	 * mode 2: CPV;
	 * mode 3: MO;
	 */
	int result;
	switch (mode){
		case 1:
			if ((value == 0) || (value==180)){
				result = 0;
			} else {
				result = 1;
			}
			break;
		case 2:
			if (value > 45) {
				result = 1;
			} else if (value == 45) {
				result = 0;
			} else {
				result = -1;
			}
			break;
	case 3:
		result = value > 0 ? 1 : -1;
		break;
	}
	return result;
}

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

	/* Initialize GLoBES and define chi^2 functions */
	glbInit(argv[0]);

	/*
	 * Read the prefix of file name and add suffix to two child string.
	 */
	char 	filename[32],
	 	 	Group_name[32],
		 	channel_info_dset_name[] = "/Spectrum/expr_0/channel_0\0",
		 	size_info_dset_name[] = "/Spectrum/expr_0/channel_0/Bin_ize\0",
		 	energy_info_dset_name[] = "/Spectrum/expr_0/channel_0/Bin_energy\0",
		 	expr_info_dset_name[] = "/Spectrum/expr_0\0";
	strcpy(filename, argv[1]);

	int  	TOTALsample = atof(argv[2]);

    /*
     * Declare the variables for hdf5.
     */
    hid_t	file_id,
			space_id,
			dset_id,
			group_id,
			sub_group_id,
			sub2_group_id,
			sub3_group_id;
    herr_t 	status;
    hsize_t	dims[2];

	/*
	 * Initialize experiment parameters.
	 */
	glbInitExperiment("./DUNE2021/DUNE_GLoBES.glb",&glb_experiment_list[0],&glb_num_of_exps);
	glbInitExperiment("./HK_globes/HK_combined_coarse.glb",&glb_experiment_list[0],&glb_num_of_exps);
	glb_params true_values = glbAllocParams();

	/* Set standard oscillation parameters (cf. hep-ph/0405172v5) */
	double theta12_c = 33.44;
	// double theta13_c = 8.57;
	// double theta23_c = 45;
	double sdm_c = 7.42;
	// double ldm_c = 2.514;

	double delta_12  = (35.86 - 31.27)/2;
	double delta_13  = (8.97 - 8.20)/2;
	double delta_23  = (51.80 - 39.60)/2;
	double delta_sdm = (8.04 - 6.82)/2;
	double delta_ldm = (2.598 - 2.431)/2;

	int    num_simulatiom;

	int    ew_low, ew_high;
	int    i, j, num_expriment, channel;

	file_id   = H5Fcreate(filename, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
	group_id = H5Gcreate(file_id, "/Spectrum", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
	for ( num_expriment = 0;  num_expriment < 2;  num_expriment++){
		double *bin_c_energy  =   glbGetBinCentersListPtr( num_expriment );
		double *bin_size      =   glbGetBinSizeListPtr( num_expriment );
		int    channel_max    =   glbGetNumberOfRules( num_expriment );

		expr_info_dset_name[15] = num_expriment+'0';
		sub_group_id = H5Gcreate(file_id, expr_info_dset_name, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
		for(channel=0;channel<channel_max;channel++){
			glbGetEnergyWindowBins( num_expriment, channel, &ew_low, &ew_high);
			int energy_window = ew_high - ew_low;
			hsize_t dims_s[1] = {energy_window};

			double dset_channel_eng[energy_window];
			double dset_channel_size[energy_window];

			j = 0;
			for ( i = ew_low; i <= ew_high; i++) {
				dset_channel_eng[j]  = bin_c_energy[i];
				dset_channel_size[j] = bin_size[i];
				j += 1;
			}

			channel_info_dset_name[15]  =  num_expriment+'0';
			channel_info_dset_name[25]  =  channel+'0';
			energy_info_dset_name[15]   =  num_expriment+'0';
			energy_info_dset_name[25]   =  channel+'0';
	 	    size_info_dset_name[15]   =  num_expriment+'0';
			size_info_dset_name[25]   =  channel+'0';

			sub2_group_id =  H5Gcreate(file_id, channel_info_dset_name, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
			space_id  =  H5Screate_simple(1, dims_s, NULL);

			dset_id       =  H5Dcreate(sub2_group_id, energy_info_dset_name, H5T_NATIVE_DOUBLE, space_id, H5P_DEFAULT, H5P_DEFAULT,
						            H5P_DEFAULT);
			status        =  H5Dwrite(dset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, dset_channel_eng);
			status        =  H5Dclose(dset_id);
			dset_id       =  H5Dcreate(sub2_group_id, size_info_dset_name, H5T_NATIVE_DOUBLE, space_id, H5P_DEFAULT, H5P_DEFAULT,
						            H5P_DEFAULT);
			status        =  H5Dwrite(dset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, dset_channel_size);
			status        =  H5Dclose(dset_id);
			status        =  H5Sclose(space_id);
			status        =  H5Gclose(sub2_group_id);
		}
		status        =  H5Gclose(sub_group_id);
	}
	status        =  H5Gclose(group_id);

	double     dataset[TOTALsample][9],
			   TRUE_RATE[TOTALsample][2][4][8],
			   TRUE_RATE_PRIME[2][4][8][TOTALsample];
	for (num_simulatiom=0;num_simulatiom<TOTALsample;num_simulatiom++){

		// double theta12 = TCRandom (theta12_c, delta_12 );
		// double theta13 = TCRandom (theta13_c, delta_13 );
		// double theta23 = TCRandom (theta23_c, delta_23 );
		// double sdm     = TCRandom (    sdm_c, delta_sdm);  //eV^2
		// double ldm     = TCRandom (    ldm_c, delta_ldm);  //eV^2

		// double theta12 = (35.86+31.27)/2 + (keithRandom()-0.5)*(35.86-31.27);
		double theta13 = ( 8.97+ 8.20)/2 + (keithRandom()-0.5)*( 8.97- 8.20);
		double theta23 = (45 + 45)/2 + (keithRandom()-0.5)*(51.80-39.60);
		// double sdm = ( 8.04 + 6.82)/2 + (keithRandom()-0.5)*( 8.04- 6.82);
		double ldm = (2.598 + 2.431)/2 + (keithRandom()-0.5)*(2.598-2.431);

		double sign_ldm = 1;
		if ( keithRandom() < 0.5) sign_ldm = -1;
		ldm = sign_ldm * ldm;

		// if(keithRandom()<0.33) {theta23 = 45;}

		double deltacp = keithRandom()*360.0*2.0-180.0;
		if (deltacp<0){
			deltacp=0;
		}	else if (deltacp>360){
			deltacp=180;
		}

		//double deltacp = keithRandom() * 360.0;

		int octant = case_judger(1, theta23);
		int cpv    = case_judger(2, deltacp);
		int mo     = case_judger(3, ldm);

		glbDefineParams(true_values,theta12_c*degree,theta13*degree,theta23*degree,deltacp*degree,1e-5*sdm_c,1e-3*ldm);
		glbSetDensityParams(true_values,1.0,GLB_ALL);

		glbSetOscillationParameters(true_values);
		glbSetRates();

		num_expriment = 0;
		channel = 0;

		for ( num_expriment = 0; num_expriment < 2; num_expriment++){
			int channel_max = glbGetNumberOfRules( num_expriment );
			for ( channel = 0; channel < channel_max ; channel++){
				glbGetEnergyWindowBins( num_expriment, channel, &ew_low, &ew_high);
				double *true_rates = glbGetRuleRatePtr( num_expriment, channel );
				int count = 0;
				for (i=ew_low; i <= ew_high; i++){
					TRUE_RATE[num_simulatiom][num_expriment][channel][count] = true_rates[i];
					count += 1;
				}
			}
		}
	 	double tmp_array[9] = {theta12_c, theta13, theta23, deltacp, sdm_c, ldm, octant, cpv, mo};
		int  count = 0;
		for ( count = 0; count < 9; count++){
			dataset[num_simulatiom][count] = tmp_array[count];
		}
	}
	int k, l;
	for ( i = 0; i < TOTALsample; i++){
		for ( j = 0; j < 2; j++){
			for ( k = 0; k < 4; k++){
				for ( l = 0; l < 8; l++){
			   		TRUE_RATE_PRIME[j][k][l][i] = TRUE_RATE[i][j][k][l];
			   	}
			}
		}
	}

	char 	channel_spec_true_name[] = "/Parameter/true/expr_0/channel_0\0",
		 	expr_spec_true_name[]    = "/Parameter/true/expr_0\0",
		 	expr_spec_sim_name[]     = "/Parameter/simulation\0";

	hsize_t dims_s[2]    =  {TOTALsample, 8};
	group_id     		 =  H5Gcreate(file_id, "/Parameter", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
	sub_group_id 		 =  H5Gcreate(file_id, "/Parameter/true/", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
	for ( i = 0; i < 2; i ++){
		expr_spec_true_name[21] = i+'0';
		sub2_group_id =  H5Gcreate(file_id, expr_spec_true_name, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
		for ( j = 0; j < 4; j++){
			channel_spec_true_name[21] = i+'0';
			channel_spec_true_name[31] = j+'0';
			space_id      =  H5Screate_simple(2, dims_s, NULL);
			dset_id       =  H5Dcreate(sub2_group_id, channel_spec_true_name, H5T_NATIVE_DOUBLE, space_id, H5P_DEFAULT, H5P_DEFAULT,
				                H5P_DEFAULT);
			status        =  H5Dwrite(dset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, TRUE_RATE_PRIME[i][j]);
			status        =  H5Dclose(dset_id);
			status        =  H5Sclose(space_id);
		}
		status        =  H5Gclose(sub2_group_id);
	}
    status        =  H5Gclose(sub_group_id);

	hsize_t dims_d[2]    = 	{TOTALsample, 9};
	sub_group_id  		 =  H5Gcreate(file_id, "/Parameter/simulation/", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
	space_id      		 =  H5Screate_simple(2, dims_d, NULL);
	dset_id       		 =  H5Dcreate(sub_group_id, "/Parameter/simulation/dataset", H5T_NATIVE_DOUBLE, space_id, H5P_DEFAULT, H5P_DEFAULT,
		                        H5P_DEFAULT);
	status        		 =  H5Dwrite(dset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, dataset);
	status        		 =  H5Dclose(dset_id);
	status        		 =  H5Sclose(space_id);
	status        		 =  H5Gclose(sub_group_id);
	status        		 =  H5Gclose(group_id);

	/* Clean up */
	glbFreeParams(true_values);

	status        =  H5Fclose(file_id);
	return 0;
}


text:
/*! \file Fitting.h
 *  \brief header file for fitting functions to data.
 */

#ifndef FITTING_H
#define FITTING_H

#include <cmath>
#include <iostream>
#include <limits>
#include <gsl/gsl_math.h>
#include <gsl/gsl_interp.h>
#include <gsl/gsl_monte.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_errno.h>
#ifdef HAVE_GSL22
#include <gsl/gsl_multifit_nlinear.h>
#endif
#include <Precision.h>
#include <Function.h>
#include <Matrix.h>
#include <GMatrix.h>
#include <Interpolate.h>
#include <ExtremaRootFinding.h>

using namespace std;
namespace Math
{
    /*! \brief Nonlinear fitting using Levenberg–Marquardt algorithm.

        Least Squares fitting of function f(x) from x[lower] to x[upper] to data y[lower]-y[upper]. Returns chi2 and altered params (which must be passed with initial guess)
        one can also pass a weight matrix (which is inverse of covariance matrix of data points), error on the relative change in chi2, and also a pointer that indicates whether
        the parameters are held fixed. finally, if data binned dof reduced by 1 and can specify maximum number of iterations
        Note the function does not explicitly take into account errors on the x and y points save in the use of the weight matrix which is the inverse of the covariance matrix
        Furthermore, the covariance matrix for the parameters is determined using the confidence level assuming normaly distributed errors. If that is not the case,
        the returned covariance matrix will not be correct and the degree to which it is correct cannot be correct for. Instead a more rigours determination of the confidence
        region used to determine the covariance matrix, such as MC methods using the data as a sample of the underlying distribution and fitting resulting data numerous times
        must be done. The resulting distribution of parameters can then be fit itself using a function if necessary.
    */
    Double_t FitNonLinLS(const math_function fitfunc, const math_function *difffuncs,
        const int nparams, Double_t *params, GMatrix &covar,
        const int npoints, const Double_t x[], const Double_t y[], GMatrix *Weight=NULL,
        Double_t error=1e-3, Double_t cl=0.95,
        int *fixparam=NULL, int binned=1, int maxiterations=1000, int iestimateerror=0,
        bool icallgsl=true);

    Double_t FitNonLinLSNoGSL(const math_function fitfunc, const math_function *difffuncs,
        const int nparams, Double_t *params, GMatrix &covar,
        const int npoints, const Double_t x[], const Double_t y[], GMatrix *Weight=NULL,
        Double_t error=1e-3, Double_t cl=0.95,
        int *fixparam=NULL, int binned=1, int maxiterations=1000, int iestimateerror=0);

#ifdef HAVE_GSL22
    Double_t FitNonLinLSWithGSL(const math_function fitfunc, const math_function *difffuncs,
        const int nparams, Double_t *params, GMatrix &covar,
        const int npoints, const Double_t x[], const Double_t y[],
        GMatrix *W, Double_t error, Double_t cl,
        int *fixparam, int binned, int maxiit, int iestimateerror);
#endif

    /*! \brief Given n points of data, determine optimal bin size and return binned data and bin centers.

        This is done by dividing the data range into N bins of width Delta. Count the number of events k_i that enter the i'th bin.
        Calculate the mean and variance of the number of events as kmean=1/N *Sum k_i and kvar=1/N*Sum (k_i -kmean)^2
        Then calculate C(Delta)= (2*kmean-kvar)/Delta^2, a cost function
        Repeat while changing Delta and find Delta' that minimizes C(Delta). Return Delta'
        Based on Shimazaki H. and Shinomoto S., A method for selecting the bin size of a time histogram Neural Computation (2007) Vol. 19(6), 1503-1527 \n
        Can also provide weights for the data points.
    */
    Double_t OptimalBins(const int npoints, Double_t *points, Double_t xmin, Double_t xmax, Double_t *weights=NULL);
}

#endif

max_stars_repo_name: ICRAR/NBodylib

max_issues_repo_name: ICRAR/NBodylib

max_forks_repo_name: ICRAR/NBodylib


text:
#include <cblas.h>
#include <assert.h>
#include <unistd.h>

#include "benchmark.hpp"

void dot_cblas(const double *A,
               const size_t A_m,
               const size_t A_n,
               const double *B,
               const size_t B_m,
               const size_t B_n,
               double *C) {
  UNUSED(B_m);
  assert(A != NULL && B != NULL && C != NULL);
  assert(A_m > 0 && A_n > 0 && B_m > 0 && B_n > 0);
  assert(A_n == B_m);

  cblas_dgemm(CblasColMajor, /* Matrix data arrangement */
              CblasNoTrans,  /* Transpose A */
              CblasNoTrans,  /* Transpose B */
              A_m,           /* Number of rows in A and C */
              B_n,           /* Number of cols in B and C */
              A_n,           /* Number of cols in A */
              1.0,           /* Scaling factor for the product of A and B */
              A,             /* Matrix A */
              A_n,           /* First dimension of A */
              B,             /* Matrix B */
              B_n,           /* First dimension of B */
              1.0,           /* Scale factor for C */
              C,             /* Output */
              A_m            /* First dimension of C */
  );
}

int main() {
  // Array - Blas
  for (size_t k = 10; k < 1000; k += 10) {
    size_t m = k;
    double *A = create_random_sq_matrix(m);
    double *B = create_random_sq_matrix(m);
    double *C = (double *) malloc(sizeof(double) * m * m);

    /* CBLAS dot() */
    /* sleep(1); */
    struct timespec t = tic();
    dot_cblas(A, m, m, B, m, m, C);
    printf("matrix_size: %ld\tdot_cblas(): %fs\n", m, toc(&t));

    free(A);
    free(B);
    free(C);
  }

  return 0;
}

max_stars_repo_name: daoran/proto

max_issues_repo_name: daoran/proto

max_forks_repo_name: daoran/proto


text:
#ifndef __GSL_VECTOR_H__
#define __GSL_VECTOR_H__

#include <gsl/gsl_vector_complex_long_double.h>
#include <gsl/gsl_vector_complex_double.h>
#include <gsl/gsl_vector_complex_float.h>

#include <gsl/gsl_vector_long_double.h>
#include <gsl/gsl_vector_double.h>
#include <gsl/gsl_vector_float.h>

#include <gsl/gsl_vector_ulong.h>
#include <gsl/gsl_vector_long.h>

#include <gsl/gsl_vector_uint.h>
#include <gsl/gsl_vector_int.h>

#include <gsl/gsl_vector_ushort.h>
#include <gsl/gsl_vector_short.h>

#include <gsl/gsl_vector_uchar.h>
#include <gsl/gsl_vector_char.h>

#endif /* __GSL_VECTOR_H__ */

max_stars_repo_licenses: Net-SNMP Xnet

max_issues_repo_licenses: Net-SNMP Xnet

max_forks_repo_licenses: Net-SNMP Xnet


text:
/*-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.

  File Name : leef.c

Purpose: embedding learning on a weighted graph

Creation Date : 15-01-2017

Last Modified : Wed 15 Mar 2017 03:07:53 PM CDT

Created By : Huan Gui (huangui2@illinois.edu)

_._._._._._._._._._._._._._._._._._._._._.*/

#include <assert.h>
#include <gsl/gsl_rng.h>
#include <math.h>
#include <pthread.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sys/time.h>
#include <unistd.h>
#include <set>

#include "ransampl.h"

#define MAX_STRING 100
#define MAX_SENTENCE 1000
#define MAX_SIGMOID 8.0
#define MAX_N_NODE_TYPE 10
#define MAX_N_EDGE_TYPE 10      // maximum number of edge types

typedef double real;          // Precision of float numbers

const int mini_batch = 1;
const int sigmoid_table_size = 2000;
const long hash_size = 30000000;     // For each type of nodes, there are 3M nodes
const int node_count_inc = 1000, neighbor_initial_count = 10;
const long neg_table_size = 1e8;
const long edge_table_size = 1e9;
const real sigmoid_coeff = sigmoid_table_size / MAX_SIGMOID / 2.0;
real sample_power = 0.75;  // defined for negative sampling as the power of frequency
real prior_power = 1;  // defined for negative sampling as the power of frequency

real max_sim = 0.5;
real max_iteration = 10;

struct ClassNeighbor {
  long node_id;
  long edge_id;
};

struct ClassNode {
  char * name;
  real degree;
  real weight;  // for the authority of nodes
  int *n_neighbor;
  int * max_n_neighbor;
  ClassNeighbor **neighbors;
};

struct ClassEdge {
  long src_id;
  long dst_id;
  real weight;
  real t_weight;
};

int n_node_type, n_edge_type;

// the node type for source node and dest node
int edge_node_type[MAX_N_EDGE_TYPE][2];

char output_folder[MAX_SENTENCE], input_folder[MAX_SENTENCE];
char embed_folder[MAX_SENTENCE], weight_folder[MAX_SENTENCE];
char train_network_file[MAX_N_EDGE_TYPE][MAX_SENTENCE];
char input_fix_embedding_file[MAX_SENTENCE];
char input_node_weight_file[MAX_SENTENCE];
char output_embedding[MAX_N_NODE_TYPE][MAX_SENTENCE];
char query_file[MAX_SENTENCE];

struct ClassNode* nodes[MAX_N_NODE_TYPE];
struct ClassEdge* edges[MAX_N_EDGE_TYPE];

set::set<long> query;

int author_type, term_type, paper_type, term_paper_edge_type;
int num_threads = 20, dim = 200;
int negative_k = 5;

// create hash table and negative table for each type of nodes and edges, respectively
long *node_hash_table[MAX_N_NODE_TYPE];
long *neg_table[MAX_N_NODE_TYPE];

ransampl_ws *edge_alias[MAX_N_EDGE_TYPE];
ransampl_ws *topic_edge_alias;

// The maximum size of each types of nodes, and the current count
long max_sz_node[MAX_N_NODE_TYPE], node_cnt[MAX_N_NODE_TYPE];

// total samples of edges to update parameters, and the current samples
long total_samples = 15, current_sample, edge_cnt[MAX_N_EDGE_TYPE];
real init_alpha = 0.025, alpha;

// the embedding of each type of nodes, relative to different types of nodes
real *emb_node[MAX_N_NODE_TYPE];

real *sigmoid_table;

//specify the types of nodes, and edge names
char node_names[MAX_N_NODE_TYPE][MAX_STRING];
char edge_names[MAX_N_EDGE_TYPE][MAX_STRING];

bool output_types[MAX_N_NODE_TYPE];
bool load_types[MAX_N_NODE_TYPE];
bool node_weight[MAX_N_NODE_TYPE];
bool edge_weight[MAX_N_NODE_TYPE];
bool valid_edge[MAX_N_NODE_TYPE][MAX_N_NODE_TYPE];

real beta[MAX_N_NODE_TYPE];
real node_sample_power[MAX_N_NODE_TYPE];

const gsl_rng_type * gsl_T;
gsl_rng * gsl_r;

clock_t start;

// the hash function is the same for all types of nodes
inline unsigned long GetHash(char *input) {
  unsigned long hash = 5381;
  char c;
  while ((c = *input) != '\0') {
    hash = ((hash << 5) + hash) + c; /* hash * 33 + c */
    ++input;
  }
  return hash % hash_size;
}

// initialize the hash table for all types of nodes
void InitHashTable() {
  for (int node_type = 0; node_type < n_node_type; ++node_type) {
    node_hash_table[node_type] = (long *) malloc(hash_size * sizeof(long));
    for (long k = 0; k < hash_size; ++k) node_hash_table[node_type][k] = -1;
  }
}

void InsertHashTable(char *key, int value, int type) {
  unsigned long loc = GetHash(key);
  while (node_hash_table[type][loc] != -1) loc = (loc + 1) % hash_size;
  node_hash_table[type][loc] = value;
}

// hashing based on the name string
long SearchHashTable(char *key, int type) {
  unsigned long addr = GetHash(key);
  while (1) {
    if (node_hash_table[type][addr] == -1) return -1;
    if (!strcmp(key, nodes[type][node_hash_table[type][addr]].name))
      return node_hash_table[type][addr];
    addr = (addr + 1) % hash_size;
  }
  return -1;
}

void AddNeighbor(int src_type, int dst_type, long src_id, long dst_id,
    int edge_type, long edge_id) {
  assert(valid_edge[src_type][dst_type]);
  ClassNode &node = nodes[src_type][src_id];
  int n_neighbor = node.n_neighbor[dst_type];

  if (n_neighbor + 2 >= node.max_n_neighbor[dst_type]) {
    node.max_n_neighbor[dst_type] *= 2;
    ClassNeighbor *tmp = (struct ClassNeighbor*) realloc(
        node.neighbors[dst_type],
        node.max_n_neighbor[dst_type] * sizeof(struct ClassNeighbor));
    if (tmp != NULL) {
      node.neighbors[dst_type] = tmp;
    } else {
      printf("Fail to reallocate\n");
      exit(1);
    }
  }
  node.neighbors[dst_type][n_neighbor].node_id = dst_id;
  node.neighbors[dst_type][n_neighbor].edge_id = edge_id;
  node.n_neighbor[dst_type] += 1;
}

// Add a node to the vertex set
int AddNode(char *name, int node_type) {
  int length = strlen(name) + 1;
  if (length > MAX_STRING) length = MAX_STRING;

  ClassNode & node = nodes[node_type][node_cnt[node_type]];
  // insert the node, and initial the value of name, node_type, and degree;
  node.name = (char *) calloc(length, sizeof(char));
  node.degree = 0;
  node.weight = 1;

  node.n_neighbor = (int *) calloc(n_node_type, sizeof(int));
  node.max_n_neighbor = (int *) calloc(n_node_type, sizeof(int));
  node.neighbors = (ClassNeighbor **) calloc(n_node_type, sizeof(ClassNeighbor *));

  for (int i = 0; i < n_node_type; ++i) {
    if (valid_edge[node_type][i]) {
      node.max_n_neighbor[i] = neighbor_initial_count;
      node.neighbors[i] = (ClassNeighbor *) calloc(neighbor_initial_count,
          sizeof(ClassNeighbor));
    }
  }

  if (length < MAX_STRING) {
    strcpy(node.name, name);
  } else {
    strncpy(node.name, name, MAX_STRING - 1);
    node.name[MAX_STRING - 1] = '\0';
  }

  ++node_cnt[node_type];
  if (node_cnt[node_type] + 2 >= max_sz_node[node_type]) {
    max_sz_node[node_type] += node_count_inc;
    struct ClassNode *tmp = (struct ClassNode *) realloc(nodes[node_type],
        max_sz_node[node_type] * sizeof(struct ClassNode));
    if (tmp != NULL) {
      nodes[node_type] = tmp;
    } else {
      printf("fail to reallocate\n");
    }
  }
  InsertHashTable(name, node_cnt[node_type] - 1, node_type);
  return node_cnt[node_type] - 1;
}

void UpdateEdgeWeight(int edge_type) {
  long src_id, dst_id;
  int src_type = edge_node_type[edge_type][0];
  int dst_type = edge_node_type[edge_type][1];

  for (long k = 0; k < edge_cnt[edge_type]; ++k) {
    src_id = edges[edge_type][k].src_id;
    dst_id = edges[edge_type][k].dst_id;

    nodes[src_type][src_id].degree -= edges[edge_type][k].weight;
    nodes[dst_type][dst_id].degree -= edges[edge_type][k].weight;

    edges[edge_type][k].weight *=
      nodes[src_type][src_id].weight * nodes[dst_type][dst_id].weight;

    nodes[src_type][src_id].degree += edges[edge_type][k].weight;
    nodes[dst_type][dst_id].degree += edges[edge_type][k].weight;
  }
  //  printf("Finish updating edge weight of edge %s \n", edge_names[edge_type]);
}

void ReadNodeWeight(int node_type) {
  FILE *fin;
  char field[MAX_STRING], str[10000];
  real weight;
  long node_id, n_lines = 0;

  sprintf(input_node_weight_file,
      "%s/%s.weights", weight_folder, node_names[node_type]);
  //  printf("%s\n", input_node_weight_file);

  fin = fopen(input_node_weight_file, "rb");
  while (fgets(str, sizeof(str), fin)) ++n_lines;
  fclose(fin);

  fin = fopen(input_node_weight_file, "rb");
  for (long k = 0; k < n_lines; ++k) {
    fscanf(fin, "%s %lf", field, &weight);
    if (prior_power != 1) {
      weight = pow(weight, prior_power);
    }
    node_id = SearchHashTable(field, node_type);
    if (node_id == -1) {
      printf("Invalid node name %s of type %s \n",
          field, node_names[node_type]);
      exit(1);
    }
    nodes[node_type][node_id].weight = weight;
  }
  //printf("Finish updating the weight for nodes %s \n", node_names[node_type]);
}

/* Read edges of different types in the network */
void ReadEdgeData(int edge_type) {
  FILE *fin;
  real weight;

  long src_id, dst_id;
  int src_type = edge_node_type[edge_type][0];
  int dst_type = edge_node_type[edge_type][1];
  char str[MAX_SENTENCE], src_name[MAX_STRING], dst_name[MAX_STRING];

  // Get how many number of edges we have
  fin = fopen(train_network_file[edge_type], "rb");

  printf("input: %s \n", train_network_file[edge_type]);
  edge_cnt[edge_type] = 0;
  while (fgets(str, sizeof(str), fin)) ++edge_cnt[edge_type];
  fclose(fin);
  printf("Number of edges (%s): %ld          \n", edge_names[edge_type], edge_cnt[edge_type]);

  edges[edge_type] = (ClassEdge *) malloc(
      edge_cnt[edge_type] * sizeof(ClassEdge));

  fin = fopen(train_network_file[edge_type], "rb");

  for (long k = 0; k < edge_cnt[edge_type]; ++k) {
    fscanf(fin, "%s %s %f", src_name, dst_name, &weight);
    assert(weight > 0);

    src_id = SearchHashTable(src_name, src_type);
    if (src_id == -1) {
      src_id = AddNode(src_name, src_type);
    }
    dst_id = SearchHashTable(dst_name, dst_type);
    if (dst_id == -1) {
      dst_id = AddNode(dst_name, dst_type);
    }
    edges[edge_type][k].src_id = src_id;
    edges[edge_type][k].dst_id = dst_id;
    edges[edge_type][k].weight = weight;

    if (src_type == author_type || src_type == term_type) {
      AddNeighbor(src_type, dst_type, src_id, dst_id, edge_type, edge_id);
      AddNeighbor(dst_type, src_type, dst_id, src_id, edge_type, edge_id);
    }
    // update node degree
    nodes[src_type][src_id].degree += weight;
    nodes[dst_type][dst_id].degree += weight;

    if (k % 50000 == 0) {
      printf("Loading edges : %.3lf%%%c", k / (real)(edge_cnt[edge_type] + 1) * 100, 13);
      fflush(stdout);
    }
  }
  fclose(fin);
}

/* Read the fixed embedding of terms */
void ReadFixEmbedding(int node_type) {
  FILE *fin;
  char str[(dim + 1) * MAX_STRING + 1000];
  long node_id;
  long fix_node_count = node_cnt[node_type], count = 0;
  char *token;

  sprintf(input_fix_embedding_file,
      "%s/%s.embeddings", embed_folder, node_names[node_type]);

  if (access(input_fix_embedding_file, F_OK) == -1) {
    printf("Can't read the file %s \n", input_fix_embedding_file);
    exit(1);
  }

  // Get how many number of edges we have
  fin = fopen(input_fix_embedding_file, "rb");

  // printf("debug: finish init\n");
  if (fgets(str, sizeof(str), fin) == NULL) {
    printf("Can't read the file %s \n", input_fix_embedding_file);
    exit(1);
  }
  while (fgets(str, sizeof(str), fin)) {
    // printf("debug: parse %s\n", str);
    if (count >= fix_node_count) break;
    ++count;

    token = strtok(str, " ");
    node_id = SearchHashTable(token, node_type);
    if (node_id == -1) {
      continue;
    }
    for (int k = 0; k < dim; ++k) {
      token = strtok(NULL, " ");
      if (token == NULL) {
        printf("Error: embedding length dismatch!\n");
        exit(1);
      }
      emb_node[node_type][node_id * dim + k] = atof(token) / 20;
    }
  }
  fclose(fin);
  printf("Embedding of node %s imported!\n", node_names[node_type]);
}

void InitNodeVector(int node_type, real *emb_node) {
  long offset = 0;
  for (long i = 0; i < node_cnt[node_type]; ++i) {
    for (int k = 0; k < dim; ++k) {
      emb_node[offset++] = (gsl_rng_uniform(gsl_r) - 0.5) / dim;
    }
  }
}

/* Initialize the vertex embedding and the context embedding */
void InitVector(){
  printf("Initializing node vector ...                       \n %c", 13);
  fflush(stdout);
  // initialize the target and context embedding of each type of node
  for (int node_type = 0; node_type < n_node_type; ++node_type) {
    posix_memalign((void **)& emb_node[node_type], 128,
        node_cnt[node_type] * dim * sizeof(real));
    if (emb_node[node_type] == NULL) {
      printf("Error: memory allocation failed\n");
      exit(1);
    }
    InitNodeVector(node_type, emb_node[node_type]);
  }
}

// Given a node type, sample a negative node
void InitNegTable(int node_type) {
  real sum = 0, cur_sum = 0, por = 0;
  long vid = 0;
  printf("Initializing negative table for %s ... \t\t\t\t          %c", node_names[node_type], 13);

  fflush(stdout);
  neg_table[node_type] = (long *)malloc(neg_table_size * sizeof(long));

  for (long k = 0; k < node_cnt[node_type]; ++k) {
    sum += pow(nodes[node_type][k].degree,
        node_sample_power[node_type]);
  }

  for (long k = 0; k < neg_table_size; ++k) {
    if ((real)(k + 1) / neg_table_size > por) {
      cur_sum += pow(nodes[node_type][vid].degree,
          node_sample_power[node_type]);
      por = cur_sum / sum;
      ++vid;
    }
    neg_table[node_type][k] = vid - 1;
  }
}

void InitEdgeAlias(int edge_type) {
  printf("Initializing alias table for edge %s ...                   %c",
      edge_names[edge_type], 13);
  real *prob = (real *) calloc(edge_cnt[edge_type], sizeof(real));
  real total_weight = 0;
  for (int i = 0; i < edge_cnt[edge_type]; ++i) {
    total_weight += edges[edge_type][i].weight;
  }
  for (int i = 0; i < edge_cnt[edge_type]; ++i) {
    prob[i] = edges[edge_type][i].weight / total_weight;
  }
  edge_alias[edge_type] = ransampl_alloc(edge_cnt[edge_type]);
  ransampl_set(edge_alias[edge_type], prob);
}

/* Fastly compute sigmoid function */
void InitSigmoidTable() {
  real x;
  printf("Initializing sigmoid table ...                     %c", 13);
  fflush(stdout);
  sigmoid_table = (real *)malloc((sigmoid_table_size + 1) * sizeof(real));
  for (int k = 0; k < sigmoid_table_size; ++k) {
    x = 2 * (real)(MAX_SIGMOID * k) / sigmoid_table_size - MAX_SIGMOID;
    sigmoid_table[k] = 1 / (1 + exp(-x));
  }
}

inline real FastSigmoid(real x) {
  if (x > MAX_SIGMOID) return 1;
  else if (x < -MAX_SIGMOID) return 0;
  return sigmoid_table[(int)((x + MAX_SIGMOID) * sigmoid_coeff)];
}

inline long SampleNegativeNode(int node_type, gsl_rng * & gsl_r_local) {
  //  return gsl_rng_uniform_int(gsl_r_local, node_cnt[node_type]);
  return neg_table[node_type][
    gsl_rng_uniform_int(gsl_r_local, neg_table_size)];
}

inline long SampleEdge(int edge_type, gsl_rng * & gsl_r_local) {
  return ransampl_draw(edge_alias[edge_type],
      gsl_rng_uniform(gsl_r_local),
      gsl_rng_uniform(gsl_r_local));
}

inline void Update(int src_type, int dst_type, long src_id,
    long end_id, long *neg_node_ids, real *context_update) {
  long target_id, target_offset, context_offset;
  real f = 0, g = 0, label = 1;
  real decay_target = alpha * beta[dst_type];
  real decay_context = alpha * beta[src_type];

  memset(context_update, 0, dim * sizeof(real));

  context_offset = src_id * dim;
  for (int i = 0; i <= negative_k; ++i) {
    if (i == negative_k) {
      target_id = end_id;
      label = 1;
    } else {
      target_id = neg_node_ids[i];
      if (target_id == end_id) continue;
      label = 0;
    }
    target_offset = target_id * dim;
    f = 0;
    for (int j = 0; j < dim; j ++) {
      f += emb_node[dst_type][target_offset + j] *
        emb_node[src_type][context_offset + j];
    }
    g = (label - FastSigmoid(f)) * alpha;
    for (int j = 0; j < dim; ++j) {
      context_update[j] += g * emb_node[dst_type][target_offset + j]
        - decay_context * emb_node[src_type][context_offset + j];
      emb_node[dst_type][target_offset + j] +=
        g * emb_node[src_type][context_offset + j]
        - decay_target * emb_node[dst_type][target_offset + j];
    }
  }

  for (int j = 0; j < dim; ++j) {
    emb_node[src_type][context_offset + j] += context_update[j];
  }

}

void *TrainThread(void *id) {
  long count = 0, last_count = 0;
  long limit = total_samples / num_threads + 2;
  int src_type, dst_type;
  long edge_id, src_id, dst_id;
  long  *neg_node_ids;

  neg_node_ids = (long *) malloc(negative_k * sizeof(long));

  gsl_rng * gsl_r_local = gsl_rng_alloc(gsl_T);
  gsl_rng_set(gsl_r_local, (unsigned long int) id);

  real *err_update = (real *) malloc(dim * sizeof(real));

  real ratio;
  clock_t now;

  while (1)
  {
    // decide when to exist the learning
    if (count > limit) {
      current_sample += count - last_count;
      break ;
    }

    // update output information
    if (count - last_count == 10000) {
      current_sample += count - last_count;
      last_count = count;
      ratio = current_sample / (real)(total_samples + 1);
      if (count % 100000 == 0) {
        now=clock();
        printf("%cAlpha: %f  Progress: %.2lf%%  Events/thread/sec: %.2fk",
            13, alpha, ratio * 100,
            current_sample / ((real)(now - start + 1) / (real)CLOCKS_PER_SEC * 1000));
        fflush(stdout);
      }
      alpha = init_alpha * (1 - ratio);
    }
    ++count;

    for (int edge_type = 0; edge_type < n_edge_type; ++edge_type) {
      src_type = edge_node_type[edge_type][0];
      dst_type = edge_node_type[edge_type][1];

      edge_id = SampleEdge(edge_type, gsl_r_local);
      src_id = edges[edge_type][edge_id].src_id;
      dst_id = edges[edge_type][edge_id].dst_id;

      //printf("%d %d %d %s %s \n", edge_type, src_type, dst_type, nodes[src_type][src_id].name, nodes[dst_type][dst_id].name);

      for (int i = 0; i < negative_k; ++i) {
        neg_node_ids[i] = SampleNegativeNode(dst_type, gsl_r_local);
      }
      Update(src_type, dst_type, src_id, dst_id, neg_node_ids, err_update);
    }
  }
  free(neg_node_ids);
  free(err_update);
  pthread_exit(NULL);
}

// output the embeddings of nodes and the context embedding
void Output() {
  // for each type of nodes, output the embeddings
  char file_name[MAX_STRING];
  FILE *fo;
  long offset;

  // output the embeddings of nodes
  for (int node_type = 0; node_type < n_node_type; ++node_type) {
    if (!output_types[node_type]) continue;
    printf("storing the embedding of nodes %s\n", node_names[node_type]);
    sprintf(file_name, "%s/%s.embeddings", output_folder, node_names[node_type]);
    fo = fopen(file_name, "wb");
    fprintf(fo, "%ld %d\n", node_cnt[node_type], dim);
    offset = 0;
    for (long i = 0; i < node_cnt[node_type]; ++i) {
      fprintf(fo, "%s ", nodes[node_type][i].name);
      for (int k = 0; k < dim; ++k) {
        fprintf(fo, "%lf ", emb_node[node_type][offset++]);
      }
      fprintf(fo, "\n");
    }
    fclose(fo);
  }
}

// query string is delimited by #
void readQuery(char *qString, set::set<long> query) {
  query.clear();
  char delim[] = "#";
  char *token = strtok(name, delim);
  long node_id;
  while (token != NULL) {
    node_id = SearchHashTable(token, term_type);
    if (node_id == -1) {
      node_id = AddNode(token, map_node_type);
    }
    query.insert(node_id);
    token = strtok(NULL, delim);
  }
}

void hitPapers(set::set<long> &query, set::set<long> &paper_set) {
  paper_set.clear();
  for (set::set<long>::iterator it = query.begin(); it != query.end(); ++it) {
    long term_id = *it;
    ClassNode &node = nodes[term_type][term_id];
    for (int j = 0; j < node.n_neighbor[paper_type]; ++j) {
      paper_set.insert(node.neighbors[paper_type][j]);
    }
  }
}

void nodeExpansion(int src_type, int dst_type, set::set<long> &src_set, set::set<long> &dst_set) {
  dst_set.clear();
  for (set::set<long>::iterator it = src_set.begin(); it != src_set.end(); ++it) {
    long src_id = *it;
    ClassNode &node = nodes[src_type][src_id];
    for (int j = 0; j < node.n_neighbor[dst_type]; ++j) {
      dst_set.insert(node.neighbors[dst_type][j]);
    }
  }
}

void queryExpansion(set::set<long> &paper_set, set::set<long> &query, set::set<long> &new_query) {
  new_query.clear();
  long t_n_edges = 0;
  std::map<long, real> topic_edges;
  for (set::set<long>::iterator it = paper_set.begin(), it != paper_set.end(); ++it) {
    long paper_id = *it;
    class &node = *it;
    for (int j = 0; j < node.n_neighbor[term_type]; ++j) {
      long edge_id = node.neighbors[term_type][j].edge_id;
      real weight = edges[term_paper_edge_type][edge_id].weight;
      topic_edges[edge_id] = weight;
    }
  }

  long * sample_topic_edges = (long *) malloc(topic_edges.size() * sizeof(long));
  real * sample_prob = (real *) malloc(topic_edges.size() * sizeof(real));

  real sum_weight = 0;
  int k = 0;
  for (std::map<long, real>::iterator it = topic_edges.begin(); it != topic_edges.end(); ++it) {
    sample_topic_edges[k] = it->first;
    sample_prob[k] = it->second;
    sum_weight += sample_prob[k];
    k ++;
  }
  topic_edge_alias = ransampl_alloc(topic_edges.size());
  ransampl_set(topic_edge_alias, sample_prob);
}

set::set<long> expandDocument(set::set<long> query) {
  long * new_query;
  set::set<long> author_set, paper_set;
  set::set<long> new_query;
  while (true) {
    hitPapers(query, paper_set);
    nodeExpansion(paper_type, author_type, paper_set, author_set);
    nodeExpansion(author_type, paper_type, author_set, paper_set);
    queryExpansion(paper_set, query, new_query);
  }
  hitPapers(query, paper_set);
  return paper_set;
}

void TrainModel() {
  gsl_rng_env_setup();
  gsl_T = gsl_rng_rand48;
  gsl_r = gsl_rng_alloc(gsl_T);
  gsl_rng_set(gsl_r, 314159265);

  InitHashTable();
  printf("HashTable Initialized! \n");

  for (int edge_type = 0; edge_type < n_edge_type; ++edge_type) {
    printf("Start to read edges of %s\n", edge_names[edge_type]);
    ReadEdgeData(edge_type);
  }

  printf("Node types: %d \n", n_node_type);
  for (int i = 0; i < n_node_type; i ++) {
    printf("%s ", node_names[i]);
  }
  printf("\n");

  for (int node_type = 0; node_type < n_node_type; ++node_type) {
    if (node_weight[node_type]) {
      ReadNodeWeight(node_type);
    }
  }

  printf("Edge types: %d \n", n_edge_type);
  for (int i = 0; i < n_edge_type; i ++) {
    printf("%s ", edge_names[i]);
  }
  printf("\n");

  printf("--------------------------------                            \n");
  for (int k = 0; k < n_node_type; ++k) {
    printf("Number of nodes (%s): %ld          \n", node_names[k], node_cnt[k]);
  }
  // map each ege into corresponding
  total_samples = (long) (1000000 * total_samples);
  printf("--------------------------------\n");
  printf("Dimension: %d\n", dim);
  printf("Initial alpha: %.3lf\n", alpha);
  printf("Thread : %d \n", num_threads);
  printf("--------------------------------\n");

  printf("Load node weight : \n");
  for (int node_type = 0; node_type < n_node_type; ++node_type) {
    printf("\t%s : %s \n", node_names[node_type],
        node_weight[node_type] ? "Yes":"No");
  }
  printf("\n");

  printf("Prior power : %f \n\n", prior_power);

  printf("Update edge weight : \n");
  for (int edge_type = 0; edge_type < n_edge_type; ++edge_type) {
    if (edge_weight[edge_type] && (node_weight[edge_node_type[edge_type][0]]
          || node_weight[edge_node_type[edge_type][1]])) {
      printf("\t%s : %s \n", edge_names[edge_type], "Yes");
      UpdateEdgeWeight(edge_type);
    }
    else {
      printf("\t%s : %s \n", edge_names[edge_type], "No");
    }
  }
  printf("\n");

  printf("Output node types:");
  for (int m = 0; m < n_node_type; ++m){
    if (output_types[m]) {
      printf(" %s", node_names[m]);
    }
  }
  printf("\n\n");

  InitVector();

  InitSigmoidTable();
  printf("Finish initialization!                \n");

  FILE *fin = fopen(query_file, "rb");
  char c;
  int nQuery = 0;
  while (fgets(str, sizeof(str), fin)) ++nQuery;
  fclose();

  fin = fopen(query_file, "rb");
  set::set<long> papers, query;
  char query_string[MAX_SENTENCE];
  for (int q = 0; q < nQuery; q++) {
    int pos = 0;
    do {
      c = fgetc(fin);
      if (c == EOF || c == '\n') break;
      else if (pos < MAX_SENTENCE - 1) query_string[pos++] = (char) c;
    }
    query_string[pos] = 0;
    readQuery(query_string, query);
    papers.clear();
    papers = expandDocument(query);

  }

  struct timeval t1, t2;
  gettimeofday(&t1, NULL);
  start = clock();
  pthread_t *pt = (pthread_t *)malloc((num_threads) * sizeof(pthread_t));

  for (long a = 0; a < num_threads; ++a) {
    pthread_create(&pt[a], NULL, TrainThread, (void *)a);
  }
  for (int a = 0; a < num_threads; ++a) {
    pthread_join(pt[a], NULL);
  }

  gettimeofday(&t2, NULL);
  printf("\nTotal optimization time: %.2lf minitues\n", (real)(t2.tv_sec - t1.tv_sec) / 60);
  Output();
}

int ArgPos(char *str, int argc, char **argv) {
  for (int a = 0; a < argc; ++a) {
    if (!strcmp(str, argv[a])) {
      if (a == argc - 1) {
        printf("Argument missing for %s\n", str);
        exit(1);
      }
      return a;
    }
  }
  return -1;
}

int main(int argc, char **argv) {
  int i, j, k;
  char str[MAX_SENTENCE], source[MAX_SENTENCE], dest[MAX_SENTENCE];
  strcpy(output_folder, "output");
  strcpy(input_folder, "data");
  InitSigmoidTable();

  if (argc == 1) {
    printf("\t leef: Local Embedding for Expert Finding\n");
    printf("Options:");
    printf("\tPlease specify the parameter file!\n");
    printf("\t-para <file> \n");
    printf("\t\tProvide the parameter file on edges\n");
    printf("\t-size <int>\n");
    printf("\t\tSet dimension of vertex embeddings; default is 200\n");
    printf("\t-scoring <int>\n");
    printf("\t\tScoring function, 0: complete element-wise product,"
        "1: complete pairwise inner product, 2: pairwise inner product with central type\n");
    printf("\t-iter <int>\n");
    printf("\t\tSet the number of training samples as the number of iterations of number of edges\n");
    printf("\t-tthread <int>\n");
    printf("\t\tUse <int> threads for tensor updating(default 1)\n");
    printf("\t-wthread <int>\n");
    printf("\t\tUse <int> threads for word2vec updating(default 1)\n");
    printf("\t-negative <int>\n");
    printf("\t\tSample <int> negative samples for word2vec updating(default 1)\n");
    printf("\t-alpha <float>\n");
    printf("\t\tSet the starting learning rate; default is 0.025\n");
    printf("\t-input input folder ...\n");
    printf("\t\tSet the folder data which read data\n");
    printf("\t-output output folder ...\n");
    printf("\t\tSet the folder to which output data\n");
    return 0;
  }

  // global parameters shared across all the edges
  if ((i = ArgPos((char *) "-size", argc, argv)) > 0)
    dim = atoi(argv[i + 1]);
  if ((i = ArgPos((char *) "-iter", argc, argv)) > 0)
    total_samples = atoi(argv[i + 1]);
  if ((i = ArgPos((char *) "-alpha", argc, argv)) > 0)
    init_alpha = atof(argv[i + 1]);
  if ((i = ArgPos((char *) "-thread", argc, argv)) > 0) {
    num_threads = atoi(argv[i + 1]);
  }
  if ((i = ArgPos((char *) "-negative", argc, argv)) > 0) {
    negative_k = atoi(argv[i + 1]);
  }

  for (i = 0; i < n_node_type; ++i) {
    for (j = 0; j < n_node_type; ++j) {
      valid_edge[i][j] = false;
    }
  }

  if ((i = ArgPos((char *) "-n_node", argc, argv)) > 0) {
    n_node_type = atoi(argv[i + 1]);
  }

  for (int i = 0; i < n_node_type; ++i) {
    nodes[i] =(struct ClassNode *) calloc(node_count_inc,
        sizeof(struct ClassNode));
    max_sz_node[i] = node_count_inc;
    node_sample_power[i] = sample_power;
  }

  if ((i = ArgPos((char *) "-nodes", argc, argv)) > 0) {
    for (j = 0; j < n_node_type; ++j) {
      strcpy(node_names[j], argv[i + j + 1]);
      if (!strcmp(node_names[j], "author")) {
        author_type = j;
      } else if (!strcmp(node_names[j], "term")) {
        term_type = j;
      } else if (!strcmp(node_names[j], "paper")) {
        paper_type = j;
      }
    }
  }

  if ((i = ArgPos((char *) "-n_edge", argc, argv)) > 0) {
    n_edge_type = atoi(argv[i + 1]);
  }

  if ((i = ArgPos((char *) "-output_folder", argc, argv)) > 0)
    strcpy(output_folder, argv[i + 1]);
  if ((i = ArgPos((char *) "-input_folder", argc, argv)) > 0)
    strcpy(input_folder, argv[i + 1]);
  if ((i = ArgPos((char *) "-embed_folder", argc, argv)) > 0)
    strcpy(embed_folder, argv[i + 1]);
  if ((i = ArgPos((char *) "-weight_folder", argc, argv)) > 0)
    strcpy(weight_folder, argv[i + 1]);

  if ((i = ArgPos((char *) "-edges", argc, argv)) > 0) {
    for (j = 0; j < n_edge_type; ++j) {
      strcpy(edge_names[j], argv[i + j + 1]);
      sprintf(train_network_file[j], "%s/%s.net", input_folder, edge_names[j]);
      printf("Input %d: %s\n", j, train_network_file[j]);
      strcpy(edge_names[j], argv[i + 1 + j]);
      strcpy(str, edge_names[j]);
      int pch = strchr(str, '_') - str;
      str[pch] = '\0';
      strcpy(source, str);
      strcpy(dest, &str[pch + 1]);
      for (k = 0; k < n_node_type; ++k) {
        if (strcmp(source, node_names[k]) == 0) {
          edge_node_type[j][0] = k;
        } else if (strcmp(dest, node_names[k]) == 0) {
          edge_node_type[j][1] = k;
        }
      }
      if (edge_node_type[j][0] == term_type && edge_node_type[j][1] == paper_type) {
        term_paper_edge_type = j;
      }
      valid_edge[edge_node_type[j][0]][edge_node_type[j][1]] = true;
      valid_edge[edge_node_type[j][1]][edge_node_type[j][0]] = true;
    }
  }

  // the regularization parameter for the embeddings
  if ((i = ArgPos((char *) "-beta", argc, argv)) >= 0) {
    for (int j = 0; j < n_node_type; ++j) {
      beta[j] = atof(argv[i + 1 + j]);
    }
  }

  if ((i = ArgPos((char *) "-prior_power", argc, argv)) >= 0) {
    prior_power = atof(argv[i + 1]);
  }

  if ((i = ArgPos((char *) "-sample_power", argc, argv)) >= 0) {
    for (int j = 0; j < n_node_type; ++j) {
      node_sample_power[j] = atof(argv[i + 1 + j]);
    }
  }

  if ((i = ArgPos((char *) "-choose_output", argc, argv)) >= 0) {
    for (int j = 0; j < n_node_type; ++j) {
      output_types[j] = (argv[i + 1][j] == '1');
    }
  }

  if ((i = ArgPos((char *) "-load_embed", argc, argv)) >= 0) {
    printf("Loading embeddings is disabled..\n");
    for (int j = 0; j < n_node_type; ++j) {
      if (argv[i + 1][j] == '1') {
        load_types[j] = true;
      }
    }
  }

  if ((i = ArgPos((char *) "-load_weight", argc, argv)) >= 0) {
    for (int j = 0; j < n_node_type; ++j) {
      if (argv[i + 1][j] == '1') {
        node_weight[j] = true;
      }
    }
  }

  if ((i = ArgPos((char *) "-edge_weight", argc, argv)) >= 0) {
    for (int j = 0; j < n_node_type; ++j) {
      if (argv[i + 1][j] == '1') {
        edge_weight[j] = true;
      } else {
        edge_weight[j] = false;
      }
    }
  }

  if ((i = ArgPos((char *) "-query_file", argc, argv)) >= 0) {
    strcpy(query_file, argv[i + 1]);
  }

  alpha = init_alpha;
  TrainModel();
  return 0;
}

max_stars_repo_name: VincentGaoHJ/taxogen-implementation

max_issues_repo_name: VincentGaoHJ/taxogen-implementation

max_forks_repo_name: VincentGaoHJ/taxogen-implementation


text:
/* 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 Foundation, Inc., 59 Temple Place, Suite
   330, Boston, MA 02111-1307 USA

   From Robert M. Ziff, "Four-tap shift-register-sequence
   random-number generators," Computers in Physics 12(4), Jul/Aug
   1998, pp 385-392.  A generalized feedback shift-register (GFSR)
   is basically an xor-sum of particular past lagged values.  A
   four-tap register looks like:
      ra[nd] = ra[nd-A] ^ ra[nd-B] ^ ra[nd-C] ^ ra[nd-D]

   Ziff notes that "it is now widely known" that two-tap registers
   have serious flaws, the most obvious one being the three-point
   correlation that comes from the defn of the generator.  Nice
   mathematical properties can be derived for GFSR's, and numerics
   bears out the claim that 4-tap GFSR's with appropriately chosen
   offsets are as random as can be measured, using the author's test.

   This implementation uses the values suggested the the author's
   example on p392, but altered to fit the GSL framework.  The "state"
   is 2^14 longs, or 64Kbytes; 2^14 is the smallest power of two that
   is larger than D, the largest offset.  We really only need a state
   with the last D values, but by going to a power of two, we can do a
   masking operation instead of a modulo, and this is presumably
   faster, though I haven't actually tried it.  The article actually
   suggested a short/fast hack:

   #define RandomInteger (++nd, ra[nd&M]=ra[(nd-A)&M]\
                          ^ra[(nd-B)&M]^ra[(nd-C)&M]^ra[(nd-D)&M])

   so that (as long as you've defined nd,ra[M+1]), then you ca do things
   like: 'if (RandomInteger < p) {...}'.

   Note that n&M varies from 0 to M, *including* M, so that the
   array has to be of size M+1.  Since M+1 is a power of two, n&M
   is a potentially quicker implementation of the equivalent n%(M+1).

   This implementation copyright (C) 1998 James Theiler, based on
   the example mt.c in the GSL, as implemented by Brian Gough.
*/

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

static inline unsigned long int gfsr4_get (void *vstate);
static double gfsr4_get_double (void *vstate);
static void gfsr4_set (void *state, unsigned long int s);

/* Magic numbers */
#define A 471
#define B 1586
#define C 6988
#define D 9689
#define M 16383 /* = 2^14-1 */
/* #define M 0x0003fff */

typedef struct
  {
    int nd;
    unsigned long ra[M+1];
  }
gfsr4_state_t;

static inline unsigned long
gfsr4_get (void *vstate)
{
  gfsr4_state_t *state = (gfsr4_state_t *) vstate;

  state->nd = ((state->nd)+1)&M;
  return state->ra[(state->nd)] =
      state->ra[((state->nd)+(M+1-A))&M]^
      state->ra[((state->nd)+(M+1-B))&M]^
      state->ra[((state->nd)+(M+1-C))&M]^
      state->ra[((state->nd)+(M+1-D))&M];

}

static double
gfsr4_get_double (void * vstate)
{
  return gfsr4_get (vstate) / 4294967296.0 ;
}

static void
gfsr4_set (void *vstate, unsigned long int s)
{
  gfsr4_state_t *state = (gfsr4_state_t *) vstate;
  int i, j;
  /* Masks for turning on the diagonal bit and turning off the
     leftmost bits */
  unsigned long int msb = 0x80000000UL;
  unsigned long int mask = 0xffffffffUL;

  if (s == 0)
    s = 4357;   /* the default seed is 4357 */

  /* We use the congruence s_{n+1} = (69069*s_n) mod 2^32 to
     initialize the state. This works because ANSI-C unsigned long
     integer arithmetic is automatically modulo 2^32 (or a higher
     power of two), so we can safely ignore overflow. */

#define LCG(n) ((69069 * n) & 0xffffffffUL)

  /* Brian Gough suggests this to avoid low-order bit correlations */
  for (i = 0; i <= M; i++)
    {
      unsigned long t = 0 ;
      unsigned long bit = msb ;
      for (j = 0; j < 32; j++)
        {
          s = LCG(s) ;
          if (s & msb)
            t |= bit ;
          bit >>= 1 ;
        }
      state->ra[i] = t ;
    }

  /* Perform the "orthogonalization" of the matrix */
  /* Based on the orthogonalization used in r250, as suggested initially
   * by Kirkpatrick and Stoll, and pointed out to me by Brian Gough
   */

  /* BJG: note that this orthogonalisation doesn't have any effect
     here because the the initial 6695 elements do not participate in
     the calculation.  For practical purposes this orthogonalisation
     is somewhat irrelevant, because the probability of the original
     sequence being degenerate should be exponentially small. */

  for (i=0; i<32; ++i) {
      int k=7+i*3;
      state->ra[k] &= mask;     /* Turn off bits left of the diagonal */
      state->ra[k] |= msb;      /* Turn on the diagonal bit           */
      mask >>= 1;
      msb >>= 1;
  }

  state->nd = i;
}

static const gsl_rng_type gfsr4_type =
{"gfsr4",                       /* name */
 0xffffffffUL,                  /* RAND_MAX  */
 0,                             /* RAND_MIN  */
 sizeof (gfsr4_state_t),
 &gfsr4_set,
 &gfsr4_get,
 &gfsr4_get_double};

const gsl_rng_type *gsl_rng_gfsr4 = &gfsr4_type;


text:
#ifndef STRUCTURES_H
#define STRUCTURES_H

#include <string>
#include <vector>

#include <gsl/gsl_matrix.h>

#define NSYSTEMS 5

#define T 3

#define UI 5
#define UJ 1

struct grain {
	//Euler angles in radians
	double phi1, phi, phi2;
	gsl_matrix *rot_mat;
	std::vector<double> crss;
	grain() {}
	grain(double p1, double p, double p2);
	~grain();
};

struct vect {
	double x, y, z;
};

struct slip_system {
	std::string label;
	//planes
	vect n;
	//directions
	vect b;
	double crss;
	//deformation matrix
	gsl_matrix* def_mat=0;

	//symetry part of deformation matrix
	gsl_matrix* symp_def=0;

	//non symetry part of deformation matrix
	gsl_matrix* asymp_def=0;

	std::vector<double> activity;

	slip_system(const std::string &l, const vect &in, const vect &ib, double r, int max_def);
	~slip_system();
};

struct slip_system_matrix {
    //m(ij)=bi*nj

    //slip systems matrix
    gsl_matrix *mat=0;

    //inverse slip systems matrix
    gsl_matrix *matinv=0;

    gsl_matrix *shear_stress=0;

    //combination of slip systems
    std::vector<int> comb;
    slip_system_matrix();
    ~slip_system_matrix();
};

struct handler {
    gsl_matrix *hand;
	handler();
	~handler();

};

struct U_handler {
    gsl_matrix *U_hand;
	U_handler();
	~U_handler();
};

#endif // STRUCTURES_H

max_stars_repo_name: bainit/hextexsim

max_issues_repo_name: bainit/hextexsim

max_forks_repo_name: bainit/hextexsim


text:
#pragma once

#include <chrono>
#include <fmt/printf.h>
#include <gsl/gsl-lite.hpp>

namespace angonoka::cli {
/**
    Pretty-print values with as much info as possible.

    @var value Value to be pretty-printed
*/
template <typename T> struct verbose {
    T value;
};
template <typename T> verbose(T) -> verbose<T>;

namespace detail {
    /**
        Helper function to print verbose durations.

        @param name Duration (s, m, d, etc)

        @return Formatter function
    */
    template <typename T>
    constexpr auto verbose_duration(gsl::czstring name)
    {
        return [=](auto& total, auto& ctx, bool& needs_space) {
            const auto dur = std::chrono::floor<T>(total);
            if (dur == dur.zero()) return;
            total -= dur;
            if (needs_space) fmt::format_to(ctx.out(), " ");
            needs_space = true;
            const auto ticks = dur.count();
            fmt::format_to(ctx.out(), "{}{}", ticks, name);
        };
    }
} // namespace detail
} // namespace angonoka::cli

namespace fmt {
using angonoka::cli::verbose;

/**
    User-defined formatter for std::chrono durations.
*/
template <typename... Ts>
struct fmt::formatter<verbose<std::chrono::duration<Ts...>>> {
    using value_type = verbose<std::chrono::duration<Ts...>>;
    static constexpr auto parse(format_parse_context& ctx)
    {
        return ctx.end();
    }

    template <typename FormatContext>
    constexpr auto format(const value_type& obj, FormatContext& ctx)
    {
        using angonoka::cli::detail::verbose_duration;
        using namespace std::chrono;
        if (obj.value == obj.value.zero())
            return format_to(ctx.out(), "0s");
        bool needs_space = false;
        auto total = duration_cast<seconds>(obj.value);
        [&](auto&&... fns) {
            (fns(total, ctx, needs_space), ...);
        }(verbose_duration<months>("mo"),
          verbose_duration<days>("d"),
          verbose_duration<hours>("h"),
          verbose_duration<minutes>("m"),
          verbose_duration<seconds>("s"));
        return ctx.out();
    }
};
} // namespace fmt

max_stars_repo_name: coffee-lord/angonoka

max_issues_repo_name: coffee-lord/angonoka

max_forks_repo_name: coffee-lord/angonoka


text:
/*
 # Copyright (C) 2009-2013 Thierry Sousbie
 # University of Tokyo / CNRS, 2009
 #
 # This file is part of porject DisPerSE
 #
 #  Author          : Thierry Sousbie
 #  Contact         : tsousbie@gmail.com
 #
 #  Licenses        : This file is 'dual-licensed', you have to choose one
 #                    of the two licenses below to apply.
 #
 #                    CeCILL-C
 #                    The CeCILL-C license is close to the GNU LGPL.
 #                    ( http://www.cecill.info/licences/Licence_CeCILL-C_V1-en.html )
 #
 #                or  CeCILL v2.0
 #                    The CeCILL license is compatible with the GNU GPL.
 #                    ( http://www.cecill.info/licences/Licence_CeCILL_V2-en.html )
 #
 #  This software is governed either by the CeCILL or the CeCILL-C license
 #  under French law and abiding by the rules of distribution of free software.
 #  You can  use, modify and or redistribute the software under the terms of
 #  the CeCILL or CeCILL-C licenses as circulated by CEA, CNRS and INRIA
 #  at the following URL : "http://www.cecill.info".
 #
 #  As a counterpart to the access to the source code and  rights to copy,
 #  modify and redistribute granted by the license, users are provided only
 #  with a limited warranty  and the software's author,  the holder of the
 #  economic rights,  and the successive licensors  have only  limited
 #  liability.
 #
 #  In this respect, the user's attention is drawn to the risks associated
 #  with loading,  using,  modifying and/or developing or reproducing the
 #  software by the user in light of its specific status of free software,
 #  that may mean  that it is complicated to manipulate,  and  that  also
 #  therefore means  that it is reserved for developers  and  experienced
 #  professionals having in-depth computer knowledge. Users are therefore
 #  encouraged to load and test the software's suitability as regards their
 #  requirements in conditions enabling the security of their systems and/or
 #  data to be ensured and,  more generally, to use and operate it in the
 #  same conditions as regards security.
 #
 #  The fact that you are presently reading this means that you have had
 #  knowledge of the CeCILL and CeCILL-C licenses and that you accept its terms.
*/
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "dtfi.h"

#ifdef HAVE_GSL
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

// mu is the average value
unsigned int randomPoisson(double mu)
{
  static int first_call = 1;
  static const gsl_rng_type *T;
  static gsl_rng * r;

  if (first_call) {
    first_call=0;
    gsl_rng_env_setup();
    T = gsl_rng_default;
    r = gsl_rng_alloc(T);
  }

  return gsl_ran_poisson (r, mu);
}

void randomSimplexCoords(double **coord, int ndims, double *result, double *P)
{
  static int first_call = 1;
  static const gsl_rng_type *T;
  static gsl_rng * r;
  int i,j,k;
  double u;
  //double P[50];

  if (first_call) {
    first_call=0;
    gsl_rng_env_setup();
    T = gsl_rng_default;
    r = gsl_rng_alloc (T);
  }

  result[0]=gsl_rng_uniform (r);

  if (ndims==3)
    {
      u=gsl_rng_uniform (r);
      if (u<result[0])
	{
	  result[1]=result[0];
	  result[0]=u;
	}
      else result[1]=u;

      u=gsl_rng_uniform (r);
      if (u<result[0])
	{
	  result[2]=result[1];
	  result[1]=result[0];
	  result[0]=u;
	}
      else if (u<result[1])
	{
	  result[2]=result[1];
	  result[1]=u;
	}
      else result[2]=u;

    }
  else
    {
      for (i=1;i<ndims;i++)
	{
	  u = gsl_rng_uniform (r);

	  for (j=0;j<i;j++)
	    if (u<result[j])
	      {
		result[i]=result[j];
		result[j]=u;
		break;
	      }

	  if (j==i) result[i]=u;
	  else {
	    j++;
	    for (;j<i;j++)
	      for (k=j+1;k<=i;k++)
		{
		  if (result[j]>result[k])
		    {
		      double tmp = result[k];
		      result[k]=result[j];
		      result[j]=tmp;
		    }
		}
	  }
	}
    }

  P[0]=result[0];
  for (i=1;i<ndims;i++) P[i]= result[i]-result[i-1];
  P[ndims]=1.-result[ndims-1];

  for (i=0;i<ndims;i++) result[i]=P[0]*coord[0][i];

  for (i=1;i<ndims+1;i++)
    for (j=0;j<ndims;j++)
      result[j]+=P[i]*coord[i][j];

}

void genDTFIInSimplex(double **coord, double *Proba, int ngen, int ndims, double **result_p)
{

  if (*result_p==NULL) {
    *result_p=(double*)malloc(sizeof(double)*ngen*ndims);
  }

  double *result = *result_p;
  double P[ndims];
  int i,n;
  double keepProba;
  double pmax_inv=0;

  static int first_call = 1;
  static const gsl_rng_type *T;
  static gsl_rng * r;

  //if (ngen>100) printf("negen = %d\n",ngen);

  if (first_call) {
    first_call=0;
    gsl_rng_env_setup();
    T = gsl_rng_default;
    r = gsl_rng_alloc (T);
  }

  if (Proba!=NULL) {
    pmax_inv=0;
    for (i=0;i<ndims+1;i++)
      if (Proba[i]>pmax_inv) pmax_inv=Proba[i];
    pmax_inv=(double)1./pmax_inv;
  }

  for (n=0;n<ngen;)
    {
      randomSimplexCoords(coord,ndims,&(result[n*ndims]),P);
      if (Proba!=NULL) {
	keepProba=Proba[0]*pmax_inv*P[0];
	for (i=1;i<ndims+1;i++) keepProba+=Proba[i]*pmax_inv*P[i];

	if (gsl_rng_uniform (r)<=keepProba) n++;
      } else n++;
    }

}

#else

unsigned int randomPoisson(double mu)
{
  fprintf(stderr,"ERROR: Function 'randomPoisson' needs GSL.\n");
  fprintf(stderr,"Please, recompile with GSL.\n");

  return -1;
}

void randomSimplexCoords(double **coord, int ndims, double *result, double *P)
{
  fprintf(stderr,"ERROR: Function 'randomSimplexCoords' needs GSL.\n");
  fprintf(stderr,"Please, recompile with GSL.\n");

  return -1;
}

void genDTFIInSimplex(double **coord, double *Proba, int ngen, int ndims, double **result_p)
{
  fprintf(stderr,"ERROR: Function 'genDTFIInSimplex' needs GSL.\n");
  fprintf(stderr,"Please, recompile with GSL.\n");
}

#endif