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

/* Author:  G. Jungman
 */

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_complex_math.h>
#include "gsl_eigen.h"

gsl_matrix *
create_hilbert_matrix(int size)
{
  int i, j;
  gsl_matrix * m = gsl_matrix_alloc(size, size);
  for(i=0; i<size; i++) {
    for(j=0; j<size; j++) {
      gsl_matrix_set(m, i, j, 1.0/(i+j+1.0));
    }
  }
  return m;
}

gsl_matrix *
create_random_symm_matrix(int size)
{
  int i, j;
  unsigned long k = 1;
  gsl_matrix * m = gsl_matrix_alloc(size, size);
  for(i=0; i<size; i++) {
    for(j=i; j<size; j++) {
      double x;
      k = (69069 * k + 1) & 0xffffffffUL;
      x = k / 4294967296.0;
      gsl_matrix_set(m, i, j, x);
      gsl_matrix_set(m, j, i, x);
    }
  }
  return m;
}

gsl_matrix_complex *
create_random_herm_matrix(int size)
{
  int i, j;
  unsigned long k = 1;
  gsl_matrix_complex * m = gsl_matrix_complex_alloc(size, size);
  for(i=0; i<size; i++) {
    for(j=i; j<size; j++) {
      gsl_complex z;
      k = (69069 * k + 1) & 0xffffffffUL;
      GSL_REAL(z) = k / 4294967296.0;
      k = (69069 * k + 1) & 0xffffffffUL;
      GSL_IMAG(z) = (i == j) ? 0 : k / 4294967296.0;
      gsl_matrix_complex_set(m, i, j, z);
      gsl_matrix_complex_set(m, j, i, gsl_complex_conjugate(z));
    }
  }
  return m;
}

void
test_eigen_results (size_t N, const gsl_matrix * m, const gsl_vector * eval,
                    const gsl_matrix * evec, const char * desc,
                    const char * desc2)
{
  size_t i,j;

  gsl_vector * x = gsl_vector_alloc(N);
  gsl_vector * y = gsl_vector_alloc(N);

  /* check eigenvalues */

  for (i = 0; i < N; i++)
    {
      double ei = gsl_vector_get (eval, i);
      gsl_vector_const_view vi = gsl_matrix_const_column(evec, i);
      gsl_vector_memcpy(x, &vi.vector);
      /* compute y = m x (should = lambda v) */
      gsl_blas_dgemv (CblasNoTrans, 1.0, m, x, 0.0, y);
      for (j = 0; j < N; j++)
        {
          double xj = gsl_vector_get (x, j);
          double yj = gsl_vector_get (y, j);
          gsl_test_rel(yj, ei * xj,  1e8 * GSL_DBL_EPSILON,
                       "%s, eigenvalue(%d,%d), %s", desc, i, j, desc2);
        }
    }

  /* check eigenvectors are orthonormal */

  for (i = 0; i < N; i++)
    {
      gsl_vector_const_view vi = gsl_matrix_const_column(evec, i);
      double nrm_v = gsl_blas_dnrm2(&vi.vector);
      gsl_test_rel (nrm_v, 1.0, N * GSL_DBL_EPSILON, "%s, normalized(%d), %s",
                    desc, i, desc2);
    }

  for (i = 0; i < N; i++)
    {
      gsl_vector_const_view vi = gsl_matrix_const_column(evec, i);
      for (j = i + 1; j < N; j++)
        {
          gsl_vector_const_view vj = gsl_matrix_const_column(evec, j);
          double vivj;
          gsl_blas_ddot (&vi.vector, &vj.vector, &vivj);
          gsl_test_abs (vivj, 0.0, N * GSL_DBL_EPSILON,
                        "%s, orthogonal(%d,%d), %s", desc, i, j, desc2);
        }
    }

  gsl_vector_free(x);
  gsl_vector_free(y);
}

void
test_eigenvalues (size_t N, const gsl_vector *eval, const gsl_vector * eval2,
                  const char * desc, const char * desc2)
{
  size_t i;
  for (i = 0; i < N; i++)
    {
      double ei = gsl_vector_get (eval, i);
      double e2i = gsl_vector_get (eval2, i);
      gsl_test_rel(ei, e2i, GSL_DBL_EPSILON, "%s, direct eigenvalue(%d), %s",
                   desc, i, desc2);
    }
}

void
test_eigen_complex_results (size_t N, const gsl_matrix_complex * m,
                            const gsl_vector * eval,
                            const gsl_matrix_complex * evec,
                            const char * desc,
                            const char * desc2)
{
  size_t i,j;

  gsl_vector_complex * x = gsl_vector_complex_alloc(N);
  gsl_vector_complex * y = gsl_vector_complex_alloc(N);

  /* check eigenvalues */

  for (i = 0; i < N; i++)
    {
      double ei = gsl_vector_get (eval, i);
      gsl_vector_complex_const_view vi = gsl_matrix_complex_const_column(evec, i);
      gsl_vector_complex_memcpy(x, &vi.vector);
      /* compute y = m x (should = lambda v) */
      gsl_blas_zgemv (CblasNoTrans, GSL_COMPLEX_ONE, m, x,
                      GSL_COMPLEX_ZERO, y);
      for (j = 0; j < N; j++)
        {
          gsl_complex xj = gsl_vector_complex_get (x, j);
          gsl_complex yj = gsl_vector_complex_get (y, j);
          gsl_test_rel(GSL_REAL(yj), ei * GSL_REAL(xj), 1e8*GSL_DBL_EPSILON,
                       "%s, eigenvalue(%d,%d), real, %s", desc, i, j, desc2);
          gsl_test_rel(GSL_IMAG(yj), ei * GSL_IMAG(xj), 1e8*GSL_DBL_EPSILON,
                       "%s, eigenvalue(%d,%d), imag, %s", desc, i, j, desc2);
        }
    }

  /* check eigenvectors are orthonormal */

  for (i = 0; i < N; i++)
    {
      gsl_vector_complex_const_view vi = gsl_matrix_complex_const_column(evec, i);
      double nrm_v = gsl_blas_dznrm2(&vi.vector);
      gsl_test_rel (nrm_v, 1.0, N * GSL_DBL_EPSILON, "%s, normalized(%d), %s",
                    desc, i, desc2);
    }

  for (i = 0; i < N; i++)
    {
      gsl_vector_complex_const_view vi = gsl_matrix_complex_const_column(evec, i);
      for (j = i + 1; j < N; j++)
        {
          gsl_vector_complex_const_view vj
            = gsl_matrix_complex_const_column(evec, j);
          gsl_complex vivj;
          gsl_blas_zdotc (&vi.vector, &vj.vector, &vivj);
          gsl_test_abs (gsl_complex_abs(vivj), 0.0, N * GSL_DBL_EPSILON,
                        "%s, orthogonal(%d,%d), %s", desc, i, j, desc2);
        }
    }

  gsl_vector_complex_free(x);
  gsl_vector_complex_free(y);
}

void
test_eigen_symm(const char * desc, const gsl_matrix * m)
{
  size_t N = m->size1;

  gsl_matrix * A = gsl_matrix_alloc(N, N);
  gsl_matrix * evec = gsl_matrix_alloc(N, N);
  gsl_vector * eval = gsl_vector_alloc(N);
  gsl_vector * eval2 = gsl_vector_alloc(N);

  gsl_eigen_symm_workspace * w1 = gsl_eigen_symm_alloc (N);
  gsl_eigen_symmv_workspace * w2 = gsl_eigen_symmv_alloc (N);

  gsl_matrix_memcpy(A, m);
  gsl_eigen_symmv(A, eval, evec, w2);
  test_eigen_results (N, m, eval, evec, desc, "unsorted");

  gsl_matrix_memcpy(A, m);
  gsl_eigen_symm(A, eval2, w1);
  test_eigenvalues (N, eval, eval2, desc, "unsorted");

  gsl_eigen_symmv_sort (eval, evec, GSL_EIGEN_SORT_VAL_ASC);
  test_eigen_results (N, m, eval, evec, desc, "val/asc");

  gsl_eigen_symmv_sort (eval, evec, GSL_EIGEN_SORT_VAL_DESC);
  test_eigen_results (N, m, eval, evec, desc, "val/desc");

  gsl_eigen_symmv_sort (eval, evec, GSL_EIGEN_SORT_ABS_ASC);
  test_eigen_results (N, m, eval, evec, desc, "abs/asc");

  gsl_eigen_symmv_sort (eval, evec, GSL_EIGEN_SORT_VAL_DESC);
  test_eigen_results (N, m, eval, evec, desc, "abs/desc");

  gsl_eigen_symm_free (w1);
  gsl_eigen_symmv_free (w2);

  gsl_matrix_free(A);
  gsl_matrix_free(evec);
  gsl_vector_free(eval);
  gsl_vector_free(eval2);
}

void
test_eigen_herm(const char * desc, const gsl_matrix_complex * m)
{
  size_t N = m->size1;

  gsl_matrix_complex * A = gsl_matrix_complex_alloc(N, N);
  gsl_matrix_complex * evec = gsl_matrix_complex_alloc(N, N);
  gsl_vector * eval = gsl_vector_alloc(N);
  gsl_vector * eval2 = gsl_vector_alloc(N);

  gsl_eigen_herm_workspace * w1 = gsl_eigen_herm_alloc (N);
  gsl_eigen_hermv_workspace * w2 = gsl_eigen_hermv_alloc (N);

  gsl_matrix_complex_memcpy(A, m);
  gsl_eigen_hermv(A, eval, evec, w2);
  test_eigen_complex_results (N, m, eval, evec, desc, "unsorted");

  gsl_matrix_complex_memcpy(A, m);
  gsl_eigen_herm(A, eval2, w1);
  test_eigenvalues (N, eval, eval2, desc, "unsorted");

  gsl_eigen_hermv_sort (eval, evec, GSL_EIGEN_SORT_VAL_ASC);
  test_eigen_complex_results (N, m, eval, evec, desc, "val/asc");

  gsl_eigen_hermv_sort (eval, evec, GSL_EIGEN_SORT_VAL_DESC);
  test_eigen_complex_results (N, m, eval, evec, desc, "val/desc");

  gsl_eigen_hermv_sort (eval, evec, GSL_EIGEN_SORT_ABS_ASC);
  test_eigen_complex_results (N, m, eval, evec, desc, "abs/asc");

  gsl_eigen_hermv_sort (eval, evec, GSL_EIGEN_SORT_VAL_DESC);
  test_eigen_complex_results (N, m, eval, evec, desc, "abs/desc");

  gsl_eigen_herm_free (w1);
  gsl_eigen_hermv_free (w2);

  gsl_matrix_complex_free(A);
  gsl_matrix_complex_free(evec);
  gsl_vector_free(eval);
  gsl_vector_free(eval2);
}

void
test_eigen_jacobi(const char * desc, const gsl_matrix * m)
{
  size_t N = m->size1;
  unsigned int nrot;

  gsl_matrix * A = gsl_matrix_alloc(N, N);
  gsl_matrix * evec = gsl_matrix_alloc(N, N);
  gsl_vector * eval = gsl_vector_alloc(N);

  gsl_matrix_memcpy(A, m);
  gsl_eigen_jacobi(A, eval, evec, 1000, &nrot);
  gsl_eigen_symmv_sort(eval, evec, GSL_EIGEN_SORT_VAL_ASC);

  test_eigen_results (N, m, eval, evec, desc, "");

  gsl_matrix_free(A);
  gsl_matrix_free(evec);
  gsl_vector_free(eval);
}

int test_invert_jacobi(void)
{
  int s = 0;
  int i, j;
  gsl_matrix * hminv = gsl_matrix_alloc(10, 10);
  gsl_matrix * id    = gsl_matrix_alloc(10, 10);

  /* 10x10 Hilbert matrix */
  gsl_matrix * hm = create_hilbert_matrix(10);
  gsl_eigen_invert_jacobi(hm, hminv, 1000);

  /* gsl_linalg_matmult(hm, hminv, id); */
  gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, hm, hminv, 0.0, id);

  for(i=0; i<10; i++) {
    for(j=0; j<10; j++) {
      double delta_ij = ( i == j ? 1.0 : 0.0 );
      double id_ij    = gsl_matrix_get(id, i, j);
      int rs = ( fabs(id_ij - delta_ij) > 5.0e-3 );
      s += rs;
    }
  }

  gsl_matrix_free(hm);
  gsl_matrix_free(hminv);
  gsl_matrix_free(id);

  return s;
}

int main()
{
  gsl_matrix *rs10 = create_random_symm_matrix (10);
  gsl_matrix_complex *rh10 = create_random_herm_matrix (10);

  double r[] =  { 0,  0, -1,  0,
                  0,  1,  0,  1,
                 -1,  0,  0,  0,
                  0,  1,  0,  0 };

  double c[] =  { 0,0,  0,0, -1,0,  0,0,
                  0,0,  1,0,  0,0,  1,0,
                 -1,0,  0,0,  0,0,  0,0,
                  0,0,  1,0,  0,0,  0,0 };

  gsl_matrix_view s4 = gsl_matrix_view_array (r, 4, 4);
  gsl_matrix_complex_view h4 = gsl_matrix_complex_view_array (c, 4, 4);

  gsl_ieee_env_setup ();

  test_eigen_symm("symm(4)", &s4.matrix);
  test_eigen_herm("herm(4)", &h4.matrix);

  test_eigen_symm("symm(10)", rs10);
  test_eigen_herm("herm(10)", rh10);

  /* gsl_matrix *h5 = create_hilbert_matrix (5); */
  /* test_eigen_jacobi("hilbert(5)", h5); */
  /* test_invert_jacobi(); */
  /* gsl_matrix_free (h5); */

  gsl_matrix_complex_free (rh10);
  gsl_matrix_free (rs10);

  exit (gsl_test_summary());
}


text:
#ifndef PROBABILITY_MODEL_H
#define PROBABILITY_MODEL_H

#include <vector>
#include <gsl/gsl_sf_gamma.h>
#include "data_input.h"
using namespace std;

class probability_model {

public:
	probability_model( unsigned long int end, unsigned long int diff);
    virtual ~probability_model();
    virtual double calculate_log_likelihood( const unsigned long int & position1, const unsigned long int & position2 ) = 0;
    virtual double calculate_log_likelihood( const vector< double > & ) = 0;
    unsigned long int get_diff() {return m_diff;}
    virtual void print_data() = 0;
    //virtual double calculate_total_log_likelihood() = 0;
    void get_cumulative_sufficient_data(const unsigned long int & time, vector< double > & sufficient_data);
    virtual void get_number_of_observations(const vector< double > & sufficient_stats, double & number_of_observations) = 0; // calculate the number of observations given a set of sufficient statistics
    unsigned long int m_end;
    unsigned long int m_diff;
    data_input * m_data;
};

probability_model::probability_model( unsigned long int end, unsigned long int diff ):m_end( end ), m_diff( diff){

}

probability_model::~probability_model(){

}

void probability_model::get_cumulative_sufficient_data(const unsigned long int & time, vector< double > & sufficient_data){
    m_data->get_cumulative_data_row(time, sufficient_data);
}

#endif

ext: h


text:
# ifndef CELL_CYCLE_H
# define CELL_CYCLE_H

#include <vector>
#include <deque>
#include <string>
#include <memory>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <math.h>

using namespace std;

#define CONSTANT_PHASE 0
#define EXPONENTIAL_PHASE 1

#endif

//////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////
class cell_cycle{//Class containing behaviour related to cell_cycle.
////////////
//////////

		private:

		int cell_state; // current state of the cycle
		double time_state; // current time in the state of the cycle
		double time_cycle; // total time in this cell_cycle
		double currentPhaseDuration; // time for the next phase
		vector<double> state_durations; // array with the cell states durations
		vector<int> state_types; // array with the cell phases types (exponential, uniform, ... )
		int divisions; // number of divisions so far
		int num_states; // total number of cell phases included
		bool divide_material; // flag to see if the species in the cell are dividied after cell division
		gsl_rng * rng; // allocator for the rng generator

		public:

		cell_cycle(){ // default constructor is a constructor with empty info
			time_state = 0;
			time_cycle = 0;
			divisions = 0;
			cell_state = 0;
			num_states = 0;
			rng= gsl_rng_alloc (gsl_rng_mt19937);
			gsl_rng_set (rng,::time(NULL)*getpid());
		}

		void addCellPhase(double duration, int type_phase){
			cout<<"Adding Cell Phase of duration"<< duration<<" \n";
			state_durations.push_back(duration);
			state_types.push_back(type_phase);
			num_states++;
			if (num_states == 1){// if this is the first state set, also use it for the current cell cycle time
				getPhaseDuration();
			}
			cout<<"Cell Phase added\n";
		}

		void ResetCellCycle(){
			time_state = 0;
			time_cycle = 0;
			cell_state = 0;
			getPhaseDuration();
			divisions = 0;
			divide_material = false;
		}

		// cell_cycle(int nstates,float duration, int type){
		// 	// The default constructor generates n cycle phases with same duration of certain type
		// 	time_state = 0;
		// 	time_cycle = 0;
		// 	divisions = 0;
		// 	cell_state = 0;
		// 	num_states = nstates;
		// 	state_durations = new double[num_states];
		// 	state_types = new int[num_states];
		// 	for(int i=0;i<nstates;i++){
		// 		state_durations[i] = duration;
		// 		state_types[i] = type;
		// 	}
		// }

		void setPhaseDuration(int phase, double duration, int type_phase){
			state_durations[phase] = duration;
			state_types[phase] = type_phase;
		}

		double getPhaseDuration(){
		// compute the value of the time for the next phase and store it in NextPhaseTime to avoid recomputing
			if (state_types[cell_state] == CONSTANT_PHASE){
				currentPhaseDuration = state_durations[cell_state];
			}
			else if (state_types[cell_state] == EXPONENTIAL_PHASE){
				currentPhaseDuration = state_durations[cell_state] * log(1.0/gsl_rng_uniform_pos(rng));
			}

			// cout<<"Next Phase duration: "<<currentPhaseDuration<<'\n';
			return currentPhaseDuration;
		}

		double getTimeForNextCellPhase(){
			return currentPhaseDuration - time_state;
		}

		int AdvanceNextPhase(){
			cell_state++;
			if (cell_state==num_states){// if the cell reaches the final state
				cell_state = 0;
				time_cycle = 0;
			}
			time_state = 0;
			getPhaseDuration();
			if (cell_state == 0){ // if we are dividing
				divisions += 1;
				divide_material = true;
			}
			else{
				divide_material = false;
			}
			return cell_state;
		}

		bool MaterialIsDivided(){
			return divide_material;
		}

		void updateCelltime(double time){
			time_cycle += time;
			time_state += time;
		}

		void setCellCyclePositionTime(int phase, double time){
			cell_state = phase;
			time_state = time;
			time_cycle = time;
		}

		void printCellProperties(){
			cout<<"total states "<<num_states<<'\n';
			cout<<"of duration ";
			for (auto d: state_durations){ cout << d << ' ';};
			cout<<"\n of type";
			for (auto d: state_types){ cout << d << ' ';};
			cout<<'\n';
		}

		int getCellPhase(){
			return cell_state;
		}

		double getPhaseTime(){
			return time_state;
		}

		void PrintState(){
			cout<< "cell state: "<<cell_state<< " time_state: "<<time_state<<'\n';
			cout<< "currentPhaseDuration: "<<currentPhaseDuration<<'\n';
			printCellProperties();
			cout<< "time cell has divided: "<<divisions<<'\n';
			cout<< "cell is ready for division: "<< divide_material<<'\n';
		}

};


text:
#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_sum.h>

#define N 20

int
main (void)
{
  double t[N];
  double sum_accel, err;
  double sum = 0;
  int n;

  gsl_sum_levin_u_workspace * w
    = gsl_sum_levin_u_alloc (N);

  const double zeta_2 = M_PI * M_PI / 6.0;

  /* terms for zeta(2) = \sum_{n=1}^{\infty} 1/n^2 */

  for (n = 0; n < N; n++)
    {
      double np1 = n + 1.0;
      t[n] = 1.0 / (np1 * np1);
      sum += t[n];
    }

  gsl_sum_levin_u_accel (t, N, w, &sum_accel, &err);

  printf ("term-by-term sum = % .16f using %d terms\n",
          sum, N);

  printf ("term-by-term sum = % .16f using %zu terms\n",
          w->sum_plain, w->terms_used);

  printf ("exact value      = % .16f\n", zeta_2);
  printf ("accelerated sum  = % .16f using %zu terms\n",
          sum_accel, w->terms_used);

  printf ("estimated error  = % .16f\n", err);
  printf ("actual error     = % .16f\n",
          sum_accel - zeta_2);

  gsl_sum_levin_u_free (w);
  return 0;
}

ext: c

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#pragma once

#include <memory>
#include <string>

#include <gsl/gsl>

#include "chainerx/backend.h"
#include "chainerx/device.h"
#include "chainerx/kernel_registry.h"

namespace chainerx {
namespace native {

class NativeDevice;
class NativeBackend;

namespace native_internal {

// Creates a device instance.
// This function is meant to be used from the backend class. Never use it for other purpose.
// This is defined in internal namespace in order to make it a friend of NativeDevice
// class.
NativeDevice* CreateDevice(NativeBackend& backend, int index);

}  // namespace native_internal

class NativeBackend : public Backend {
public:
    static constexpr const char* kDefaultName = "native";

    using Backend::Backend;

    std::string GetName() const override;

    int GetDeviceCount() const override;

    bool SupportsTransfer(Device& src_device, Device& dst_device) override;

    static KernelRegistry& GetGlobalKernelRegistry() {
        static gsl::owner<KernelRegistry*> global_kernel_registry = new KernelRegistry{};
        return *global_kernel_registry;
    }

protected:
    KernelRegistry& GetParentKernelRegistry() override { return GetGlobalKernelRegistry(); }

private:
    std::unique_ptr<Device> CreateDevice(int index) override;
};

}  // namespace native
}  // namespace chainerx

max_stars_repo_name: tkerola/chainer

max_issues_repo_name: tkerola/chainer

max_forks_repo_name: tkerola/chainer


text:
#ifndef _fim_h_included_
#define _fim_h_included_

#include <assert.h>
#include <stdlib.h>
#include <algorithm>
#include <cctype>
#include <fstream>
#include <functional>
#include <iostream>
#include <map>
#include <numeric>
#include <set>
#include <sstream>
#include <vector>
#include <sstream>
#include <string>
#include <boost/config.hpp>
#include <boost/foreach.hpp>
#include <boost/random.hpp>
#include <boost/shared_ptr.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>

#include <boost/numeric/ublas/io.hpp>

#include <jsc/bioinfo/gene_anno.hpp>
#include <jsc/util/log.hpp>

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

#include "splicing_graph.h"
#include "read.h"
#include "linalg.h"
#include "accessible_read_starts.h"

using namespace std;
using namespace boost;
using namespace boost::numeric;
using namespace jsc::bioinfo;
using namespace jsc::util;

class FIM {
	public:
		FIM () {
			ofim_call_count = 0;
		}

		double get_ofim_call_count() const {
			return ofim_call_count;
		};

		static double estimate_mle_variance(
				ublas::matrix<double> const & I) {
			return estimate_mle_variance_by_diag(I);
		}

		static double estimate_mle_variance_by_diag(
				ublas::matrix<double> const & I) {
			unsigned long K = I.size1() + 1;
			double sum = 0;
			for (unsigned long p = 0; p < K - 1; ++p) {
				sum += 1.0 / I(p,p);
			}
			return sum;
		}

		static double estimate_mle_variance_by_inv(
				ublas::matrix<double> const & I) {
			assert(I.size1() == I.size2());
			L_(debug) << I;
			ublas::matrix<double> inv =
				ublas::zero_matrix<double>(I.size1(), I.size1());
			linalg::invert_matrix(I, inv);

			L_(debug) << inv;

			unsigned long K = I.size1() + 1;
			double sum = 0;
			for (unsigned long p = 0; p < K - 1; ++p) {
				for (unsigned long q = 0; q < K - 1; ++q) {
					sum += inv(p,q);
				}
				sum += inv(p,p);
			}

			return sum;
		}

		static void print_fim_for_R(ostream & os, ublas::matrix<double> I) {
			unsigned long K = I.size1() + 1;
			for (unsigned long p = 0; p < K - 1; ++p) {
				for (unsigned long q = 0; q < K - 1; ++q) {
					os << I(p,q) << "\t";
				}
				os << endl;
			}
		}

		static void print_fim_diag_for_R(ostream & os, ublas::matrix<double> I) {
			unsigned long K = I.size1() + 1;
			os << "(";
			for (unsigned long p = 0; p < K - 1; ++p) {
				os << I(p,p) << ",";
			}
			os << ")";
			os << endl;
		}

		ublas::matrix<double> bruteforce_fim(shared_ptr<Read> readp,
				vector<double> const & iso_probs) {
			L_(debug) << "FIM::bruteforce_fim";
			reset_ofim_call_count();
			shared_ptr<AccessibleReadStarts> ars = readp->get_ARS();
			unsigned long K = ars->get_num_isoforms();
			ublas::matrix<double> I = ublas::zero_matrix<double>(K-1, K-1);

			for (unsigned long k = 0; k < K; ++k) {
				if (iso_probs[k] == 0) {
					continue;
				}
				ublas::matrix<int> sign;
				L_(debug) << "Possible Isoform #" << k
					<< ": Generating reads with starts from 0 to "
					<< ars->get_iso_ARS_total_length(k)
					<< " in ARS";
				for (unsigned long a = 0;
						a < ars->get_iso_ARS_total_length(k);
						++a) {
					// generate a read starting from a in ARS
					readp->generate_read(k, ars->ARStart2IsoStart(k, a));
					double G = readp->prob_generated_by_iso(k);
					double log_scaler = log(G) + log(iso_probs[k]);
					ublas::matrix<double> m = ofim(readp, iso_probs, sign);
					for (unsigned long p = 0; p < K - 1; ++p) {
						for (unsigned long q = 0; q < K - 1; ++q) {
							if (sign(p,q) > 0) {
								m(p,q) = exp(m(p,q) + log_scaler);
							} else if (sign(p,q) < 0) {
								m(p,q) = - exp(m(p,q) + log_scaler);
							} else {
								// m(p,q) must be zero, no need to perform scaling
							}
							assert(!isinf(m(p,q)) && !isnan(m(p,q)));
						}
					}
					I += m;
				}
			}
			return I;
		}

		ublas::matrix<double> fast_fim(shared_ptr<Read> readp,
				vector<double> const & iso_probs) {
			L_(debug) << "FIM::fast_fim";
			reset_ofim_call_count();
			shared_ptr<AccessibleReadStarts> ars = readp->get_ARS();
			unsigned long K = ars->get_num_isoforms();
			ublas::matrix<double> I = ublas::zero_matrix<double>(K-1, K-1);
			for (unsigned long k = 0; k < K; ++k) {
				if (iso_probs[k] == 0) {
					continue;
				}
				ublas::matrix<int> sign;
				L_(debug) << "Possible Isoform #" << k
					<< ": Generating reads with starts from 0 to "
					<< ars->get_iso_ARS_total_length(k)
					<< " in ARS";

				unsigned long a = 0;
				while (a < ars->get_iso_ARS_total_length(k)) {
					// generate a read starting from a
					readp->generate_read(k, ars->ARStart2IsoStart(k, a));

					unsigned long N_overlapping_exons = readp->exon_indices.size();
					unsigned long N_eq_samples = min(
							ars->exon_lengths[readp->exon_indices[0]]
							- readp->start_at_first_exon,
							ars->exon_lengths[readp->exon_indices[N_overlapping_exons - 1]]
							- readp->end_at_last_exon + 1);

					while ( (a + N_eq_samples - 1 >= ars->get_iso_ARS_total_length(k)) || (ars->ARStart2IsoStart(k, a + N_eq_samples - 1) != (ars->ARStart2IsoStart(k, a) + N_eq_samples - 1)) ){
						N_eq_samples--;
					}

					double G = readp->prob_generated_by_iso(k);
					ublas::matrix<double> m = ofim(readp, iso_probs, sign);
					bool zero_scaler = ((G == 0) || (iso_probs[k] == 0));
					double log_scaler = log(N_eq_samples)
						+ log(G)
						+ log(iso_probs[k]);
					for (unsigned long p = 0; p < K - 1; ++p) {
						for (unsigned long q = 0; q < K - 1; ++q) {
							if (zero_scaler) {
								m(p,q) = 0;
							} else if (sign(p,q) > 0) {
								m(p,q) = exp(m(p,q) + log_scaler);
							} else if (sign(p,q) < 0) {
								m(p,q) = - exp(m(p,q) + log_scaler);
							} else {
								// m(p,q) must be zero, no need to perform scaling
							}
							assert(!isinf(m(p,q)) && !isnan(m(p,q)));
						}
					}
					I += m;
					a += N_eq_samples;
				}
			}
			return I;
		}

		ublas::matrix<double> faster_fim(shared_ptr<Read> readp,
				vector<double> const & iso_probs) {
			L_(debug) << "FIM::faster_fim";
			reset_ofim_call_count();
			shared_ptr<AccessibleReadStarts> ars = readp->get_ARS();
			unsigned long K = ars->get_num_isoforms();
			ublas::matrix<double> I = ublas::zero_matrix<double>(K-1, K-1);
			vector<interval_list<long > > covered_sample_starts;
			for (unsigned long k = 0; k < K; ++k) {
				interval_list<long> il;
				covered_sample_starts.push_back(il);
			}
			for (unsigned long k = 0; k < K; ++k) {
				if (iso_probs[k] == 0) {
					continue;
				}

				unsigned long a = covered_sample_starts[k].find_min_uncovered(0);
				L_(debug) << "Possible Isoform #" << k
					<< ": Generating reads with starts from " << a << " to "
					<< ars->get_iso_ARS_total_length(k)
					<< " in ARS";
				while (a < ars->get_iso_ARS_total_length(k)) {
					// generate a read starting from a
					readp->generate_read(k, ars->ARStart2IsoStart(k, a));

					unsigned long N_overlapping_exons = readp->exon_indices.size();
					unsigned long N_eq_samples = min(
							ars->exon_lengths[readp->exon_indices[0]]
							- readp->start_at_first_exon,
							ars->exon_lengths[readp->exon_indices[N_overlapping_exons - 1]]
							- readp->end_at_last_exon + 1);
					while ( (a + N_eq_samples - 1 >= ars->get_iso_ARS_total_length(k)) || (ars->ARStart2IsoStart(k, a + N_eq_samples - 1) != (ars->ARStart2IsoStart(k, a) + N_eq_samples - 1)) ){
						N_eq_samples--;
					}

					ublas::matrix<int> sign;
					ublas::matrix<double> m = ofim(readp, iso_probs, sign);

					double G = readp->prob_generated_by_iso(k);
					double scaler = G * N_eq_samples * iso_probs[k];

					if (covered_sample_starts[k].compute_overlap(a, a + N_eq_samples)
							> 0) {
						L_(debug2) << "Cov for iso # " << k
							<< ": "	<< covered_sample_starts[k];
						L_(error) << "Overlap > 0!";
					}
					covered_sample_starts[k].add_interval(a, a + N_eq_samples);
					L_(debug2) << "Adding [" << a
						<< "," << (a + N_eq_samples)
						<< ") to possible isoform #" << k;
					for (unsigned long j = 0; j < K; ++j) {
						if (j == k) continue;
						long first_matching_idx =
							is_connected_exons_compatible_with_isoform(
									readp->exon_indices,
									ars->iso_exon_indices[j]);
						if (first_matching_idx >= 0) {
							double G_j = readp->prob_generated_by_iso(j);
							scaler += G_j * N_eq_samples * iso_probs[j];
							unsigned long a2 = readp->start_at_first_exon
								+ ars->iso_exon_total_lengths[j][first_matching_idx]
								- ars->exon_lengths[readp->exon_indices[0]];
							a2 = ars->IsoStart2ARStart(j, a2);

							if (covered_sample_starts[j].compute_overlap(a2, a2 + N_eq_samples)
									> 0) {
								L_(debug2) << "Cov for iso # " << j
									<< ": "	<< covered_sample_starts[j];
								L_(error) << "Overlap > 0!";
							}
							covered_sample_starts[j].add_interval(a2, a2 + N_eq_samples);
							L_(debug2) << "Adding [" << a2
								<< "," << (a2 + N_eq_samples)
								<< ") to possible isoform #" << j;
						}
					}

					double log_scaler = log(scaler);
					for (unsigned long p = 0; p < K - 1; ++p) {
						for (unsigned long q = 0; q < K - 1; ++q) {
							if (sign(p,q) > 0) {
								m(p,q) = exp(m(p,q) + log_scaler);
							} else if (sign(p,q) < 0) {
								m(p,q) = - exp(m(p,q) + log_scaler);
							} else {
								// m(p,q) must be zero, no need to perform scaling
							}
							assert(!isinf(m(p,q)) && !isnan(m(p,q)));
						}
					}
					I += m;
					a = covered_sample_starts[k].find_min_uncovered(0);
				}
			}
			return I;
		}
	protected:
		void reset_ofim_call_count() {
			ofim_call_count = 0;
		};

		// \mathfrak{I}_{s}^{(m)}(\Theta)
		ublas::matrix<double> ofim(shared_ptr<const Read> readp,
				vector<double> const & iso_probs,
				ublas::matrix<int> & sign) {
			++ofim_call_count;
			shared_ptr<AccessibleReadStarts> ars = readp->get_ARS();
			unsigned long K = ars->get_num_isoforms();
			ublas::matrix<double> m = ublas::zero_matrix<double>(K-1,K-1);
			sign = ublas::zero_matrix<int>(K-1, K-1);
			vector<double> v_delta_G(K, 0);
			for (unsigned long k = 0; k < K; ++k) {
				if (readp->is_compatible_with_iso(k)) {
					v_delta_G[k] = readp->prob_generated_by_iso(k);
				}
			}

			double sum = 0;
			for (unsigned long k = 0; k < K; ++k) {
				if (v_delta_G[k] > 0 && iso_probs[k] > 0) {
					sum += iso_probs[k] * v_delta_G[k];
				}
			}

			assert(sum > 0);

			double log_prod = log(sum) + log(sum);
			for (unsigned long p = 0; p < K - 1; ++p) {
				for (unsigned long q = 0; q < K - 1; ++q) {
					m(p,q) = (v_delta_G[p] - v_delta_G[K - 1])*
						(v_delta_G[q] - v_delta_G[K - 1]);
					if (m(p,q) != 0) {
						sign(p,q) = m(p,q) > 0 ? 1 : -1;
						if (m(p,q) < 0) {
							m(p,q) = - m(p,q);
						}
						m(p,q) = log(m(p,q)) - log_prod;
					}
				}
			}

			if (readp->exon_indices[0] == 0
					&& readp->start_at_first_exon == 2730) {
				L_(info) << *readp << " " << m(0,0) << " " << sign(0,0);
			}

			return m;
		}

	private:
		double ofim_call_count;
};

#endif

max_stars_repo_name: gersteinlab/LESSeq

max_issues_repo_name: gersteinlab/LESSeq

max_forks_repo_name: gersteinlab/LESSeq


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

#include <config.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_permute_vector.h>
#include <gsl/gsl_linalg.h>

#include "apply_givens.c"

/* The purpose of this package is to speed up QR-decomposition for
   large matrices.  Because QR-decomposition is column oriented, but
   GSL uses a row-oriented matrix format, there can considerable
   speedup obtained by computing the LQ-decomposition of the
   transposed matrix instead.  This package provides LQ-decomposition
   and related algorithms.  */

/* Factorise a general N x M matrix A into
 *
 *   P A = L Q
 *
 * where Q is orthogonal (M x M) and L is lower triangular (N x M).
 * When A is rank deficient, r = rank(A) < n, then the permutation is
 * used to ensure that the lower n - r columns of L are zero and the first
 * l rows of Q form an orthonormal basis for the rows of A.
 *
 * Q is stored as a packed set of Householder transformations in the
 * strict upper triangular part of the input matrix.
 *
 * L is stored in the diagonal and lower triangle of the input matrix.
 *
 * P: column j of P is column k of the identity matrix, where k =
 * permutation->data[j]
 *
 * The full matrix for Q can be obtained as the product
 *
 *       Q = Q_k .. Q_2 Q_1
 *
 * where k = MIN(M,N) and
 *
 *       Q_i = (I - tau_i * v_i * v_i')
 *
 * and where v_i is a Householder vector
 *
 *       v_i = [1, m(i,i+1), m(i,i+2), ... , m(i,M)]
 *
 * This storage scheme is the same as in LAPACK.  See LAPACK's
 * dgeqpf.f for details.
 *
 */

int
gsl_linalg_PTLQ_decomp (gsl_matrix * A, gsl_vector * tau, gsl_permutation * p, int *signum, gsl_vector * norm)
{
  const size_t N = A->size1;
  const size_t M = A->size2;

  if (tau->size != GSL_MIN (M, N))
    {
      GSL_ERROR ("size of tau must be MIN(M,N)", GSL_EBADLEN);
    }
  else if (p->size != N)
    {
      GSL_ERROR ("permutation size must be N", GSL_EBADLEN);
    }
  else if (norm->size != N)
    {
      GSL_ERROR ("norm size must be N", GSL_EBADLEN);
    }
  else
    {
      size_t i;

      *signum = 1;

      gsl_permutation_init (p); /* set to identity */

      /* Compute column norms and store in workspace */

      for (i = 0; i < N; i++)
        {
          gsl_vector_view c = gsl_matrix_row (A, i);
          double x = gsl_blas_dnrm2 (&c.vector);
          gsl_vector_set (norm, i, x);
        }

      for (i = 0; i < GSL_MIN (M, N); i++)
        {
          /* Bring the column of largest norm into the pivot position */

          double max_norm = gsl_vector_get(norm, i);
          size_t j, kmax = i;

          for (j = i + 1; j < N; j++)
            {
              double x = gsl_vector_get (norm, j);

              if (x > max_norm)
                {
                  max_norm = x;
                  kmax = j;
                }
            }

          if (kmax != i)
            {
              gsl_matrix_swap_rows (A, i, kmax);
              gsl_permutation_swap (p, i, kmax);
              gsl_vector_swap_elements(norm,i,kmax);

              (*signum) = -(*signum);
            }

          /* Compute the Householder transformation to reduce the j-th
             column of the matrix to a multiple of the j-th unit vector */

          {
            gsl_vector_view c = gsl_matrix_subrow (A, i, i, M - i);
            double tau_i = gsl_linalg_householder_transform (&c.vector);

            gsl_vector_set (tau, i, tau_i);

            /* Apply the transformation to the remaining columns */

            if (i + 1 < N)
              {
                gsl_matrix_view m = gsl_matrix_submatrix (A, i +1, i, N - (i+1), M - i);
                gsl_linalg_householder_mh (tau_i, &c.vector, &m.matrix);
              }
          }

          /* Update the norms of the remaining columns too */

          if (i + 1 < M)
            {
              for (j = i + 1; j < N; j++)
                {
                  double x = gsl_vector_get (norm, j);

                  if (x > 0.0)
                    {
                      double y = 0;
                      double temp= gsl_matrix_get (A, j, i) / x;

                      if (fabs (temp) >= 1)
                        y = 0.0;
                      else
                        y = x * sqrt (1 - temp * temp);

                      /* recompute norm to prevent loss of accuracy */

                      if (fabs (y / x) < sqrt (20.0) * GSL_SQRT_DBL_EPSILON)
                        {
                          gsl_vector_view c = gsl_matrix_subrow (A, j, i + 1, M - (i + 1));
                          y = gsl_blas_dnrm2 (&c.vector);
                        }

                      gsl_vector_set (norm, j, y);
                    }
                }
            }
        }

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_PTLQ_decomp2 (const gsl_matrix * A, gsl_matrix * q, gsl_matrix * r, gsl_vector * tau, gsl_permutation * p, int *signum, gsl_vector * norm)
{
  const size_t N = A->size1;
  const size_t M = A->size2;

  if (q->size1 != M || q->size2 !=M)
    {
      GSL_ERROR ("q must be M x M", GSL_EBADLEN);
    }
  else if (r->size1 != N || r->size2 !=M)
    {
      GSL_ERROR ("r must be N x M", GSL_EBADLEN);
    }
  else if (tau->size != GSL_MIN (M, N))
    {
      GSL_ERROR ("size of tau must be MIN(M,N)", GSL_EBADLEN);
    }
  else if (p->size != N)
    {
      GSL_ERROR ("permutation size must be N", GSL_EBADLEN);
    }
  else if (norm->size != N)
    {
      GSL_ERROR ("norm size must be N", GSL_EBADLEN);
    }

  gsl_matrix_memcpy (r, A);

  gsl_linalg_PTLQ_decomp (r, tau, p, signum, norm);

  /* FIXME:  aliased arguments depends on behavior of unpack routine! */

  gsl_linalg_LQ_unpack (r, tau, q, r);

  return GSL_SUCCESS;
}

/* Solves the system x^T A = b^T using the P^T L Q  factorisation,

   z^T L = b^T Q^T

   x = P z;

   to obtain x. Based on SLATEC code. */

int
gsl_linalg_PTLQ_solve_T (const gsl_matrix * QR,
                         const gsl_vector * tau,
                         const gsl_permutation * p,
                         const gsl_vector * b,
                         gsl_vector * x)
{
  if (QR->size1 != QR->size2)
    {
      GSL_ERROR ("QR matrix must be square", GSL_ENOTSQR);
    }
  else if (QR->size2 != p->size)
    {
      GSL_ERROR ("matrix size must match permutation size", GSL_EBADLEN);
    }
  else if (QR->size2 != b->size)
    {
      GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
    }
  else if (QR->size1 != x->size)
    {
      GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
    }
  else
    {
      gsl_vector_memcpy (x, b);

      gsl_linalg_PTLQ_svx_T (QR, tau, p, x);

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_PTLQ_svx_T (const gsl_matrix * LQ,
                       const gsl_vector * tau,
                       const gsl_permutation * p,
                       gsl_vector * x)
{
  if (LQ->size1 != LQ->size2)
    {
      GSL_ERROR ("LQ matrix must be square", GSL_ENOTSQR);
    }
  else if (LQ->size2 != p->size)
    {
      GSL_ERROR ("matrix size must match permutation size", GSL_EBADLEN);
    }
  else if (LQ->size1 != x->size)
    {
      GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
    }
  else
    {
      /* compute sol = b^T Q^T */

      gsl_linalg_LQ_vecQT (LQ, tau, x);

      /* Solve  L^T x = sol, storing x inplace in sol */

      gsl_blas_dtrsv (CblasLower, CblasTrans, CblasNonUnit, LQ, x);

      gsl_permute_vector_inverse (p, x);

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_PTLQ_LQsolve_T (const gsl_matrix * Q, const gsl_matrix * L,
                           const gsl_permutation * p,
                           const gsl_vector * b,
                           gsl_vector * x)
{
  if (Q->size1 != Q->size2 || L->size1 != L->size2)
    {
      return GSL_ENOTSQR;
    }
  else if (Q->size1 != p->size || Q->size1 != L->size1
           || Q->size1 != b->size)
    {
      return GSL_EBADLEN;
    }
  else
    {
      /* compute b' = Q b */

      gsl_blas_dgemv (CblasNoTrans, 1.0, Q, b, 0.0, x);

      /* Solve L^T x = b', storing x inplace */

      gsl_blas_dtrsv (CblasLower, CblasTrans, CblasNonUnit, L, x);

      /* Apply permutation to solution in place */

      gsl_permute_vector_inverse (p, x);

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_PTLQ_Lsolve_T (const gsl_matrix * LQ,
                        const gsl_permutation * p,
                        const gsl_vector * b,
                        gsl_vector * x)
{
  if (LQ->size1 != LQ->size2)
    {
      GSL_ERROR ("LQ matrix must be square", GSL_ENOTSQR);
    }
  else if (LQ->size1 != b->size)
    {
      GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
    }
  else if (LQ->size2 != x->size)
    {
      GSL_ERROR ("matrix size must match x size", GSL_EBADLEN);
    }
  else if (p->size != x->size)
    {
      GSL_ERROR ("permutation size must match x size", GSL_EBADLEN);
    }
  else
    {
      /* Copy x <- b */

      gsl_vector_memcpy (x, b);

      /* Solve L^T x = b, storing x inplace */

      gsl_blas_dtrsv (CblasLower, CblasTrans, CblasNonUnit, LQ, x);

      gsl_permute_vector_inverse (p, x);

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_PTLQ_Lsvx_T (const gsl_matrix * LQ,
                        const gsl_permutation * p,
                        gsl_vector * x)
{
  if (LQ->size1 != LQ->size2)
    {
      GSL_ERROR ("LQ matrix must be square", GSL_ENOTSQR);
    }
  else if (LQ->size2 != x->size)
    {
      GSL_ERROR ("matrix size must match x size", GSL_EBADLEN);
    }
  else if (p->size != x->size)
    {
      GSL_ERROR ("permutation size must match x size", GSL_EBADLEN);
    }
  else
    {
      /* Solve L^T x = b, storing x inplace */

      gsl_blas_dtrsv (CblasLower, CblasTrans, CblasNonUnit, LQ, x);

      gsl_permute_vector_inverse (p, x);

      return GSL_SUCCESS;
    }
}

/* Update a P^T L Q factorisation for P A= L Q ,  A' =  A +  v u^T,
                                                 PA' = PA + Pv u^T

 * P^T L' Q' = P^T LQ + v u^T
 *       = P^T (L + (P v) u^T Q^T) Q
 *       = P^T (L + (P v) w^T) Q
 *
 * where w = Q^T u.
 *
 * Algorithm from Golub and Van Loan, "Matrix Computations", Section
 * 12.5 (Updating Matrix Factorizations, Rank-One Changes)
 */

int
gsl_linalg_PTLQ_update (gsl_matrix * Q, gsl_matrix * L,
                        const gsl_permutation * p,
                        const gsl_vector * v, gsl_vector * w)
{
  if (Q->size1 != Q->size2 || L->size1 != L->size2)
    {
      return GSL_ENOTSQR;
    }
  else if (L->size1 != Q->size2 || v->size != Q->size2 || w->size != Q->size2)
    {
      return GSL_EBADLEN;
    }
  else
    {
      size_t j, k;
      const size_t N = Q->size1;
      const size_t M = Q->size2;
      double w0;

      /* Apply Given's rotations to reduce w to (|w|, 0, 0, ... , 0)

         J_1^T .... J_(n-1)^T w = +/- |w| e_1

         simultaneously applied to L,  H = J_1^T ... J^T_(n-1) L
         so that H is upper Hessenberg.  (12.5.2) */

      for (k = M - 1; k > 0; k--)
        {
          double c, s;
          double wk = gsl_vector_get (w, k);
          double wkm1 = gsl_vector_get (w, k - 1);

          gsl_linalg_givens (wkm1, wk, &c, &s);
          gsl_linalg_givens_gv (w, k - 1, k, c, s);
          apply_givens_lq (M, N, Q, L, k - 1, k, c, s);
        }

      w0 = gsl_vector_get (w, 0);

      /* Add in v w^T  (Equation 12.5.3) */

      for (j = 0; j < N; j++)
        {
          double lj0 = gsl_matrix_get (L, j, 0);
          size_t p_j = gsl_permutation_get (p, j);
          double vj = gsl_vector_get (v, p_j);
          gsl_matrix_set (L, j, 0, lj0 + w0 * vj);
        }

      /* Apply Givens transformations L' = G_(n-1)^T ... G_1^T H
         Equation 12.5.4 */

      for (k = 1; k < N; k++)
        {
          double c, s;
          double diag = gsl_matrix_get (L, k - 1, k - 1);
          double offdiag = gsl_matrix_get (L, k - 1, k );

          gsl_linalg_givens (diag, offdiag, &c, &s);
          apply_givens_lq (M, N, Q, L, k - 1, k, c, s);
        }

      return GSL_SUCCESS;
    }
}

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


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

#ifndef __GSL_MATRIX_FLOAT_H__
#define __GSL_MATRIX_FLOAT_H__

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

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_inline.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_vector_float.h>

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

__BEGIN_DECLS

typedef struct
{
  size_t size1;
  size_t size2;
  size_t tda;
  float * data;
  gsl_block_float * block;
  int owner;
} gsl_matrix_float;

typedef struct
{
  gsl_matrix_float matrix;
} _gsl_matrix_float_view;

typedef _gsl_matrix_float_view gsl_matrix_float_view;

typedef struct
{
  gsl_matrix_float matrix;
} _gsl_matrix_float_const_view;

typedef const _gsl_matrix_float_const_view gsl_matrix_float_const_view;

/* Allocation */

GSL_FUN gsl_matrix_float *
gsl_matrix_float_alloc (const size_t n1, const size_t n2);

GSL_FUN gsl_matrix_float *
gsl_matrix_float_calloc (const size_t n1, const size_t n2);

GSL_FUN gsl_matrix_float *
gsl_matrix_float_alloc_from_block (gsl_block_float * b,
                                   const size_t offset,
                                   const size_t n1,
                                   const size_t n2,
                                   const size_t d2);

GSL_FUN gsl_matrix_float *
gsl_matrix_float_alloc_from_matrix (gsl_matrix_float * m,
                                    const size_t k1,
                                    const size_t k2,
                                    const size_t n1,
                                    const size_t n2);

GSL_FUN gsl_vector_float *
gsl_vector_float_alloc_row_from_matrix (gsl_matrix_float * m,
                                        const size_t i);

GSL_FUN gsl_vector_float *
gsl_vector_float_alloc_col_from_matrix (gsl_matrix_float * m,
                                        const size_t j);

GSL_FUN void gsl_matrix_float_free (gsl_matrix_float * m);

/* Views */

GSL_FUN _gsl_matrix_float_view
gsl_matrix_float_submatrix (gsl_matrix_float * m,
                            const size_t i, const size_t j,
                            const size_t n1, const size_t n2);

GSL_FUN _gsl_vector_float_view
gsl_matrix_float_row (gsl_matrix_float * m, const size_t i);

GSL_FUN _gsl_vector_float_view
gsl_matrix_float_column (gsl_matrix_float * m, const size_t j);

GSL_FUN _gsl_vector_float_view
gsl_matrix_float_diagonal (gsl_matrix_float * m);

GSL_FUN _gsl_vector_float_view
gsl_matrix_float_subdiagonal (gsl_matrix_float * m, const size_t k);

GSL_FUN _gsl_vector_float_view
gsl_matrix_float_superdiagonal (gsl_matrix_float * m, const size_t k);

GSL_FUN _gsl_vector_float_view
gsl_matrix_float_subrow (gsl_matrix_float * m, const size_t i,
                         const size_t offset, const size_t n);

GSL_FUN _gsl_vector_float_view
gsl_matrix_float_subcolumn (gsl_matrix_float * m, const size_t j,
                            const size_t offset, const size_t n);

GSL_FUN _gsl_matrix_float_view
gsl_matrix_float_view_array (float * base,
                             const size_t n1,
                             const size_t n2);

GSL_FUN _gsl_matrix_float_view
gsl_matrix_float_view_array_with_tda (float * base,
                                      const size_t n1,
                                      const size_t n2,
                                      const size_t tda);

GSL_FUN _gsl_matrix_float_view
gsl_matrix_float_view_vector (gsl_vector_float * v,
                              const size_t n1,
                              const size_t n2);

GSL_FUN _gsl_matrix_float_view
gsl_matrix_float_view_vector_with_tda (gsl_vector_float * v,
                                       const size_t n1,
                                       const size_t n2,
                                       const size_t tda);

GSL_FUN _gsl_matrix_float_const_view
gsl_matrix_float_const_submatrix (const gsl_matrix_float * m,
                                  const size_t i, const size_t j,
                                  const size_t n1, const size_t n2);

GSL_FUN _gsl_vector_float_const_view
gsl_matrix_float_const_row (const gsl_matrix_float * m,
                            const size_t i);

GSL_FUN _gsl_vector_float_const_view
gsl_matrix_float_const_column (const gsl_matrix_float * m,
                               const size_t j);

GSL_FUN _gsl_vector_float_const_view
gsl_matrix_float_const_diagonal (const gsl_matrix_float * m);

GSL_FUN _gsl_vector_float_const_view
gsl_matrix_float_const_subdiagonal (const gsl_matrix_float * m,
                                    const size_t k);

GSL_FUN _gsl_vector_float_const_view
gsl_matrix_float_const_superdiagonal (const gsl_matrix_float * m,
                                      const size_t k);

GSL_FUN _gsl_vector_float_const_view
gsl_matrix_float_const_subrow (const gsl_matrix_float * m, const size_t i,
                               const size_t offset, const size_t n);

GSL_FUN _gsl_vector_float_const_view
gsl_matrix_float_const_subcolumn (const gsl_matrix_float * m, const size_t j,
                                  const size_t offset, const size_t n);

GSL_FUN _gsl_matrix_float_const_view
gsl_matrix_float_const_view_array (const float * base,
                                   const size_t n1,
                                   const size_t n2);

GSL_FUN _gsl_matrix_float_const_view
gsl_matrix_float_const_view_array_with_tda (const float * base,
                                            const size_t n1,
                                            const size_t n2,
                                            const size_t tda);

GSL_FUN _gsl_matrix_float_const_view
gsl_matrix_float_const_view_vector (const gsl_vector_float * v,
                                    const size_t n1,
                                    const size_t n2);

GSL_FUN _gsl_matrix_float_const_view
gsl_matrix_float_const_view_vector_with_tda (const gsl_vector_float * v,
                                             const size_t n1,
                                             const size_t n2,
                                             const size_t tda);

/* Operations */

GSL_FUN void gsl_matrix_float_set_zero (gsl_matrix_float * m);
GSL_FUN void gsl_matrix_float_set_identity (gsl_matrix_float * m);
GSL_FUN void gsl_matrix_float_set_all (gsl_matrix_float * m, float x);

GSL_FUN int gsl_matrix_float_fread (FILE * stream, gsl_matrix_float * m) ;
GSL_FUN int gsl_matrix_float_fwrite (FILE * stream, const gsl_matrix_float * m) ;
GSL_FUN int gsl_matrix_float_fscanf (FILE * stream, gsl_matrix_float * m);
GSL_FUN int gsl_matrix_float_fprintf (FILE * stream, const gsl_matrix_float * m, const char * format);

GSL_FUN int gsl_matrix_float_memcpy(gsl_matrix_float * dest, const gsl_matrix_float * src);
GSL_FUN int gsl_matrix_float_swap(gsl_matrix_float * m1, gsl_matrix_float * m2);

GSL_FUN int gsl_matrix_float_swap_rows(gsl_matrix_float * m, const size_t i, const size_t j);
GSL_FUN int gsl_matrix_float_swap_columns(gsl_matrix_float * m, const size_t i, const size_t j);
GSL_FUN int gsl_matrix_float_swap_rowcol(gsl_matrix_float * m, const size_t i, const size_t j);
GSL_FUN int gsl_matrix_float_transpose (gsl_matrix_float * m);
GSL_FUN int gsl_matrix_float_transpose_memcpy (gsl_matrix_float * dest, const gsl_matrix_float * src);

GSL_FUN float gsl_matrix_float_max (const gsl_matrix_float * m);
GSL_FUN float gsl_matrix_float_min (const gsl_matrix_float * m);
GSL_FUN void gsl_matrix_float_minmax (const gsl_matrix_float * m, float * min_out, float * max_out);

GSL_FUN void gsl_matrix_float_max_index (const gsl_matrix_float * m, size_t * imax, size_t *jmax);
GSL_FUN void gsl_matrix_float_min_index (const gsl_matrix_float * m, size_t * imin, size_t *jmin);
GSL_FUN void gsl_matrix_float_minmax_index (const gsl_matrix_float * m, size_t * imin, size_t * jmin, size_t * imax, size_t * jmax);

GSL_FUN int gsl_matrix_float_isnull (const gsl_matrix_float * m);
GSL_FUN int gsl_matrix_float_ispos (const gsl_matrix_float * m);
GSL_FUN int gsl_matrix_float_isneg (const gsl_matrix_float * m);
GSL_FUN int gsl_matrix_float_isnonneg (const gsl_matrix_float * m);

GSL_FUN int gsl_matrix_float_add (gsl_matrix_float * a, const gsl_matrix_float * b);
GSL_FUN int gsl_matrix_float_sub (gsl_matrix_float * a, const gsl_matrix_float * b);
GSL_FUN int gsl_matrix_float_mul_elements (gsl_matrix_float * a, const gsl_matrix_float * b);
GSL_FUN int gsl_matrix_float_div_elements (gsl_matrix_float * a, const gsl_matrix_float * b);
GSL_FUN int gsl_matrix_float_scale (gsl_matrix_float * a, const double x);
GSL_FUN int gsl_matrix_float_add_constant (gsl_matrix_float * a, const double x);
GSL_FUN int gsl_matrix_float_add_diagonal (gsl_matrix_float * a, const double x);

/***********************************************************************/
/* The functions below are obsolete                                    */
/***********************************************************************/
GSL_FUN int gsl_matrix_float_get_row(gsl_vector_float * v, const gsl_matrix_float * m, const size_t i);
GSL_FUN int gsl_matrix_float_get_col(gsl_vector_float * v, const gsl_matrix_float * m, const size_t j);
GSL_FUN int gsl_matrix_float_set_row(gsl_matrix_float * m, const size_t i, const gsl_vector_float * v);
GSL_FUN int gsl_matrix_float_set_col(gsl_matrix_float * m, const size_t j, const gsl_vector_float * v);
/***********************************************************************/

/* inline functions if you are using GCC */

GSL_FUN INLINE_DECL float   gsl_matrix_float_get(const gsl_matrix_float * m, const size_t i, const size_t j);
GSL_FUN INLINE_DECL void    gsl_matrix_float_set(gsl_matrix_float * m, const size_t i, const size_t j, const float x);
GSL_FUN INLINE_DECL float * gsl_matrix_float_ptr(gsl_matrix_float * m, const size_t i, const size_t j);
GSL_FUN INLINE_DECL const float * gsl_matrix_float_const_ptr(const gsl_matrix_float * m, const size_t i, const size_t j);

#ifdef HAVE_INLINE
INLINE_FUN
float
gsl_matrix_float_get(const gsl_matrix_float * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_VAL("first index out of range", GSL_EINVAL, 0) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_VAL("second index out of range", GSL_EINVAL, 0) ;
        }
    }
#endif
  return m->data[i * m->tda + j] ;
}

INLINE_FUN
void
gsl_matrix_float_set(gsl_matrix_float * m, const size_t i, const size_t j, const float x)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_VOID("first index out of range", GSL_EINVAL) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_VOID("second index out of range", GSL_EINVAL) ;
        }
    }
#endif
  m->data[i * m->tda + j] = x ;
}

INLINE_FUN
float *
gsl_matrix_float_ptr(gsl_matrix_float * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
        }
    }
#endif
  return (float *) (m->data + (i * m->tda + j)) ;
}

INLINE_FUN
const float *
gsl_matrix_float_const_ptr(const gsl_matrix_float * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
        }
    }
#endif
  return (const float *) (m->data + (i * m->tda + j)) ;
}

#endif

__END_DECLS

#endif /* __GSL_MATRIX_FLOAT_H__ */

max_stars_repo_name: berkus/music-cs

max_issues_repo_name: berkus/music-cs

max_forks_repo_name: berkus/music-cs


text:
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_interp.h>

static double fun(double x) {
    return x + cos(x * x);
}

int main(void) {
    const double a = 1.0;
    const double b = 10.0;
    const int steps = 10;
    double xi, yi, x[100], y[100], dx;
    FILE *input, *output;
    int i;
    const gsl_interp_type *ts[] = {
            gsl_interp_linear,
            gsl_interp_polynomial,
            gsl_interp_cspline,
            gsl_interp_cspline_periodic,
            gsl_interp_akima,
            gsl_interp_akima_periodic,
            gsl_interp_steffen,
    };
    char *tsn[] = {
            "gsl_interp_linear.txt",
            "gsl_interp_polynomial.txt",
            "gsl_interp_cspline.txt",
            "gsl_interp_cspline_periodic.txt",
            "gsl_interp_akima.txt",
            "gsl_interp_akima_periodic.txt",
            "gsl_interp_steffen.txt",
    };
    int tsc = sizeof(ts) / sizeof(const gsl_interp_type *);
    input = fopen("wartosci.txt", "w");
    dx = (b - a) / (double) steps;
    for (i = 0; i <= steps; ++i) {
        x[i] = a + (double) i * dx + 0.5 * sin((double) i * dx);
        y[i] = fun(x[i]);
        fprintf(input, "%g %g\n", x[i], y[i]);
    }
    for (i = 0; i < tsc; i++) {
        output = fopen(tsn[i], "w");
        gsl_interp_accel *acc = gsl_interp_accel_alloc();
        gsl_spline *spline = gsl_spline_alloc(ts[i], steps + 1);
        gsl_spline_init(spline, x, y, steps + 1);
        for (xi = a; xi <= b; xi += 0.01) {
            yi = gsl_spline_eval(spline, xi, acc);
            fprintf(output, "%g %g\n", xi, yi);
        }
        gsl_spline_free(spline);
        gsl_interp_accel_free(acc);
    }
    return EXIT_SUCCESS;
}

max_stars_repo_name: mistyfiky/agh-mownit

max_issues_repo_name: mistyfiky/agh-mownit

max_forks_repo_name: mistyfiky/agh-mownit


text:
/**
* @license Apache-2.0
*
* Copyright (c) 2018 The Stdlib Authors.
*
* 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.
*/

#include <cblas.h>
#include <stdio.h>

int main() {
	double A[ 6 ] = { 1.0, 2.0, 1.0, -3.0, 4.0, -1.0 };
	double B[ 6 ] = { 1.0, 2.0, 1.0, -3.0, 4.0, -1.0 };
	double C[ 9 ] = { 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5 };
	int i = 0;

	cblas_dgemm( CblasColMajor, CblasNoTrans, CblasTrans, 3, 3, 2, 1, A, 3, B, 3, 2, C, 3 );

	for( i = 0; i < 9; i++ ) {
		printf( "%lf ", C[ i ] );
	}
	printf( "\n" );
}

max_stars_repo_name: ghalimi/stdlib

max_issues_repo_name: ghalimi/stdlib

max_forks_repo_name: sthagen/stdlib


text:
#ifndef CWANNIER_EVCACHE_H
#define CWANNIER_EVCACHE_H

#include <stdbool.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include "input.h"
#include "ecache.h"
#include "submesh.h"

typedef struct {
    int na;
    int nb;
    int nc;
    int num_bands;
    int G_order[3];
    int G_neg[3];
    UEInputFn UEfn;
    gsl_vector **energies;
    gsl_matrix_complex **evecs;
} EvecCache;

EvecCache* init_EvecCache(int na, int nb, int nc, int num_bands, int G_order[3], int G_neg[3], UEInputFn UEfn);

void free_EvecCache(EvecCache *evCache);

void EvecCache_MinMaxVals(EvecCache *evCache, double *emin, double *emax);

EnergyCache* copy_to_Ecache(EvecCache *evCache);

#endif //CWANNIER_EVCACHE_H

max_stars_repo_name: tflovorn/ctetra

max_issues_repo_name: tflovorn/ctetra

max_forks_repo_name: tflovorn/ctetra


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

#include <iostream>
#include <fstream>
#include <string>
#include <cmath>
#include <vector>
#include <algorithm>
#include "Rand.h"
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_sf_bessel.h>

using namespace std;

const int         XMAX = 600;
const long double EXPXMAX = exp((double)XMAX);
const long double EXP     = exp(1.0);

double Beta(double t, double x, double Itot);
double Gaussian(double x, double mean, double var);
long double Poisson(int lambda, int k);

void Statistics(const vector<double>& data, double& ave, double& squareave, double& var);

void Statistics(const vector<double>& data, double& max, double& min, double& ave, double& var, double& error);
void Statistics(const vector<int>& data, double& max, double& min, double& ave, double& var, double& error);
void Statistics(const vector<int>& data, double& xmax, double& xmin, double& xave, double& xvar, double& xerror, double& x2ave);

//void Statistics(const vector<int>& data, double& xmax, double& xmin, double& xave, double& xvar, double& xerror, double& x2ave, double& x1inverseave);

void Statistics(const vector<int>& data, double& Xmin, double& Xmax,
		double& X1inverseave, //  <1/(X+1)>
		double& logXave,      //  < logX >
		double& Xsrave,       //  < X^0.5 >
		double& Xave,         //  < X >
		double& X2ave );      //  < x^2 >

void moment(const vector<int>& data, double &ave, double &adev, double &sdev, double &var, double &skew, double &curt);

void StatisticalDispersion(vector<int>& X, double& H7, double& H8, double& H9, double& H10);

double GetXsecondmin(const vector<double>& data, double Xmin);

void Statistics(const vector<double>& data,
		double& min, double& Q1, double& Q2, double& Q3, double& max,
		double& ave, double& var, double& error,
		ofstream& fout1, ofstream& fout2);

void Binning(const vector<double>& data, int B, vector<double>& Binsize, vector<double>& Q);
int GetBinIndex(const double& x, const vector<double>& Q);

#endif

max_stars_repo_name: spxuw/RFIM

max_issues_repo_name: spxuw/RFIM

max_forks_repo_name: spxuw/RFIM


text:
#include <mpi.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sys/types.h>
#include <unistd.h>
#include <math.h>
#include <gsl/gsl_math.h>

#ifdef INVARIANCETEST
#define DO_NOT_REDEFINE_MPI_COMM_WORLD
#endif

#include "allvars.h"
#include "proto.h"

/*! \file main.c
 *  \brief start of the program
 */
/*!
 *  This function initializes the MPI communication packages, and sets
 *  cpu-time counters to 0. Then begrun() is called, which sets up
 *  the simulation either from IC's or from restart files.  Finally,
 *  run() is started, the main simulation loop, which iterates over
 *  the timesteps.
 */
int main(int argc, char **argv)
{
  int i;

  MPI_Init(&argc, &argv);
  MPI_Comm_rank(MPI_COMM_WORLD, &ThisTask);
  MPI_Comm_size(MPI_COMM_WORLD, &NTask);

#ifdef INVARIANCETEST
  World_ThisTask = ThisTask;
  World_NTask = NTask;

  if(World_NTask != (INVARIANCETEST_SIZE1 + INVARIANCETEST_SIZE2))
    {
      printf("wrong number of procs for invariance-test\n");
      MPI_Finalize();		/* clean up & finalize MPI */
      return 0;
    }

  if(World_ThisTask < INVARIANCETEST_SIZE1)
    Color = 0;
  else
    Color = 1;

  MPI_Comm_split(MPI_COMM_WORLD, Color, 0, &MPI_CommLocal);
  MPI_Comm_rank(MPI_CommLocal, &ThisTask);
  MPI_Comm_size(MPI_CommLocal, &NTask);
#endif

  for(PTask = 0; NTask > (1 << PTask); PTask++);

  if(argc < 2)
    {
      if(ThisTask == 0)
	{
	  printf("Parameters are missing.\n");
	  printf("Call with <ParameterFile> [<RestartFlag>] [<RestartSnapNum>]\n");
	}
      endrun(0);
    }

  strcpy(ParameterFile, argv[1]);

  if(argc >= 3)
    RestartFlag = atoi(argv[2]);
  else
    RestartFlag = 0;

  if(argc >= 4)
    RestartSnapNum = atoi(argv[3]);
  else
    RestartSnapNum = -1;

  /* initialize CPU-time/Wallclock-time measurement */
  for(i = 0; i < CPU_PARTS; i++)
    All.CPU_Sum[i] = CPU_Step[i] = 0;

  CPUThisRun = 0;
  WallclockTime = second();

  begrun();			/* set-up run  */

  run();			/* main simulation loop */

  MPI_Finalize();		/* clean up & finalize MPI */

  return 0;
}

/* ----------------------------------------------------------------------
   The rest of this file contains documentation for compiling and
   running the code, in a format appropriate for 'doxygen'.
   ----------------------------------------------------------------------
 */

/*! \mainpage Reference documentation of GADGET-2

\author Volker Springel \n
        Max-Planck-Institute for Astrophysics \n
        Garching, Germany \n
        volker@mpa-garching.mpg.de \n

\n

\section install Compilation

GADGET-2 needs the following non-standard libraries for compilation:

- \b MPI - the Message Passing Interface (version 1.0 or higher). Many
  vendor supplied versions exist, in addition to excellent open source
  implementations, e.g.  MPICH
  (http://www-unix.mcs.anl.gov/mpi/mpich/) or LAM
  (http://www.lam-mpi.org/).
- \b GSL - the GNU scientific library. This open-source package can be
  obtained at http://www.gnu.org/software/gsl, for example. GADGET
  only needs this library for a few very simple cosmological
  integrations at start-up.
- \b FFTW - the Fastest Fourier Transform in the West. This
  open-source package can be obtained at http://www.fftw.org,
  for example. It is only needed for simulations that use the TreePM algorithm.

Note that FFTW needs to be compiled with parallel support enabled.
This can be achieved by passing the option --enable-mpi to the
configure script. When at it, you might as well add
--enable-type-prefix to obtain the libraries in both a single
and double precision version. If this has not been done, you should
set the option \b NOTYPEPREFIX_FFTW in GADGET's \ref Gadget-Makefile
"Makefile".

Note that if any of the above libraries is not installed in standard
locations on your system, the \ref Gadget-Makefile "Makefile" provided
with the code may need slight adjustments. Similarly, compiler
options, particularly with respect to optimisations, may need
adjustment to the C-compiler that is used. Finally, the \ref
Gadget-Makefile "Makefile" contains a number of compile-time options
that need to be set appropriately for the type of simulation that is
simulated.

The provided makefile is compatible with GNU-make, i.e. typing \b make
or \b gmake should then build the executable P-Gadget2.  If your
site does not have GNU-make, get it, or write your own makefile.

\section howtorun Running the code

In order to start the code, a \ref parameterfile "parameterfile" needs
to be specified.

*/

/*! \page parameterfile  Parameterfile of GADGET-2

- \b ViscositySourceScaling \n Parameter for dynamic viscosity (e.g. 2.5)

- \b ViscosityDecayLength \n Parameter for dynamic viscosity (e.g. 2.0)

- \b ViscosityAlphaMin \n  Parameter for dynamic viscosity (e.g. 0.01)

- \b ViscosityAlphaMin \n  value for w of equation of state (e.g. -0.4)

- \b DarkEnergyFile \n  file containing z and w(z) (e.g. wz_0.83_RP.txt)

- \b VelIniScale \n  factor to rescale v_ini (e.g. 1.10)

*/

/*! \page Gadget-Makefile  Makefile of GADGET-2

There are a number of new types of simulations that can be run with
GADGET-2, most importantly some that use different variants of the
TreePM algorithm. A schematic overview of these simulation-types is
given in Table 1.

The TreePM algorithm is switched on by passing the desired PM
mesh-size at compile time via the makefile to the code. The relevant
parameter is \b PMGRID, see below. Using an explicit force split, the
long-range force is then computed with Fourier techniques, while the
short-range force is done with the tree. Because the tree needs only
be walked locally, a speed-up can arise, particularly for near to
homogeneous particle distributions, but not only restricted to
them. Periodic and non-periodic boundary conditions are
implemented. In the latter case, the code will internally compute FFTs
of size \b 2*PMGRID in order to allow for the required zero-padding.

Pure SPH simulations can also be run in periodic boxes whose
x-dimension is an integer multiple of the dimensions in the y- and
z-directions. This can be useful for hydrodynamic test calculations,
e.g. shock tubes in 3D.

Many aspects of these new features are controlled with compile-time
options in the makefile rather than by the parameterfile. This was
done in order to allow the generation of highly optimised binaries by
the compiler, even when the underlying source allows for many
different ways to run the code. Unfortunately, this technique has the
disadvantage that different simulations may require different binaries
of GADGET. If several simulations are run concurrently, there is
hence the danger that a simulation is started/resumed with the `wrong'
binary. Note that while GADGET checks the plausibility of some of
the most important code options, this is not done for all of them. To
minimise the risk of using the wrong code for a simulation, my
recommendation is therefore to produce a separate executable for each
simulation that is run. For example, a good strategy is to make a copy
of the whole code together with its makefile in the output directory
of each simulation run, and then to use this copy to compile the code
and to run the simulation.

The makefile contains a dummy list of all available compile-time code
options, with most of them commented out by default. To activate a
certain feature, the corresponding parameter should be commented in,
and given the desired value, where appropriate.  At the beginning of
the makefile, there is also a brief explanation of each of the
options.

Important Note:Whenever you change one of the makefile options
described below, a full recompilation of the code is necessary. To
guarantee that this is done, you should give the command make
clean, which will erase all object files, followed by make.

\section secmake1 Options that describe the physics of the simulation

- \b PERIODIC \n Set this if you want to have periodic boundary
   conditions.

- \b COOLING \n This enables radiative cooling and heating. It also
   enables an external UV background which is read from a file called
   TREECOOL.

- \b SFR \n This enables star formation using an effective multiphase
   models. This option requires cooling.

- \b METALS \n This activates the tracking of metal enrichment in gas
   and stars. Note that cooling by metal-lines is not yet included in
   the code.

- \b STELLARAGE \n This stores the formation times of each star
   particle, and includes this information in snapshot files.

- \b WINDS \n This activates galactic winds. Requires star formation.

- \b ISOTROPICWINDS \n This makes the wind isotropic. If not set, the
   wind is spawned in an axial way. Requires winds to be activated.

- \b NOGRAVITY \n This switches off gravity. Makes only sense for pure
   SPH simulations in non-expanding space.

- \b LONG_X/Y/Z \n These options can be used together with PERIODIC and
     NOGRAVITY only.  When set, the options define numerical factors that
     can be used to distorts the periodic simulation cube into a
     parallelepiped of arbitrary aspect ratio. This can be useful for
     idealized SPH tests.

- \b TWODIMS \n This effectively switches of one dimension in SPH,
     i.e. the code follows only 2d hydrodynamics in the xy-, yz-, or
     xz-plane. This only works with NOGRAVITY, and if all coordinates of
     the third axis are exactly equal. Can be useful for idealized SPH
     tests.

- \b SPH_BND_PARTICLES \n If this is set, particles with a particle-ID
     equal to zero do not receive any SPH acceleration. This can be useful
     for idealized SPH tests, where these particles represent fixed
     "walls".

\section secmake2 Options that affect SPH

- \b NOFIXEDMASSINKERNEL \n If set, the number of SPH particles in the
   kernel is kept constant instead of the mass contained in the
   kernel.

- \b NOGRADHSML \n If activated, an SPH equation of motion without
   \f$\nabla h\f$-correction factors is used.  Note, for the default
   `entropy'-formulation of SPH (Springel and Hernquist, 2002, MNRAS,
   333, 649), the switches NOFIXEDMASSINKERNEL and
   NOGRADHSML should not be set.

- \b NOVISCOSITYLIMITER \n If this is set, there is no explicit upper
   limit on the viscosity.  In the default version, this limiter will
   try to protect against possible particle `reflections', which may
   in principle occur for poor timestepping in the presence of strong
   shocks.

\section secmake3 Options for the numerical algorithms

- \b PMGRID=128 \n This enables the TreePM method, i.e. the long-range
   force is computed with a PM-algorithm, and the short range force
   with the tree. The parameter has to be set to the size of the mesh
   that should be used, e.g.~64, 96, 128, etc. The mesh dimensions
   need not necessarily be a power of two, but the FFT is fastest for
   such a choice.  Note: If the simulation is not in a periodic box,
   then a FFT method for vacuum boundaries is employed, using a mesh
   with dimension twice that specified by PMGRID.

- \b PLACEHIGHRESREGION=1+8 \n If this option is set (will only work
   together with \b PMGRID, then the long range force is computed in
   two stages: One Fourier-grid is used to cover the whole simulation
   volume, allowing the computation of the large-scale force.  A
   second Fourier mesh is placed on the region occupied by
   `high-resolution' particles, allowing the computation of an
   intermediate-scale force. Finally, the force on very small scales
   is computed by the tree. This procedure can be useful for
   `zoom-simulations', where the majority of particles (the high-res
   particles) are occupying only a small fraction of the volume. To
   activate this option, the parameter needs to be set to an integer
   that encodes the particle types that make up the high-res particles
   in the form of a bit mask. For example, if types 0, 1, and 4 are
   the high-res particles, then the parameter should be set to
   PLACEHIGHRESREGION=1+2+16, i.e. to the sum
   \f$2^0+2^1+2^4\f$. The spatial region covered by the high-res grid
   is determined automatically from the initial conditions. Note: If a
   periodic box is used, the high-res zone may not intersect the box
   boundaries.

- \b ENLARGEREGION=1.1 \n The spatial region covered by the high-res
   zone normally has a fixed size during the simulation, which
   initially is set to the smallest region that encompasses all
   high-res particles. Normally, the simulation will be interrupted,
   if high-res particles leave this region in the course of the
   run. However, by setting this parameter to a value larger than one,
   the high-res region can be expanded.  For example, setting it to
   1.4 will enlarge its side-length by 40% (it remains centred on the
   high-res particles). Hence, with such a setting, the high-res
   region may expand or move by a limited amount. If in addition \b
   SYNCHRONIZATION is activated, then the code will be able to
   continue even if high-res particles leave the initial high-res
   grid. In this case, the code will update the size and position of
   the grid that is placed onto the high-resolution region
   automatically. To prevent that this potentially happens every
   single PM step, one should nevertheless assign a value slightly
   larger than 1 to \b ENLARGEREGION.

- \b DOUBLEPRECISION \n This makes the code store and compute internal
   particle data in double precision. Note that output files are
   nevertheless written by converting the values that are saved to
   single precision.

- \b NOTREERND \n If this is not set, the tree construction will
   succeed even when there are a few particles at identical
   locations. This is done by `rerouting' particles once the node-size
   has fallen below \f$10^{-3}\f$ of the softening length. When this
   option is activated, this will be suppressed and the tree
   construction will always fail if there are particles at extremely
   close or identical coordinates.

- \b NOSTOP_WHEN_BELOW_MINTIMESTEP \n If this is
   activated, the code will not terminate when the timestep falls
   below the value of \b MinSizeTimestep specified in the
   parameterfile. This is useful for runs where one wants to enforce a
   constant timestep for all particles. This can be done by activating
   this option, and by setting \b MinSizeTimestep} and \b
   MaxSizeTimestep to an equal value.

- \b NOPMSTEPADJUSTMENT \n When this is set, the long-range timestep
   for the PM force computation is always determined by \b
   MaxSizeTimeStep.  Otherwise, it is set to the minimum of \b
   MaxSizeTimeStep and the timestep obtained for the maximum
   long-range force with an effective softening scale equal to the PM
   smoothing-scale.

\section secmake4 Architecture specific options

- \b T3E \n The code assumes that \b sizeof(int)=4 holds. A few
   machines (like Cray T3E) have \b sizeof(int)=8. In this case, set
   the T3E flag.

- \b NOTYPEPREFIX_FFTW \n If this is set, the fftw-header/libraries
   are accessed without type prefix (adopting whatever was chosen as
   default at compile-time of fftw). Otherwise, the type prefix 'd'
   for double-precision is used.

\section secmake5 Input options

- \b MOREPARAMS \n Activate this to allow a set of additional
   parameters in the parameterfile which control the star formation
   and feedback sector. This option must be activated when star
   formation is switched on.

\section secmake6 Output options

- \b OUTPUTPOTENTIAL \n This will force the code to compute
   gravitational potentials for all particles each time a snapshot
   file is generated. These values are then included in the snapshot
   files. Note that the computation of the values of the potential
   costs additional time.

- \b OUTPUTACCELERATION \n This will include the physical acceleration
   of each particle in snapshot files.

- \b OUTPUTCHANGEOFENTROPY \n This will include the rate of change of
   entropy of gas particles in snapshot files.

- \b OUTPUTTIMESTEP \n This will include the timesteps actually taken
   by each particle in the snapshot files.

\section secmake7 Miscellaneous options

- \b PEANOHILBERT \n This is a tuning option. When set, the code will
   bring the particles into Peano-Hilbert order after each domain
   decomposition. This improves cache utilisation and performance.

- \b WALLCLOCK \n If set, a wallclock timer is used by the code to
   measure internal time consumption (see cpu-log file).  Otherwise, a
   timer that measures consumed processor ticks is used.

- \b FORCETEST=0.01 \n This can be set to check the force accuracy of
   the code, and is only included as a debugging option. The option
   needs to be set to a number between 0 and 1 (e.g. 0.01), which
   specifies the fraction of randomly chosen particles for which at
   each timestep forces by direct summation are computed. The normal
   tree-forces and the `correct' direct summation forces are then
   collected in a file \b forcetest.txt for later inspection. Note
   that the simulation itself is unaffected by this option, but it
   will of course run much(!)  slower, particularly if
   FORCETEST*NumPart*NumPart \f$>>\f$ \b NumPart. Note: Particle IDs
   must be set to numbers \f$>=1\f$ for this option to work.

*/

max_stars_repo_name: egpbos/egp

max_issues_repo_name: egpbos/egp

max_forks_repo_name: egpbos/egp


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

double f (double x, void *params);
double f (double x, void *params)
{
    double n = *(double *) params;
    double f = exp (-x) * pow (x, n);

    return f;
}
void integralGen (int nmin, int nmax, double vals[]);
void integralGen (int nmin, int nmax, double vals[])
{
    double result, error;
    double a = 0., b = 1.;
    double abserr = 1.e-9, relerr = 1.e-9;
    double n;
    size_t np = 1000;
    gsl_integration_workspace *w = gsl_integration_workspace_alloc (np);

    gsl_function F;

    F.function = &f;

    for (int i = nmin; i <= nmax; i++)
    {
        n = (double) i;
        F.params = &n;

        gsl_integration_qag (&F, a, b, abserr, relerr, np, GSL_INTEG_GAUSS15,
            w, &result, &error);
        vals[i] = result;

    }

    gsl_integration_workspace_free (w);
}


text:
/* eigen/jacobi.c
 *
 * Copyright (C) 2004 Brian Gough, Gerard Jungman
 *
 * 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.
 */

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_eigen.h>

/* Algorithm 8.4.3 - Cyclic Jacobi.  Golub & Van Loan, Matrix Computations */

static inline double
symschur2 (gsl_matrix * A, size_t p, size_t q, double *c, double *s)
{
  double Apq = gsl_matrix_get (A, p, q);

  if (Apq != 0.0)
    {
      double App = gsl_matrix_get (A, p, p);
      double Aqq = gsl_matrix_get (A, q, q);
      double tau = (Aqq - App) / (2.0 * Apq);
      double t, c1;

      if (tau >= 0.0)
        {
          t = 1.0 / (tau + hypot (1.0, tau));
        }
      else
        {
          t = -1.0 / (-tau + hypot (1.0, tau));
        }

      c1 = 1.0 / hypot (1.0, t);

      *c = c1;
      *s = t * c1;
    }
  else
    {
      *c = 1.0;
      *s = 0.0;
    }

  /* reduction in off(A) is 2*(A_pq)^2 */

  return fabs (Apq);
}

inline static void
apply_jacobi_L (gsl_matrix * A, size_t p, size_t q, double c, double s)
{
  size_t j;
  const size_t N = A->size2;

  /* Apply rotation to matrix A,  A' = J^T A */

  for (j = 0; j < N; j++)
    {
      double Apj = gsl_matrix_get (A, p, j);
      double Aqj = gsl_matrix_get (A, q, j);
      gsl_matrix_set (A, p, j, Apj * c - Aqj * s);
      gsl_matrix_set (A, q, j, Apj * s + Aqj * c);
    }
}

inline static void
apply_jacobi_R (gsl_matrix * A, size_t p, size_t q, double c, double s)
{
  size_t i;
  const size_t M = A->size1;

  /* Apply rotation to matrix A,  A' = A J */

  for (i = 0; i < M; i++)
    {
      double Aip = gsl_matrix_get (A, i, p);
      double Aiq = gsl_matrix_get (A, i, q);
      gsl_matrix_set (A, i, p, Aip * c - Aiq * s);
      gsl_matrix_set (A, i, q, Aip * s + Aiq * c);
    }
}

inline static double
norm (gsl_matrix * A)
{
  size_t i, j, M = A->size1, N = A->size2;
  double sum = 0.0, scale = 0.0, ssq = 1.0;

  for (i = 0; i < M; i++)
    {
      for (j = 0; j < N; j++)
        {
          double Aij = gsl_matrix_get (A, i, j);

          if (Aij != 0.0)
            {
              double ax = fabs (Aij);

              if (scale < ax)
                {
                  ssq = 1.0 + ssq * (scale / ax) * (scale / ax);
                  scale = ax;
                }
              else
                {
                  ssq += (ax / scale) * (ax / scale);
                }
            }

        }
    }

  sum = scale * sqrt (ssq);

  return sum;
}

int
gsl_eigen_jacobi (gsl_matrix * a,
                  gsl_vector * eval,
                  gsl_matrix * evec, unsigned int max_rot, unsigned int *nrot)
{
  size_t i, p, q;
  const size_t M = a->size1, N = a->size2;
  double red, redsum = 0.0;

  if (M != N)
    {
      GSL_ERROR ("eigenproblem requires square matrix", GSL_ENOTSQR);
    }
  else if (M != evec->size1 || M != evec->size2)
    {
      GSL_ERROR ("eigenvector matrix must match input matrix", GSL_EBADLEN);
    }
  else if (M != eval->size)
    {
      GSL_ERROR ("eigenvalue vector must match input matrix", GSL_EBADLEN);
    }

  gsl_vector_set_zero (eval);
  gsl_matrix_set_identity (evec);

  for (i = 0; i < max_rot; i++)
    {
      double nrm = norm (a);

      if (nrm == 0.0)
        break;

      for (p = 0; p < N; p++)
        {
          for (q = p + 1; q < N; q++)
            {
              double c, s;

              red = symschur2 (a, p, q, &c, &s);
              redsum += red;

              /* Compute A <- J^T A J */
              apply_jacobi_L (a, p, q, c, s);
              apply_jacobi_R (a, p, q, c, s);

              /* Compute V <- V J */
              apply_jacobi_R (evec, p, q, c, s);
            }
        }
    }

  *nrot = i;

  for (p = 0; p < N; p++)
    {
      double ep = gsl_matrix_get (a, p, p);
      gsl_vector_set (eval, p, ep);
    }

  if (i == max_rot)
    {
      return GSL_EMAXITER;
    }

  return GSL_SUCCESS;
}

int
gsl_eigen_invert_jacobi (const gsl_matrix * a,
                         gsl_matrix * ainv, unsigned int max_rot)
{
  if (a->size1 != a->size2 || ainv->size1 != ainv->size2)
    {
      GSL_ERROR("jacobi method requires square matrix", GSL_ENOTSQR);
    }
  else if (a->size1 != ainv->size2)
    {
     GSL_ERROR ("inverse matrix must match input matrix", GSL_EBADLEN);
    }

  {
    const size_t n = a->size2;
    size_t i,j,k;
    unsigned int nrot = 0;
    int status;

    gsl_vector * eval = gsl_vector_alloc(n);
    gsl_matrix * evec = gsl_matrix_alloc(n, n);
    gsl_matrix * tmp = gsl_matrix_alloc(n, n);

    gsl_matrix_memcpy (tmp, a);

    status = gsl_eigen_jacobi(tmp, eval, evec, max_rot, &nrot);

    for(i=0; i<n; i++)
      {
        for(j=0; j<n; j++)
          {
            double ainv_ij = 0.0;

            for(k = 0; k<n; k++)
              {
                double f = 1.0 / gsl_vector_get(eval, k);
                double vik = gsl_matrix_get (evec, i, k);
                double vjk = gsl_matrix_get (evec, j, k);
                ainv_ij += vik * vjk * f;
              }
            gsl_matrix_set (ainv, i, j, ainv_ij);
          }
      }

    gsl_vector_free(eval);
    gsl_matrix_free(evec);
    gsl_matrix_free(tmp);

    if (status)
      {
        return status;
      }
    else
      {
        return GSL_SUCCESS;
      }
  }
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
/* specfunc/gsl_sf_dilog.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/* Author:  G. Jungman */

#ifndef __GSL_SF_DILOG_H__
#define __GSL_SF_DILOG_H__

#include <gsl/gsl_sf_result.h>

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

__BEGIN_DECLS

/* Real part of DiLogarithm(x), for real argument.
 * In Lewin's notation, this is Li_2(x).
 *
 *   Li_2(x) = - Re[ Integrate[ Log[1-s] / s, {s, 0, x}] ]
 *
 * Note that Im[Li_2(x)] = { 0 for x <= 1, -Pi*log(x) for x > 1 }
 */
int     gsl_sf_dilog_e(const double x, gsl_sf_result * result);
double     gsl_sf_dilog(const double x);

/* DiLogarithm(z), for complex argument z = r Exp[i theta].
 */
int gsl_sf_complex_dilog_e(const double r, double theta, gsl_sf_result * result_re, gsl_sf_result * result_im);

__END_DECLS

#endif /* __GSL_SF_DILOG_H__ */


text:
#include <assert.h>
#include <time.h>
#include <stdlib.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_rng.h>
#include "utils/utils.h"
#include "utils/compress.h"
#include "rule.h"

#define DIRICHLET 1
#define D_CONST .4

void populate_buf(uint64_t grule_size, uint8_t* rule_array, char* rule_buf)
{
  /* Populate buffer with newly created rule */
  for (uint64_t v = 0; v < grule_size; v++) {
    sprintf(&rule_buf[v], "%"PRIu8, rule_array[v]);
  }
  rule_buf[grule_size] = '\0';
}

/**
 * Symmetrize the rule by setting all the states and their symmetries to having
 * the same output.
 */
void symmetrize_rule(uint64_t grule_size,
                     uint8_t rule_array[grule_size],
                     int states, int horizon)
{
  int side = 2 * horizon + 1;
  int neigh_size = side * side - 1;

  uint32_t position_180;
  uint32_t position_90;
  uint32_t position_270;
  uint32_t position_vflip;
  uint32_t position_hflip;
  uint32_t position_dflip;
  uint32_t position_aflip;

  /* grule_size can be very big, this array is better on the heap */
  /* Keep track of already seen positions with book-keeping */
  uint8_t* book_keep = calloc(grule_size, sizeof(uint8_t));

  int pos;

  for (uint64_t i = 0; i < grule_size; ++i) {
    /* Skip already seen positions when looping through the rule */
    if (book_keep[i] == 1) {
      continue;
    }

    position_180 = 0;
    position_90 = 0;
    position_270 = 0;
    position_vflip = 0;
    position_hflip = 0;
    position_dflip = 0;
    position_aflip = 0;

    /* Create the representation of the symmetrized position by swapping the
       states in its number representation. */
    for (int p = 0; p < neigh_size + 1; ++p) {
      /* 180° rotation */
      position_180 += (uint32_t)ipow(states, p) *
        ((i / (uint32_t)ipow(states, neigh_size  - p)) % states);

      /* 90° rotation */
      pos = (neigh_size - side + 1 - (side * (p%side)) + p/side);
      position_90 += (uint32_t)ipow(states, p) *
        ((i / (uint32_t)ipow(states, pos)) % states);

      /* 270° rotation */
      position_270 += (uint32_t)ipow(states, p) *
        ((i / (uint32_t)ipow(states, neigh_size  - pos)) % states);

      /* Vertical flip */
      pos = (side * (p / side)) + (side - 1 - (p % 3));
      position_vflip += (uint32_t)ipow(states, p) *
        ((i / (uint32_t)ipow(states, pos)) % states);

      /* Horizontal flip */
      if (p/side < side/2) {
        pos = (p - side + neigh_size + 1)%(neigh_size + 1);
      } else if (p/side > side/2) {
        pos = (p + side)%(neigh_size + 1);
      } else {
        pos = p;
      }
      position_hflip += (uint32_t)ipow(states, p) *
        ((i / (uint32_t)ipow(states, pos)) % states);

      /* Diagonal flip */
      pos = neigh_size - ((p % side) * side  + (p / side));
      position_dflip += (uint32_t)ipow(states, p) *
        ((i / (uint32_t)ipow(states, pos)) % states);

      /* Antidiagonal flip */
      pos = ((p % side) * side  + (p / side));
      position_aflip += (uint32_t)ipow(states, p) *
        ((i / (uint32_t)ipow(states, pos)) % states);

    }

    /* Add all seen positions to the book to not process them again */
    book_keep[i] = 1;
    book_keep[position_180] = 1;
    book_keep[position_90] = 1;
    book_keep[position_270] = 1;
    book_keep[position_vflip] = 1;
    book_keep[position_hflip] = 1;
    book_keep[position_dflip] = 1;
    book_keep[position_aflip] = 1;

    rule_array[position_180] = rule_array[i];
    rule_array[position_90] = rule_array[i];
    rule_array[position_270] = rule_array[i];
    rule_array[position_vflip] = rule_array[i];
    rule_array[position_hflip] = rule_array[i];
    rule_array[position_dflip] = rule_array[i];
    rule_array[position_aflip] = rule_array[i];
  }
  free(book_keep);
}

/**
 * @brief Build a rule from the provided command-line arguments.
 * @param grule_size The size of the rule.
 * @param rule_array The array holding the rule.
 * @param rule_buf A mirror buffer containing a string representation of the
 *                 rule.
 * @param rule_arg The string from which the rule is read.
 * @param states Number of expected states in the rule.
 */
void build_rule_from_args(uint64_t grule_size,
                          uint8_t rule_array[grule_size],
                          char rule_buf[grule_size + 1],
                          char* rule_arg, int states)
{
  /* Rule is given in base-(states - 1) format */
  for (uint64_t s = 0 ; s < grule_size; ++s) {
    rule_array[s] = rule_arg[s] - '0';
    rule_buf[s] = rule_arg[s];
    assert(rule_array[s] < states);
  }
  rule_buf[grule_size] = '\0';
}

/**
 * Simple double comparison function used for sorting in the rule generation
 * function.
 */
int comp (const void * elem1, const void * elem2)
{
  float f = *((float*)elem1);
  float s = *((float*)elem2);
  return (f > s) - (f < s);
}

/**
 * Generate a random general rule. The rule is written to rule_array and a
 * string  representation is also saved in rule_buf.
 */
void generate_general_rule(uint64_t grule_size,
                           uint8_t rule_array[grule_size],
                           char rule_buf[grule_size + 1],
                           int states, int horizon)
{
  int inc;

  /* Choose lambda parameter at random as well as the proportion of
     transitions to other states */
  double alphas[states];
  double theta[states], lambda[states - 1], rand_num;

  /* This method samples the transition probability simplex according to a
     Dirichlet distribution with parameter D_CONST */
  if (DIRICHLET == 1) {
    const gsl_rng_type * T;

    gsl_rng * r;
    gsl_rng_env_setup();

    T = gsl_rng_default;
    r = gsl_rng_alloc(T);
    time_t t;
    gsl_rng_set(r, (unsigned long)time(&t));

    for (int i = 0; i < states; ++i) {
      alphas[i] = D_CONST;
    }
    gsl_ran_dirichlet(r, states, alphas, theta);
    lambda[0] = theta[0];
    for (int i = 1; i < states - 1; ++i) {
      lambda[i] = theta[i] + lambda[i - 1];
    }

    gsl_rng_free(r);
  }
  /* Second method that samples the simplex uniformly */
  else {
    for (int i = 0; i < states - 1; ++i) {
      lambda[i] = ((double)rand() / (double)((unsigned)RAND_MAX + 1));
    }
    qsort(lambda, sizeof(lambda)/sizeof(*lambda), sizeof(*lambda), comp);
  }

  for (uint64_t v = 0; v < grule_size; v++) {
    /* Assign the rule to the first state that passes the test */
    rand_num = (double)rand() / (double)((unsigned)RAND_MAX + 1);
    inc = 0;
    while (lambda[inc] < rand_num && inc < states - 1) {
      ++inc;
    }
    rule_array[v] = (uint8_t)inc;
  }

  symmetrize_rule(grule_size, rule_array, states, horizon);

  populate_buf(grule_size, rule_array, rule_buf);
}

void cross_breed(uint64_t grule_size, uint8_t* parent_rule_A,
                 uint8_t* parent_rule_B, uint8_t* child,
                 char rule_buf[grule_size], double rate,
                 int horizon, int states)
{
  for (uint64_t i = 0; i < grule_size; ++i) {
    double rand_num = (double)rand() / (double)((unsigned)RAND_MAX + 1);
    child[i] = (rate > rand_num) ? parent_rule_A[i]: parent_rule_B[i];
  }

  perturb_rule(grule_size, child, rule_buf, states, horizon, 0.05);
}

void perturb_rule(uint64_t grule_size,
                  uint8_t rule_array[grule_size],
                  char rule_buf[grule_size + 1],
                  int states, int horizon, double rate)
{
  for (uint64_t v = 0; v < grule_size; ++v) {
    /* Perturb transisition outcome with probability rate */
    double rand_num = (double)rand() / (double)((unsigned)RAND_MAX + 1);
    if (rand_num < rate) {
      rule_array[v] += (1 + rand() % (states - 1));
      rule_array[v] %= states;
    }
  }

  symmetrize_rule(grule_size, rule_array, states, horizon);
  populate_buf(grule_size, rule_array, rule_buf);
}

void make_map(struct Options2D* opts, char* rule_buf, int step)
{
  FILE* dic_file;
  char* fname;

  asprintf(&fname, "%s/map/%lu.map", opts->data_dir_name, hash(rule_buf));
  dic_file = fopen(fname, "w+");

  if (dic_file == NULL) {
    fprintf(stderr, "Error while opening map file %s.", fname);
    exit(1);
  }

  fprintf(dic_file, "%s", rule_buf);
  fclose(dic_file);
  free(fname);

  sprintf(rule_buf, "%lu", hash(rule_buf));
  printf("%i: Rule %s\n", step, rule_buf);
}

void generate_totalistic_rule(uint64_t rule_size, uint8_t rule_array[rule_size],
                              char rule_buf[rule_size + 1], int states)
{
  unsigned long rule_number = 0UL;

  for (uint64_t s = 0 ; s < rule_size; ++s) {
    rule_array[s] = rand() % states;
    rule_number += rule_array[s] * ipow(states, s);
    if (states >= 3) {
      rule_buf[s] = '0' + rule_array[s];
    }
  }

  if (states == 2) {
    sprintf(rule_buf, "%lu", rule_number);
  }
}


text:
#ifndef _UI_HELPERS_H_
#define _UI_HELPERS_H_

#define OEMRESOURCE

#include <algorithm>
#include <functional>
#include <memory>
#include <optional>
#include <span>
#include <string>
#include <type_traits>
#include <unordered_map>
#include <vector>

#include <ppl.h>

#include <gsl/gsl>
#include <range/v3/all.hpp>

// Included before windows.h, because pfc.h includes winsock2.h
#include "../pfc/pfc.h"

#include <windows.h>
#include <WindowsX.h>
#include <SHLWAPI.H>
#include <vssym32.h>
#include <uxtheme.h>
#include <shldisp.h>
#include <ShlObj.h>
#include <gdiplus.h>
#include <Usp10.h>
#include <CommonControls.h>
#include <intsafe.h>

#include <wil/resource.h>

#include "../mmh/stdafx.h"

#ifndef RECT_CX
#define RECT_CX(rc) ((rc).right - (rc).left)
#endif

#ifndef RECT_CY
#define RECT_CY(rc) ((rc).bottom - (rc).top)
#endif

#include "literals.h"

#include "handle.h"
#include "win32_helpers.h"
#include "ole.h"

#include "dialog.h"
#include "dpi.h"
#include "container_window.h"
#include "message_hook.h"
#include "trackbar.h"
#include "solid_fill.h"
#include "theming.h"

#include "gdi.h"
#include "text_drawing.h"
#include "list_view/list_view.h"
#include "drag_image.h"
#include "info_box.h"
#include "menu.h"

#include "ole/data_object.h"
#include "ole/enum_format_etc.h"

#endif