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

/* newton.c -- newton root finding algorithm

   This is the classical Newton-Raphson iteration.

   x[i+1] = x[i] - f(x[i])/f'(x[i])

*/

#include <config.h>

#include <stddef.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <float.h>

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

#include "roots.h"

typedef struct
  {
    double f, df;
  }
newton_state_t;

static int newton_init (void * vstate, gsl_function_fdf * fdf, double * root);
static int newton_iterate (void * vstate, gsl_function_fdf * fdf, double * root);

static int
newton_init (void * vstate, gsl_function_fdf * fdf, double * root)
{
  newton_state_t * state = (newton_state_t *) vstate;

  const double x = *root ;

  state->f = GSL_FN_FDF_EVAL_F (fdf, x);
  state->df = GSL_FN_FDF_EVAL_DF (fdf, x) ;

  return GSL_SUCCESS;

}

static int
newton_iterate (void * vstate, gsl_function_fdf * fdf, double * root)
{
  newton_state_t * state = (newton_state_t *) vstate;

  double root_new, f_new, df_new;

  if (state->df == 0.0)
    {
      GSL_ERROR("derivative is zero", GSL_EZERODIV);
    }

  root_new = *root - (state->f / state->df);

  *root = root_new ;

  GSL_FN_FDF_EVAL_F_DF(fdf, root_new, &f_new, &df_new);

  state->f = f_new ;
  state->df = df_new ;

  if (!finite(f_new))
    {
      GSL_ERROR ("function not continuous", GSL_EBADFUNC);
    }

  if (!finite (df_new))
    {
      GSL_ERROR ("function not differentiable", GSL_EBADFUNC);
    }

  return GSL_SUCCESS;
}

static const gsl_root_fdfsolver_type newton_type =
{"newton",				/* name */
 sizeof (newton_state_t),
 &newton_init,
 &newton_iterate};

const gsl_root_fdfsolver_type  * gsl_root_fdfsolver_newton = &newton_type;


text:
/*
  File Name: grad.c
  Last Modification Date:	3/9/94	10:03:45
  Current Version: grad.c	1.7
  %File Creation Date: 26/10/91
  Author: Ramesh Gopinath  <ramesh@dsp.rice.edu>

  Copyright: All software, documentation, and related files in this distribution
  are Copyright (c) 1993  Rice University

  Permission is granted for use and non-profit distribution providing that this
  notice be clearly maintained. The right to distribute any portion for profit
  or as part of any commercial product is specifically reserved for the author.
*/

/*
#################################################################################
grad.mex4:
function [g,f,h] = grad(k,M,N,fs); Computes the gradient of stopband energy of the
                   linear-phase prototype filter in a cosine-modulated filter bank
		   with respect to lattice parameters of the J lattices.

      Input:
            k : Vector of denormalized lattice parameters for the J lattices of
	        CMFB
            M : M-channel, M-band etc
            N : 2Mm, length of the prototype filter
           fs : Stopband edge (as a fraction of pi)
      Output:
            f : Stopband Energy
            h : Prototype filter (second half)
#################################################################################
*/

#ifndef PROTOTYPE_DESIGN_H
#define PROTOTYPE_DESIGN_H

#include <stdio.h>
#include <assert.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_linalg.h>
#include <gsl/gsl_blas.h>

#include "common/refcount.h"
#include "common/jexception.h"

// ----- definition for class `CosineModulatedPrototypeDesign' -----
//
class CosineModulatedPrototypeDesign {
 public:
  CosineModulatedPrototypeDesign(int M = 256, int N = 3072, double fs = 1.0);
  ~CosineModulatedPrototypeDesign();

  void fcn(const double* x, double* f);
  void grad(const double* x, double* g);

  int M() const { return _M; }
  int N() const { return _N; }
  int m() const { return _m; }
  int J() const { return _J; }

  // return (one-half) of the prototype filter
  const gsl_vector* proto();

 private:
  const int					_M;	// No. of channels
  const int					_N;	// filter length 2 * _M * _m
  const int					_M2;	// 2 * _M
  const int					_m;
  const int					_Mm;
  const int					_J;
  const int					_Jm;
  const int					_oddM;
  const int					_oddm;

  double*					_sinews;
  double*					_h;
  double*					_hh;
  double*					_Ph;
  double*					_jac;

  int*						_index;

  gsl_vector*					_proto;
};

double design_f(const gsl_vector* v, void* params);
void   design_df(const gsl_vector* v, void* params, gsl_vector* df);
void   design_fdf(const gsl_vector* v, void* params, double* f, gsl_vector* df);

// ----- definition for class `PrototypeDesignBase' -----
//
class PrototypeDesignBase {
public:
  PrototypeDesignBase(int M, int m, unsigned r, double wp, int tau);
  ~PrototypeDesignBase();

protected:
  virtual void _calculateAb();
  void _calculateC();
  void _svd(gsl_matrix* U, gsl_matrix* V, gsl_vector* S, gsl_vector* workSpace);
  gsl_matrix* _nullSpace(const gsl_matrix* A, double tolerance);
  gsl_vector* _pseudoInverse(const gsl_matrix* A, const gsl_vector* b, double tolerance);
  void _solveNonSingular(const gsl_matrix* H, const gsl_vector* c0, const gsl_matrix* P, double tolerance);
  void _solveSingular(gsl_matrix* Ac, gsl_vector* bc, gsl_matrix* K, gsl_vector* dp, double tolerance);
  double _condition(gsl_matrix* H, gsl_matrix* P);
  void _printMatrix(const gsl_matrix* A, FILE* fp = stdout) const;
  void _printVector(const gsl_vector* b, FILE* fp = stdout) const;

  void save(const String& fileName);

  const int					_M;	// the number of subbands
  const int					_m;
  const	double					_wp;
  const int					_R;
  const int					_D;
  const int					_L;
  	int					_tau;

  gsl_matrix*					_A;
  gsl_vector*					_b;
  gsl_matrix*					_cpA;
  gsl_vector*					_cpb;
  gsl_matrix*					_C;
  gsl_matrix*					_cpC;

  gsl_vector*					_singularVals;
  gsl_vector*					_scratch;
  gsl_vector*					_workSpace;

  gsl_vector*					_prototype;
};

// ----- definition for class `AnalysisOversampledDFTDesign' -----
//
typedef Countable AnalysisOversampledDFTDesignCountable;
class AnalysisOversampledDFTDesign : public AnalysisOversampledDFTDesignCountable, protected PrototypeDesignBase {
public:
  AnalysisOversampledDFTDesign(int M = 512, int m = 2, int r = 1, double wp = 1.0, int tau_h = -1 );
  ~AnalysisOversampledDFTDesign();

  // design analysis prototype
  // tolerance for double precision: 2.2204e-16
  // tolerance for single precision: 1.1921e-07
  const gsl_vector* design(double tolerance = 2.2204e-16);

  // calculate distortion measures
  const gsl_vector* calcError(bool doPrint = true);

  void save(const String& fileName) { PrototypeDesignBase::save(fileName); }

protected:
  virtual void _solve(double tolerance);
  double _passbandResponseError();
  double _inbandAliasingDistortion();

  gsl_vector*					_error;	// [0] the total response error
                                                        // [1] the residual aliasing distortion
                                                        // [2] objective fuction
};

typedef refcountable_ptr<AnalysisOversampledDFTDesign> AnalysisOversampledDFTDesignPtr;

// ----- definition for class `SynthesisOversampledDFTDesign' -----
//
typedef Countable SynthesisOversampledDFTDesignCountable;
class SynthesisOversampledDFTDesign : public SynthesisOversampledDFTDesignCountable, protected PrototypeDesignBase {
public:
  SynthesisOversampledDFTDesign(const gsl_vector* h, int M = 512, int m = 2, int r = 1,
				double v = 0.01, double wp = 1.0, int tau_T = -1 );
  ~SynthesisOversampledDFTDesign();

  // design synthesis prototype
  // tolerance for double precision: 2.2204e-16
  // tolerance for single precision: 1.1921e-07
  const gsl_vector* design(double tolerance = 2.2204e-16);

  // calculate distortion measures
  virtual const gsl_vector* calcError(bool doPrint = true);

  void save(const String& fileName) { PrototypeDesignBase::save(fileName); }

protected:
  void _calculateEfP();
  virtual void _solve(double tolerance);
  double _totalResponseError();
  double _residualAliasingDistortion();

        gsl_vector*				_h;
  const double					_v;

  gsl_matrix*					_E;
  gsl_matrix*					_P;
  gsl_vector*					_f;

  gsl_matrix*					_cpE;
  gsl_matrix*					_cpP;
  gsl_vector*					_cpf;
  gsl_vector*					_error;	// [0] the total response error
                                                        // [1] the residual aliasing distortion
                                                        // [2] objective fuction
};

typedef refcountable_ptr<SynthesisOversampledDFTDesign> SynthesisOversampledDFTDesignPtr;

// ----- definition for class `AnalysisNyquistMDesign' -----
//
class AnalysisNyquistMDesign : public AnalysisOversampledDFTDesign {
public:
  AnalysisNyquistMDesign(int M = 512, int m = 2, int r = 1, double wp = 1.0, int tau_h = -1 );
  ~AnalysisNyquistMDesign();

private:
  virtual void _solve(double tolerance);
  void _calculateFd();
  void _calculateKdp();

  gsl_matrix*					_F;
  gsl_vector*					_d;
  gsl_matrix*					_K;
  gsl_vector*					_dp;
};

typedef Inherit<AnalysisNyquistMDesign, AnalysisOversampledDFTDesignPtr>	AnalysisNyquistMDesignPtr;

// ----- definition for class `SynthesisNyquistMDesign' -----
//
class SynthesisNyquistMDesign : public SynthesisOversampledDFTDesign {
public:
  SynthesisNyquistMDesign(const gsl_vector* h, int M = 512, int m = 2, int r = 1,
			  double wp = 1.0, int tau_g = -1);
  ~SynthesisNyquistMDesign();

  virtual const gsl_vector* calcError(bool doPrint = true);

private:
  virtual void _solve(double tolerance);

  double _passbandResponseError();
  void _calculateHc0();
  void _calculateJcp();

  gsl_matrix*					_H;
  gsl_vector*					_c0;
  gsl_matrix*					_J;
  gsl_vector*					_cp;
};

typedef Inherit<SynthesisNyquistMDesign, SynthesisOversampledDFTDesignPtr> SynthesisNyquistMDesignPtr;

// ----- definition for class `SynthesisNyquistMDesignCompositeResponse' -----
//
class SynthesisNyquistMDesignCompositeResponse : public SynthesisNyquistMDesign {
public:
  SynthesisNyquistMDesignCompositeResponse(const gsl_vector* h, int M = 512, int m = 2, int r = 1,
					   double wp = 1.0, int tau = -1);
  ~SynthesisNyquistMDesignCompositeResponse();

private:
  virtual void _calculateAb();
  inline double _sinc(int m);

  double* 					_sincValues;
};

typedef Inherit<SynthesisNyquistMDesignCompositeResponse, SynthesisNyquistMDesignPtr> SynthesisNyquistMDesignCompositeResponsePtr;

double SynthesisNyquistMDesignCompositeResponse::_sinc(int m)
{
  int am     = abs(m);

  if (am >= 2 * _L)
    throw jindex_error("Problem: index %d > 2 * L = %d", am, 2 * _L);

  double val = _sincValues[am];
  if (val != -HUGE_VAL) return val;

  if (m == 0)
    val = _sincValues[am] = 1.0;
  else
    val = _sincValues[am] = sin(_wp * am) / (_wp * am);

  return val;
}

#endif

ext: h


text:
#include <stdio.h>
#include <gsl/gsl_fit.h>

int
main (void)
{
  int i, n = 4;
  double x[4] = { 1970, 1980, 1990, 2000 };
  double y[4] = {   12,   11,   14,   13 };
  double w[4] = {  0.1,  0.2,  0.3,  0.4 };

  double c0, c1, cov00, cov01, cov11, chisq;

  gsl_fit_wlinear (x, 1, w, 1, y, 1, n,
                   &c0, &c1, &cov00, &cov01, &cov11,
                   &chisq);

  printf ("# best fit: Y = %g + %g X\n", c0, c1);
  printf ("# covariance matrix:\n");
  printf ("# [ %g, %g\n#   %g, %g]\n",
          cov00, cov01, cov01, cov11);
  printf ("# chisq = %g\n", chisq);

  for (i = 0; i < n; i++)
    printf ("data: %g %g %g\n",
                   x[i], y[i], 1/sqrt(w[i]));

  printf ("\n");

  for (i = -30; i < 130; i++)
    {
      double xf = x[0] + (i/100.0) * (x[n-1] - x[0]);
      double yf, yf_err;

      gsl_fit_linear_est (xf,
                          c0, c1,
                          cov00, cov01, cov11,
                          &yf, &yf_err);

      printf ("fit: %g %g\n", xf, yf);
      printf ("hi : %g %g\n", xf, yf + yf_err);
      printf ("lo : %g %g\n", xf, yf - yf_err);
    }
  return 0;
}

ext: c

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


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

#pragma once

#include "iterators.h"
#include "macros.h"
#include "type_traits.h"
#include "string.h"
#include "path.h"

#include <cereal/cereal.hpp>

#include <algorithm>
#include <assert.h>
#include <cmath>
#include <gsl/gsl>
#include <iterator>
#include <map>
#include <sstream>
#include <string>
#include <vector>

namespace mira
{
    struct property_bag;

    struct property_interface
    {
        property_interface(const char* type, const char* name)
            : m_type{ type }
            , m_name{ name }
        {}

        virtual ~property_interface()
        {}

        const char* name() const
        {
            return m_name;
        }

        const char* type() const
        {
            return m_type;
        }

    private:
        const char* const m_type;
        const char* const m_name;
    };

    struct simple_property : public property_interface
    {
        virtual void from_string(const std::string& str) = 0;
        virtual void from_stream(std::istream& str) = 0;

        virtual void from_other(const property_interface& other) = 0;

        virtual std::type_index serializable_type() const = 0;

        virtual std::string to_string() const = 0;

        simple_property(const char* type, const char* name)
            : property_interface(type, name)
        {}
    };

    template<typename SerializableType>
    struct serializable_property : public simple_property
    {
        serializable_property(const char* type, const char* name)
            : simple_property(type, name)
        {}

        virtual ~serializable_property() {}

        std::type_index serializable_type() const override
        {
            return { typeid(SerializableType) };
        }

        // serialization functions
        template<class Archive>
        SerializableType save_minimal(Archive const&) const
        {
            return get_serializable();
        }

        template<class Archive>
        void load_minimal(Archive const&, SerializableType const& val)
        {
            set_serializable(val);
        }

    protected:
        virtual SerializableType get_serializable() const = 0;
        virtual void set_serializable(const SerializableType& value) = 0;
    };

    template<typename T> struct property;

    struct property_bag : public property_interface
    {
        bool set(const std::string& name, const std::string& value)
        {
            auto itr = m_properties.find(name.c_str());
            if (itr == m_properties.end())
                return false;

            itr->second->from_string(value);

            return true;
        }

        void from_other(const property_bag& other)
        {
            assert(strcmp(name(), other.name()) == 0);
            assert(other.m_properties.size() == m_properties.size());
            assert(other.m_groups.size() == m_groups.size());

            for (auto prop : m_properties)
            {
                prop.second->from_other(*other.m_properties.at(prop.first));
            }

            for (size_t i = 0; i < m_groups.size(); ++i)
            {
                assert(strcmp(m_groups[i]->name(), other.m_groups[i]->name()) == 0);
                m_groups[i]->from_other(*other.m_groups[i]);
            }
        }

        std::string get(const std::string& name) const
        {
            return m_properties.find(name.c_str())->second->to_string();
        }

        std::vector<simple_property*> properties()
        {
            std::vector<simple_property*> props;
            std::transform(m_properties.begin(), m_properties.end(), std::back_inserter(props), [](const auto& entry) {
                return entry.second;
            });
            return props;
        }

        std::vector<const simple_property*> properties() const
        {
            std::vector<const simple_property*> props;
            std::transform(m_properties.begin(), m_properties.end(), std::back_inserter(props), [](const auto& entry) {
                return entry.second;
            });
            return props;
        }

        gsl::span<property_bag*> groups()
        {
            return m_groups;
        }

        gsl::span<const property_bag*> groups() const
        {
            return { const_cast<property_bag const**>(m_groups.data()), (gsl::span<const property_bag*>::size_type)m_groups.size() };
        }

        // serialization function
        template<typename ArchiveT>
        void serialize(ArchiveT& ar)
        {
            using arch_func = void(ArchiveT & ar, property_interface * prop);

            static std::unordered_map<std::type_index, arch_func*> converters{
                { typeid(bool), &archive_prop<ArchiveT, serializable_property<bool>> },
                { typeid(int), &archive_prop<ArchiveT, serializable_property<int>> },
                { typeid(unsigned int), &archive_prop<ArchiveT, serializable_property<unsigned int>> },
                { typeid(uint32_t), &archive_prop<ArchiveT, serializable_property<uint32_t>> },
                { typeid(uint64_t), &archive_prop<ArchiveT, serializable_property<uint64_t>> },
                { typeid(size_t), &archive_prop<ArchiveT, serializable_property<size_t>> },
                { typeid(float), &archive_prop<ArchiveT, serializable_property<float>> },
                { typeid(std::string), &archive_prop<ArchiveT, serializable_property<std::string>> },
                { typeid(double), &archive_prop<ArchiveT, serializable_property<double>> },
            };

            for (auto prop : m_properties)
            {
                auto found = converters.find(prop.second->serializable_type());
                if (found == converters.end())
                    throw std::invalid_argument("unable to convert type, conversion function not found");

                found->second(ar, prop.second);
            }

            for (auto group : m_groups)
            {
                ar(cereal::make_nvp(group->name(), *group));
            }
        }

    protected:
        property_bag(const char* name)
            : property_bag{ name, nullptr }
        {}

        property_bag(const char* name, property_bag* owner)
            : property_interface{ "bag", name }
        {
            if (owner != nullptr)
            {
                owner->add(this);
            }
        }

        property_bag(const property_bag& other) = delete;
        property_bag& operator=(const property_bag& other) = delete;

    private:
        template<typename ArchiveT, typename T>
        static void archive_prop(ArchiveT& ar, property_interface* prop)
        {
            ar(cereal::make_nvp(prop->name(), *static_cast<T*>(prop)));
        }

        template<typename T>
        friend struct property;

        void add(simple_property* prop)
        {
            m_properties.insert({ prop->name(), prop });
        }

        void add(property_bag* prop)
        {
            m_groups.push_back(prop);
        }

        std::map<const char*, simple_property*, string_compare> m_properties;

        std::vector<property_bag*> m_groups;
    };

    template<typename T, bool isEnum = std::is_enum<T>::value>
    struct property_serialization_traits;

    template<typename T>
    struct property_serialization_traits<T, true>
    {
        using serialization_type = std::underlying_type_t<T>;

        static serialization_type to_serializable(const T& value)
        {
            return static_cast<serialization_type>(value);
        }

        static T from_serializable(const serialization_type& value)
        {
            return static_cast<T>(value);
        }

        static void from_stream(std::istream& stream, T& value)
        {
            serialization_type underlyingValue{};

            stream >> underlyingValue;

            assert(!stream.fail() && "Property bag failed to convert string to value");

            value = static_cast<T>(underlyingValue);
        }

        static void to_stream(std::ostream& stream, const T& value)
        {
            auto underlying = static_cast<serialization_type>(value);

            std::stringstream ss;
            stream.precision(std::numeric_limits<serialization_type>::max_digits10);
            stream << underlying;

            assert(!stream.fail() && "Property bag failed to convert value to string");
        }
    };

    template<typename T>
    struct property_serialization_traits<T, false>
    {
        using serialization_type = T;

        static serialization_type to_serializable(const T& value)
        {
            return value;
        }

        static T from_serializable(const serialization_type& value)
        {
            return value;
        }

        static void from_stream(std::istream& stream, T& value)
        {
            if (std::is_same<T, bool>())
            {
                stream >> std::boolalpha >> value;
            }
            else
            {
                stream >> value;
            }

            assert((is_string<T>() || !stream.fail()) && "Property bag failed to convert string to value");
        }

        static void to_stream(std::ostream& stream, const T& value)
        {
            if (std::is_same<T, bool>())
            {
                stream << std::boolalpha << value;
            }
            else
            {
                stream.precision(std::numeric_limits<T>::max_digits10);
                stream << value;
            }
            assert(!stream.fail() && "Property bag failed to convert value to string");
        }
    };

    template<>
    struct property_serialization_traits<path, false>
    {
        using serialization_type = std::string;

        static serialization_type to_serializable(const path& value)
        {
            return value.string();
        }

        static path from_serializable(const serialization_type& value)
        {
            return { value };
        }

        static void from_stream(std::istream& stream, path& value)
        {
            stream >> value;
        }

        static void to_stream(std::ostream& stream, const path& value)
        {
            stream << value;
        }
    };

    template<typename T>
    struct property : public serializable_property<typename property_serialization_traits<T>::serialization_type>
    {
        using traits = property_serialization_traits<T>;

        using type = std::decay_t<T>;

        property(property_bag* lookup, const char* name, const char* type, const T& defaultValue)
            : serializable_property<typename traits::serialization_type>{ type, name }
            , value{ defaultValue }
        {
            lookup->add(this);
        }

        property& operator=(const property& other)
        {
            value = other.value;
            return *this;
        }

        virtual void from_string(const std::string& str) override
        {
            std::stringstream ss{ str };
            from_stream(ss);
        }

        virtual void from_stream(std::istream& stream) override
        {
            traits::from_stream(stream, value);
        }

        virtual void from_other(const property_interface& other) override
        {
            assert(dynamic_cast<const property<T>*>(&other) && "invalid type");
            value = static_cast<const property<T>&>(other).value;
        }

        virtual std::string to_string() const override
        {
            std::stringstream ss;
            traits::to_stream(ss, value);
            return ss.str();
        }

        operator T() const
        {
            return value;
        }

        property& operator=(const T& val)
        {
            value = val;
            return *this;
        }

        T value;

    protected:
        typename traits::serialization_type get_serializable() const override
        {
            return traits::to_serializable(value);
        }

        void set_serializable(const typename traits::serialization_type& val) override
        {
            value = traits::from_serializable(val);
        }
    };

    template<typename T>
    std::ostream& operator<<(std::ostream& os, const property<T>& prop)
    {
        return os << prop.value;
    }
}

#define PROPERTY(TYPE, NAME, DEFAULT) ::mira::property<TYPE> NAME{ this, #NAME, #TYPE, DEFAULT };

#define WEIRDLY_NAMED_PROPERTY(TYPE, STRNAME, DEFAULT)                                                                 \
    ::mira::property<TYPE> CONCATENATE_MACRO(setting_, __COUNTER__){ this, STRNAME, #TYPE, DEFAULT };

#define BAG_PROPERTY(TYPE) TYPE TYPE{ this };

#define NAMED_BAG_PROPERTY(TYPE, NAME) TYPE NAME{ this, #NAME };

#define PROPERTYBAG(NAME, ...)                                                                                \
    struct NAME : public ::mira::property_bag                                                                          \
    {                                                                                \
        NAME(const NAME& other)                                                                                \
            : NAME{ other.name() }                                                                                \
        {                                                                                \
            from_other(other);                                                                                \
        }                                                                                \
        NAME& operator=(const NAME& other)                                                                             \
        {                                                                                \
            from_other(other);                                                                                \
            return *this;                                                                                \
        }                                                                                \
        NAME()                                                                                \
            : ::mira::property_bag{ #NAME }                                                                            \
        {}                                                                                \
        NAME(::mira::property_bag* owner)                                                                              \
            : ::mira::property_bag{ #NAME, owner }                                                                     \
        {}                                                                                \
        NAME(const char* name)                                                                                \
            : ::mira::property_bag{ name }                                                                             \
        {}                                                                                \
        NAME(::mira::property_bag* owner, const char* name)                                                            \
            : ::mira::property_bag{ name, owner }                                                                      \
        {}                                                                                \
        __VA_ARGS__                                                                                \
    }

max_stars_repo_name: syntheticmagus/mageslam

max_issues_repo_name: syntheticmagus/mageslam

max_forks_repo_name: syntheticmagus/mageslam


text:
/**
 * \author Sylvain Marsat, University of Maryland - NASA GSFC
 *
 * \brief C code headers for the geometric coefficients entering the response for LISA-like detectors.
 *
 */

#ifndef _LISAGEOMETRY_H
#define _LISAGEOMETRY_H

#define _XOPEN_SOURCE 500

#ifdef __GNUC__
#define UNUSED __attribute__ ((unused))
#else
#define UNUSED
#endif

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <complex.h>
#include <time.h>
#include <unistd.h>
#include <getopt.h>
#include <stdbool.h>
#include <string.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_bspline.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_min.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_complex.h>

#include "constants.h"
#include "struct.h"

#if defined(__cplusplus)
#define complex _Complex
extern "C" {
#elif 0
} /* so that editors will match preceding brace */
#endif

/*****************************************************/
/**************** TDI variables **********************/

typedef enum TDItag {
  delayO,                  /* Orbital delay */
  y12L,                    /* Constellation-only response y12L */
  y12,                     /* Complete response y12 (includes orbital delay) */
  TDIXYZ,
  TDIalphabetagamma,
  TDIAETXYZ,
  TDIAETalphabetagamma,
  TDIX,
  TDIalpha,
  TDIAXYZ,
  TDIEXYZ,
  TDITXYZ,
  TDIAalphabetagamma,
  TDIEalphabetagamma,
  TDITalphabetagamma
} TDItag;

/* Enumerator to choose what level of low-f approximation to do in the response */
typedef enum ResponseApproxtag {
  full,
  lowfL,
  lowf
} ResponseApproxtag;

/**************************************************/
/************** LISAconstellation *****************/

typedef enum LISANoiseType{
  LISA2010noise,
  LISA2017noise,
  LISAProposalnoise
} LISANoiseType;

typedef struct tagLISAconstellation {
  double OrbitOmega;
  double OrbitPhi0;
  double OrbitR;
  double ConstOmega;
  double ConstPhi0;
  double ConstL;
  LISANoiseType noise;
} LISAconstellation;
extern LISAconstellation LISAProposal;
extern LISAconstellation LISA2017;
extern LISAconstellation LISA2010;
extern LISAconstellation slowOrbitLISA;
extern LISAconstellation tinyOrbitLISA;
extern LISAconstellation fastOrbitLISA;
extern LISAconstellation bigOrbitLISA;

/**************************************************/
/**************** Prototypes **********************/

/* Function to convert string input TDI string to TDItag */
TDItag ParseTDItag(char* string);

/* Function to convert string input ResponseApprox to tag */
ResponseApproxtag ParseResponseApproxtag(char* string);

/* Function cardinal sine */
double sinc(const double x);

/* Compute Solar System Barycenter time tSSB from retarded time at the center of the LISA constellation tL */
/* NOTE: depends on the sky position given in SSB parameters */
double tSSBfromLframe(const LISAconstellation *variant, const double tL, const double lambdaSSB, const double betaSSB);
/* Compute retarded time at the center of the LISA constellation tL from Solar System Barycenter time tSSB */
double tLfromSSBframe(const LISAconstellation *variant, const double tSSB, const double lambdaSSB, const double betaSSB);
/* Convert L-frame params to SSB-frame params */
/* NOTE: no transformation of the phase -- approximant-dependence with e.g. EOBNRv2HMROM setting phiRef at fRef, and freedom in definition */
int ConvertLframeParamsToSSBframe(
  double* tSSB,
  double* lambdaSSB,
  double* betaSSB,
  double* polSSB,
  const double tL,
  const double lambdaL,
  const double betaL,
  const double psiL,
  const LISAconstellation *variant);
/* Convert SSB-frame params to L-frame params */
/* NOTE: no transformation of the phase -- approximant-dependence with e.g. EOBNRv2HMROM setting phiRef at fRef, and freedom in definition */
int ConvertSSBframeParamsToLframe(
  double* tL,
  double* lambdaL,
  double* betaL,
  double* polL,
  const double tSSB,
  const double lambdaSSB,
  const double betaSSB,
  const double psiSSB,
  const LISAconstellation *variant);

/* Function to compute, given a value of a sky position and polarization, all the complicated time-independent trigonometric coefficients entering the response */
void SetCoeffsG(const double lambda, const double beta, const double psi);

/* Functions evaluating the G_AB functions, combining the two polarization with the spherical harmonics factors */
double complex G21mode(const LISAconstellation *LISAvariant, const double f, const double t, const double complex Yfactorplus, const double complex Yfactorcross);
double complex G12mode(const LISAconstellation *LISAvariant, const double f, const double t, const double complex Yfactorplus, const double complex Yfactorcross);
double complex G32mode(const LISAconstellation *LISAvariant, const double f, const double t, const double complex Yfactorplus, const double complex Yfactorcross);
double complex G23mode(const LISAconstellation *LISAvariant, const double f, const double t, const double complex Yfactorplus, const double complex Yfactorcross);
double complex G13mode(const LISAconstellation *LISAvariant, const double f, const double t, const double complex Yfactorplus, const double complex Yfactorcross);
double complex G31mode(const LISAconstellation *LISAvariant, const double f, const double t, const double complex Yfactorplus, const double complex Yfactorcross);
int EvaluateGABmode(
  const LISAconstellation *LISAvariant,
  double complex* G12,                     /* Output for G12 */
  double complex* G21,                     /* Output for G21 */
  double complex* G23,                     /* Output for G23 */
  double complex* G32,                     /* Output for G32 */
  double complex* G31,                     /* Output for G31 */
  double complex* G13,                     /* Output for G13 */
  const double f,                          /* Frequency */
  const double t,                          /* Time */
  const double complex Yfactorplus,        /* Spin-weighted spherical harmonic factor for plus */
  const double complex Yfactorcross,       /* Spin-weighted spherical harmonic factor for cross */
  const int tagdelayR,                     /* Tag: when 1, include the phase term of the R-delay */
  const ResponseApproxtag responseapprox); /* Tag to select possible low-f approximation level in FD response */

/* Functions evaluating the Fourier-domain factors (combinations of the GAB's) for TDI observables */
/* NOTE: factors have been scaled out, in parallel of what is done for the noise function */
/* Note: in case only one channel is considered, amplitudes for channels 2 and 3 are simply set to 0 */
/* (allows minimal changes from the old structure that assumed KTV A,E,T - but probably not optimal) */
int EvaluateTDIfactor3Chan(
  const LISAconstellation *variant,    /* Description of LISA variant */
  double complex* factor1,                       /* Output for factor for TDI channel 1 */
  double complex* factor2,                       /* Output for factor for TDI channel 2 */
  double complex* factor3,                       /* Output for factor for TDI channel 3 */
  const double complex G12,                      /* Input for G12 */
  const double complex G21,                      /* Input for G21 */
  const double complex G23,                      /* Input for G23 */
  const double complex G32,                      /* Input for G32 */
  const double complex G31,                      /* Input for G31 */
  const double complex G13,                      /* Input for G13 */
  const double f,                                /* Frequency */
  const TDItag tditag,                           /* Selector for the TDI observables */
  const ResponseApproxtag responseapprox);       /* Tag to select possible low-f approximation level in FD response */
/* int EvaluateTDIfactor1Chan( */
/*   double complex* factor,                       /\* Output for factor for TDI channel *\/ */
/*   const double complex G12,                      /\* Input for G12 *\/ */
/*   const double complex G21,                      /\* Input for G21 *\/ */
/*   const double complex G23,                      /\* Input for G23 *\/ */
/*   const double complex G32,                      /\* Input for G32 *\/ */
/*   const double complex G31,                      /\* Input for G31 *\/ */
/*   const double complex G13,                      /\* Input for G13 *\/ */
/*   const double f,                                /\* Frequency *\/ */
/*   const TDItag tditag);                          /\* Selector for the TDI observable *\/ */
/* Function evaluating the Fourier-domain factors that have been scaled out of TDI observables */
/* The factors scaled out, parallel what is done for the noise functions */
/* Note: in case only one channel is considered, factors for channels 2 and 3 are simply set to 0 */
int ScaledTDIfactor3Chan(
  const LISAconstellation *variant,    /* Description of LISA variant */
  double complex* factor1,                       /* Output for factor for TDI factor 1 */
  double complex* factor2,                       /* Output for factor for TDI factor 2 */
  double complex* factor3,                       /* Output for factor for TDI factor 3 */
  const double f,                                /* Frequency */
  const TDItag tditag);                          /* Selector for the TDI observables */
/* Function restoring the factor that have been scaled out of the TDI observables */
/* NOTE: the operation is made in-place, and the input is overwritten */
int RestoreInPlaceScaledFactorTDI(
  const LISAconstellation *variant,    /* Description of LISA variant */
  ListmodesCAmpPhaseFrequencySeries* listtdi,     /* Output/Input: list of mode contributions to TDI observable */
  TDItag tditag,                                  /* Tag selecting the TDI observable */
  int nchannel);                                  /* TDI channel number */

/* Functions for the response in time domain */

/* Basic yslr observables (including orbital delay) from hplus, hcross */
double y12TD(
  const LISAconstellation *variant,    /* Description of LISA variant */
  gsl_spline* splinehp,                    /* Input spline for TD hplus */
  gsl_spline* splinehc,                    /* Input spline for TD hcross */
  gsl_interp_accel* accelhp,               /* Accelerator for hp spline */
  gsl_interp_accel* accelhc,               /* Accelerator for hc spline */
  const double t);                         /* Time */
double y21TD(
  const LISAconstellation *variant,    /* Description of LISA variant */
  gsl_spline* splinehp,                    /* Input spline for TD hplus */
  gsl_spline* splinehc,                    /* Input spline for TD hcross */
  gsl_interp_accel* accelhp,               /* Accelerator for hp spline */
  gsl_interp_accel* accelhc,               /* Accelerator for hc spline */
  const double t);                         /* Time */
double y23TD(
  const LISAconstellation *variant,    /* Description of LISA variant */
  gsl_spline* splinehp,                    /* Input spline for TD hplus */
  gsl_spline* splinehc,                    /* Input spline for TD hcross */
  gsl_interp_accel* accelhp,               /* Accelerator for hp spline */
  gsl_interp_accel* accelhc,               /* Accelerator for hc spline */
  const double t);                         /* Time */
double y32TD(
  const LISAconstellation *variant,    /* Description of LISA variant */
  gsl_spline* splinehp,                    /* Input spline for TD hplus */
  gsl_spline* splinehc,                    /* Input spline for TD hcross */
  gsl_interp_accel* accelhp,               /* Accelerator for hp spline */
  gsl_interp_accel* accelhc,               /* Accelerator for hc spline */
  const double t);                         /* Time */
double y31TD(
  const LISAconstellation *variant,    /* Description of LISA variant */
  gsl_spline* splinehp,                    /* Input spline for TD hplus */
  gsl_spline* splinehc,                    /* Input spline for TD hcross */
  gsl_interp_accel* accelhp,               /* Accelerator for hp spline */
  gsl_interp_accel* accelhc,               /* Accelerator for hc spline */
  const double t);                         /* Time */
double y13TD(
  const LISAconstellation *variant,    /* Description of LISA variant */
  gsl_spline* splinehp,                    /* Input spline for TD hplus */
  gsl_spline* splinehc,                    /* Input spline for TD hcross */
  gsl_interp_accel* accelhp,               /* Accelerator for hp spline */
  gsl_interp_accel* accelhc,               /* Accelerator for hc spline */
  const double t);                         /* Time */
/* TDI observables (including orbital delay) from hplus, hcross */
int EvaluateTDIXYZTDhphc(
  const LISAconstellation *variant,    /* Description of LISA variant */
  double* TDIX,                            /* Output: value of TDI observable X */
  double* TDIY,                            /* Output: value of TDI observable Y */
  double* TDIZ,                            /* Output: value of TDI observable Z */
  gsl_spline* splinehp,                    /* Input spline for TD hplus */
  gsl_spline* splinehc,                    /* Input spline for TD hcross */
  gsl_interp_accel* accelhp,               /* Accelerator for hp spline */
  gsl_interp_accel* accelhc,               /* Accelerator for hc spline */
  const double t);                         /* Time */
int EvaluateTDIAETXYZTDhphc(
  const LISAconstellation *variant,    /* Description of LISA variant */
  double* TDIA,                            /* Output: value of TDI observable X */
  double* TDIE,                            /* Output: value of TDI observable Y */
  double* TDIT,                            /* Output: value of TDI observable Z */
  gsl_spline* splinehp,                    /* Input spline for TD hplus */
  gsl_spline* splinehc,                    /* Input spline for TD hcross */
  gsl_interp_accel* accelhp,               /* Accelerator for hp spline */
  gsl_interp_accel* accelhc,               /* Accelerator for hc spline */
  const double t);                         /* Time */

/* Generate hO orbital-delayed for one mode contribution from amp, phase */
int Generateh22TDO(
  const LISAconstellation *variant,    /* Description of LISA variant */
  AmpPhaseTimeSeries** h22tdO,             /* Output: amp/phase time series for h22TDO */
  gsl_spline* splineamp,                   /* Input spline for TD mode amplitude */
  gsl_spline* splinephase,                 /* Input spline for TD mode phase */
  gsl_interp_accel* accelamp,              /* Accelerator for amp spline */
  gsl_interp_accel* accelphase,            /* Accelerator for phase spline */
  gsl_vector* times,                       /* Vector of times to evaluate */
  int nbptmargin);                         /* Margin set to 0 on both side to avoid problems with delays out of the domain */
/* Generate y12L from orbital-delayed h22 in amp/phase form */
int Generatey12LTD(
  const LISAconstellation *variant,    /* Description of LISA variant */
  RealTimeSeries** y12Ltd,                 /* Output: real time series for y12L */
  gsl_spline* splineamp,                   /* Input spline for TD mode amplitude */
  gsl_spline* splinephase,                 /* Input spline for TD mode phase */
  gsl_interp_accel* accelamp,              /* Accelerator for amp spline */
  gsl_interp_accel* accelphase,            /* Accelerator for phase spline */
  gsl_vector* times,                       /* Vector of times to evaluate */
  double Theta,                            /* Inclination */
  double Phi,                              /* Phase */
  int nbptmargin);                         /* Margin set to 0 on both side to avoid problems with delays out of the domain */

/* Generate TDI observables (including orbital delay) for one mode contritbution from amp, phase */
/* NOTE: developed for testing purposes, so far supports only dO and y12L */
int GenerateTDITD3Chanhlm(
  AmpPhaseTimeSeries** hlm,                /* Output: real time series for TDI channel 1 */
  RealTimeSeries** TDI2,                   /* Output: real time series for TDI channel 2 */
  RealTimeSeries** TDI3,                   /* Output: real time series for TDI channel 3 */
  gsl_spline* splineamp,                   /* Input spline for TD mode amplitude */
  gsl_spline* splinephase,                 /* Input spline for TD mode phase */
  gsl_interp_accel* accelamp,              /* Accelerator for amp spline */
  gsl_interp_accel* accelphase,            /* Accelerator for phase spline */
  gsl_vector* times,                       /* Vector of times to evaluate */
  int nbptsmargin,                         /* Margin set to 0 on both side to avoid problems with delays out of the domain */
  double theta,                            /* Inclination angle - used to convert hlm to hplus, hcross for y12L - ignored for dO */
  double phi,                              /* Observer phase - used to convert hlm to hplus, hcross for y12L - ignored for dO */
  TDItag tditag);                          /* Tag selecting the TDI observables */

/* Generate TDI observables (including orbital delay) for one mode contritbution from hplus, hcross */
int GenerateTDITD3Chanhphc(
  const LISAconstellation *variant,    /* Description of LISA variant */
  RealTimeSeries** TDI1,                   /* Output: real time series for TDI channel 1 */
  RealTimeSeries** TDI2,                   /* Output: real time series for TDI channel 2 */
  RealTimeSeries** TDI3,                   /* Output: real time series for TDI channel 3 */
  gsl_spline* splinehp,                    /* Input spline for TD hplus */
  gsl_spline* splinehc,                    /* Input spline for TD hcross */
  gsl_interp_accel* accelhp,               /* Accelerator for hp spline */
  gsl_interp_accel* accelhc,               /* Accelerator for hc spline */
  gsl_vector* times,                       /* Vector of times to evaluate */
  int nbptsmargin,                         /* Margin set to 0 on both side to avoid problems with delays out of the domain */
  TDItag tditag);                          /* Tag selecting the TDI observables */

#if 0
{ /* so that editors will match succeeding brace */
#elif defined(__cplusplus)
}
#endif

#endif /* _LISAGEOMETRY_H */

max_stars_repo_name: JohnGBaker/flare

max_issues_repo_name: JohnGBaker/flare

max_forks_repo_name: JohnGBaker/flare


text:
#pragma once

#include <gsl\gsl>
#include <winrt\Windows.Foundation.h>
#include <d3d11.h>
#include "DrawableGameComponent.h"
#include "MatrixHelper.h"

namespace Rendering
{
	class AmbientLightingMaterial;

	class AmbientLightingDemo final : public Library::DrawableGameComponent
	{
	public:
		AmbientLightingDemo(Library::Game& game, const std::shared_ptr<Library::Camera>& camera);

		bool AnimationEnabled() const;
		void SetAnimationEnabled(bool enabled);

		float AmbientLightIntensity() const;
		void SetAmbientLightIntensity(float intensity);

		virtual void Initialize() override;
		virtual void Update(const Library::GameTime& gameTime) override;
		virtual void Draw(const Library::GameTime& gameTime) override;

	private:
		inline static const float RotationRate{ DirectX::XM_PI };

		std::shared_ptr<AmbientLightingMaterial> mMaterial;
		DirectX::XMFLOAT4X4 mWorldMatrix{ Library::MatrixHelper::Identity };
		winrt::com_ptr<ID3D11Buffer> mVertexBuffer;
		winrt::com_ptr<ID3D11Buffer> mIndexBuffer;
		std::uint32_t mIndexCount{ 0 };
		float mModelRotationAngle{ 0.0f };
		bool mAnimationEnabled{ true };
		bool mUpdateMaterial{ true };
	};
}


text:
/*
 * Gsl_blas_ddot: test gsl_blas_ddot (vector . vector)
 *
 * Copyright (c) 2012 Jérémie Decock
 *
 * Required: GSL library (libgsl0-dev)
 * Usage: gcc gsl_blas_ddot.c -lgsl -lgslcblas -lm
 *    or: gcc gsl_blas_ddot.c $(pkg-config --libs gsl)
 */

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

int main (void)
{
    gsl_vector * v1 = gsl_vector_alloc(3);
    gsl_vector * v2 = gsl_vector_alloc(3);

    gsl_vector_set_all(v1, 3.0);
    gsl_vector_set_all(v2, 2.0);

    double r;
    gsl_blas_ddot(v1, v2, &r);

    printf("%f\n", r);

    gsl_vector_free(v1);
    gsl_vector_free(v2);

    return 0;
}

max_stars_repo_name: jeremiedecock/snippets

max_issues_repo_name: jeremiedecock/snippets

max_forks_repo_name: jeremiedecock/snippets


text:
#ifndef __DUMMY_H__
#define __DUMMY_H__
#include<stdio.h>
//#define NDEBUG
#include<stdint.h>
#include<omp.h>
#include<float.h>
#include<assert.h>
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <time.h>
#include <math.h>
#include<string.h>
#include<dirent.h>
#include<sys/stat.h>
#include<sys/param.h>
#include<sys/types.h>
#include<signal.h>
//GSL library
#include <gsl/gsl_randist.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_cblas.h>
#include <gsl/gsl_blas.h>
#include<gsl/gsl_matrix.h>
#include <utility>
#include <string>
#include <iostream>
#include<iterator>
#include<map>
#include<vector>
#include <cctype>
#include<algorithm>
using namespace std ;

const uint32_t loop_max=400;
uint32_t loop=0;
const double  epsilon= 0.05; //error associated with user variable !
const double  epsilon_m=    0.000001; //error associated with user variable !
const double epsilon_phi=   0.000001; //error associated with user variable !
const double epsilon_theta= 0.000001; //error associated with user variable !
const double lower_limit_theta= 0.0000000000001;
const double lower_limit_m =    0.0000000000001;
const double lower_limit_phi =    0.00000000001;
clock_t  start;
clock_t  end;
gsl_rng *rand_obj ;

double maxima_phi_i_possible=0.00000000000000000000000000000000000001;
double log_Fx_maxima_phi_i_possible=-DBL_MAX;

//double * diversity;
uint32_t *freq_single_den; //  freqency of theta_w in denominator

typedef struct {
  bool* word_caps; // keeps tag if word at this pos was
  bool* comma_present_before;
  bool* period_present_before;
  vector<string> word_list ;
  uint32_t *int_word_list ; // Its an array to words
  uint32_t  *topic_usage;
  int real_word_count ; //Number of words in document
  int true_label ;
  int sampled_label;
  uint8_t *topic_list ; //Topic associated with each word
  char* document_name ;
  bool* priority_flag ;
  bool  hasPriorityWord;
  int mygamma;
}document ;

document *d_ptr ;

//Number of topics defined in begining
int updated_doc=0 ;

//GIBBS  samples count
//
int gibbs_counter=0;
const int gibbs_sample_count =25 ;
double *M_gibbs  ;
double *THETA_gibbs  ;
double *PHI_gibbs ;
const double  log_ratio=log(1000);

#define topic_count 5  // 1 to 30 !
#define document_count 790
#define total_words_count 10000

double *M ;
double *THETA ;
double *PHI;

//File names

const char* M_file="M.bin";
const char* M_IJ_CW_file="m_ij_cw.bin";
const char* THETA_file="THETA.bin";
const char* PHI_file="PHI.bin";
const char* TOPIC_file="topics_list.bin";

//For keeping last maxima , this will  speed up  search
double *last_maxima_theta_iw;
double *last_maxima_m_ij;
double *last_maxima_phi ;

//Thread info
const  int m_thread_cnt=4;
const  int theta_thread_cnt=8;

//Priors
const double m_diag=1;
const double m_nondiag=1;
const double m_beta=1;

const double theta_k=1;
const double theta_beta=1;
/*
double phi_sigma_square=10000;
*/
double phi_sigma_square=1;
double phi_mu=0;
//FUNCTIONS
//INIT  functions
void init(); //wrapper   for  all the inits
void init_model_parameters() ; //Inits model params
void init_document_params(); //INIT  document structure
void init_variables();
void init_topic_assignment(document* current_doc);
void test_model();
void free_things();

//Topic update for each document
inline void  update_topic_usage(document* current_doc_ptr, const uint8_t  new_topic,const int old_topic  );
void topic_update(document* current_doc_ptr);
void update_TOPIC();
double get_factor12(const uint32_t word_curr,const int word_next,const uint8_t t_prime,const int t_last,const int t_next,document* doc,const uint32_t mypos);
inline double get_factor2_t_i(document* current_doc_ptr, double sum);
double get_current_sum(document* doc,const uint32_t N_d, int sz);
inline void scale_pdf(double* pdf ,const int size,double sum );
inline void make_cdf_from_pdf(double* pdf ,double* cdf , const int size );
void update_m_i_j_count_word(const uint32_t word_curr,const int word_next,const int t_prime,const int t_last,const int t_next,document* doc,const int  mypos);

//update m_i_j
void  update_M();
void  update_m_ij(const int row,const int col,double * sum_den,uint32_t* freq_den,uint32_t*  freq_num);
void get_max_f_m_ij(const uint8_t i, const uint8_t j,double*log_F_max,double* maxima,double*sum_den,uint32_t*freq_den,uint32_t* freq_num,double *maxima_possible,double *value_possible);
double evaluate_m_ij_ll(double x0,const uint8_t i,const uint8_t j,double* sum_den,uint32_t* freq_den,uint32_t* freq_num,double *max, double* val_at_maxi);
void m_compute_double_term_word_stat(const uint8_t i,double* sum_den ,uint32_t* freq_den,uint32_t* freq_num_i);

//UPDATE THETA matrix
void update_THETA();
void update_theta_iw(const uint8_t i, const uint32_t w,double* den_single_term_header,double*  sum_den,uint32_t * freq_den,const  int freq_single_num);
void   get_max_f_theta_iw(const uint8_t i, const uint32_t w,double* log_Fx,double* maxima,const uint32_t freq_single_den_w);
double evaluate_theta_iw_ll(double x0,const uint8_t i,const uint32_t w,const double den_single_term_header,double * sum_den,uint32_t * freq_den,const uint32_t freq_single_num,uint32_t* freq_num,const  uint32_t freq_single_den_w,double *maxima_theta_iw_possible,double *log_Fx_maxima_theta_iw_possible);
double get_den_theta_jw(const uint8_t j , const uint32_t w);
double compute_theta_single_term_denominator(const uint32_t  w );
void theta_compute_double_term_word_stat(const uint32_t w,double* theta_double_term_sum_den,uint32_t* theta_double_term_freq_den,uint32_t* freq_single_num_w,uint32_t* freq_double_num);

//Common codes
double inline  get_number_between_two_numbers(double low, double high);
inline double rand_0_1();
void ctrl_c_handler(int single_num);
void  write_model();
double timetaken(clock_t start,clock_t  end , const char* mesaage);
double gamma_approx_sampler(double r,double log_F_max,double log_f_r,double maxima);
void   guess_lebel(document* current_doc_ptr);
double m_diagonistics(const int i,const int j,double* sum_den,uint32_t* freq_den,uint32_t* freq_num ,double* possible_maxima ,double * value_possible,const char*msg ,double maxima=-1);

double theta_diagonistics(const int i, const int w, const double den_single_term_header,double * sum_den,uint32_t* freq_den,const uint32_t freq_single_num,uint32_t *  freq_num,const uint32_t freq_single_den_w,double *maxima1,double * val_at_max1,const char*msg,double maxima=-1) ;

int linear_scan(double*cdf , const int size,double data);
void load_model_params();
void load_topics();

//update PHI
double  evaluate_phi_i(const double x0,const uint8_t i );
void update_PHI();
void   get_max_f_phi_i(const uint8_t i, double* log_Fx,double* maxima ) ; // calculate maximum of function of log !

void print_topic_corpus();
void print_model_params();
void print_system_config();

void reset_topic_usage();

uint32_t  get_key_3D_m(int x,int y,int z,int x_size,int y_size);
void print_time();
void init_gibbs_params();
void write_model_gibbs();
void do_gibbs_expectation();

const string highlight_tag_start="<:t:>";
const string highlight_tag_end="<:/t:>";
const string comma=",";
const string period=".";
void make_word_id(string word);
std::map<string , int > word_id_map;
std::map<int , string > id_word_map;
string  get_word_from_id(int id);

int  get_id_from_word(string word);
//while we read document we will silently populate  this map for furtehr uses !
void  create_word_id_id_word_map(string word);
//Makes word  of document
string make_lowercase_string(string word);
void make_int_words(document* doc,vector<string>& word_list );
bool doesWordStartsUpper(string data);
bool isThisDocumentUnderLined(string name);

vector<string> dictionary;
const char* dict_file="/data1/rahul/projects/text_mining/branch/punctuation/chris_data/word_list.txt" ;
const char* data="/data1/rahul/projects/text_mining/branch/punctuation/chris_data/AD";
//const char* dict_file="/home/rkumar2/final_test/data/word_list.txt" ;
//const char* data="/home/rkumar2/final_test/data/AE";

void  make_Dictionary();
void inline set_d_2D_m(double *ptr,int row,int col,double value,int row_length);
double inline get_d_2D_m(double* ptr,int row,int col,int row_length);
void  set_i_3D_m(uint32_t *ptr,int x,int y,int z,int  value,int y_size,int z_size);
void inline set_i_2D_m(uint32_t *ptr,int row,int col,int value,int row_length);
uint32_t inline get_i_2D_m(uint32_t* ptr,int row,int col,int row_length);
#endif


text:
#pragma once

#include <functional>
#include <gsl/gsl_util>
#include <vector>

struct Executor;
namespace detail {
Executor*& get_current_executor_holder();
}

struct Executor {
  using WorkQueue = std::vector<std::function<void()>>;

  void add_work(std::function<void()> fn) {
    work_.emplace_back(std::move(fn));
  }

  void execute() {
    const auto cleanup                    = gsl::finally([&] { //
      detail::get_current_executor_holder() = nullptr;
    });
    detail::get_current_executor_holder() = this;

    while (not work_.empty()) {
      auto pending_work = WorkQueue{};
      swap(pending_work, work_);
      for (auto&& fn : pending_work) {
        fn();
      }
    }
  }

  WorkQueue work_;
};

inline Executor& get_current_executor() {
  return *detail::get_current_executor_holder();
}


text:
/* Copyright (c) 2014, Giuseppe Argentieri <giuseppe.argentieri@ts.infn.it>

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

/*
 * =====================================================================================
 *
 *       Filename:  Reg_sc.c
 *
 *    Description:
 *
 *        Version:  1.0
 *        Created:  21/04/2014 12:27:24
 *       Revision:  none
 *       Compiler:  gcc
 *
 *         Author:  Giuseppe Argentieri (ga), giuseppe.argentieri@ts.infn.it
 *   Organization:  Università degli Studi di Trieste
 *
 * =====================================================================================
 */

/*
   Re g_sc is made of 4 integrands: 2 with principal values and 2 without.

   The common factor function is k_func = alpha*k*exp(-k/omega_c)/tanh(b*k/2)

   The poles are in (omega1-O) and (omega1+O)

   The sign of the four fractions

   1/(k+(omega1+O)) , 1/(k+(omega1-O)) , 1/(k-(omega1-O)) , 1/(k-(omega1+O))

   are

   + , - , + , -
*/

#include <stdio.h>

#include "funcs.h"
#include <gsl/gsl_integration.h>

/*
 * ===  FUNCTION  ======================================================================
 *         Name:  k_func
 *  Description:  Common integrand factor
 * =====================================================================================
 */
double k_func ( double k, void* params )
{

	/* Copying parameters */
	struct f_params* pars = (struct f_params*) params ;
	double o_c, b, O, o_1, alpha ;
	assign_p( pars, &o_c, &b, &O, &o_1 ) ;
	alpha = pars->alpha ;

	double temp = alpha*k*exp(-k/o_c) ;
	double val = temp/tanh(b*k/2) ;

	return val;
}		/* -----  end of function k_func  ----- */

/*
 * ===  FUNCTION  ======================================================================
 *         Name:  k_func_1
 *  Description:  First integrand: k_func/(k+(omega1+O))
 * =====================================================================================
 */
double k_func_1 ( double k, void* params )
{
	/* Copying parameters */
	struct f_params* pars = (struct f_params*) params ;
	double o_c, b, O, o_1 ;
	assign_p( pars, &o_c, &b, &O, &o_1 ) ;

	/* Invoking k_func and returning value */
	double v , v1 ;
	v = k_func( k, pars ) ;
	v1 = v/(k+(o_1+O)) ;

	return v1 ;
}		/* -----  end of function k_func_1  ----- */

/*
 * ===  FUNCTION  ======================================================================
 *         Name:  first_int
 *  Description:  Integrating k_func/(k+(omega1+O))
 * =====================================================================================
 */
int first_int ( double* val, double* error, void* params )
{
	/* Copying parameters */
	struct f_params* pars = (struct f_params*) params ;
	double o_c, b, O, o_1 ;
	assign_p( pars, &o_c, &b, &O, &o_1 ) ;

	/* Allocating workspace for integration */
	gsl_integration_workspace* first_int_ws =
		gsl_integration_workspace_alloc(WS_SZ) ;

	gsl_function F ;
	F.function = &k_func_1 ;
	F.params = pars ;

	/* Integration over (0,+Infinity) */
	double res, err ;
	int status = gsl_integration_qagiu( &F, 0, 1e-9, 1e-3, WS_SZ, first_int_ws,
			&res, &err ) ;

	*val = res ; *error = err ;

	gsl_integration_workspace_free(first_int_ws) ;
	return status;
}		/* -----  end of function first_int  ----- */

/*
 *      FUNCTION
 *         Name:  k_func_2
 *  Description:  Second integrand: k_func/(k+omega1-O)
 *
 */
double k_func_2 ( double k , void* params )
{
	/* Copying parameters */
	struct f_params* pars = (struct f_params*) params ;
	double o_c, b, O, o_1 ;
	assign_p( pars, &o_c, &b, &O, &o_1 ) ;

	/* Invoking k_func and returning value */
	double v, v2 ;
	v = k_func( k, pars ) ;
	v2 = v/(k+(o_1-O)) ;

	return v2;
}		/* -----  end of function k_func_2  ----- */

/*
 *      FUNCTION
 *         Name:  second_int
 *  Description:  Integration k_func/(k+(omega1-O))
 *
 */
int second_int ( double* val, double* error, void* params )
{
	/* Copying parameters */
	struct f_params* pars = (struct f_params*) params ;
	double o_c, b, O, o_1 ;
	assign_p( pars, &o_c, &b, &O, &o_1 ) ;

	/* Allocating workspace for integration */
	gsl_integration_workspace* second_int_ws =
		gsl_integration_workspace_alloc(WS_SZ) ;

	gsl_function F;
	F.function = &k_func_2 ;
	F.params = pars ;

	double res, err ;
	int status = gsl_integration_qagiu( &F, 0, 1e-9, 1e-3, WS_SZ, second_int_ws,
			&res, &err ) ;
	*val = res ; *error = err ;

	gsl_integration_workspace_free(second_int_ws) ;

	return status;
}		/* -----  end of function second_int  ----- */

/*
 *      FUNCTION
 *         Name:  k_func_3
 *  Description:  k_func/(k-(omega1-O)) to be integrated on the tail +inf
 *
 */
double k_func_3 ( double k, void* params )
{
	/* Copying parameters */
	struct f_params* pars = (struct f_params*) params ;
	double o_c, b, O, o_1 ;
	assign_p( pars, &o_c, &b, &O, &o_1 ) ;

	double temp , val ;
	temp = k_func( k, pars ) ;
	val = temp/(k-(o_1-O)) ;

	return val;
}		/* -----  end of function k_func_3  ----- */

/*
 *      FUNCTION
 *         Name:  third_int
 *  Description:  Integration of the p.v. k_func/(k-(omega1-O))
 *
 */
int third_int ( double* val, double* err, void* params )
{
	/* Copying parameters */
	struct f_params* pars = (struct f_params*) params ;
	double o_c, b, O, o_1 ;
	assign_p( pars, &o_c, &b, &O, &o_1 ) ;

	/* Integration of the p.v. */
	gsl_function F ;
	F.function = &k_func ;
	F.params = pars ;

	gsl_integration_workspace* third_int_ws =
		gsl_integration_workspace_alloc(WS_SZ) ;

	double v1 , err1 ;
	int status1 = gsl_integration_qawc( &F, (o_1-O)/2, 3*(o_1-O)/2, o_1-O,
			1e-9, 1e-3, WS_SZ, third_int_ws, &v1, &err1 ) ;

	gsl_integration_workspace_free(third_int_ws) ;

	/* Integration to +infinity */
	gsl_function G ;
	G.function = &k_func_3 ;
	G.params = pars ;

	gsl_integration_workspace* third_int_ws_2 =
		gsl_integration_workspace_alloc(WS_SZ) ;

	double v2, err2 ;
	int status2 = gsl_integration_qagiu( &G, 3*(o_1-O)/2, 1e-9, 1e-3, WS_SZ,
			third_int_ws_2, &v2, &err2 ) ;

	gsl_integration_workspace_free(third_int_ws_2) ;

	/* Integration on ( 0, (o_1-O)/2 ) */
	gsl_integration_workspace* third_int_ws_3 =
		gsl_integration_workspace_alloc(WS_SZ) ;

	gsl_function H ;
	H.function = &k_func_3 ;
	H.params = pars ;

	double v3, err3 ;
	int status3 = gsl_integration_qag( &H, 0, (o_1-O)/2, 1e-9, 1e-3, WS_SZ, 6,
			third_int_ws_3, &v3, &err3 ) ;

	gsl_integration_workspace_free(third_int_ws_3) ;

	*val = v1 + v2 + v3 ; *err = err1 + err2 + err3 ;

	return status1 + status2 + status3 ;
}		/* -----  end of function third_int  ----- */

/*
 *      FUNCTION
 *         Name:  k_func_4
 *  Description:  k_func/(k-(omega1+O))
 *
 */
double k_func_4 ( double k, void* params )
{
	/* Copying parameters */
	struct f_params* pars = (struct f_params*) params ;
	double o_c, b, O, o_1 ;
	assign_p( pars, &o_c, &b, &O, &o_1 ) ;

	double temp , val ;
	temp = k_func( k, pars ) ;
	val = temp/(k-(o_1+O)) ;

	return val;
}		/* -----  end of function k_func_4  ----- */

/*
 *      FUNCTION
 *         Name:  fourth_integ
 *  Description:  Integration of the p.v. kfunc/(k-(omega1+O))
 *
 */
int fourth_int ( double* val, double* err, void* params )
{
	/* Copying parameters */
	struct f_params* pars = (struct f_params*) params ;
	double o_c, b, O, o_1 ;
	assign_p( pars, &o_c, &b, &O, &o_1 ) ;

	/* Integration of the p.v. */
	gsl_function F ;
	F.function = &k_func ;
	F.params = pars ;

	gsl_integration_workspace* fourth_int_ws =
		gsl_integration_workspace_alloc(WS_SZ) ;

	double v1 , err1 ;
	int status1 = gsl_integration_qawc( &F, (o_1+O)/2, 3*(o_1+O)/2, o_1+O,
			1e-9, 1e-3, WS_SZ, fourth_int_ws, &v1, &err1 ) ;

	gsl_integration_workspace_free(fourth_int_ws) ;

	/* Integration to +infinity */
	gsl_function G ;
	G.function = &k_func_4 ;
	G.params = pars ;

	gsl_integration_workspace* fourth_int_ws_2 =
		gsl_integration_workspace_alloc(WS_SZ) ;

	double v2, err2 ;
	int status2 = gsl_integration_qagiu( &G, 3*(o_1+O)/2, 1e-9, 1e-3, WS_SZ,
			fourth_int_ws_2, &v2, &err2 ) ;

	gsl_integration_workspace_free(fourth_int_ws_2) ;

	/* Integration on ( 0 , (O+o_1)/2 ) */
	gsl_integration_workspace* fourth_int_ws_3 =
		gsl_integration_workspace_alloc(WS_SZ) ;

	gsl_function K ;
	K.function = &k_func_4 ;
	K.params = pars ;
	double v3, err3 ;

	int status3 = gsl_integration_qag( &K, 0, (O+o_1)/2, 1e-9, 1e-3, WS_SZ,
			GSL_INTEG_GAUSS61, fourth_int_ws_3, &v3, &err3 ) ;

	gsl_integration_workspace_free(fourth_int_ws_3) ;

	*val = v1 + v2 + v3 ; *err = err1 + err2 + err3 ;

	return status1 + status2 + status3;
}		/* -----  end of function fourth_integ  ----- */

/*
 *      FUNCTION
 *         Name:  re_gsc
 *  Description:
 *
 */
int re_gsc ( double* result, double* error, void* params )
{
	/* Copying parameters */
	struct f_params* pars = (struct f_params*) params ;
	double o_c, b, O, o_1 ;
	assign_p( pars, &o_c, &b, &O, &o_1 ) ;

	double v1, v2, v3, v4, err1, err2, err3, err4 ;
	int status = first_int( &v1, &err1, pars ) + second_int( &v2, &err2, pars )
		+ third_int( &v3, &err3, pars ) + fourth_int( &v4, &err4, pars ) ;
	*result = -(v1 - v2 + v3 - v4)/4.0 ;
	*error = (err1 + err2 + err3 + err4)/4.0 ;

	return status;
}		/* -----  end of function re_gsc  ----- */

/*
 *      FUNCTION
 *         Name:  re_gcs
 *  Description:  Re g_cs is the same of Re g_sc but with O and o_1 exchanged.
 *  		  Therefore the signs are
 *
 *  		  + , + , - , -
 *
 */
int re_gcs ( double* result, double* error, void* params )
{
	/* Copying parameters */
	struct f_params* pars = (struct f_params*) params ;
	double o_c, b, O, o_1 ;
	assign_p( pars, &o_c, &b, &O, &o_1 ) ;

	double v1, v2, v3, v4, err1, err2, err3, err4 ;
	int status =  first_int( &v1, &err1, pars ) + second_int( &v2, &err2, pars )
		+ third_int( &v3, &err3, pars ) + fourth_int( &v4, &err4, pars ) ;

	*result = -(v1 + v2 - v3 - v4)/4 ; *error = (err1 + err2 + err3 + err4)/4 ;

	return status;
}		/* -----  end of function re_gcs  ----- */


text:
/*
 * GSLIntegrator.h
 *
 *  Created on: 6 cze 2020
 *      Author: tomhof
 */

#ifndef SRC_GSLINTEGRATOR_H_
#define SRC_GSLINTEGRATOR_H_

#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_monte.h>
#include <gsl/gsl_monte_plain.h>
#include <gsl/gsl_monte_miser.h>
#include <gsl/gsl_monte_vegas.h>

class GSLIntegrator {
private:
	size_t calls = 50000;
public:
	GSLIntegrator();
	virtual ~GSLIntegrator();
    double integrate(double (*f)(double * x_array, size_t dim, void * params), double a1, double b1, double a2, double b2, double* int_err);
    void display_results(char *title, double result, double error);

	size_t getCalls() const {
		return calls;
	}

	void setCalls(size_t calls = 50000) {
		this->calls = calls;
	}
};

#endif /* SRC_GSLINTEGRATOR_H_ */


text:
#pragma once

#include <gsl/gsl_sf_zeta.h>
#include <cmath>

double real_trigamma_cmplxarg(double x, double y, int m, int l){
    //compute the real part of the trigamma function at z = x + iy using
    //my accelerated series formula
    //argument should have positive real part (x > 0)
    //m gives the number of zeta function terms to use
    //l gives the number of terms to use in the residual series
    double result = 0;

    //if x is small, the result will be large and thus, inaccurate.  Use the
    //polygamma shift formula to give a larger value of x for the computation
    if(x < 1000){
        return ((x*x - y*y)/(x*x + y*y)/(x*x + y*y)) + real_trigamma_cmplxarg(x + 1, y, m, l);
    }
    else{
        double phase, ypow = 1.0, xpow, denom;
        //compute finite sum of Hurwitz zeta functions
        for (int i = 1; i < m; ++i){
            phase = 2*(i/2) == i ? -1.0 : 1.0;
            result += phase*(2*i - 1)*ypow*gsl_sf_hzeta(2*i, x);
            ypow = ypow*y*y;
        }
        //compute the infinite sum of residuals
        phase = 2*(m/2) == m ? 1.0 : -1.0;
        for (int n = 0; n < l; ++n){
            xpow = pow(x + n, 2*m);
            denom = ((x + n)*(x + n) + y*y)*((x + n)*(x + n) + y*y);
            result += phase*ypow*((2*m - 1)*y*y + (2*m + 1)*(x + n)*(x + n)) / xpow / denom;
        }
        return result;
    }
}


text:
/* linalg/invtri.c
 *
 * Copyright (C) 2016, 2017, 2018, 2019 Patrick Alken
 *
 * This is free software; you can redistribute it and/or modify it
 * under the terms of the GNU General Public License as published by the
 * Free Software Foundation; either version 3, or (at your option) any
 * later version.
 *
 * This source is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
 * for more details.
 *
 * 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.
 *
 * This module contains code to invert triangular matrices
 */

#include <config.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>

#include "recurse.h"

static int triangular_inverse_L2(CBLAS_UPLO_t Uplo, CBLAS_DIAG_t Diag, gsl_matrix * T);
static int triangular_inverse_L3(CBLAS_UPLO_t Uplo, CBLAS_DIAG_t Diag, gsl_matrix * T);
static int triangular_singular(const gsl_matrix * T);

int
gsl_linalg_tri_invert(CBLAS_UPLO_t Uplo, CBLAS_DIAG_t Diag, gsl_matrix * T)
{
  const size_t N = T->size1;

  if (N != T->size2)
    {
      GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
    }
  else
    {
      int status;

      status = triangular_singular(T);
      if (status)
        return status;

      return triangular_inverse_L3(Uplo, Diag, T);
    }
}

/*
triangular_inverse_L2()
  Invert a triangular matrix T

Inputs: Uplo - CblasUpper or CblasLower
        Diag - unit triangular?
        T    - on output the upper (or lower) part of T
               is replaced by its inverse

Return: success/error

Notes:
1) Based on LAPACK routine DTRTI2 using Level 2 BLAS
*/

static int
triangular_inverse_L2(CBLAS_UPLO_t Uplo, CBLAS_DIAG_t Diag, gsl_matrix * T)
{
  const size_t N = T->size1;

  if (N != T->size2)
    {
      GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
    }
  else
    {
      gsl_matrix_view m;
      gsl_vector_view v;
      size_t i;

      if (Uplo == CblasUpper)
        {
          for (i = 0; i < N; ++i)
            {
              double aii;

              if (Diag == CblasNonUnit)
                {
                  double *Tii = gsl_matrix_ptr(T, i, i);
                  *Tii = 1.0 / *Tii;
                  aii = -(*Tii);
                }
              else
                {
                  aii = -1.0;
                }

              if (i > 0)
                {
                  m = gsl_matrix_submatrix(T, 0, 0, i, i);
                  v = gsl_matrix_subcolumn(T, i, 0, i);

                  gsl_blas_dtrmv(CblasUpper, CblasNoTrans, Diag,
                                 &m.matrix, &v.vector);

                  gsl_blas_dscal(aii, &v.vector);
                }
            } /* for (i = 0; i < N; ++i) */
        }
      else
        {
          for (i = 0; i < N; ++i)
            {
              double ajj;
              size_t j = N - i - 1;

              if (Diag == CblasNonUnit)
                {
                  double *Tjj = gsl_matrix_ptr(T, j, j);
                  *Tjj = 1.0 / *Tjj;
                  ajj = -(*Tjj);
                }
              else
                {
                  ajj = -1.0;
                }

              if (j < N - 1)
                {
                  m = gsl_matrix_submatrix(T, j + 1, j + 1,
                                           N - j - 1, N - j - 1);
                  v = gsl_matrix_subcolumn(T, j, j + 1, N - j - 1);

                  gsl_blas_dtrmv(CblasLower, CblasNoTrans, Diag,
                                 &m.matrix, &v.vector);

                  gsl_blas_dscal(ajj, &v.vector);
                }
            } /* for (i = 0; i < N; ++i) */
        }

      return GSL_SUCCESS;
    }
}

/*
triangular_inverse_L3()
  Invert a triangular matrix T

Inputs: Uplo - CblasUpper or CblasLower
        Diag - unit triangular?
        T    - on output the upper (or lower) part of T
               is replaced by its inverse

Return: success/error

Notes:
1) Based on ReLAPACK routine DTRTRI using Level 3 BLAS
*/

static int
triangular_inverse_L3(CBLAS_UPLO_t Uplo, CBLAS_DIAG_t Diag, gsl_matrix * T)
{
  const size_t N = T->size1;

  if (N != T->size2)
    {
      GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
    }
  else if (N <= CROSSOVER_INVTRI)
    {
      /* use Level 2 BLAS code */
      return triangular_inverse_L2(Uplo, Diag, T);
    }
  else
    {
      /*
       * partition matrix:
       *
       * T11 T12
       * T21 T22
       *
       * where T11 is N1-by-N1
       */
      int status;
      const size_t N1 = GSL_LINALG_SPLIT(N);
      const size_t N2 = N - N1;
      gsl_matrix_view T11 = gsl_matrix_submatrix(T, 0, 0, N1, N1);
      gsl_matrix_view T12 = gsl_matrix_submatrix(T, 0, N1, N1, N2);
      gsl_matrix_view T21 = gsl_matrix_submatrix(T, N1, 0, N2, N1);
      gsl_matrix_view T22 = gsl_matrix_submatrix(T, N1, N1, N2, N2);

      /* recursion on T11 */
      status = triangular_inverse_L3(Uplo, Diag, &T11.matrix);
      if (status)
        return status;

      if (Uplo == CblasLower)
        {
          /* T21 = - T21 * T11 */
          gsl_blas_dtrmm(CblasRight, Uplo, CblasNoTrans, Diag, -1.0, &T11.matrix, &T21.matrix);

          /* T21 = T22 * T21^{-1} */
          gsl_blas_dtrsm(CblasLeft, Uplo, CblasNoTrans, Diag, 1.0, &T22.matrix, &T21.matrix);
        }
      else
        {
          /* T12 = - T11 * T12 */
          gsl_blas_dtrmm(CblasLeft, Uplo, CblasNoTrans, Diag, -1.0, &T11.matrix, &T12.matrix);

          /* T12 = T12 * T22^{-1} */
          gsl_blas_dtrsm(CblasRight, Uplo, CblasNoTrans, Diag, 1.0, &T22.matrix, &T12.matrix);
        }

      /* recursion on T22 */
      status = triangular_inverse_L3(Uplo, Diag, &T22.matrix);
      if (status)
        return status;

      return GSL_SUCCESS;
    }
}

static int
triangular_singular(const gsl_matrix * T)
{
  size_t i;

  for (i = 0; i < T->size1; ++i)
    {
      double Tii = gsl_matrix_get(T, i, i);
      if (Tii == 0.0)
        return GSL_ESING;
    }

  return GSL_SUCCESS;
}

#ifndef GSL_DISABLE_DEPRECATED

int
gsl_linalg_tri_upper_invert(gsl_matrix * T)
{
  int status = triangular_singular(T);
  if (status)
    return status;

  return triangular_inverse_L3(CblasUpper, CblasNonUnit, T);
}

int
gsl_linalg_tri_lower_invert(gsl_matrix * T)
{
  int status = triangular_singular(T);
  if (status)
    return status;

  return triangular_inverse_L3(CblasLower, CblasNonUnit, T);
}

int
gsl_linalg_tri_upper_unit_invert(gsl_matrix * T)
{
  int status = triangular_singular(T);
  if (status)
    return status;

  return triangular_inverse_L3(CblasUpper, CblasUnit, T);
}

int
gsl_linalg_tri_lower_unit_invert(gsl_matrix * T)
{
  int status = triangular_singular(T);
  if (status)
    return status;

  return triangular_inverse_L3(CblasLower, CblasUnit, T);
}

#endif

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
/*
* C functions for the NIG marginal distribution.
*
* This file is part of Fieldosophy, a toolkit for random fields.
*
* Copyright (C) 2021 Anders Gunnar Felix Hildeman <fieldosophySPDEC@gmail.com>
*
* This Source Code is subject to the terms of the BSD 3-Clause License.
* If a copy of the license was not distributed with this file, you can obtain one at https://opensource.org/licenses/BSD-3-Clause.
*
*
* Author: Anders Gunnar Felix Hildeman
* Date: 2020-04
*/

#include <stdio.h>
#include <string.h>
#include <math.h>
#include <omp.h>
#include <gsl/gsl_sf_bessel.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>

extern "C"
{

    // Make sure that alpha tilde is not too large
    void NIG_filterParams( double *pAlphaTilde, double *pBetaTilde)
    {
        if (*pAlphaTilde > 100.0d)
        {
            *pBetaTilde = *pBetaTilde * ( 100.0d / (*pAlphaTilde) );
            *pAlphaTilde = 100.0d;
        }

        return;
    }

    // Modified Bessel function of the second kind
    double Kv(const int v, const double x)
    {
        double lOut = 0.0d;
        int status;

        gsl_sf_result result;

        switch(v)
        {
        case 0:
            status = gsl_sf_bessel_K0_e(x, &result);
            break;
        case 1:
            status = gsl_sf_bessel_K1_e(x, &result);
            break;
        default:
            status = gsl_sf_bessel_Kn_e(v, x, &result);
            break;
        }

        // If some kind of error
        if (status != 0)
        {
            //printf("Error code: %d in Kv", status);
            lOut = GSL_NAN;
        }
        else
        {
            lOut = result.val;
        }

        return lOut;
    }

    // Logarithm of modified Bessel function of second kind
    double lKv(const int v, const double x)
    {
        double lOut = Kv( v, x );
        // Logarithmize
        lOut = log( lOut );

        return lOut;
    }

    // logarithm of PDF for NIG distribution
    unsigned int NIG_lpdf( double * const pX, const unsigned int pN, const double alphaTilde, const double betaTilde, const double mu, const double delta )
    {
        // turn off gsl error handler
        gsl_set_error_handler_off();

        // Loop through each element in vector
        #pragma omp parallel for
        for ( unsigned int iter = 0; iter < pN; iter++)
        {
            // Get normalized data value
            const double lCur = (pX[iter] - mu) / delta;
            // Compute q function
            const double lQ = sqrt( (1.0d + lCur * lCur) );
            // Compute bessel function
            const double lLogK = lKv(1, alphaTilde * lQ);

            // compute log pdf
            double lCurOut = lLogK;
            lCurOut -= log(lQ);
            lCurOut += betaTilde * lCur;
            // Add constants
            lCurOut += log( alphaTilde / M_PI ) + sqrt( alphaTilde * alphaTilde - betaTilde * betaTilde );
            lCurOut -= log(delta);

            // Store full log pdf
            pX[iter] = lCurOut;
        }

        return pN;
    }

    // Compute log-likelihood
    double NIG_logLikelihood( double * lOut, const double * const pX, const unsigned int pN, const double alphaTilde, const double betaTilde, const double mu, const double delta )
    {
        // Copy content
        memcpy( lOut, pX, sizeof(double) * pN );

        // Compute log-pdf
        NIG_lpdf( lOut, pN, alphaTilde, betaTilde, mu, delta );
        // Sum
        double lSum = 0.0d;
        for (unsigned int iter = 0; iter < pN; iter++)
        {
            lSum += lOut[iter];
        }

        return lSum;
    }

    // EM-algorithm for finding MLE of NIG
    unsigned int NIG_MLE( double * const pLogLik, double * const pX, const unsigned int pN,
        double * const pAlphaTilde, double * const pBetaTilde, double * const pMu, double * const pDelta,
        const double pXBar, const unsigned int pM, const double pTol )
    {

        // turn off gsl error handler
        gsl_set_error_handler_off();

        // Init
        double lAlphaTilde = *pAlphaTilde;
        double lBetaTilde = *pBetaTilde;
        double lMu = *pMu;
        double lDelta = *pDelta;
        // Filter tilde parameters
        NIG_filterParams( &lAlphaTilde, &lBetaTilde);

        // Preallocate arrays
        double * lTempArray1 = (double*) malloc(sizeof(double) * pN);
        double * lTempArray2 = (double*) malloc(sizeof(double) * pN);
        if (lTempArray1 == NULL || lTempArray2 == NULL)
            return 1;

        // Compute log-likelihood
        *pLogLik = NIG_logLikelihood( lTempArray1, pX, pN, lAlphaTilde, lBetaTilde, lMu, lDelta );

        // If likelihood is bad
        if (isnan(*pLogLik))
            return 1;

        // Iterate EM-algorithm
        for (unsigned int iter = 0; iter < pM; iter++)
        {
            // --- E-step ---

            double lSBar = 0.0d;
            double lWBar = 0.0d;
            #pragma omp parallel for reduction(+:lSBar) reduction(+:lWBar)
            for (unsigned int iter2 = 0; iter2 < pN; iter2++)
            {
                const double lPhiSqrt = sqrt( 1.0d + ( (pX[iter2]-lMu)/lDelta ) * ( (pX[iter2]-lMu)/lDelta )  );
                // compute besselk function values
                const double lK0 = Kv( 0, lAlphaTilde * lPhiSqrt );
                const double lK1 = Kv( 1, lAlphaTilde * lPhiSqrt );
                const double lK2 = Kv( 2, lAlphaTilde * lPhiSqrt );
                // Compute pseudo values
                lTempArray1[iter2] = ( lDelta * lPhiSqrt * lK0 ) / ( (lAlphaTilde * lK1) / lDelta ); // S
                lTempArray2[iter2] = ( (lAlphaTilde * lK2) / lDelta ) / ( lDelta * lPhiSqrt * lK1 ); // W
                lSBar += lTempArray1[iter2];
                lWBar += lTempArray2[iter2];
            }
            lSBar /= ((double)pN);
            lWBar /= ((double)pN);

            // --- M-step ---

            // Compute delta
            lDelta = sqrt( 1.0d / ( lWBar - (1.0d / lSBar) ) );
            double lGamma = lDelta / lSBar;
            // Compute beta
            double lBeta = 0.0d;
            for (unsigned int iter2 = 0; iter2 < pN; iter2++)
            {
                lBeta += pX[iter2] * lTempArray2[iter2];
            }
            lBeta = lBeta - pXBar * lWBar * ((double)pN);
            lBeta /= ( (double)pN - lSBar * lWBar * ((double)pN) );
            lBetaTilde = lBeta * lDelta;
            // Compute alpha
            double lAlpha = sqrt( lGamma * lGamma + lBeta * lBeta );
            lAlphaTilde = lAlpha * lDelta;
            // Filter tilde values
            NIG_filterParams( &lAlphaTilde, &lBetaTilde);
            lGamma = sqrt( lAlphaTilde * lAlphaTilde - lBetaTilde * lBetaTilde ) / lDelta;
            // Compute mu
            lMu = pXBar - lBetaTilde / lGamma;

            // Compute loglik
            const double llLikNew = NIG_logLikelihood( lTempArray1, pX, pN, lAlphaTilde, lBetaTilde, lMu, lDelta );

            // If likelihood is bad
            if (isnan(llLikNew))
                return(iter+2);    // Mark that likelihood went bad during iterations

            // Compute likelihood difference
            const double lLLikDiff = llLikNew - *pLogLik;

            // If new likelihood is higher
            if (lLLikDiff > 0)
            {
                // Set new log-likelihood as old one
                *pLogLik = llLikNew;
                // Update output parameters
                *pAlphaTilde = lAlphaTilde;
                *pBetaTilde = lBetaTilde;
                *pMu = lMu;
                *pDelta = lDelta;

                // If should break due to tolerance
                if ( lLLikDiff < pTol )
                    break;
            }

        }   // end of EM-iterations

        // Free allocated arrays
        free(lTempArray1);
        free(lTempArray2);

         // Return no problem
         return( 0 );
    }

}

max_stars_repo_name: andyGFHill/fieldosophy

max_issues_repo_name: andyGFHill/fieldosophy

max_forks_repo_name: andyGFHill/fieldosophy


text:
#include <stdio.h>
#include <stdarg.h>
#include <string.h>
#include <math.h>
#include <gbpLib.h>
#include <gbpRNG.h>
#include <gbpMCMC.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_fit.h>
#include <gsl/gsl_interp.h>

void init_MCMC_arrays(MCMC_info *MCMC) {
    int           i_DS;
    MCMC_DS_info *current_DS;
    MCMC_DS_info *next_DS;

    if(MCMC->n_M == NULL) {
        SID_log("Initializing MCMC target arrays for %d dataset(s)...", SID_LOG_OPEN, MCMC->n_DS);

        // Initialize chain arrays
        MCMC->n_M       = (int *)SID_malloc(sizeof(int) * MCMC->n_DS);
        MCMC->M_new     = (double **)SID_malloc(sizeof(double *) * MCMC->n_DS);
        MCMC->M_last    = (double **)SID_malloc(sizeof(double *) * MCMC->n_DS);
        MCMC->n_M_total = 0;
        i_DS            = 0;
        current_DS      = MCMC->DS;
        while(current_DS != NULL) {
            next_DS            = current_DS->next;
            MCMC->n_M[i_DS]    = current_DS->n_M;
            MCMC->M_new[i_DS]  = (double *)SID_malloc(sizeof(double) * MCMC->n_M[i_DS]);
            MCMC->M_last[i_DS] = (double *)SID_malloc(sizeof(double) * MCMC->n_M[i_DS]);
            MCMC->n_M_total += MCMC->n_M[i_DS];
            current_DS = next_DS;
            i_DS++;
        }

        // Initialize results arrays
        MCMC->P_min   = (double *)SID_malloc(sizeof(double) * MCMC->n_P);
        MCMC->P_max   = (double *)SID_malloc(sizeof(double) * MCMC->n_P);
        MCMC->P_avg   = (double *)SID_malloc(sizeof(double) * MCMC->n_P);
        MCMC->dP_avg  = (double *)SID_malloc(sizeof(double) * MCMC->n_P);
        MCMC->P_best  = (double *)SID_malloc(sizeof(double) * MCMC->n_P);
        MCMC->P_peak  = (double *)SID_malloc(sizeof(double) * MCMC->n_P);
        MCMC->P_lo_68 = (double *)SID_malloc(sizeof(double) * MCMC->n_P);
        MCMC->P_hi_68 = (double *)SID_malloc(sizeof(double) * MCMC->n_P);
        MCMC->P_lo_95 = (double *)SID_malloc(sizeof(double) * MCMC->n_P);
        MCMC->P_hi_95 = (double *)SID_malloc(sizeof(double) * MCMC->n_P);

        MCMC->ln_likelihood_DS      = (double *)SID_malloc(sizeof(double) * MCMC->n_DS);
        MCMC->ln_likelihood_DS_best = (double *)SID_malloc(sizeof(double) * MCMC->n_DS);
        MCMC->ln_likelihood_DS_peak = (double *)SID_malloc(sizeof(double) * MCMC->n_DS);
        MCMC->n_DoF_DS              = (int *)SID_malloc(sizeof(int) * MCMC->n_DS);
        MCMC->n_DoF_DS_best         = (int *)SID_malloc(sizeof(int) * MCMC->n_DS);
        MCMC->n_DoF_DS_peak         = (int *)SID_malloc(sizeof(int) * MCMC->n_DS);

        // If MCMC_MODE_MINIMIZE_IO is on, we need to allocate buffers to store the chain
        if(SID_CHECK_BITFIELD_SWITCH(MCMC->mode, MCMC_MODE_MINIMIZE_IO)) {
            MCMC->flag_success_buffer      = (char *)SID_malloc(sizeof(char) * MCMC->n_iterations * MCMC->n_avg);
            MCMC->ln_likelihood_new_buffer = (double *)SID_malloc(sizeof(double) * MCMC->n_iterations * MCMC->n_avg);
            MCMC->P_new_buffer             = (double *)SID_malloc(sizeof(double) * MCMC->n_iterations * MCMC->n_avg * MCMC->n_P);
            MCMC->M_new_buffer             = (double *)SID_malloc(sizeof(double) * MCMC->n_iterations * MCMC->n_avg * MCMC->n_M_total);
        }

        SID_log("Done.", SID_LOG_CLOSE);
    }
}

max_stars_repo_name: gbpoole/gbpCode

max_issues_repo_name: gbpoole/gbpCode

max_forks_repo_name: gbpoole/gbpCode