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

/* Author:  G. Jungman */

/* Calculate the matrix exponential, following
 * Moler + Van Loan, SIAM Rev. 20, 801 (1978).
 */

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_mode.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_blas.h>

#include "gsl_linalg.h"

/* store one of the suggested choices for the
 * Taylor series / square  method from Moler + VanLoan
 */
struct moler_vanloan_optimal_suggestion
{
  int k;
  int j;
};
typedef  struct moler_vanloan_optimal_suggestion  mvl_suggestion_t;

/* table from Moler and Van Loan
 * mvl_tab[gsl_mode_t][matrix_norm_group]
 */
static mvl_suggestion_t mvl_tab[3][6] =
{
  /* double precision */
  {
    { 5, 1 }, { 5, 4 }, { 7, 5 }, { 9, 7 }, { 10, 10 }, { 8, 14 }
  },

  /* single precision */
  {
    { 2, 1 }, { 4, 0 }, { 7, 1 }, { 6, 5 }, { 5, 9 }, { 7, 11 }
  },

  /* approx precision */
  {
    { 1, 0 }, { 3, 0 }, { 5, 1 }, { 4, 5 }, { 4, 8 }, { 2, 11 }
  }
};

inline
static double
sup_norm(const gsl_matrix * A)
{
  double min, max;
  gsl_matrix_minmax(A, &min, &max);
  return GSL_MAX_DBL(fabs(min), fabs(max));
}

static
mvl_suggestion_t
obtain_suggestion(const gsl_matrix * A, gsl_mode_t mode)
{
  const unsigned int mode_prec = GSL_MODE_PREC(mode);
  const double norm_A = sup_norm(A);
  if(norm_A < 0.01) return mvl_tab[mode_prec][0];
  else if(norm_A < 0.1) return mvl_tab[mode_prec][1];
  else if(norm_A < 1.0) return mvl_tab[mode_prec][2];
  else if(norm_A < 10.0) return mvl_tab[mode_prec][3];
  else if(norm_A < 100.0) return mvl_tab[mode_prec][4];
  else if(norm_A < 1000.0) return mvl_tab[mode_prec][5];
  else
  {
    /* outside the table we simply increase the number
     * of squarings, bringing the reduced matrix into
     * the range of the table; this is obviously suboptimal,
     * but that is the price paid for not having those extra
     * table entries
     */
    const double extra = log(1.01*norm_A/1000.0) / M_LN2;
    const int extra_i = (unsigned int) ceil(extra);
    mvl_suggestion_t s = mvl_tab[mode][5];
    s.j += extra_i;
    return s;
  }
}

/* use series representation to calculate matrix exponential;
 * this is used for small matrices; we use the sup_norm
 * to measure the size of the terms in the expansion
 */
static void
matrix_exp_series(
  const gsl_matrix * B,
  gsl_matrix * eB,
  int number_of_terms
  )
{
  int count;
  gsl_matrix * temp = gsl_matrix_calloc(B->size1, B->size2);

  /* init the Horner polynomial evaluation,
   * eB = 1 + B/number_of_terms; we use
   * eB to collect the partial results
   */
  gsl_matrix_memcpy(eB, B);
  gsl_matrix_scale(eB, 1.0/number_of_terms);
  gsl_matrix_add_diagonal(eB, 1.0);
  for(count = number_of_terms-1; count >= 1; --count)
  {
    /*  mult_temp = 1 + B eB / count  */
    gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, B, eB, 0.0, temp);
    gsl_matrix_scale(temp, 1.0/count);
    gsl_matrix_add_diagonal(temp, 1.0);

    /*  transfer partial result out of temp */
    gsl_matrix_memcpy(eB, temp);
  }

  /* now eB holds the full result; we're done */
  gsl_matrix_free(temp);
}

int
gsl_linalg_exponential_ss(
  const gsl_matrix * A,
  gsl_matrix * eA,
  gsl_mode_t mode
  )
{
  if(A->size1 != A->size2)
  {
    GSL_ERROR("cannot exponentiate a non-square matrix", GSL_ENOTSQR);
  }
  else if(A->size1 != eA->size1 || A->size2 != eA->size2)
  {
    GSL_ERROR("exponential of matrix must have same dimension as matrix", GSL_EBADLEN);
  }
  else
  {
    int i;
    const mvl_suggestion_t sugg = obtain_suggestion(A, mode);
    const double divisor = exp(M_LN2 * sugg.j);

    gsl_matrix * reduced_A = gsl_matrix_alloc(A->size1, A->size2);

    /*  decrease A by the calculated divisor  */
    gsl_matrix_memcpy(reduced_A, A);
    gsl_matrix_scale(reduced_A, 1.0/divisor);

    /*  calculate exp of reduced matrix; store in eA as temp  */
    matrix_exp_series(reduced_A, eA, sugg.k);

    /*  square repeatedly; use reduced_A for scratch */
    for(i = 0; i < sugg.j; ++i)
    {
      gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, eA, eA, 0.0, reduced_A);
      gsl_matrix_memcpy(eA, reduced_A);
    }

    gsl_matrix_free(reduced_A);

    return GSL_SUCCESS;
  }
}


text:
/* Cholesky Decomposition
 *
 * Copyright (C) 2000  Thomas Walter
 *
 * 3 May 2000: Modified for GSL by Brian Gough
 *
 * 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 2, 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.
 */

/*
 * Cholesky decomposition of a symmetrix positive definite matrix.
 * This is useful to solve the matrix arising in
 *    periodic cubic splines
 *    approximating splines
 *
 * This algorthm does:
 *   A = L * L'
 * with
 *   L  := lower left triangle matrix
 *   L' := the transposed form of L.
 *
 */

#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>

int
gsl_linalg_cholesky_decomp (gsl_matrix * A)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if (M != N)
    {
      GSL_ERROR("cholesky decomposition requires square matrix", GSL_ENOTSQR);
    }
  else
    {
      size_t i,j,k;
      int status = 0;

      /* Do the first 2 rows explicitly.  It is simple, and faster.  And
       * one can return if the matrix has only 1 or 2 rows.
       */

      double A_00 = gsl_matrix_get (A, 0, 0);

      double L_00 = sqrt(A_00);

      if (A_00 <= 0)
        {
          status = GSL_EDOM ;
        }

      gsl_matrix_set (A, 0, 0, L_00);

      if (M > 1)
        {
          double A_10 = gsl_matrix_get (A, 1, 0);
          double A_11 = gsl_matrix_get (A, 1, 1);

          double L_10 = A_10 / L_00;
          double diag = A_11 - L_10 * L_10;
          double L_11 = sqrt(diag);

          if (diag <= 0)
            {
              status = GSL_EDOM;
            }

          gsl_matrix_set (A, 1, 0, L_10);
          gsl_matrix_set (A, 1, 1, L_11);
        }

      for (k = 2; k < M; k++)
        {
          double A_kk = gsl_matrix_get (A, k, k);

          for (i = 0; i < k; i++)
            {
              double sum = 0;

              double A_ki = gsl_matrix_get (A, k, i);
              double A_ii = gsl_matrix_get (A, i, i);

              gsl_vector_view ci = gsl_matrix_row (A, i);
              gsl_vector_view ck = gsl_matrix_row (A, k);

              if (i > 0) {
                gsl_vector_view di = gsl_vector_subvector(&ci.vector, 0, i);
                gsl_vector_view dk = gsl_vector_subvector(&ck.vector, 0, i);

                gsl_blas_ddot (&di.vector, &dk.vector, &sum);
              }

              A_ki = (A_ki - sum) / A_ii;
              gsl_matrix_set (A, k, i, A_ki);
            }

          {
            gsl_vector_view ck = gsl_matrix_row (A, k);
            gsl_vector_view dk = gsl_vector_subvector (&ck.vector, 0, k);

            double sum = gsl_blas_dnrm2 (&dk.vector);
            double diag = A_kk - sum * sum;

            double L_kk = sqrt(diag);

            if (diag <= 0)
              {
                status = GSL_EDOM;
              }

            gsl_matrix_set (A, k, k, L_kk);
          }
        }

      /* Now copy the transposed lower triangle to the upper triangle,
       * the diagonal is common.
       */

      for (i = 1; i < M; i++)
        {
          for (j = 0; j < i; j++)
            {
              double A_ij = gsl_matrix_get (A, i, j);
              gsl_matrix_set (A, j, i, A_ij);
            }
        }

      if (status == GSL_EDOM)
        {
          GSL_ERROR ("matrix must be positive definite", GSL_EDOM);
        }

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_cholesky_solve (const gsl_matrix * LLT,
                           const gsl_vector * b,
                           gsl_vector * x)
{
  if (LLT->size1 != LLT->size2)
    {
      GSL_ERROR ("cholesky matrix must be square", GSL_ENOTSQR);
    }
  else if (LLT->size1 != b->size)
    {
      GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
    }
  else if (LLT->size2 != x->size)
    {
      GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
    }
  else
    {
      /* Copy x <- b */

      gsl_vector_memcpy (x, b);

      /* Solve for c using forward-substitution, L c = b */

      gsl_blas_dtrsv (CblasLower, CblasNoTrans, CblasNonUnit, LLT, x);

      /* Perform back-substitution, U x = c */

      gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, LLT, x);

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_cholesky_svx (const gsl_matrix * LLT,
                         gsl_vector * x)
{
  if (LLT->size1 != LLT->size2)
    {
      GSL_ERROR ("cholesky matrix must be square", GSL_ENOTSQR);
    }
  else if (LLT->size2 != x->size)
    {
      GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
    }
  else
    {
      /* Solve for c using forward-substitution, L c = b */

      gsl_blas_dtrsv (CblasLower, CblasNoTrans, CblasNonUnit, LLT, x);

      /* Perform back-substitution, U x = c */

      gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, LLT, x);

      return GSL_SUCCESS;
    }
}

ext: c


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

/* Calculate the exponential integral function by invoking the GNU scientific
   library */
/* Only intended to be used with R; not the most elegant integration but it
   may do for now */

/* One argument, the value (double) at which to evaluate Ei(x) */

double x, ei; /* argument to exponential integral, value */
char *program_name; /* name program is invoked under, for errors */

int main(int argc, char* argv[]) {
  void usage(void);	/* Warn users about proper usage */

  program_name = argv[0];
  if (argc != 2) {
    usage();
  }
  x = atof(&argv[1][0]);
  ei = gsl_sf_expint_E1(x);
  printf("%.18e\n",ei);
  return(0);
}

void usage(void) {
  (void) fprintf(stderr, "Usage is %s [floating-point argument]\n", program_name);
  exit(8);
}

max_stars_repo_name: jguerber/spatialwarnings

max_stars_repo_licenses: Unlicense MIT

max_issues_repo_name: jguerber/spatialwarnings

max_issues_repo_licenses: Unlicense MIT

max_forks_repo_name: jguerber/spatialwarnings

max_forks_repo_licenses: Unlicense MIT


text:
/*
 * Copyright 2021 The DAPHNE Consortium
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

#pragma once

#include <runtime/local/context/DaphneContext.h>
#include <runtime/local/datastructures/DataObjectFactory.h>
#include <runtime/local/datastructures/DenseMatrix.h>

#include <cblas.h>
#include <lapacke.h>

#include <cassert>
#include <cstddef>

// ****************************************************************************
// Struct for partial template specialization
// ****************************************************************************

template<class DTRes, class DTLhs, class DTRhs>
struct Solve {
    static void apply(DTRes *& res, const DTLhs * lhs, const DTRhs * rhs, DCTX(ctx)) = delete;
};

// ****************************************************************************
// Convenience function
// ****************************************************************************

template<class DTRes, class DTLhs, class DTRhs>
void solve(DTRes *& res, const DTLhs * lhs, const DTRhs * rhs, DCTX(ctx)) {
    Solve<DTRes, DTLhs, DTRhs>::apply(res, lhs, rhs, ctx);
}

// ****************************************************************************
// (Partial) template specializations for different data/value types
// ****************************************************************************

// ----------------------------------------------------------------------------
// DenseMatrix <- DenseMatrix, DenseMatrix
// ----------------------------------------------------------------------------

template<>
struct Solve<DenseMatrix<float>, DenseMatrix<float>, DenseMatrix<float>> {
    static void apply(DenseMatrix<float> *& res, const DenseMatrix<float> * lhs, const DenseMatrix<float> * rhs, DCTX(ctx)) {
        const auto nr1 = static_cast<int>(lhs->getNumRows());
        const auto nc1 = static_cast<int>(lhs->getNumCols());
        const auto nc2 = static_cast<int>(rhs->getNumCols());
        assert((nr1 == static_cast<int>(rhs->getNumRows())) && "#rows of lhs and #rows of rhs must be the same");
        assert((nr1 == nc1) && "#rows and #cols of lhs must be the same");
        assert((static_cast<int>(lhs->getRowSkip()) == nc1) && "#cols of lhs must match row skip");
        assert((nc2==1) && "#cols of rhs must be 1");

        if(res == nullptr)
            res = DataObjectFactory::create<DenseMatrix<float>>(nr1, nc2, false);

        // solve system of equations via LU decomposition
        int ipiv[nr1];       // permutation indexes
        float work[nr1*nc1]; // LU factorization of gesv
        memcpy(work, lhs->getValues(), nr1*nc1*sizeof(float));         //for in-place A
        memcpy(res->getValues(), rhs->getValues(), nr1*sizeof(float)); //for in-place b-out
        [[maybe_unused]] int info = LAPACKE_sgesv(LAPACK_ROW_MAJOR, nr1, 1, work, nc1, ipiv, res->getValues(), 1);
        assert((info<=0) && "A factor Ui is exactly singular, so the solution could not be computed");
    }
};

template<>
struct Solve<DenseMatrix<double>, DenseMatrix<double>, DenseMatrix<double>> {
    static void apply(DenseMatrix<double> *& res, const DenseMatrix<double> * lhs, const DenseMatrix<double> * rhs, DCTX(ctx)) {
        const auto nr1 = static_cast<int>(lhs->getNumRows());
        const auto nc1 = static_cast<int>(lhs->getNumCols());
        const auto nc2 = static_cast<int>(rhs->getNumCols());
        assert((nr1 == static_cast<int>(rhs->getNumRows())) && "#rows of lhs and #rows of rhs must be the same");
        assert((nr1 == nc1) && "#rows and #cols of lhs must be the same");
        assert((static_cast<int>(lhs->getRowSkip()) == nc1) && "#cols of lhs must match row skip");
        assert((nc2==1) && "#cols of rhs must be 1");

        if(res == nullptr)
            res = DataObjectFactory::create<DenseMatrix<double>>(nr1, nc2, false);

        // solve system of equations via LU decomposition
        int ipiv[nr1];       // permutation indexes
        double work[nr1*nc1]; // LU factorization of gesv
        memcpy(work, lhs->getValues(), nr1*nc1*sizeof(double));         //for in-place A
        memcpy(res->getValues(), rhs->getValues(), nr1*sizeof(double)); //for in-place b-out
        [[maybe_unused]] int info = LAPACKE_dgesv(LAPACK_ROW_MAJOR, nr1, 1, work, nc1, ipiv, res->getValues(), 1);
        assert((info<=0) && "A factor Ui is exactly singular, so the solution could not be computed");
    }
};

ext: h

max_stars_repo_name: daphne-eu/daphne

max_issues_repo_name: daphne-eu/daphne

max_forks_repo_name: daphne-eu/daphne


text:
#ifndef _EM_IO_H_
#define _EM_IO_H_ 1

#include <petsc.h>

struct EMContext;

PetscErrorCode read_mdl(const char *, EMContext *);
PetscErrorCode read_emd(const char *, EMContext *);

PetscErrorCode save_mesh(EMContext *, const char *, int);
PetscErrorCode save_rsp(EMContext *, const char *);

#endif

max_stars_repo_name: emfem/emfem

max_issues_repo_name: emfem/emfem

max_forks_repo_name: emfem/emfem


text:
#pragma once

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

namespace Library
{
	class ProxyModel;
}

namespace Rendering
{
	class PointLightMaterial;

	class PointLightDemo final : public Library::DrawableGameComponent
	{
	public:
		PointLightDemo(Library::Game& game, const std::shared_ptr<Library::Camera>& camera);
		PointLightDemo(const PointLightDemo&) = delete;
		PointLightDemo(PointLightDemo&&) = default;
		PointLightDemo& operator=(const PointLightDemo&) = default;
		PointLightDemo& operator=(PointLightDemo&&) = default;
		~PointLightDemo();

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

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

		float PointLightIntensity() const;
		void SetPointLightIntensity(float intensity);

		const DirectX::XMFLOAT3& LightPosition() const;
		const DirectX::XMVECTOR LightPositionVector() const;
		void SetLightPosition(const DirectX::XMFLOAT3& position);
		void SetLightPosition(DirectX::FXMVECTOR position);

		float LightRadius() const;
		void SetLightRadius(float radius);

		float SpecularIntensity() const;
		void SetSpecularIntensity(float intensity);

		float SpecularPower() const;
		void SetSpecularPower(float power);

		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<PointLightMaterial> mMaterial;
		DirectX::XMFLOAT4X4 mWorldMatrix{ Library::MatrixHelper::Identity };
		winrt::com_ptr<ID3D11Buffer> mVertexBuffer;
		winrt::com_ptr<ID3D11Buffer> mIndexBuffer;
		std::uint32_t mIndexCount{ 0 };
		Library::PointLight mPointLight;
		std::unique_ptr<Library::ProxyModel> mProxyModel;
		float mModelRotationAngle{ 0.0f };
		bool mAnimationEnabled{ true };
		bool mUpdateMaterial{ true };
	};
}


text:
#ifndef GBNET_MODELBASE
#define GBNET_MODELBASE

#include <algorithm>
#include <cmath>
#include <csignal>
#include <time.h>

#include <gsl/gsl_rstat.h>

#include "GraphBase.h"

namespace gbn
{
    typedef std::pair<std::string, double> varid_stat_pair_t;
    typedef std::vector<varid_stat_pair_t> posterior_stat_vector_t;
    typedef std::vector<std::pair<std::string, double>> gelman_rubin_vector_t;

    class ModelBase
    {
        private:
        protected:
            int current_signal = 0;
        public:
            network_t network;
            evidence_dict_t evidence = evidence_dict_t();
            prior_active_tf_set_t active_tf_set = prior_active_tf_set_t();
            unsigned int n_graphs = 3;
            bool noise_listen_children = false;

            std::vector<GraphBase *> graphs;

            ModelBase(unsigned int = 3);
            virtual ~ModelBase();

            gelman_rubin_vector_t get_gelman_rubin();
            double get_max_gelman_rubin();

            void sample(unsigned int = 50000, unsigned int = 30);
            void sample_n(unsigned int);
            void burn_stats();

            void print_stats();

            posterior_stat_vector_t get_posterior_means(std::string);
            posterior_stat_vector_t get_posterior_sdevs(std::string);

            void set_signal(int);
    };
}

#endif

max_stars_repo_name: umbibio/gbnet

max_issues_repo_name: umbibio/gbnet

max_forks_repo_name: umbibio/gbnet


text:
#pragma once

#include <halley/maths/vector2.h>
#include <halley/time/halleytime.h>
#include <halley/core/graphics/sprite/sprite.h>
#include <gsl/gsl>
#include <thread>
#include <mutex>
#include <atomic>

namespace Halley
{
	class RenderContext;
	class TextureDescriptor;
	class VideoAPI;
	class AudioAPI;
	class StreamingAudioClip;
	class IAudioHandle;
	class TextureRenderTarget;

	enum class MoviePlayerState
	{
		Uninitialised,
		Loading,
		Paused,
		StartingToPlay,
		Playing,
		Finished
	};

	enum class MoviePlayerStreamType
	{
		Unknown,
		Video,
		Audio,
		Other
	};

	class MoviePlayerStream
	{
	public:
		MoviePlayerStreamType type = MoviePlayerStreamType::Unknown;
		bool playing = true;
		bool eof = false;
	};

	struct PendingFrame
	{
		std::shared_ptr<Texture> texture;
		Time time;
	};

	struct MoviePlayerAliveFlag
	{
		bool isAlive = true;
		mutable std::mutex mutex;
	};

	class MoviePlayer
	{
	public:
		MoviePlayer(VideoAPI& video, AudioAPI& audio);
		virtual ~MoviePlayer();

		void play();
		void pause();
		void reset();
		void stop();

		virtual bool hasError() const;

		virtual void update(Time t);
		void render(Resources& resources, RenderContext& rc);
		Sprite getSprite(Resources& resources);

		MoviePlayerState getState() const;
		Vector2i getSize() const;

	protected:
		virtual void requestVideoFrame() = 0;
		virtual void requestAudioFrame() = 0;
		virtual void onReset();
		virtual void onStartPlay();
		virtual void waitForVideoInfo();
		virtual bool needsYV12Conversion() const;
		virtual bool shouldRecycleTextures() const;
		virtual void onDoneUsingTexture(std::shared_ptr<Texture> texture);
		virtual Rect4i getCropRect() const;

		void setVideoSize(Vector2i size);

		VideoAPI& getVideoAPI() const;
		AudioAPI& getAudioAPI() const;

		void onVideoFrameAvailable(Time time, TextureDescriptor&& descriptor);
		void onVideoFrameAvailable(Time time, std::shared_ptr<Texture> texture);
		void onAudioFrameAvailable(Time time, gsl::span<const short> samples);
		void onAudioFrameAvailable(Time time, gsl::span<const float> samples);

		std::shared_ptr<MoviePlayerAliveFlag> getAliveFlag() const;

		std::vector<MoviePlayerStream> streams;
		std::list<std::shared_ptr<Texture>> recycleTexture;
		std::list<PendingFrame> pendingFrames;
		int maxVideoFrames;
		int maxAudioSamples;

	private:
		VideoAPI& video;
		AudioAPI& audio;

		MoviePlayerState state = MoviePlayerState::Uninitialised;

		Vector2i videoSize;
		std::shared_ptr<Texture> currentTexture;
		std::shared_ptr<TextureRenderTarget> renderTarget;
		std::shared_ptr<Texture> renderTexture;

		std::shared_ptr<StreamingAudioClip> streamingClip;
		std::shared_ptr<IAudioHandle> audioHandle;

		std::atomic<bool> threadRunning;
		std::atomic<bool> threadAborted;
		std::thread workerThread;
		std::shared_ptr<MoviePlayerAliveFlag> aliveFlag;

		Time time = 0;

		void startThread();
		void stopThread();
		void threadEntry();

		virtual bool useCustomThreads() const;
		virtual void stopCustomThreads();

		bool needsMoreVideoFrames() const;
		bool needsMoreAudioFrames() const;
	};
}

max_stars_repo_name: sunhay/halley

max_issues_repo_name: sunhay/halley

max_forks_repo_name: sunhay/halley


text:
/* This file is auto-generated from weight_functions.h.src */
#ifndef DOUBLE_PREC
#define DOUBLE_PREC
#endif
// # -*- mode: c -*-
#pragma once

#include "defs.h"

#include <stdint.h>

#include <gsl/gsl_matrix_double.h>
#include <gsl/gsl_linalg.h>

// Info about a particle pair that we will pass to the weight function
typedef struct
{
    int64_t nprojbins;
    int64_t nsbins;
    double sqr_s;
    double *supp_sqr; //proper way to have array in struct?

} proj_struct_double;

//typedef double (*weight_func_t_double)(const pair_struct_double*);

//////////////////////////////////
// Projection functions
//////////////////////////////////

static inline void tophat_double(const proj_struct_double *proj, double *u){

    int ins = -1;
    for(int p=0;p<proj->nsbins;p++){
        u[p] = 0;
        if (proj->sqr_s >= proj->supp_sqr[p] && proj->sqr_s < proj->supp_sqr[p+1]){
            ins = p;
        }
    }
    if (ins>=0 && ins<proj->nprojbins){
        u[ins] = 1.0;
    }
}

//////////////////////////////////
// Utility functions
//////////////////////////////////

static inline void compute_amplitudes(int nprojbins, int nd1, int nd2, int nr1, int nr2,
            double *dd, double *dr, double *rd, double *rr, double *qq, double *amps){

    printf("Computing amps\n");
    printf("qq:\n");
    for(int i=0;i<nprojbins;i++){
        for (int j=0; j<nprojbins; j++){
            printf(" %f", qq[i*nprojbins+j]);
        }
        printf("\n");
    }

    // Computer numerator of estimator
    double numerator[nprojbins];
    double qqnorm[nprojbins*nprojbins];
    for (int i=0; i<nprojbins; i++){
        double ddnorm = dd[i]/((double)nd1*(double)nd2);
        double drnorm = dr[i]/((double)nd1*(double)nr2);
        double rdnorm = rd[i]/((double)nr1*(double)nd2);
        double rrnorm = rr[i]/((double)nr1*(double)nr2);
        numerator[i] = ddnorm - drnorm - rdnorm + rrnorm;
        for (int j=0; j<nprojbins; j++){
            printf("%f", qq[i*nprojbins+j]);
            qqnorm[i*nprojbins+j] = qq[i*nprojbins+j]/((double)nr1*(double)nr2);
        }
    }

    printf("qqnorm:\n");
    for(int i=0;i<nprojbins;i++){
        for (int j=0; j<nprojbins; j++){
            printf(" %f", qqnorm[i*nprojbins+j]);
        }
        printf("\n");
    }

	int s;
	// Define all the used matrices
	gsl_matrix *qq_mat = gsl_matrix_alloc(nprojbins, nprojbins);
	gsl_matrix *qq_mat_inv = gsl_matrix_alloc(nprojbins, nprojbins);
	gsl_permutation *perm = gsl_permutation_alloc(nprojbins);
	// Fill the matrix m
	for (int i=0; i<nprojbins; i++){
        for (int j=0; j<nprojbins; j++){
            gsl_matrix_set(qq_mat, i, j, qqnorm[i*nprojbins+j]);
        }
    }
	// Make LU decomposition of matrix m
	gsl_linalg_LU_decomp(qq_mat, perm, &s);
	// Invert the matrix m
	gsl_linalg_LU_invert(qq_mat, perm, qq_mat_inv);

    printf("qqinv:\n");
    for(int i=0;i<nprojbins;i++){
        for (int j=0; j<nprojbins; j++){
            printf(" %f", gsl_matrix_get(qq_mat_inv, i, j));
        }
        printf("\n");
    }

    // Take inner product of qqinv * numerator, get amplitude vector
    // TODO: double check
    for (int i=0; i<nprojbins; i++){
        double aval = 0;
        for (int j=0; j<nprojbins; j++){
            aval += gsl_matrix_get(qq_mat_inv, i, j) * numerator[j];
        }
        amps[i] = aval;
    }
}

static inline void evaluate_xi(int nprojbins, double *amps, int nsvals, double *svals,
                                                int nsbins, double *sbins, double *xi){

    // will need to generalize, projbins won't always be rela   ted to sbins
    double supp_sqr[nsbins];
    //plus 1 because one more edge than number of bins
    for (int i=0; i<nsbins+1; i++){
        supp_sqr[i] = pow(sbins[i], 2);
    }

    // ns: number of s values at which to evaluate xi
    for (int i=0; i<nsvals; i++){
        //get basis function u for given value of s
        double u[nprojbins];
        double sqr_s = pow(svals[i], 2);
        proj_struct_double projdata = {.nprojbins=nprojbins, .nsbins=nsbins, .sqr_s=sqr_s, .supp_sqr=supp_sqr};
        tophat_double(&projdata, u);

        //multiply u by the amplitudes to get xi in that s bin (xi is vector of length ns_toeval)
        double xival = 0;
        for (int j=0; j<nprojbins; j++){
            xival += amps[j]*u[j];

        }
        xi[i] = xival;
    }

}

/* Gives a pointer to the weight function for the given weighting method
 * and instruction set.
 */
//static inline weight_func_t_double get_weight_func_by_method_double(const weight_method_t method){
//    switch(method){
//        case PAIR_PRODUCT:
//            return &pair_product_double;
//        default:
//        case NONE:
//            return NULL;
//    }
//}

max_stars_repo_name: kstoreyf/foo

max_issues_repo_name: kstoreyf/foo

max_forks_repo_name: kstoreyf/foo


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

// uncomment to check that local allocs actually get freed
// #define AA_ALLOC_STACK_MAX 0

#include <stdlib.h>
#include <cblas.h>

#include "amino.h"
#include "amino/mat.h"
#include "amino/mat_internal.h"

#define VEC_LEN(X) ((int)(X->len))
#define MAT_ROWS(X) ((int)(X->rows))
#define MAT_COLS(X) ((int)(X->cols))

/******************/
/* Error Handling */
/******************/

static void
s_err_default( const char *message )
{
    fprintf(stderr, "AMINO ERROR: %s\n",message);
    abort();
    exit(EXIT_FAILURE);
}

static aa_la_err_fun  *s_err_fun = s_err_default;

void
aa_la_set_err( aa_la_err_fun *fun )
{
    s_err_fun = fun;
}

AA_API void
aa_la_err( const char *message ) {
    s_err_fun(message);
}

AA_API void
aa_la_fail_size( size_t a, size_t b )
{
        const size_t size = 256;
        char buf[size];
        snprintf(buf,size, "Mismatched sizes: %lu != %lu\n", a, b);
        aa_la_err(buf);
}

/****************/
/* Construction */
/****************/

AA_API void
aa_dvec_view( struct aa_dvec *vec, size_t len, double *data, size_t inc )
{
    *vec = AA_DVEC_INIT(len,data,inc);
}

AA_API void
aa_dvec_slice( const struct aa_dvec *src,
               size_t start,
               size_t stop,
               size_t step,
               struct aa_dvec *dst )
{
    if( stop > src->len || stop < start ) {
        aa_la_err("Slice out-of-bounds\n");
    }
    *dst = AA_DVEC_INIT( (stop - start) / step,
                         src->data + src->inc*start,
                         src->inc*step );
}

AA_API void
aa_dmat_row_vec( const struct aa_dmat *src, size_t row, struct aa_dvec *dst )
{
    if( row >= src->rows ) {
        aa_la_err("Row vector out-of-bounds\n");
    }
    aa_dvec_view(dst,src->cols, src->data + row, src->ld);
}

AA_API void
aa_dmat_col_vec( const struct aa_dmat *src, size_t col, struct aa_dvec *dst )
{
    if( col >= src->cols ) {
        aa_la_err("Row vector out-of-bounds\n");
    }
    aa_dvec_view(dst, src->rows, src->data + col*src->ld, 1);
}

AA_API void
aa_dmat_diag_vec( const struct aa_dmat *src, struct aa_dvec *dst )
{
    aa_dvec_view(dst, AA_MIN(src->rows, src->cols), src->data, src->ld + 1);
}

AA_API void
aa_dmat_view( struct aa_dmat *mat, size_t rows, size_t cols,
              double *data, size_t ld )
{
    *mat = AA_DMAT_INIT(rows,cols,data,ld);
}

AA_API void
aa_dmat_view_block(  struct aa_dmat *dst,
                     const struct aa_dmat *src,
                     size_t row_start, size_t col_start,
                     size_t rows, size_t cols )
{
    aa_dmat_block( src,
                   row_start, col_start,
                   row_start + rows, col_start + cols,
                   dst );
}

AA_API void
aa_dmat_block( const struct aa_dmat *src,
               size_t row_start, size_t col_start,
               size_t row_end, size_t col_end,
               struct aa_dmat *dst )
{
    size_t m = src->rows, n = src->cols;
    if( row_start >= row_end ||
        row_start >= m ||
        row_end > m ||
        col_start >= col_end ||
        col_start >= n ||
        col_end > n )
    {
        aa_la_err("Block out-of-bounds\n");
    }

    aa_dmat_view( dst,
                  row_end - row_start, col_end - col_start,
                  src->data + col_start*src->ld + row_start, src->ld );
}

AA_API struct aa_dvec *
aa_dvec_malloc( size_t len ) {
    struct aa_dvec *r;
    const size_t s_desc = sizeof(*r);
    const size_t s_elem = sizeof(r->data[0]);
    const size_t pad =  s_elem - (s_desc % s_elem);
    const size_t size = s_desc + pad + len * s_elem;

    const size_t off = s_desc + pad;

    char *ptr = (char*)malloc( size );

    r = (struct aa_dvec*)ptr;
    aa_dvec_view(r, len, (double*)(ptr+off), 1);
    return r;
}

AA_API struct aa_dvec *
aa_dvec_alloc( struct aa_mem_region *reg, size_t len ) {
    struct aa_dvec *r;
    const size_t s_desc = sizeof(*r);
    const size_t s_elem = sizeof(r->data[0]);
    const size_t pad =  s_elem - (s_desc % s_elem);
    char *ptr = (char*)aa_mem_region_alloc(reg, s_desc + pad + len * s_elem );

    r = (struct aa_dvec*)ptr;
    aa_dvec_view(r, len, (double*)(ptr+s_desc+pad), 1);
    return r;
}

AA_API struct aa_dvec *
aa_dvec_dup( struct aa_mem_region *reg, const struct aa_dvec *src)
{
    struct aa_dvec *dst = aa_dvec_alloc(reg,src->len);
    aa_dvec_copy(src,dst);
    return dst;
}

AA_API struct aa_dmat *
aa_dmat_dup( struct aa_mem_region *reg, const struct aa_dmat *src)
{
    struct aa_dmat *dst = aa_dmat_alloc(reg,src->rows,src->cols);
    aa_dmat_copy(src,dst);
    return dst;
}

AA_API struct aa_dmat *
aa_dmat_malloc( size_t rows, size_t cols )
{
    struct aa_dmat *r;
    const size_t s_desc = sizeof(*r);
    const size_t s_elem = sizeof(r->data[0]);
    const size_t pad =  s_elem - (s_desc % s_elem);
    char *ptr = (char*)malloc( s_desc + pad + rows*cols * s_elem );

    r = (struct aa_dmat*)ptr;
    aa_dmat_view(r, rows, cols, (double*)(ptr+s_desc+pad), rows);

    return r;
}

AA_API struct aa_dmat *
aa_dmat_alloc( struct aa_mem_region *reg, size_t rows, size_t cols )
{
    struct aa_dmat *r;
    const size_t s_desc = sizeof(*r);
    const size_t s_elem = sizeof(r->data[0]);
    const size_t pad =  s_elem - (s_desc % s_elem);
    size_t size = s_desc + pad + rows*cols * s_elem ;

    char *ptr = (char*)aa_mem_region_alloc(reg, size);

    r = (struct aa_dmat*)ptr;
    aa_dmat_view( r, rows, cols, (double*)(ptr+s_desc+pad), rows );

    return r;
}

AA_API void
aa_dmat_set( struct aa_dmat *A, double alpha, double beta )
{
    int mi   = (int)(A->rows);
    int ni   = (int)(A->cols);
    int ldai = (int)(A->ld);

    dlaset_( "G", &mi, &ni,
             &alpha, &beta,
             A->data, &ldai );

}

AA_API void
aa_dvec_set( struct aa_dvec *vec, double alpha )
{
    double *end = vec->data + vec->len*vec->inc;
    for( double *x = vec->data; x < end; x += vec->inc ) {
        *x = alpha;
    }
}

void
aa_dvec_zero( struct aa_dvec *vec )
{
    aa_dvec_set(vec,0);
}

AA_API void
aa_dmat_zero( struct aa_dmat *mat )
{
    aa_dmat_set(mat, 0, 0);
}

/* Level 1 BLAS */
AA_API void
aa_dvec_swap( struct aa_dvec *x, struct aa_dvec *y )
{
    aa_la_check_size(x->len, y->len);
    cblas_dswap( VEC_LEN(x), AA_VEC_ARGS(x), AA_VEC_ARGS(y) );
}

AA_API void
aa_dvec_scal( double a, struct aa_dvec *x )
{
    cblas_dscal( VEC_LEN(x), a, AA_VEC_ARGS(x) );
}

static void s_inc( size_t n, double alpha, double *x, size_t inc ) {
    for( double *end = x + n*inc; x < end; x += inc ) {
        *x += alpha;
    }
}

AA_API void
aa_dvec_inc( double alpha, struct aa_dvec *x )
{
    s_inc( x->len, alpha, x->data, x->inc );
}

AA_API void
aa_dvec_copy( const struct aa_dvec *x, struct aa_dvec *y )
{
    aa_la_check_size(x->len, y->len);
    cblas_dcopy( VEC_LEN(x), AA_VEC_ARGS(x), AA_VEC_ARGS(y) );
}

AA_API void
aa_dvec_axpy( double a, const struct aa_dvec *x, struct aa_dvec *y )
{
    aa_la_check_size(x->len, y->len);
    cblas_daxpy( VEC_LEN(x), a, AA_VEC_ARGS(x), AA_VEC_ARGS(y) );
}

AA_API double
aa_dvec_dot( const struct aa_dvec *x, struct aa_dvec *y )
{
    aa_la_check_size(x->len, y->len);
    return cblas_ddot( VEC_LEN(x), AA_VEC_ARGS(x), AA_VEC_ARGS(y) );
}

AA_API double
aa_dvec_nrm2( const struct aa_dvec *x )
{
    return cblas_dnrm2( VEC_LEN(x), AA_VEC_ARGS(x) );
}

/* Level 2 BLAS */
AA_API void
aa_dmat_gemv( CBLAS_TRANSPOSE trans,
              double alpha, const struct aa_dmat *A,
              const struct aa_dvec *x,
              double beta, struct aa_dvec *y )
{
    if( CblasTrans == trans ) {
        aa_la_check_size( A->rows, x->len );
        aa_la_check_size( A->cols, y->len );
    } else {
        aa_la_check_size( A->rows, y->len );
        aa_la_check_size( A->cols, x->len );
    }

    cblas_dgemv( CblasColMajor, trans,
                 MAT_ROWS(A), MAT_COLS(A),
                 alpha, AA_MAT_ARGS(A),
                 AA_VEC_ARGS(x),
                 beta, AA_VEC_ARGS(y) );

}

/* Level 3 BLAS */
AA_API void
aa_dmat_gemm( CBLAS_TRANSPOSE transA, CBLAS_TRANSPOSE transB,
              double alpha, const struct aa_dmat *A,
              const struct aa_dmat *B,
              double beta, struct aa_dmat *C )
{
    aa_la_check_size( A->rows, C->rows );
    aa_la_check_size( A->cols, B->rows );
    aa_la_check_size( B->cols, C->cols );

    assert( A->rows <= A->ld );
    assert( B->rows <= B->ld );
    assert( C->rows <= C->ld );

    cblas_dgemm( CblasColMajor,
                 transA, transB,
                 MAT_ROWS(A), MAT_COLS(B), MAT_COLS(A),
                 alpha, AA_MAT_ARGS(A),
                 AA_MAT_ARGS(B),
                 beta, AA_MAT_ARGS(C) );
}

/* LAPACK */

AA_API void
aa_dmat_lacpy( const char uplo[1],
               const struct aa_dmat *A,
               struct aa_dmat *B )
{

    aa_la_check_size(A->rows,B->rows);
    aa_la_check_size(A->cols,B->cols);
    int mi   = (int)(A->rows);
    int ni   = (int)(A->cols);
    int ldai = (int)(A->ld);
    int ldbi = (int)(B->ld);

    dlacpy_(uplo, &mi, &ni,
            A->data, &ldai,
            B->data, &ldbi);
}

AA_API void
aa_dmat_copy(  const struct aa_dmat *A, struct aa_dmat *B)
{
    aa_dmat_lacpy("G",A,B);
}

/* Matrix Functions */

static double s_ssd( size_t n,
                     double a,
                     double *x, size_t incx,
                     double *y, size_t incy )
{
    for( double *e = x + n*incx; x < e; x += incx, y+=incy ) {
        double d = *x - *y;
        a += d*d;
    }
    return a;
}

AA_API double
aa_dvec_ssd( const struct aa_dvec *x, const struct aa_dvec *y)
{
    aa_la_check_size( x->len, y->len );
    return s_ssd( x->len, 0,
                  x->data, x->inc,
                  y->data, y->inc );
}

AA_API double
aa_dmat_ssd( const struct aa_dmat *A, const struct aa_dmat *B)
{
    size_t m = A->rows, n = A->cols;
    aa_la_check_size( m, B->rows );
    aa_la_check_size( n, B->cols );
    double a = 0;
    for( double *Ac=A->data, *Bc=B->data, *ec=A->data + n*A->ld;
         Ac < ec;
         Ac+=A->ld, Bc+=B->ld )
    {
        a = s_ssd( m, a, Ac, 1, Bc, 1 );
    }
    return a;
}

AA_API double
aa_dmat_nrm2( const struct aa_dmat *A )
{
    int m=(int)A->rows, n=(int)A->cols, ld=(int)A->ld;
    return dlange_("F", &m, &n, A->data, &ld, NULL );
}

AA_API void
aa_dmat_scal( struct aa_dmat *x, double alpha )
{
    int m = (int)x->rows, n=(int)x->cols, ld=(int)x->ld;
    double cfrom=1;
    int info;
    dlascl_("G", NULL, NULL,
            &cfrom, &alpha,
            &m, &n, x->data, &ld,
            &info);
}

AA_API void
aa_dmat_inc( struct aa_dmat *A, double alpha )
{
    size_t m=A->rows, n=A->cols, ld=A->ld;
    double *x = A->data;

    if( m == n ) {
        s_inc(m*n, alpha, x, 1);
    } else {
        for( double *e = x + n*ld; x < e; x+=ld ) {
            s_inc(m, alpha, x, 1);
        }
    }

}

AA_API void
aa_dmat_axpy( double alpha, const struct aa_dmat *X, struct aa_dmat *Y)
{
    size_t m=X->rows, n=X->cols, ldX=X->ld, ldY=Y->ld;
    double *x=X->data, *y=Y->data;

    aa_la_check_size(m, Y->rows);
    aa_la_check_size(n, Y->cols);

    if( m == ldX && m == ldY ) {
        cblas_daxpy( (int)(m*n), alpha, x,1, y,1 );
    } else {
        int mi = (int)m;
        for( double *e = x + n*ldX; x < e; x+=ldX, y+=ldY ) {
            cblas_daxpy( mi, alpha, x,1, y,1 );
        }
    }
}

void
aa_dmat_trans( const struct aa_dmat *A, struct aa_dmat *B)
{
    const size_t m   = A->rows;
    const size_t n   = A->cols;
    const size_t lda = A->ld;
    const size_t ldb = B->ld;

    aa_la_check_size(m, B->cols);
    aa_la_check_size(n, B->rows);

    for ( double *Acol = A->data, *Brow=B->data, *Be=B->data + n;
          Brow < Be;
          Brow++, Acol+=lda )
    {
        cblas_dcopy( (int)m, Acol, 1, Brow, (int)ldb );
    }
}

AA_API int
aa_dmat_inv( struct aa_dmat *A )
{
    aa_la_check_size(A->rows,A->cols);
    int info;

    int *ipiv = (int*)
        aa_mem_region_local_alloc(sizeof(int)*A->rows);

    // LU-factor
    info = aa_cla_dgetrf( MAT_ROWS(A), MAT_COLS(A), AA_MAT_ARGS(A), ipiv );

    int lwork = -1;
    while(1) {
        double *work = (double*)
            aa_mem_region_local_tmpalloc( sizeof(double)*
                                      (size_t)(lwork < 0 ? 1 : lwork) );
        aa_cla_dgetri( MAT_ROWS(A), AA_MAT_ARGS(A), ipiv, work, lwork );
        if( lwork > 0 ) break;
        assert( -1 == lwork );
        lwork = (int)work[0];
    }

    aa_mem_region_local_pop(ipiv);

    return info;
}

int s_svd_helper ( const struct aa_dmat *A,
                   struct aa_mem_region *reg,
                   size_t *pkmin, size_t *pkmax,
                   double **pU, double **pVt, double **pS )
{

    size_t m = A->rows;
    size_t n = A->cols;

    // find  min/max dimensions
    if( m < n ) {
        *pkmin = m;
        *pkmax = n;
    } else {
        *pkmin = n;
        *pkmax = m;
    }

    // This method uses the SVD
    double *W  = (double*)aa_mem_region_alloc( reg, sizeof(double) * (m*m + n*n + *pkmin) );
    *pU  = W;        // size m*m
    *pVt = *pU + m*m; // size n*n
    *pS  = *pVt + n*n; // size min(m,n)

    // A = U S V^T
    aa_la_d_svd( m,n,
                 A->data, A->ld,
                 *pU, m, *pS, *pVt, n );

    return 0;
}

int
aa_dmat_pinv( const struct aa_dmat *A, double tol, struct aa_dmat *As )
{

    aa_la_check_size(A->rows, As->cols);
    aa_la_check_size(A->cols, As->rows);

    struct aa_mem_region *reg =  aa_mem_region_local_get();
    void *ptrtop = aa_mem_region_ptr(reg);

    size_t kmin, kmax;
    double *U, *Vt, *S;
    s_svd_helper( A, reg,
                  &kmin, &kmax, &U, &Vt, &S );

    size_t m = A->rows;
    const int mi = (int)(A->rows);
    const int ni = (int)(A->cols);

    if( tol < 0 ) {
        tol = (double)kmax * S[0] * DBL_EPSILON;
    }

    // \sum 1/s_i * v_i * u_i^T
    aa_dmat_zero(As);
    double *Asd = As->data;
    for( size_t i = 0; i < kmin && S[i] > tol; i ++ ) {
        cblas_dger( CblasColMajor, ni, mi, 1/S[i],
                    Vt + i, ni,
                    U + m*i, 1,
                    Asd, ni
            );
    }

    aa_mem_region_pop( reg, ptrtop );

    return 0;

    /* struct aa_mem_region *reg = aa_mem_region_local_get(); */
    /* struct aa_dmat *Ap = aa_dmat_alloc(reg,m,n); */

    /* // B = AA^T */
    /* double *B; */
    /* int ldb; */
    /* if( m <= n ) { */
    /*     /\* Use A_star as workspace when it's big enough *\/ */
    /*     B = As->data; */
    /*     ldb = (int)(As->ld); */
    /* } else { */
    /*     B = AA_MEM_REGION_NEW_N( reg, double, m*m ); */
    /*     ldb = (int)m; */
    /* } */

    /* aa_la_dlacpy( "A", A, Ap ); */

    /* // B is symmetric.  Only compute the upper half. */
    /* cblas_dsyrk( CblasColMajor, CblasUpper, CblasNoTrans, */
    /*              MAT_ROWS(Ap), MAT_COLS(Ap), */
    /*              1, AA_MAT_ARGS(Ap), */
    /*              0, B, ldb ); */

    /* // B += kI */
    /* /\* for( size_t i = 0; i < m*(size_t)ldb; i += ((size_t)ldb+1) ) *\/ */
    /* /\*     B[i] += k; *\/ */

    /* /\* Solve via Cholesky Decomp *\/ */
    /* /\* B^T (A^*)^T = A    and     B = B^T (Hermitian) *\/ */

    /* aa_cla_dposv( 'U', (int)m, (int)n, */
    /*               B, ldb, */
    /*               AA_MAT_ARGS(Ap) ); */

    /* aa_dmat_trans( Ap, As ); */

    /* aa_mem_region_pop( reg, Ap ); */

}

int
aa_dmat_dpinv( const struct aa_dmat *A, double k, struct aa_dmat *As)
{
    aa_la_check_size(A->rows, As->cols);
    aa_la_check_size(A->cols, As->rows);

    // TODO: Try DGESV for non-positive-definite matrices

    size_t m = A->rows;
    size_t n = A->cols;

    struct aa_mem_region *reg = aa_mem_region_local_get();
    void *ptrtop = aa_mem_region_ptr(reg);
    int r = -1;

    /* As = inv(A'*A - k*I) * A' A'*inv(A*A' - k*I) */
    if( m <= n ) {
        /*
         * (A^*) = A^T*inv(A*A^T - k*I)
         * (A^*) = A^T*inv(B)
         * (A^*)*B = A^T
         * B^T (A^*)^T = A    and     B = B^T (Hermitian)
         *
         */

        /* Use A_star as workspace for B */
        double *B = As->data;
        int ldb = (int)(As->ld);

        /* Ap = A */
        struct aa_dmat *Ap = aa_dmat_alloc(reg,m,n);
        aa_dmat_lacpy( "A", A, Ap );

        // B is symmetric.  Only compute the upper half.
        // B = A * A'
        cblas_dsyrk( CblasColMajor, CblasUpper, CblasNoTrans,
                     MAT_ROWS(Ap), MAT_COLS(Ap),
                     1, AA_MAT_ARGS(Ap),
                     0, B, ldb );

        /* B += kI */
        for( double *x = B, *e = B+(int)m*ldb; x < e; x+=ldb+1 )
            *x += k;

        /* B^T (A^*)^T = A    and     B = B^T (Hermitian) */

        /* Solve via Cholesky (positive definite) */
        r = aa_cla_dposv( 'U', (int)m, (int)n,
                           B, ldb,
                           AA_MAT_ARGS(Ap) );

        /* Solve via LU */
        /* int *ipiv = AA_MEM_REGION_NEW_N(reg,int,m); */
        /* r = aa_la_d_sysv( "U", m, n, */
        /*                   B, (size_t)ldb, */
        /*                   ipiv, */
        /*                   Ap->data, Ap->ld ); */

        aa_dmat_trans( Ap, As );

    } else {

        /*
         * (A^*) = inv(A^T*A - k*I) * A^T
         * (A^*) = inv(B) * A^T
         * B*(A^*) = A^T
         */

        /* Use A_star as workspace for B */
        double *B = AA_MEM_REGION_NEW_N(reg,double,n*n);
        int ldb = (int)(n);

        /* As = A^T */
        aa_dmat_trans( A, As );

        // B is symmetric.  Only compute the upper half.
        // B = A' * A = As * As'
        cblas_dsyrk( CblasColMajor, CblasUpper, CblasNoTrans,
                     MAT_ROWS(As), MAT_COLS(As),
                     1, AA_MAT_ARGS(As),
                     0, B, ldb );

        /* B += kI */
        //for( size_t i = 0; i < n*(size_t)ldb; i += ((size_t)ldb+1) )
        for( double *x = B, *e = B+n*n; x < e; x+=ldb+1 )
            *x += k;

        /* B * (A^*)^T = A    and     B = B^T (Hermitian) */

        /* Solve via Cholesky (positive definite) */
        r = aa_cla_dposv( 'U', (int)n, (int)m,
                           B, ldb,
                           AA_MAT_ARGS(As) );

        /* Solve via LU */
        /* int *ipiv = AA_MEM_REGION_NEW_N(reg,int,n); */
        /* r = aa_la_d_sysv( "U", n, m, */
        /*                   B, (size_t)ldb, */
        /*                   ipiv, */
        /*                   As->data, As->ld ); */

    }

    aa_mem_region_pop( reg, ptrtop );

    return r;
}

int
aa_dmat_dzdpinv(  const struct aa_dmat *A, double s_min, struct aa_dmat *As)
{

    aa_la_check_size(A->rows, As->cols);
    aa_la_check_size(A->cols, As->rows);

    struct aa_mem_region *reg =  aa_mem_region_local_get();
    void *ptrtop = aa_mem_region_ptr(reg);

    size_t kmin, kmax;
    double *U, *Vt, *S;
    s_svd_helper( A, reg,
                  &kmin, &kmax, &U, &Vt, &S );

    size_t m = A->rows;
    const int mi = (int)(A->rows);
    const int ni = (int)(A->cols);

    // \sum s_i/(s_i**2+k) * v_i * u_i^T
    aa_dmat_zero(As);
    double *Asd = As->data;
    size_t i = 0;

    // Undamped parts
    for( ; i < kmin && S[i] >= s_min; i ++ ) {
        cblas_dger( CblasColMajor, ni, mi, 1/S[i],
                    Vt + i, ni,
                    U + m*i, 1,
                    Asd, ni
            );
    }

    // Damped parts
    double s2 = s_min*s_min;
    for( ; i < kmin; i ++ ) {
        cblas_dger( CblasColMajor, ni, mi, S[i]/s2,
                    Vt + i, ni,
                    U + m*i, 1,
                    Asd, ni
            );
    }

    aa_mem_region_pop( reg, ptrtop );

    return 0;

}

max_stars_repo_name: dyalab/amino

max_issues_repo_name: dyalab/amino

max_forks_repo_name: dyalab/amino


text:
#include <gsl/gsl_math.h>
#include "gsl_cblas.h"
#include "cblas.h"

float
cblas_sdsdot (const int N, const float alpha, const float *X, const int incX,
	      const float *Y, const int incY)
{
#define INIT_VAL  alpha
#define ACC_TYPE  double
#define BASE float
#include "source_dot_r.h"
#undef ACC_TYPE
#undef BASE
#undef INIT_VAL
}


text:
/*
 * Copyright (c) 2016-2021 lymastee, All rights reserved.
 * Contact: lymastee@hotmail.com
 *
 * This file is part of the gslib project.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 * SOFTWARE.
 */

#pragma once

#ifndef image_a41c29de_bd9b_4bbd_9012_6dc956f524b3_h
#define image_a41c29de_bd9b_4bbd_9012_6dc956f524b3_h

#include <gslib/type.h>
#include <gslib/string.h>
#include <ariel/dirty.h>

__ariel_begin__

class image
{
public:
    enum image_format
    {
        fmt_gray,                   /* gray8 */
        fmt_rgba,                   /* rgba8888 */
    };

    friend class imageio;

public:
    image();
    ~image() { destroy(); }
    bool is_valid() const;
    image_format get_format() const { return _format; }
    int get_depth() const { return _depth; }
    bool create(image_format fmt, int w, int h);
    void destroy();
    void init(const color& cr);
    void enable_alpha_channel(bool b) { _is_alpha_channel_valid = b; }
    bool has_alpha() const { return _is_alpha_channel_valid; }
    int get_width() const { return _width; }
    int get_height() const { return _height; }
    int get_bytes_per_line() const { return _bytes_per_line; }
    int get_size() const { return _color_bytes; }
    void set_xdpi(int dpi) { _xdpi = dpi; }
    void set_ydpi(int dpi) { _ydpi = dpi; }
    int get_xdpi() const { return _xdpi; }
    int get_ydpi() const { return _ydpi; }
    byte* get_data(int x, int y) const;
    byte* get_bits() const { return _data; }
    bool load(const string& filepath);
    bool load(const gchar* filepath, int len) { return load(string(filepath, len)); }
    bool save(const string& filepath) const;
    void clear(byte b, const rect* rc = nullptr);
    void clear(const color& cr, const rect* rc = nullptr);
    void copy(const image& img);
    void copy(const image& img, int x, int y, int cx, int cy, int sx, int sy);

protected:
    image_format        _format;
    int                 _width;
    int                 _height;
    int                 _depth;
    int                 _color_bytes;
    int                 _bytes_per_line;
    int                 _xdpi;
    int                 _ydpi;
    byte*               _data;
    bool                _is_alpha_channel_valid;
};

__ariel_end__

#endif

max_stars_repo_name: lymastee/gslib

max_issues_repo_name: lymastee/gslib

max_forks_repo_name: lymastee/gslib


text:
/**
 * @file      MeshHelpers.h
 *
 * @brief     This header contains mesh related helper functions
 *
 * @author    Stefan Reinhold
 * @copyright Copyright (C) 2018 Stefan Reinhold  -- All Rights Reserved.
 *            You may use, distribute and modify this code under the terms of
 *            the AFL 3.0 license; see LICENSE for full license details.
 */

#pragma once

#include "Mesh.h"

#include <Eigen/Core>
#include <gsl/gsl>
#include <igl/cotmatrix.h>
#include <igl/per_vertex_normals.h>

namespace CortidQCT {

namespace Internal {

template <class T = float>
using VertexMatrix = Eigen::Matrix<T, Eigen::Dynamic, 3>;
template <class T = float> using NormalMatrix = VertexMatrix<T>;
using FacetMatrix =
    Eigen::Matrix<typename Mesh<float>::Index, Eigen::Dynamic, 3>;
using LabelVector =
    Eigen::Matrix<typename Mesh<float>::Label, Eigen::Dynamic, 1>;
template <class T = float> using LaplacianMatrix = Eigen::SparseMatrix<T>;

/// Returns the Nx3 vertex matrix of the mesh
template <class T> inline VertexMatrix<T> vertexMatrix(Mesh<T> const &mesh) {
  return mesh
      .withUnsafeVertexPointer([&mesh](auto const *ptr) {
        return Eigen::Map<Eigen::Matrix<T, 3, Eigen::Dynamic> const>{
            ptr, 3, gsl::narrow<Eigen::Index>(mesh.vertexCount())};
      })
      .transpose();
}

/// Returns the Mx3 index matrix of the mesh
template <class T> inline FacetMatrix facetMatrix(Mesh<T> const &mesh) {
  using Index = typename Mesh<T>::Index;
  return mesh
      .withUnsafeIndexPointer([&mesh](auto const *ptr) {
        return Eigen::Map<Eigen::Matrix<Index, 3, Eigen::Dynamic> const>{
            ptr, 3, gsl::narrow<Eigen::Index>(mesh.triangleCount())};
      })
      .transpose();
}

/// Returns a N-vector containing the per vertex labels
template <class T> inline LabelVector labelVector(Mesh<T> const &mesh) {
  using Label = typename Mesh<T>::Label;
  return mesh.withUnsafeLabelPointer([&mesh](auto const *ptr) {
    return Eigen::Map<Eigen::Matrix<Label, Eigen::Dynamic, 1> const>{
        ptr, gsl::narrow<Eigen::Index>(mesh.vertexCount())};
  });
}

/// Returns a Nx3 matrix with per-vertex normals
template <class DerivedV, class DerivedF>
inline NormalMatrix<typename DerivedV::Scalar>
perVertexNormalMatrix(Eigen::MatrixBase<DerivedV> const &V,
                      Eigen::MatrixBase<DerivedF> const &F) {
  // output matrix
  NormalMatrix<typename DerivedV::Scalar> normals;

  igl::per_vertex_normals(V, F, igl::PER_VERTEX_NORMALS_WEIGHTING_TYPE_ANGLE,
                          normals);

  return normals;
}

/// Returns a Nx3 matrix with per-vertex normals
template <class T>
inline NormalMatrix<T> perVertexNormalMatrix(Mesh<T> const &mesh) {
  return perVertexNormalMatrix(vertexMatrix(mesh), facetMatrix(mesh));
}

/// Returns the NxN sparse laplacian matrix (using cotangent weights)
template <class DerivedV, class DerivedF>
inline LaplacianMatrix<typename DerivedV::Scalar>
laplacianMatrix(Eigen::MatrixBase<DerivedV> const &V,
                Eigen::MatrixBase<DerivedF> const &F) {
  // output matrix
  Eigen::SparseMatrix<typename DerivedV::Scalar> laplacian;
  Expects(V.rows() > 0 && V.cols() == 3);
  Expects(F.rows() > 0 && F.cols() == 3);
  igl::cotmatrix(V, F, laplacian);

  return laplacian;
}

/// Returns the NxN sparse laplacian matrix (using cotangent weights)
template <class T>
inline LaplacianMatrix<T> laplacianMatrix(Mesh<T> const &mesh) {
  return laplacianMatrix(vertexMatrix(mesh), facetMatrix(mesh));
}

} // namespace Internal
} // namespace CortidQCT

max_stars_repo_name: ithron/CortidQCT

max_issues_repo_name: ithron/CortidQCT

max_forks_repo_name: ithron/CortidQCT


text:
/**
 *
 * @file core_stsmlq_corner.c
 *
 *  PLASMA core_blas kernel
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Hatem Ltaief
 * @author Mathieu Faverge
 * @author Azzam Haidar
 * @date 2010-11-15
 * @generated s Tue Jan  7 11:44:48 2014
 *
 **/
#include <lapacke.h>
#include "common.h"
#undef COMPLEX
#define REAL

/***************************************************************************//**
 *
 * @ingroup CORE_float
 *
 *  CORE_stsmlq_corner: see CORE_stsmlq
 *
 * This kernel applies left and right transformations as depicted below:
 * |I -VTV'| * | A1  A2 | * |I - VT'V'|
 *             | A2' A3 |
 * where A1 and A3 are symmetric matrices.
 * Only the lower part is referenced.
 * This is an adhoc implementation, can be further optimized...
 *
 *******************************************************************************
 *
 * @param[in] m1
 *         The number of rows of the tile A1. m1 >= 0.
 *
 * @param[in] n1
 *         The number of columns of the tile A1. n1 >= 0.
 *
 * @param[in] m2
 *         The number of rows of the tile A2. m2 >= 0.
 *
 * @param[in] n2
 *         The number of columns of the tile A2. n2 >= 0.
 *
 * @param[in] m3
 *         The number of rows of the tile A3. m3 >= 0.
 *
 * @param[in] n3
 *         The number of columns of the tile A3. n3 >= 0.
 *
 * @param[in] k
 *         The number of elementary reflectors whose product defines
 *         the matrix Q.
 *
 * @param[in] ib
 *         The inner-blocking size.  ib >= 0.
 *
 * @param[in] nb
 *         The blocking size.  nb >= 0.
 *
 * @param[in,out] A1
 *         On entry, the m1-by-n1 tile A1.
 *         On exit, A1 is overwritten by the application of Q.
 *
 * @param[in] lda1
 *         The leading dimension of the array A1. lda1 >= max(1,m1).
 *
 * @param[in,out] A2
 *         On entry, the m2-by-n2 tile A2.
 *         On exit, A2 is overwritten by the application of Q.
 *
 * @param[in] lda2
 *         The leading dimension of the tile A2. lda2 >= max(1,m2).
 *
 * @param[in,out] A3
 *         On entry, the m3-by-n3 tile A3.
 *
 * @param[in] lda3
 *         The leading dimension of the tile A3. lda3 >= max(1,m3).
 *
 * @param[in] V
 *         The i-th row must contain the vector which defines the
 *         elementary reflector H(i), for i = 1,2,...,k, as returned by
 *         CORE_STSLQT in the first k rows of its array argument V.
 *
 * @param[in] ldv
 *         The leading dimension of the array V. ldv >= max(1,K).
 *
 * @param[in] T
 *         The IB-by-n1 triangular factor T of the block reflector.
 *         T is upper triangular by block (economic storage);
 *         The rest of the array is not referenced.
 *
 * @param[in] ldt
 *         The leading dimension of the array T. ldt >= IB.
 *
 * @param[out] WORK
 *         Workspace array of size
 *             LDWORK-by-m1 if side == PlasmaLeft
 *             LDWORK-by-IB if side == PlasmaRight
 *
 * @param[in] ldwork
 *         The leading dimension of the array WORK.
 *             LDWORK >= max(1,IB) if side == PlasmaLeft
 *             LDWORK >= max(1,n1) if side == PlasmaRight
 *
 *******************************************************************************
 *
 * @return
 *          \retval PLASMA_SUCCESS successful exit
 *          \retval <0 if -i, the i-th argument had an illegal value
 *
 ******************************************************************************/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_stsmlq_corner = PCORE_stsmlq_corner
#define CORE_stsmlq_corner PCORE_stsmlq_corner
#define CORE_stsmlq PCORE_stsmlq
int  CORE_stsmlq(PLASMA_enum side, PLASMA_enum trans,
                 int m1, int n1, int m2, int n2, int K, int IB,
                 float *A1, int lda1,
                 float *A2, int lda2,
                 const float *V, int ldv,
                 const float *T, int ldt,
                 float *WORK, int LDWORK);
#endif
int CORE_stsmlq_corner( int m1, int n1, int m2, int n2, int m3, int n3,
                        int k, int ib, int nb,
                        float *A1, int lda1,
                        float *A2, int lda2,
                        float *A3, int lda3,
                        const float *V, int ldv,
                        const float *T, int ldt,
                        float *WORK, int ldwork)
{
    PLASMA_enum side;
    PLASMA_enum trans;
    int i, j;

    if ( m1 != n1 ) {
        coreblas_error(1, "Illegal value of M1, N1");
        return -1;
    }

    /* Rebuild the symmetric block: WORK <- A1 */
    for (i = 0; i < m1; i++)
        for (j = i; j < n1; j++){
            *(WORK + i + j*ldwork) = *(A1 + i + j*lda1);
            if (j > i){
                *(WORK + j + i*ldwork) =  ( *(WORK + i + j*ldwork) );
            }
        }

    /*  Copy the transpose of A2: WORK+nb*ldwork <- A2' */
    for (j = 0; j < n2; j++)
        for (i = 0; i < m2; i++){
            *(WORK + j + (i + nb) * ldwork) = ( *(A2 + i + j*lda2) );
        }

    side = PlasmaRight;
    trans = PlasmaTrans;

    /*  Right application on |A1 A2| */
    CORE_stsmlq(side, trans, m1, n1, m2, n2, k, ib,
                WORK, ldwork, A2, lda2,
                V, ldv, T, ldt,
                WORK+3*nb*ldwork, ldwork);

    /*  Rebuild the symmetric block: WORK+2*nb*ldwork <- A3 */
    for (i = 0; i < m3; i++)
        for (j = i; j < n3; j++){
            *(WORK + i + (j + 2*nb) * ldwork) = *(A3 + i + j*lda3);
            if (j > i){
                *(WORK + j + (i + 2*nb) * ldwork) =   ( *(WORK + i + (j + 2*nb) * ldwork) );
            }
        }

    /*  Right application on | A2' A3 | */
    CORE_stsmlq(side, trans, n2, m2, m3, n3, k, ib,
                WORK+nb*ldwork, ldwork, WORK+2*nb*ldwork, ldwork,
                V, ldv, T, ldt,
                WORK + 3*nb*ldwork, ldwork);

    side = PlasmaLeft;
    trans = PlasmaNoTrans;

    /*  Left application on | A1  | */
    /*                      | A2' | */
    CORE_stsmlq(side, trans, m1, n1, n2, m2, k, ib,
                WORK, ldwork, WORK+nb*ldwork, ldwork,
                V, ldv, T, ldt,
                WORK + 3*nb*ldwork, ldwork);

    /*  Copy back the final result to the upper part of A1 */
    /*  A1 = WORK */
    for (i = 0; i < m1; i++)
        for (j = i; j < n1; j++)
            *(A1 + i + j*lda1) = *(WORK + i + j*ldwork);

    /*  Left application on | A2 | */
    /*                      | A3 | */
    CORE_stsmlq(side, trans, m2, n2, m3, n3, k, ib,
                A2, lda2, WORK+2*nb*ldwork, ldwork,
                V, ldv, T, ldt,
                WORK + 3*nb*ldwork, ldwork);

    /*  Copy back the final result to the upper part of A3 */
    /*  A3 = WORK+2*nb*ldwork */
    for (i = 0; i < m3; i++)
        for (j = i; j < n3; j++)
            *(A3 + i + j*lda3) = *(WORK + i + (j+ 2*nb) * ldwork);

    return PLASMA_SUCCESS;
}


text:
//
// Created by Harold on 2020/10/7.
//

#ifndef M_MATH_M_HIST_H
#define M_MATH_M_HIST_H

#include <gsl/gsl_histogram.h>
#include <gsl/gsl_histogram2d.h>
#include <vector>
#include <utility>

namespace M_MATH {
    template<typename T>
    class Histogram1D {
    public:
        // [min, max), it: <x>
        template<typename IT>
        Histogram1D(IT begin, IT end, T min, T max, size_t N);
        ~Histogram1D() { gsl_histogram_free(h); gsl_histogram_pdf_free(p); };
        // counts in every bin
        std::vector<size_t> counts();
        // bins with [lower, upper] boundaries
        std::vector<std::pair<T, T>> bins();
        // probability distribution, r in [0, 1]
        double pdf_sample(double r);
        // standard deviation of the histogrammed variable
        double sigma() { return gsl_histogram_sigma(h); }
        // mean of the histogrammed variable
        double mean() { return gsl_histogram_mean(h); }
        // sum of all bin values
        double sum() { return gsl_histogram_sum(h); }
        // find x in which bin
        size_t find(double x) { size_t idx; gsl_histogram_find(h, x, &idx); return idx; }

    private:
        gsl_histogram *h;
        gsl_histogram_pdf *p;
    };

    template<typename T>
    class Histogram2D {
    public:
        // [min, max), it: <x, y>
        template<typename IT>
        Histogram2D(IT begin, IT end, T xmin, T xmax, T ymin, T ymax, size_t xdim, size_t ydim);
        ~Histogram2D() { gsl_histogram2d_free(h); gsl_histogram2d_pdf_free(p); };
        // counts in every bin
        std::vector<size_t> counts();
        // bins with [lower, upper] boundaries, <xrange, yrange>
        std::pair<std::vector<std::pair<T, T>>, std::vector<std::pair<T, T>>> bins();
        // probability distribution <x, y>, r1, r2 in [0, 1]
        std::pair<double, double> pdf_sample(double r1, double r2);
        // standard deviation of the histogrammed x and y variables
        std::pair<double, double> sigma() { return std::make_pair(gsl_histogram2d_xsigma(h), gsl_histogram2d_ysigma(h)); }
        // mean of the histogrammed x and y variables
        std::pair<double, double> mean() { return std::make_pair(gsl_histogram2d_xmean(h), gsl_histogram2d_ymean(h)); }
        // sum of all bin values
        double sum() { return gsl_histogram2d_sum(h); }
        // covariance of the histogrammed x and y variables
        double cov() { return gsl_histogram2d_cov(h); }
        // find x in which bin, <x, y>
        std::pair<size_t, size_t> find(double x, double y) { size_t xidx, yidx; gsl_histogram2d_find(h, x, y, &xidx, &yidx); return std::make_pair(xidx, yidx); }
    private:
        gsl_histogram2d *h;
        gsl_histogram2d_pdf *p;
    };

    template<typename T>
    template<typename IT>
    Histogram1D<T>::Histogram1D(IT const begin, IT const end, T min, T max, size_t N) : h(nullptr), p(nullptr) {
        h = gsl_histogram_alloc(N);
        gsl_histogram_set_ranges_uniform(h, min, max);
        for (auto it = begin; it != end; ++it)
            gsl_histogram_increment(h, *it);
    }

    template<typename T>
    std::vector<size_t> Histogram1D<T>::counts() {
        std::vector<size_t> res(h->n);
        for (auto i = 0; i < h->n; ++i)
            res[i] = gsl_histogram_get(h, i);
        return res;
    }

    template<typename T>
    std::vector<std::pair<T, T>> Histogram1D<T>::bins() {
        std::vector<std::pair<T, T>> res(h->n);
        double lower, upper;
        for (auto i = 0; i < h->n; ++i) {
            gsl_histogram_get_range(h, i, &lower, &upper);
            res[i] = std::make_pair(lower, upper);
        }
        return res;
    }

    template<typename T>
    double Histogram1D<T>::pdf_sample(double r) {
        if (p == nullptr) {
            p = gsl_histogram_pdf_alloc(h->n);
            gsl_histogram_pdf_init(p, h);
        }
        return gsl_histogram_pdf_sample(p, r);
    }

    template<typename T>
    template<typename IT>
    Histogram2D<T>::Histogram2D(IT begin, IT end,
                                T xmin, T xmax,
                                T ymin, T ymax,
                                size_t xdim, size_t ydim) : h(nullptr), p(nullptr) {
        h = gsl_histogram2d_alloc(xdim, ydim);
        gsl_histogram2d_set_ranges_uniform(h, xmin, xmax, ymin, ymax);
        for (auto it = begin; it != end; ++it)
            gsl_histogram2d_increment(h, (*it).x, (*it).y);
    }

    template<typename T>
    std::vector<size_t> Histogram2D<T>::counts() {
        std::vector<size_t> res((h->nx) * (h->ny));
        for (auto i = 0; i < h->nx; ++i)
            for (auto j = 0; j < h->ny; ++j)
                res[j * (h->nx) + i] = gsl_histogram2d_get(h, i, j);
        return res;
    }

    template<typename T>
    std::pair<std::vector<std::pair<T, T>>, std::vector<std::pair<T, T>>> Histogram2D<T>::bins() {
        std::vector<std::pair<T, T>> xrange(h->nx), yrange(h->ny);
        double lower, upper;
        for (auto i = 0; i < h->nx; ++i) {
            gsl_histogram2d_get_xrange(h, i, &lower, &upper);
            xrange[i] = std::make_pair(lower, upper);
        }
        for (auto i = 0; i < h->ny; ++i) {
            gsl_histogram2d_get_yrange(h, i, &lower, &upper);
            yrange[i] = std::make_pair(lower, upper);
        }
        return std::make_pair(xrange, yrange);
    }

    template<typename T>
    std::pair<double, double> Histogram2D<T>::pdf_sample(double r1, double r2) {
        if (p == nullptr) {
            p = gsl_histogram2d_pdf_alloc(h->nx, h->ny);
            gsl_histogram2d_pdf_init(p, h);
        }
        double x, y;
        gsl_histogram2d_pdf_sample(p, r1, r2, &x, &y);
        return std::make_pair(x, y);
    }
}

#endif //M_MATH_M_HIST_H


text:
/* include */
#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"
#include <stdarg.h>

#define degree 0.0174
double min (int n, ...) {
	/*
	n是參數個數，後面才是參數本身
	*/
	int 	i;
	double 	min_num = 1e20;
	double 	input;
	va_list vl;
	va_start(vl,n);
	for ( i = 0 ; i < n ; i++ ) {
		input = va_arg(vl,double);
		min_num = min_num > input ? input : min_num;
	}
	va_end(vl);
	return min_num;
 }

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

/* 定義3 sigma range 的Prior (For NO) */
double prior_3sigma_NO(const glb_params in, void* user_data) {
	glb_params 		central_values 	= 	glbAllocParams();
	glb_params 		input_errors 	= 	glbAllocParams();
	glb_projection 	p 				= 	glbAllocProjection();

	glbGetCentralValues(central_values);
	glbGetInputErrors(input_errors);
	glbGetProjection(p);
	int i;
	double 	pv = 0.0,
		  	fit_theta12,
		   	fit_theta13,
		   	fit_theta23,
		   	fit_deltacp,
		   	fit_ldm,
		   	fit_sdm;

    /* 取得參數目前Fit Value */
	fit_theta12 = glbGetOscParams(in,0);
	fit_theta13 = glbGetOscParams(in,1);
	fit_theta23 = glbGetOscParams(in,2);
	fit_deltacp = glbGetOscParams(in,3);
	fit_sdm     = glbGetOscParams(in,4);
	fit_ldm     = glbGetOscParams(in,5);

    /* 判斷參數是否要引入Prior */
	if (glbGetProjectionFlag(p,0)==GLB_FREE){
		if (fit_theta12  > 35.86 *degree || fit_theta12 < 31.27 *degree){
		    pv += 1e20;
		}
	}
	if (glbGetProjectionFlag(p,1)==GLB_FREE){
		if (fit_theta13  > 8.97 *degree || fit_theta13 < 8.20 *degree){
		    pv += 1e20;
		}
	}
	if (glbGetProjectionFlag(p,2)==GLB_FREE){
		if (fit_theta23 > 51.80 *degree  || fit_theta23 < 39.60 *degree ){
		    pv += 1e20;
		}
	}
	if (glbGetProjectionFlag(p,3)==GLB_FREE){
		if (fit_deltacp == 0 *degree || fit_deltacp == 180 *degree || fit_deltacp < 0 *degree|| fit_deltacp > 360 *degree){
		    pv += 1e20;
		}
	}
	if (glbGetProjectionFlag(p,4)==GLB_FREE){
		if (fit_sdm > 8.04 *1e-5  || fit_sdm  < 6.82 *1e-5){
		    pv += 1e20;
		}
	}
	if (glbGetProjectionFlag(p,5)==GLB_FREE){
		if (fit_ldm  > 2.598 *1e-3 || fit_ldm  < 2.431 *1e-3){
		    pv += 1e20;
		}
	}

	glbFreeParams(central_values);
	glbFreeParams(input_errors);
	glbFreeProjection(p);
	return pv;
}

 /* 定義3 sigma range 的Prior (For IO) */
double prior_3sigma_IO(const glb_params in, void* user_data)
{
	glb_params 		central_values 	= 	glbAllocParams();
	glb_params 		input_errors 	= 	glbAllocParams();
	glb_projection 	p 				= 	glbAllocProjection();

	glbGetCentralValues(central_values);
	glbGetInputErrors(input_errors);
	glbGetProjection(p);
	int i;
	double 	pv = 0.0,
		   	fit_theta12,
		   	fit_theta13,
		   	fit_theta23,
		   	fit_deltacp,
		   	fit_ldm,
		   	fit_sdm;

	/* 取得參數目前Fit Value */
	fit_theta12 = glbGetOscParams(in,0);
	fit_theta13 = glbGetOscParams(in,1);
	fit_theta23 = glbGetOscParams(in,2);
	fit_deltacp = glbGetOscParams(in,3);
	fit_sdm     = glbGetOscParams(in,4);
	fit_ldm     = glbGetOscParams(in,5);

	/* 判斷參數是否要引入Prior */
	if(glbGetProjectionFlag(p,0)==GLB_FREE){
		if(fit_theta12  > 35.87 *degree || fit_theta12 < 31.27 *degree){
		    pv += 1e20;
		}
	}
	if(glbGetProjectionFlag(p,1)==GLB_FREE){
		if(fit_theta13  > 8.98 *degree || fit_theta13 < 8.24 *degree){
		    pv += 1e20;
		}
	}
	if(glbGetProjectionFlag(p,2)==GLB_FREE){
		if(fit_theta23 > 52.00 *degree  || fit_theta23 < 39.90 *degree ){
		    pv += 1e20;
		}
	}
	if(glbGetProjectionFlag(p,3)==GLB_FREE){
		if(fit_deltacp == 0 *degree || fit_deltacp == 180 *degree ){
		    pv += 1e20;
		}
	}
	if(glbGetProjectionFlag(p,4)==GLB_FREE){
		if(fit_sdm > 8.04 *1e-5  || fit_sdm  < 6.82 *1e-5){
		    pv += 1e20;
		}
	}
	if(glbGetProjectionFlag(p,5)==GLB_FREE){
		if(fit_ldm  < -2.583 *1e-3 || fit_ldm  > -2.412 *1e-3){
		    pv += 1e20;
		}
	}

	glbFreeParams(central_values);
	glbFreeParams(input_errors);
	glbFreeProjection(p);
	return pv;
}

/* Poisson亂數生成器 */
int random_poisson(double mu) {
	const 	gsl_rng_type * T;
	gsl_rng * r;
	int 	test;
	int 	i;
	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;
	r = gsl_rng_alloc (T);
	gsl_rng_set(r, mySeed);
	unsigned int k = gsl_ran_poisson (r, mu);
	gsl_rng_free (r);
	return k;
}

/* 定義Poisson Likelihood Function */
inline double poisson_likelihood(double true_rate, double fit_rate) {
	double res = 0 ;
	true_rate = true_rate == 0 ? (true_rate+1e-9) : true_rate;
	fit_rate = fit_rate == 0 ? (fit_rate+1e-9) : fit_rate;

	res = fit_rate - true_rate;
	if (fit_rate <= 0.0) {
		res = 1e100;
	} else if (true_rate > 0) {
		res += true_rate * log(true_rate/fit_rate);
	}

	return 2.0 * res;
}

/* 對True Value Spectrum 做Poisson Fluctuation */
void do_poisson_fluctuation(glb_params test_values, int mode_expr, double* dataset, int* dset_info) {

	glbSetOscillationParameters(test_values);
	glbSetRates();

	int		ew_low,
			ew_high,
			i,
			ladder = 0;
	int    	count_dpf = 0;
	int 	num_channel,
			num_bins,
			cumu_bins;
	if (mode_expr == 0){
		/* Deprecated
		double *ve_dune      = glbGetRuleRatePtr(0, 0);
		double *vebar_dune   = glbGetRuleRatePtr(0, 1);
		double *vu_dune      = glbGetRuleRatePtr(0, 2);
		double *vubar_dune   = glbGetRuleRatePtr(0, 3);
		*/

		num_channel = 4;
		cumu_bins	= 0;
		for (int channel = 0; channel < num_channel; channel++){
			double *target = glbGetRuleRatePtr(mode_expr, channel);
			glbGetEnergyWindowBins(mode_expr, channel, &ew_low, &ew_high);
			num_bins = ew_high - ew_low + 1;
			count_dpf = 0;
			for (i=ew_low; i <= ew_high; i++) {
				*(dataset + cumu_bins + count_dpf) = random_poisson(target[i]);
				count_dpf += 1;
			}
			cumu_bins += num_bins;
			*(dset_info + channel) = num_bins;
		}
	} else if (mode_expr == 1){
		/* Deprecated
		double *ve_t2hk      = glbGetRuleRatePtr(1, 0);
		double *vu_t2hk      = glbGetRuleRatePtr(1, 1);
		double *vebar_t2hk   = glbGetRuleRatePtr(1, 2);
		double *vubar_t2hk   = glbGetRuleRatePtr(1, 3);
		*/
		num_channel = 4;
		cumu_bins	= 0;
		for (int channel = 0; channel < num_channel; channel++){
			double *target = glbGetRuleRatePtr(mode_expr, channel);
			glbGetEnergyWindowBins(mode_expr, channel, &ew_low, &ew_high);
			num_bins = ew_high - ew_low + 1;
			count_dpf = 0;
			for (i=ew_low; i <= ew_high; i++) {
				*(dataset + cumu_bins + count_dpf) = random_poisson(target[i]);
				count_dpf += 1;
			}
			cumu_bins += num_bins;
			*(dset_info + channel) = num_bins;
		}
	} else if (mode_expr == -1) {
		/* Deprecated
		double *ve_dune      = glbGetRuleRatePtr(0, 0);
		double *vebar_dune   = glbGetRuleRatePtr(0, 1);
		double *vu_dune      = glbGetRuleRatePtr(0, 2);
		double *vubar_dune   = glbGetRuleRatePtr(0, 3);
		double *ve_t2hk      = glbGetRuleRatePtr(1, 0);
		double *vu_t2hk      = glbGetRuleRatePtr(1, 1);
		double *vebar_t2hk   = glbGetRuleRatePtr(1, 2);
		double *vubar_t2hk   = glbGetRuleRatePtr(1, 3);
		*/
		int num_expr	= 2;
		num_channel 	= 4;
		cumu_bins		= 0;
		int expr;
		for ( expr = 0; expr < num_expr; expr++) {
			for (int channel = 0; channel < num_channel; channel++){
				double *target = glbGetRuleRatePtr(expr, channel);
				glbGetEnergyWindowBins(expr, channel, &ew_low, &ew_high);
				num_bins = ew_high - ew_low + 1;
				count_dpf = 0;
				for (i=ew_low; i <= ew_high; i++) {
					*(dataset + cumu_bins + count_dpf) = random_poisson(target[i]);
					count_dpf += 1;
				}
				cumu_bins += num_bins;
				*(dset_info + channel) = num_bins;
			}
		}
	} else {
		printf("Please inpuit a correct experiment protocol.");
	}
}

/* 定義 Chi Square */
double chi2_poisson(int exp, int rule, int np, double *x, double *errors, void* user_data) {
	double 	*signal_fit_rate = glbGetSignalFitRatePtr(exp, rule);
	double 	*bg_fit_rate     = glbGetBGFitRatePtr(exp, rule);
	double 	fit_rate;
	double 	chi2 = 0.0;
	int 	i;
	int 	ew_low,
			ew_high;
	glbGetEnergyWindowBins(exp, rule, &ew_low, &ew_high);
	int 	sum_y = 0,
			index,
			pos,
			count_cp;
	int 	data_info_dune[4] = { 66, 66, 66, 66}; /* {(y_axis_length)} */
	int 	data_info_t2hk[4] = { 8, 12, 8, 12};
	int 	data_info_all[8] = { 66, 66, 66, 66, 8, 12, 8, 12};

	switch (exp) {
		case 0:
			for ( pos = 0; pos < rule; pos ++){
				sum_y += data_info_dune[pos];
			}

			for (i=ew_low; i <= ew_high; i++) {
				fit_rate = signal_fit_rate[i] + bg_fit_rate[i];
				chi2 += poisson_likelihood( *( (double*) user_data + sum_y + i), fit_rate);
			}
			break;
		case 1:
			for ( pos = 0; pos < rule; pos ++){
				sum_y += data_info_t2hk[pos];
			}

			for (i=ew_low; i <= ew_high; i++) {
				fit_rate = signal_fit_rate[i] + bg_fit_rate[i];
				chi2 += poisson_likelihood( *( (double*) user_data + sum_y + i), fit_rate);
			}
			break;
	}
    return chi2;
}

/* 定義 Test Statistic (Delta Chi-Square) */
//參數:{CP : [0, 180, 1], MO : [1, -1] , CPV Hypothesis的deltacp}
// 0 : CPC at deltacp=0; 180 : CPC at deltacp=180; 1 : CPV at deltacp
//選定實驗EXP:[0,1,GLB_ALL]

double delta_chi2 (int CP, int MO, double deltacp, int EXP) {
	double  a = 0,
			b = 0;
	double  chi_0_NO,
			chi_0_IO,
			chi_pi_NO,
			chi_pi_IO,
			chi_cpv_NO,
			chi_cpv_IO;

	int data_info_dune[6] = {4, 264, 66, 66, 66, 66}; /* {(y_axis_length)} */
	int data_info_t2hk[6] = {4, 40, 8, 12, 8, 12};
	int data_info_all[10] = {8, 304, 66, 66, 66, 66, 8, 12, 8, 12};
	double *darray;
	int sum_y  = 0,
		ladder = 0,
		index_sum,
		i_idx, j_idx;
	int LENGTH = EXP != -1 ? 4 : 8;
	int dset_y_info[LENGTH];

	/* 定義global fit參數(Normal Ordering, NuFIT 5.0, 2020) */
	double theta12_N = 33.44;
	double theta13_N = 8.57;
	double theta23_N = 49;
	double sdm_N = 7.42;
	double ldm_N = 2.514;

	/* 定義global fit參數(Inverse Ordering, NuFIT 5.0, 2020) */
	double theta12_I = 33.45;
	double theta13_I = 8.61;
	double theta23_I = 49.3;
	double sdm_I = 7.42;
	double ldm_I = -2.497;

	switch (EXP){
		case 0:
			for ( index_sum = 2; index_sum < 6; index_sum++ ){
				sum_y += data_info_dune[index_sum];
			}
			darray = malloc( data_info_dune[0] * sum_y * sizeof(double));
			break;
		case 1:
			for ( index_sum = 2; index_sum < 6; index_sum++ ){
				sum_y += data_info_t2hk[index_sum];
			}
			darray = malloc( data_info_t2hk[0] * sum_y * sizeof(double));
			break;
		case -1:
			for ( index_sum = 2; index_sum < 10; index_sum++ ){
				sum_y += data_info_all[index_sum];
			}
			darray = malloc( data_info_all[0] * sum_y * sizeof(double));
			break;
	}

	//根據CP、MO的假設，生成Poisson Sample，計算其test statistic
	/* 定義glb_params */
	glb_params test_values_cpc_0_NO  = glbAllocParams();
	glb_params test_values_cpc_0_IO  = glbAllocParams();
	glb_params test_values_cpc_pi_NO = glbAllocParams();
	glb_params test_values_cpc_pi_IO = glbAllocParams();
	glb_params test_values_cpv_NO    = glbAllocParams();
	glb_params test_values_cpv_IO    = glbAllocParams();
	glb_params input_errors = glbAllocParams();
	glb_params minimum = glbAllocParams(); //////////

	/* 定義test_values_cpc_0_NO */
	glbDefineParams(test_values_cpc_0_NO, theta12_N*degree, theta13_N*degree, theta23_N*degree, 0*degree, 1e-5*sdm_N, 1e-3*ldm_N);
	glbSetDensityParams(test_values_cpc_0_NO,1.0,GLB_ALL);

	/* 定義test_values_cpc_0_IO */
	glbDefineParams(test_values_cpc_0_IO, theta12_I*degree, theta13_I*degree, theta23_I*degree, 0*degree, 1e-5*sdm_I, 1e-3*ldm_I);
	glbSetDensityParams(test_values_cpc_0_IO,1.0,GLB_ALL);

	/* 定義test_values_cpc_pi_NO */
	glbDefineParams(test_values_cpc_pi_NO, theta12_N*degree, theta13_N*degree, theta23_N*degree, 180*degree, 1e-5*sdm_N, 1e-3*ldm_N);
	glbSetDensityParams(test_values_cpc_pi_NO,1.0,GLB_ALL);

	/* 定義test_values_cpc_pi_IO */
	glbDefineParams(test_values_cpc_pi_IO, theta12_I*degree, theta13_I*degree, theta23_I*degree, 180*degree, 1e-5*sdm_I, 1e-3*ldm_I);
	glbSetDensityParams(test_values_cpc_pi_IO,1.0,GLB_ALL);

	/* 定義test_values_cpv_NO */                                                         //deltacp為Input的值
	glbDefineParams(test_values_cpv_NO, theta12_N*degree, theta13_N*degree, theta23_N*degree, deltacp*degree, 1e-5*sdm_N, 1e-3*ldm_N);
	glbSetDensityParams(test_values_cpv_NO,1.0,GLB_ALL);

	/* 定義test_values_cpv_IO */                                                         //deltacp為Input的值
	glbDefineParams(test_values_cpv_IO, theta12_I*degree, theta13_I*degree, theta23_I*degree, deltacp*degree, 1e-5*sdm_I, 1e-3*ldm_I);
	glbSetDensityParams(test_values_cpv_IO,1.0,GLB_ALL);

	/* 設定Projection */
	glb_projection projection_cp_fixed = glbAllocProjection();
	glb_projection projection_cp_free  = glbAllocProjection();

	//GLB_FIXED/GLB_FREE                      theta12    theta13  theta23    deltacp     m21        m31
	glbDefineProjection(projection_cp_fixed, GLB_FIXED, GLB_FREE, GLB_FREE, GLB_FIXED, GLB_FIXED, GLB_FREE);//deltacp theta12 m21 不動，其他可變
	glbSetDensityProjectionFlag(projection_cp_fixed,GLB_FIXED,GLB_ALL);//matter density不變

	//GLB_FIXED/GLB_FREE                      theta12    theta13  theta23    deltacp     m21        m31
	glbDefineProjection(projection_cp_free,  GLB_FIXED, GLB_FREE, GLB_FREE, GLB_FREE, GLB_FIXED, GLB_FREE);// theta12 m21 不動，其他可變
	glbSetDensityProjectionFlag(projection_cp_free,GLB_FIXED,GLB_ALL);//matter density不變

	/* 關閉系統誤差 */
	glbSwitchSystematics(GLB_ALL,GLB_ALL,GLB_OFF);

	/* 設定Input_errors */
	glbDefineParams(input_errors,0,0,0,0,0,0);
	glbSetDensityParams(input_errors,0,GLB_ALL);
	glbSetInputErrors(input_errors);

	/* 根據MO,CP的假設，生成 Poisson Spectrum */
	switch (MO) {
		case 1:
	  		if (CP == 0){ //CPC at deltacp=0
				printf("生成deltacp = 0, NO 的Poisson Spectrum \n");
				/* 根據CPC_0_NO的假設，生成Poisson True Spectrum */
				do_poisson_fluctuation(test_values_cpc_0_NO, EXP, darray, dset_y_info);
			}
	  		if (CP == 180){ //CPC at deltacp=180
				printf("生成deltacp = 180, NO 的Poisson Spectrum \n");
				/* 根據CPC_0_IO的假設，生成Poisson True Spectrum */
		  		do_poisson_fluctuation(test_values_cpc_pi_NO, EXP, darray, dset_y_info);
			}
	  		if (CP == 1){ //CPV
				printf("生成deltacp = input value, NO 的Poisson Spectrum \n");
				/* 根據CPV_NO的假設，生成Poisson True Spectrum */
				do_poisson_fluctuation(test_values_cpv_NO, EXP, darray, dset_y_info);
			}
			break;
		case -1:
	  		if (CP == 0){ //CPC at deltacp=0
				printf("生成deltacp = 0, IO 的Poisson Spectrum \n");
				/* 根據CPC_0_NO的假設，生成Poisson True Spectrum */
		  		do_poisson_fluctuation(test_values_cpc_0_IO, EXP, darray, dset_y_info);
			}
	  		if (CP == 180){ //CPC at deltacp=180
				printf("生成deltacp = 180, IO 的Poisson Spectrum \n");
				/* 根據CPC_0_IO的假設，生成Poisson True Spectrum */
		  		do_poisson_fluctuation(test_values_cpc_pi_IO, EXP, darray, dset_y_info);
			}
	  		if (CP == 1){ //CPV
				printf("生成deltacp = input value, IO 的Poisson Spectrum \n");
				/* 根據CPV_IO的假設，生成Poisson True Spectrum */
		  		do_poisson_fluctuation(test_values_cpv_IO, EXP, darray, dset_y_info);
			}
			break;
	}

	/* 計算CPC Hypothesis (4種情況)*/

	/* 設定Prior (3 sigma range, Normal Ordering)*/
	glbRegisterPriorFunction(prior_3sigma_NO,NULL,NULL,NULL);

	/* 計算Chi square under cpc_0_NO */
	glbSetProjection(projection_cp_fixed); //設定Projection deltacp_Fixed
	glbSetOscillationParameters(test_values_cpc_0_NO);
	glbSetRates();
	glbDefineChiFunction(&chi2_poisson, 0, "chi2_poisson", darray);
	glbSetChiFunction(GLB_ALL, GLB_ALL, GLB_OFF, "chi2_poisson", NULL);
	glbSetCentralValues(test_values_cpc_0_NO);
	chi_0_NO = glbChiNP(test_values_cpc_0_NO, minimum ,EXP);

	/* 計算Chi square under cpc_pi_NO */
	glbSetProjection(projection_cp_fixed); //設定Projection deltacp_Fixed
	glbSetOscillationParameters(test_values_cpc_pi_NO);
	glbSetRates();
	glbDefineChiFunction(&chi2_poisson, 0, "chi2_poisson", darray);
	glbSetChiFunction(GLB_ALL, GLB_ALL, GLB_OFF, "chi2_poisson", NULL);
	glbSetCentralValues(test_values_cpc_pi_NO);
	chi_pi_NO = glbChiNP(test_values_cpc_pi_NO, minimum ,EXP);

	/* 設定Prior (3 sigma range, Inverse Ordering)*/
	glbRegisterPriorFunction(prior_3sigma_IO,NULL,NULL,NULL);

	/* 計算Chi square under cpc_0_IO */
	glbSetProjection(projection_cp_fixed); //設定Projection deltacp_Fixed
	glbSetOscillationParameters(test_values_cpc_0_IO);
	glbSetRates();
	glbDefineChiFunction(&chi2_poisson, 0, "chi2_poisson", darray);
	glbSetChiFunction(GLB_ALL, GLB_ALL, GLB_OFF, "chi2_poisson", NULL);
	glbSetCentralValues(test_values_cpc_0_IO);
	chi_0_IO = glbChiNP(test_values_cpc_0_IO, minimum ,EXP);

	/* 計算Chi square under cpc_pi_IO */
	glbSetProjection(projection_cp_fixed); //設定Projection deltacp_Fixed
	glbSetOscillationParameters(test_values_cpc_pi_IO);
	glbSetRates();
	glbDefineChiFunction(&chi2_poisson, 0, "chi2_poisson", darray);
	glbSetChiFunction(GLB_ALL, GLB_ALL, GLB_OFF, "chi2_poisson", NULL);
	glbSetCentralValues(test_values_cpc_pi_IO);
	chi_pi_IO = glbChiNP(test_values_cpc_pi_IO, minimum ,EXP);

	/* 取 chi_0_NO , chi_pi_NO, chi_0_IO , chi_pi_IO  四者之最小值 */
	a = min(4, chi_0_NO, chi_pi_NO, chi_0_IO, chi_pi_IO);

	/* 計算CPV Hypothesis (2種情況)*/

	/* 設定Prior (3 sigma range, Normal Ordering)*/
	glbRegisterPriorFunction(prior_3sigma_NO, NULL, NULL, NULL);

	/* 設定Projection (deltacp_Free)*/
	glbSetProjection(projection_cp_free);

	/* 計算Chi square under cpv_NO */
	glbSetOscillationParameters(test_values_cpv_NO);
	glbSetRates();
	glbDefineChiFunction(&chi2_poisson, 0, "chi2_poisson", darray);
	glbSetChiFunction(GLB_ALL, GLB_ALL, GLB_OFF, "chi2_poisson", NULL);
	glbSetCentralValues(test_values_cpv_NO);
	chi_cpv_NO = glbChiNP(test_values_cpv_NO, minimum, EXP);

	/* 設定Prior (3 sigma range, Inverse Ordering)*/
	glbRegisterPriorFunction(prior_3sigma_IO, NULL, NULL, NULL);

	/* 設定Projection (deltacp_Free)*/
	glbSetProjection(projection_cp_free);

	/* 計算Chi square under cpv_IO */
	glbSetOscillationParameters(test_values_cpv_IO);
	glbSetRates();
	glbDefineChiFunction(&chi2_poisson, 0, "chi2_poisson", darray);
	glbSetChiFunction(GLB_ALL, GLB_ALL, GLB_OFF, "chi2_poisson", NULL);
	glbSetCentralValues(test_values_cpv_IO);
	chi_cpv_IO = glbChiNP(test_values_cpv_IO, minimum, EXP);

	/* 取 chi_cpv_NO , chi_cpv_IO  兩者之最小值 */
	b = min(2, chi_cpv_NO, chi_cpv_IO);

	/* 輸出Delta Chi square */
	printf("a = %g, b = %g  \n",a,b);
	printf("a - b = %g  \n",a-b);
	free(darray);
	return a-b;
}

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

	char   filename[32];
	strcpy(filename, argv[1]);
	int    TOTALsample = atof(argv[2]);
	double angle       = atof(argv[3]);
	int    expr        = atof(argv[4]);

	/*
     * Declare the variables for hdf5.
     */
	hid_t		file_id,
				space_id,
				dset_id,
				group_id,
				sub_group_id;
	herr_t		status;
	hsize_t		dims_q0[2] = {TOTALsample, 4},
				dims_q1[1] = {TOTALsample};

	glbInit(argv[0]);
	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);

	printf("File name: %s\n",filename);

	double 	q0,  q1;
	double 	Q0[TOTALsample][4],
			Q1[TOTALsample];
	for (int num_simulatiom = 0; num_simulatiom < TOTALsample; num_simulatiom++) {
		q0 = delta_chi2(0 , 1 , angle, expr);
		Q0[num_simulatiom][0] = q0;

		q0 = delta_chi2(180, 1 , angle, expr);
		Q0[num_simulatiom][1] = q0;

		q0 = delta_chi2(0, -1 , angle, expr);
		Q0[num_simulatiom][2] = q0;

		q0 = delta_chi2(180, -1 , angle, expr);
		Q0[num_simulatiom][3] = q0;

		q1 = delta_chi2(1 , 1, angle, expr);
		Q1[num_simulatiom] = q1;
	}

	file_id   	= 	H5Fcreate(filename, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
	space_id  	=  	H5Screate_simple(2, dims_q0, NULL);
	dset_id    	=  	H5Dcreate(file_id, "q0", H5T_NATIVE_DOUBLE, space_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
	status     	=  	H5Dwrite(dset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, Q0);
	status    	=  	H5Dclose(dset_id);
	status    	=  	H5Sclose(space_id);

	space_id  	=  	H5Screate_simple(1, dims_q1, NULL);
	dset_id    	=  	H5Dcreate(file_id, "q1", H5T_NATIVE_DOUBLE, space_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
	status     	=  	H5Dwrite(dset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, Q1);
	status    	=  	H5Dclose(dset_id);
	status    	=  	H5Sclose(space_id);
	status     	=  	H5Fclose(file_id);
	return 0;
}


text:
/* -*- Mode: C++; tab-width: 4; c-basic-offset: 4; indent-tabs-mode: nil -*- */
#pragma once

#include <event.h>
#include <memcached/engine_error.h>
#include <memcached/types.h>
#include <platform/socket.h>
#include <subdoc/operations.h>
#include <gsl/gsl>

#include <mutex>
#include <queue>
#include <unordered_map>
#include <vector>

/** \file
 * The main memcached header holding commonly used data
 * structures and function prototypes.
 */

/** Maximum length of a key. */
#define KEY_MAX_LENGTH 250

#define MAX_SENDBUF_SIZE (256 * 1024 * 1024)

/* Maximum length of config which can be validated */
#define CONFIG_VALIDATE_MAX_LENGTH (64 * 1024)

/* Maximum IOCTL get/set key and payload (body) length */
#define IOCTL_KEY_LENGTH 128
#define IOCTL_VAL_LENGTH 128

#define MAX_VERBOSITY_LEVEL 2

class Cookie;
class Connection;
struct thread_stats;

void initialize_buckets();
void cleanup_buckets();

void associate_initial_bucket(Connection& connection);

/*
 * Functions such as the libevent-related calls that need to do cross-thread
 * communication in multithreaded mode (rather than actually doing the work
 * in the current thread) are called via "dispatch_" frontends, which are
 * also #define-d to directly call the underlying code in singlethreaded mode.
 */
void worker_threads_init();

void threads_shutdown();
void threads_cleanup();

/**
 * Create a socketpair and make it non-blocking
 *
 * @param sockets Where to store the sockets
 * @return true if success, false otherwise (and the error reason logged)
 */
bool create_nonblocking_socketpair(std::array<SOCKET, 2>& sockets);

class ListeningPort;
void dispatch_conn_new(SOCKET sfd, std::shared_ptr<ListeningPort>& interface);

void threadlocal_stats_reset(std::vector<thread_stats>& thread_stats);

void notify_io_complete(gsl::not_null<const void*> cookie,
                        ENGINE_ERROR_CODE status);
void safe_close(SOCKET sfd);
int add_conn_to_pending_io_list(Connection* c,
                                Cookie* cookie,
                                ENGINE_ERROR_CODE status);
const char* get_server_version();
bool is_memcached_shutting_down();

/**
 * Connection-related functions
 */

/**
 * Increments topkeys count for the key specified within the command context
 * provided by the cookie.
 */
void update_topkeys(const Cookie& cookie);

struct ServerApi;
ServerApi* get_server_api();

void shutdown_server();
bool associate_bucket(Connection& connection, const char* name);
void disassociate_bucket(Connection& connection);

void disable_listen();
bool is_listen_disabled();

/**
 * The executor pool used to pick up the result for requests spawn by the
 * client io threads and dispatched over to a background thread (in order
 * to allow for out of order replies).
 */
namespace cb {
class ExecutorPool;
}
extern std::unique_ptr<cb::ExecutorPool> executorPool;

void iterate_all_connections(std::function<void(Connection&)> callback);

void start_stdin_listener(std::function<void()> function);


text:
/* Output the current topic and POV assignments in text in text form
   to stdout, optionally iteratively maximizing the assignments first
   (to find a high-probability assignments). Even if performing
   maximization, the current assignments should be post-burn-in for
   best results. */

#include <assert.h>
#include <gsl/gsl_rng.h>
#include <inttypes.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>

#include "index.h"
#include "parse_mmaps.h"
#include "probability.h"
#include "sample.h"

int main(int argc, char **argv) {
  if (argc != 4) {
    printf("Usage: %s mmap_directory maximization_iterations threads\n",
           argv[0]);
    exit(1);
  }
  int maximization_iterations = atoi(argv[2]);
  int num_threads = atoi(argv[3]);
  struct mmap_info mmap_info = open_mmaps_readonly(argv[1]);
  int64_t count_revisions;
  struct revision_assignment* revision_assignments;
  get_revision_assignment_array(&mmap_info, &count_revisions, &revision_assignments);

  if (maximization_iterations > 0) {
    struct sample_threads sample_threads;
    initialize_threads(&sample_threads, num_threads, &mmap_info);
    for (int i = 0; i < maximization_iterations; ++i) {
      resample_maximize(&sample_threads);
    }
    destroy_threads(&sample_threads);
  }
  for (int64_t revision_num = 0; revision_num < count_revisions; ++revision_num) {
    if (revision_assignments[revision_num].pov < 0
	|| revision_assignments[revision_num].topic < 0) {
      continue;
    }
    printf("%" PRId64 " %d %d\n",
	   revision_num,
	   revision_assignments[revision_num].topic,
	   revision_assignments[revision_num].pov);
  }
  close_mmaps(mmap_info);
}

max_stars_repo_name: allenlavoie/topic-pov

max_issues_repo_name: allenlavoie/topic-pov

max_forks_repo_name: allenlavoie/topic-pov