text:
#pragma once

#include <Strategic/Solver.h>
#include <Sudoku/Board.h>
#include <Sudoku/Location_Utilities.h>
#include <Sudoku/Options.h>
#include <gsl/gsl>

#include <string>

#include <iomanip> // setw(), setfill()
#include <sstream>
#include <utility>

namespace Sudoku
{
// Solver-function declarations:
template<int N>
void test_solver_unique(Board<Options<elem_size<N>>, N>&);
template<int N>
void test_solver_exclusive(Board<Options<elem_size<N>>, N>&);
template<int N>
Board<Value, N> getResult(Board<Options<elem_size<N>>, N> const&) noexcept;

//====--------------------------------------------------------------------====//
class Console
{
	struct delimiter
	{
		char space;
		char empty;
		std::string row_block;
		std::string col_block;
		std::string newl;
		std::string block_cross;
	};

	const delimiter space{' ', ' ', "", "", "", ""};
	const delimiter space2{' ', ' ', "", "", "\n", ""};
	const delimiter display{' ', ' ', "-", "|", "\n", "o"};
	// const delimiter csv;
	// const delimiter xml;

public:
	Console() noexcept;
	explicit Console(delimiter);
	//~Console() = default;

	template<int N>
	std::stringstream print_row(const Board<int, N>&, gsl::index row_id) const;
	template<int N, int E>
	std::stringstream
		print_row(const Board<Options<E>, N>&, gsl::index row_id) const;
	template<int N>
	std::stringstream print_board(const Board<int, N>&) const;
	template<int N, int E>
	std::stringstream print_board(const Board<Options<E>, N>&) const;

private:
	const delimiter d;
	// bool Format::find_option(const Board<std::set<int>>&,Location,int value);
	// format elem
	// format col-block section
	// format row-block = separator line
};

//====--------------------------------------------------------------------====//
// Member-functions:

inline Console::Console() noexcept : d(display)
{
	// empty constructor
}

inline Console::Console(delimiter del) : d(std::move(del))
{
	// empty constructor
}

namespace impl
{
	[[nodiscard]] constexpr int pow(int const base, int exp)
	{
		if (exp == 0)
		{
			return 1;
		}
		return base * pow(base, --exp);
	}

	[[nodiscard]] constexpr int charsize(int value, int length)
	{
		constexpr int decimal{10};
		if (value < pow(decimal, length))
		{
			return length;
		}
		++length;
		return charsize(value, length);
	}

	[[nodiscard]] constexpr int charsize(int value)
	{
		constexpr int decimal{10};
		if (value < decimal)
		{
			return 1;
		}
		return charsize(value, 2);
	}
} // namespace impl

template<int N>
std::stringstream
	Console::print_row(const Board<int, N>& input, gsl::index row_id) const
{
	std::stringstream stream;
	constexpr int chars = impl::charsize(elem_size<N>) + 1;

	stream << d.col_block << std::setfill(d.space);
	for (gsl::index i = 0; i < elem_size<N>; ++i)
	{
		if (input[row_id][i] == 0) // no value
		{
			stream << std::setw(chars) << d.space;
		}
		else
		{
			stream << std::setw(chars) << input[row_id][i];
		}
		if ((i + 1) % base_size<N> == 0)
		{
			stream << std::setw(2) << d.col_block;
		}
	}
	return stream;
}

template<int N>
std::stringstream Console::print_board(const Board<int, N>& input) const
{
	std::stringstream stream;
	std::stringstream temp;
	constexpr int chars = impl::charsize(elem_size<N>) + 1;

	// opening bar
	temp << d.block_cross;
	for (gsl::index j = 0; j < base_size<N>; ++j)
	{
		temp << std::setfill(d.row_block[0])
			 << std::setw(chars * base_size<N> + 2) << d.block_cross;
	}
	std::string bar;
	temp >> bar;
	stream << bar << '\n';

	// loop rows
	for (gsl::index i = 0; i < elem_size<N>; ++i)
	{
		stream << print_row(input, i).str() << d.newl;
		if ((i + 1) % base_size<N> == 0)
		{
			stream << bar << '\n';
		}
	}
	return stream;
}

template<int N, int E>
std::stringstream Console::print_board(const Board<Options<E>, N>& input) const
{
	static_assert(E == N * N);
	constexpr gsl::index block_size = elem_size<N> + base_size<N> + 2;
	constexpr gsl::index row_length = base_size<N> * block_size;
	/*
	9   9   9
	o-----------------------------------------o
	| 123 123 123 | 123 123 123 | 123 123 123 |
	| 456 456 456 | 456 456 456 | 456 456 456 | 9
	| 789 789
	|
	| 123 123
	| ...
	o-----------------------------------------o
	|
	*/
	std::stringstream stream;
	std::stringstream n0; // horizontal block separator
	std::stringstream n4; // vertical element separator
	n0 << std::setfill(d.row_block[0]) << d.block_cross << std::setw(row_length)
	   << d.block_cross << d.newl;
	n4 << std::setfill(' ') << std::setw(block_size) << d.col_block;

	stream << std::setfill(d.empty) << std::setw(base_size<N>);

	stream << '\n' << n0.str();

	for (gsl::index row{0}; row < elem_size<N>; ++row)
	{
		stream << print_row(input, row).str();

		if ((row + 1) % base_size<N> == 0)
		{
			stream << n0.str();
		}
		else
		{
			stream << '|';
			for (gsl::index index{}; index < base_size<N>; ++index)
			{
				stream << n4.str();
			}
			stream << '\n';
		}
	}
	return stream;
}

template<int N, int E>
std::stringstream Console::print_row(
	const Board<Options<E>, N>& input, gsl::index row_id) const
{
	std::stringstream stream;

	gsl::index X{1};
	for (gsl::index k{0}; k < base_size<N>; ++k)
	{
		stream << d.col_block << d.space;
		for (gsl::index col{0}; col < elem_size<N>; ++col)
		{
			for (gsl::index i{X}; i < X + base_size<N>; ++i)
			{
				if (input[row_id][col].test(to_Value<N>(i)))
				{
					stream << i;
				}
				else
				{
					stream << d.empty;
				}
			}
			if ((col + 1) % base_size<N> == 0)
			{
				stream << std::setfill(d.space) << std::setw(2) << d.col_block
					   << d.space;
			}
			else
			{
				stream << d.space;
			}
		}
		stream << d.newl;
		X += base_size<N>;
	}
	return stream;
}

//====--------------------------------------------------------------------====//
// Solver function applications:
template<int N> // NOLINTNEXTLINE(runtime/references)
void test_solver_unique(Board<Options<elem_size<N>>, N>& board)
{
	int found{1};
	while (found > 0)
	{
		found = 0;
		for (gsl::index i = 0; i < elem_size<N>; ++i)
		{
			found += unique_in_section(board, board.row(i));
		}
		for (gsl::index i = 0; i < elem_size<N>; ++i)
		{
			found += unique_in_section(board, board.col(i));
		}
		for (gsl::index i = 0; i < elem_size<N>; ++i)
		{
			found += unique_in_section(board, board.block(i));
		}
	}
}

template<int N> // NOLINTNEXTLINE(runtime/references)
void test_solver_exclusive(Board<Options<elem_size<N>>, N>& board)
{
	int found{1};
	while (found > 0)
	{
		found = 0;
		for (gsl::index i = 0; i < elem_size<N>; ++i)
		{
			found += section_exclusive(board, board.row(i));
		}
		for (gsl::index i = 0; i < elem_size<N>; ++i)
		{
			found += section_exclusive(board, board.col(i));
		}
		for (gsl::index i = 0; i < elem_size<N>; ++i)
		{
			found += section_exclusive(board, board.block(i));
		}
	}
}

template<int N>
Board<Value, N>
	getResult(const Board<Options<elem_size<N>>, N>& options) noexcept
{
	Board<Value, N> result{};
	for (gsl::index i = 0; i < full_size<N>; ++i)
	{
		if (options[Location<N>(i)].is_answer())
		{
			result[Location<N>(i)] = get_answer(options[Location<N>(i)]);
		}
	}
	return result;
}

} // namespace Sudoku

max_stars_repo_name: Farwaykorse/fwkSudoku

max_issues_repo_name: Farwaykorse/fwkSudoku

max_forks_repo_name: Farwaykorse/fwkSudoku


text:
#ifndef KABSCH_RMSD
#define KABSCH_RMSD

#include <valarray>

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

namespace kabsch
{

template <class M>
double rmsd(
        const M P,
        const M Q)
{
    double rmsd {0.0};
    const unsigned int D {3};
    const unsigned int size = P.size();
    const unsigned int N = size / D;

    for(unsigned int i = 0; i < size; ++i) {
        rmsd += (P[i] - Q[i])*(P[i] - Q[i]);
    }

    return sqrt(rmsd/N);
}

template <class M>
M centroid(M coordinates)
{

    double x {0};
    double y {0};
    double z {0};
    unsigned int size = coordinates.size();
    unsigned int n_atoms = size / 3;

    unsigned int i = 0;
    while(i<size)
    {
        x += coordinates[i++];
        y += coordinates[i++];
        z += coordinates[i++];
    }

    x /= n_atoms;
    y /= n_atoms;
    z /= n_atoms;

    return M {x, y, z};
}

template <class Matrix>
Matrix multiply(Matrix A, Matrix B,
    const int M,
    const int N,
    const int K)
{
    double one = 1.0;

    Matrix C(M*N);

    cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
        M, N, K, one,
        begin(A), K,
        begin(B), N, 1.0,
        begin(C), N);

    return C;
}

template <class Matrix>
Matrix transpose_multiply(Matrix A, Matrix B,
    const int M,
    const int N,
    const int K)
{
    double one = 1.0;

    Matrix C(M*N);

    cblas_dgemm(CblasRowMajor, CblasTrans, CblasNoTrans,
        M, N, K, one,
        begin(A), M,
        begin(B), N, one,
        begin(C), N);

    return C;
}

template <class Matrix>
std::tuple<Matrix, Matrix, Matrix> matrix_svd(Matrix A, int rows, int cols)
{
    // lapack_int LAPACKE_dgesvd( int matrix_layout, char jobu, char jobvt,
    //     lapack_int m, lapack_int n,
    //     double* a, lapack_int lda,
    //     double* s, double* u, lapack_int ldu,
    //     double* vt, lapack_int ldvt,
    //     double* superb );

    Matrix U(cols*rows);
    Matrix S(rows);
    Matrix VT(cols*rows);
    Matrix superb(cols*rows);
    int info;

    info = LAPACKE_dgesvd(LAPACK_ROW_MAJOR, 'A', 'A',
        rows, cols,
        begin(A), rows,
        begin(S),
        begin(U), rows,
        begin(VT), rows,
        begin(superb));

    // TODO check if SVD is unconverged
    if(info > 0)
    {
        // Failed
    }

    return make_tuple(U, S, VT);
}

template <class M>
double determinant3x3(M A)
{
    // determinant of a square 3x3 matrix
    double det = A[0]*A[4]*A[8]
        +A[1]*A[5]*A[6]
        +A[2]*A[3]*A[7]
        -A[2]*A[4]*A[6]
        -A[1]*A[3]*A[8]
        -A[0]*A[5]*A[7];

    return det;
}

template <class M, class T>
M kabsch(M P, M Q, const T n_atoms)
{
    // const unsigned int L = P.size();
    // const unsigned int D = 3;

    M U;
    M V;
    M S;

    M C = transpose_multiply(P, Q, 3, 3, n_atoms);

    tie(U, S, V) = matrix_svd(C, 3, 3);

    // Getting the sign of the det(U)*(V) to decide whether we need to correct
    // our rotation matrix to ensure a right-handed coordinate system.
    if(determinant3x3(U)*determinant3x3(V) < 0.0)
    {
        // TODO More numpy'ish way to do this?
        U[std::slice( 2, 3, 3 )] = {-U[3*0+2], -U[3*1+2], -U[3*2+2]};
    }

    M rotation = multiply(U, V, 3, 3, 3);

    return rotation;
}

template <class M, class T>
M kabsch_rotate(
        const M P,
        const M Q,
        const T n_atoms)
{
    M U = kabsch(P, Q, n_atoms);
    M product = multiply(P, U, n_atoms, 3, 3);
    return product;
}

template <class M, class T>
double kabsch_rmsd(
        const M P,
        const M Q,
        const T n_atoms)
{
    M P_rotated = kabsch_rotate(P, Q, n_atoms);
    return rmsd(P_rotated, Q);
}

} // namespace rmsd

#endif

ext: h

max_stars_repo_name: kjappelbaum/rmsd

max_issues_repo_name: kjappelbaum/rmsd

max_forks_repo_name: kjappelbaum/rmsd


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

#include <config.h>
#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_histogram.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_ieee_utils.h>

void test1d (void);
void test2d (void);
void test1d_resample (void);
void test2d_resample (void);
void test1d_trap (void);
void test2d_trap (void);

int
main (void)
{
  test1d();
  test2d();
  test1d_resample();
  test2d_resample();
  test1d_trap();
  test2d_trap();

  exit (gsl_test_summary ());
}

ext: c

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
#pragma once

#include <gsl/span>

#include "aas/aas_math.h"

namespace aas {

	class Skeleton;
	struct SkelAnim;
	struct SkelBoneState;

	struct AnimPlayerStreamHeader {
		int field0;
		int field1;
		int field2;
		int8_t pad[3];
	};

	struct AnimPlayerStreamBone
	{
		DX::XMFLOAT4 prevScale;
		DX::XMFLOAT4 scale;
		DX::XMFLOAT4 prevRotation;
		DX::XMFLOAT4 rotation;
		DX::XMFLOAT4 prevTranslation;
		DX::XMFLOAT4 translation;
		int16_t scaleFrame;
		int16_t scaleNextFrame;
		int16_t rotationFrame;
		int16_t rotationNextFrame;
		int16_t translationFrame;
		int16_t translationNextFrame;
		int field6C;
	};

	struct AnimPlayerStream {
		const Skeleton* skeleton;
		const SkelAnim* animation;
		int streamIdx;
		float scaleFactor;
		float translationFactor;
		float currentFrame;
		void* keyframePtr;
		int boneCount;
		AnimPlayerStreamBone bones[1]; //set according to skafile bone count

		AnimPlayerStream(AnimPlayerStream&) = delete;
		AnimPlayerStream(AnimPlayerStream&&) = delete;
		AnimPlayerStream& operator=(AnimPlayerStream&) = delete;
		AnimPlayerStream& operator=(AnimPlayerStream&&) = delete;

		void SetFrame(float frame);
		void Initialize(const SkelAnim *animation, int streamIdx);
		void GetBoneState(gsl::span<SkelBoneState> boneStateOut);
		float GetCurrentFrame() const;
	};

}

max_stars_repo_name: edoipi/TemplePlus

max_issues_repo_name: edoipi/TemplePlus

max_forks_repo_name: edoipi/TemplePlus


text:
#ifndef POINTCLOUD_H_8O5EBVHZ
#define POINTCLOUD_H_8O5EBVHZ

#include <Eigen/Core>
#include <algorithm>
#include <gsl/gsl>
#include <sens_loc/math/coordinate.h>
#include <sens_loc/math/image.h>

namespace sens_loc::math {

/// A pointcloud is a set of camera coordinates.
template <typename Real>
using pointcloud   = std::vector<camera_coord<Real>>;
using pointcloud_t = pointcloud<float>;

/// A set of pixel coordinates.
template <typename Real>
using imagepoints   = std::vector<pixel_coord<Real>>;
using imagepoints_t = imagepoints<float>;

/// Calculate the point-wise euclidean distance between each point in \c c0
/// and \c c1.
/// The result is a row-vector with the distance for point 'i' at position 'i'.
/// \tparam one of \c pointcloud or \c imagepoints.
/// \pre c0.size() == c1.size() => the same number of points
/// \post each element in the result is non-negative.
template <typename Points>
std::vector<typename Points::value_type::real>
pointwise_distance(const Points& c0, const Points& c1) noexcept {
    Expects(c0.size() == c1.size());

    std::vector<typename Points::value_type::real> distances;
    if (c0.empty())
        return distances;

    distances.reserve(c0.size());
    for (std::size_t i = 0; i < c0.size(); ++i)
        distances.emplace_back((c0[i] - c1[i]).norm());

    Ensures(distances.size() == c0.size());
    Ensures(distances.size() == c1.size());
    Ensures(*std::min_element(distances.begin(), distances.end()) >= 0.0F);

    return distances;
}

/// A pose is an affine transformation in 3 dimensions. Internally this is
/// processed with homogeneous coordinates.
/// \sa pointcloud_t
using pose_t = Eigen::Matrix4f;

/// Calculate the relative pose to get from \c from to \c to.
///
/// Let \c O be the origin of a coordinate system, \c from and \c to are both
/// absolute poses in this coordinate system.
/// \c relative_pose calculates the relative pose to get from \c from to \c to.
/// \pre both poses need to be valid
/// \pre \c from must be invertable => pose is valid
inline pose_t relative_pose(const pose_t& from, const pose_t& to) noexcept {
    return from.inverse() * to;
}

/// Translate and rotate all points in \c points with the transformation \c p.
inline pointcloud_t operator*(const pose_t&       p,
                              const pointcloud_t& points) noexcept {
    pointcloud_t result;
    result.reserve(points.size());

    for (const auto& pt : points) {
        Eigen::Vector4f tr = p * Eigen::Vector4f{pt.X(), pt.Y(), pt.Z(), 1.0F};
        result.emplace_back(tr.x(), tr.y(), tr.z());
    }
    return result;
}

}  // namespace sens_loc::math

#endif /* end of include guard: POINTCLOUD_H_8O5EBVHZ */

max_stars_repo_name: JonasToth/depth-conversions

max_issues_repo_name: JonasToth/depth-conversions

max_forks_repo_name: JonasToth/depth-conversions


text:
#ifndef config_e9763ac5_66db_439d_99b9_4c84c1b41e14_h
#define config_e9763ac5_66db_439d_99b9_4c84c1b41e14_h

#include <gslib\config.h>

#ifdef __cplusplus

#define __rathen_begin__    namespace gs { namespace rathen {
#define __rathen_end__      }; };

#endif

/* rathen source postfix & binary postfix */
#define rathen_spf  _t(".rs")
#define rathen_bpf  _t(".rb")

#endif

max_stars_repo_name: lymastee/gslib

max_issues_repo_name: lymastee/gslib

max_forks_repo_name: lymastee/gslib


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

#include <float.h>
#include <stdio.h>
#include <gsl/gsl_ieee_utils.h>

int
main ()
{
  int i ;
  double x = 10 ;
  float f = 1.0/3.0 ;
  double d = 1.0/3.0 ;

  double fd = f ; /* promote from float to double */

  gsl_ieee_env_setup() ;

  printf("     float 1/3 = ") ; gsl_ieee_printf_float(&f) ; printf("\n") ;
  printf("promoted float = ") ; gsl_ieee_printf_double(&fd) ; printf("\n") ;
  printf("    double 1/3 = ") ; gsl_ieee_printf_double(&d) ; printf("\n") ;

  for (i=0;i<10;i++) {
    x = 0.5 *(x + 2.0/x) ;
    printf("%.18g ",x) ;
    gsl_ieee_printf_double(&x) ;
    printf("\n") ;
  }

  x = 10 * x * GSL_DBL_MAX ;

  printf("%.18g ",x) ; gsl_ieee_printf_double(&x) ; printf("\n") ;

  x = x / 0 ; ;

  printf("%.18g ",x) ; gsl_ieee_printf_double(&x) ; printf("\n") ;

  f = -1.0/3.0 ;

  while (f < 0) {
    f = f / 10 ;
    printf("%.18g ",f) ; gsl_ieee_printf_float(&f) ; printf("\n") ;
  }

}

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim


text:
/* DomSelfAdaptFWD.c

Forward-in-time simulation of an adaptive alleles with dominance and selfing, with linked neutral fragment.

This is the BATCH version - to be used on cluster machine to produce many simulation replicates at once.

Burn-in program: generates neutral background diversity for other simulations to use. See README for more information.

Simulation uses routines found with the GNU Scientific Library (GSL)
(http://www.gnu.org/software/gsl/)
Since GSL is distributed under the GNU General Public License
(http://www.gnu.org/copyleft/gpl.html), you must download it
separately from this file.

*/

/* Preprocessor statements */
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <stddef.h>
#include <string.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <sys/stat.h>
#include <sys/types.h>

#define INITTS 1000

/* Function prototypes */
void Wait();
unsigned int bitswitch(unsigned int x);
void generation_BI(unsigned int N, double self, double R, unsigned int *nums, unsigned int **neutin, unsigned int **neutout, unsigned int npoly, double *polypos, const gsl_rng *r);
void rec_sort(double *index, unsigned int nrow);
void addpoly(unsigned int N, unsigned int **neutin, double *polyloc, unsigned int *npoly, double theta, const gsl_rng *r);
void reassign2(unsigned int **neutin, unsigned int **neutout, double *posin, double *posout, unsigned int npoly, unsigned int N, unsigned int aft);
void popprint(unsigned int **neutin, double *posin, unsigned int npoly, unsigned int N, unsigned int ei, unsigned int suffix);
void polyprint(unsigned int **neutin, double *posin, unsigned int npoly, unsigned int N);
void ptrim(unsigned int size, unsigned int **neutin, unsigned int **neutout, double *posin, double *posout, unsigned int npoly, unsigned int *npolyT, double *avpi);
unsigned int uniqueH(unsigned int **neutin, unsigned int npoly, unsigned int N);

void Wait(){
	printf("Press Enter to Continue");
	while( getchar() != '\n' );
	printf("\n");
}

/* Bit-switching routine */
unsigned int bitswitch(unsigned int x){
	unsigned int ret;
	switch(x)
	{
		case 0:
			ret = 1;
			break;
		case 1:
			ret = 0;
			break;
		default:	/* If none of these cases chosen, exit with error message */
			fprintf(stderr,"Error in bitswitch: input is not 0 or 1.\n");
			exit(1);
            break;
	}
	return ret;
}

/* Sorting rec events after choosing them */
void rec_sort(double *index, unsigned int nrow){
	unsigned int j, i;		/* Sorting indices */
	double temp0;		/* For swapping */

	for(j = 1; j < nrow; j++){
		i = j;
		while( (i > 0) && ( *(index + (i - 1) ) > *(index + i) )){
			/* Swapping entries */
			temp0 = *(index + (i - 1));
			*(index + (i - 1)) =  *(index + (i));
			*(index + (i)) = temp0;
			i--;
		}
	}

}

/* Reproduction routine */
void generation_BI(unsigned int N, double self, double R, unsigned int *nums, unsigned int **neutin, unsigned int **neutout, unsigned int npoly, double *polypos, const gsl_rng *r){
	unsigned int i, j, k;		/* Pop counter, neutral marker counter, rec event counter */
	unsigned int isself = 0;	/* Is this a selfing reproduction? */
	unsigned int choose1 = 0;	/* Chromosome to be selected */
	unsigned int choose2 = 0;	/* Chromosome to be selected (2nd with outcrossing) */
	unsigned int wc1 = 0;		/* Which chrom (parent 1) */
	unsigned int wc2 = 0;		/* Which chrom (parent 2) */
	unsigned int index1 = 0;	/* Index of chosen sample 1 */
	unsigned int index2 = 0;	/* Index of chosen sample 2 */
	unsigned int nself = 0;		/* Number of selfing events */
	unsigned int selfix = 0;	/* Selfing index */
	unsigned int recix1 = 0;	/* Rec index chr 1 */
	unsigned int recix2 = 0;	/* Rec index chr 2 */
	unsigned int nrec1 = 0;		/* Number rec events chr 1 */
	unsigned int nrec2 = 0;		/* Number rec events chr 2 */

	/* First deciding number of selfing events; creating vector of events + choosing */
	nself = gsl_ran_binomial(r, self, N);
	unsigned int *selfev = calloc(nself,sizeof(unsigned int));
	gsl_ran_choose(r, selfev, nself, nums, N, sizeof(unsigned int));

	for(i = 0; i < N; i++){		/* Regenerating population */

    	/* Choosing first parent, and relevant chromosomes */
    	/* ADJUSTED SO PURELY NEUTRAL SELECTION ONLY */
		choose1 = gsl_rng_uniform_int(r,N);
		wc1 = gsl_ran_bernoulli(r,0.5);
		wc2 = gsl_ran_bernoulli(r,0.5);
		index1 = 2*choose1 + wc1;

		/* Copying over selected sites, first checking if selfing occurs */
		isself = 0;
		if(nself != 0){
			if(*(selfev + selfix) == i){
				isself = 1;
				selfix++;
			}
		}
		if(isself == 1){
			index2 = 2*choose1 + wc2;
		}else if(isself == 0){
			choose2 = gsl_rng_uniform_int(r,N);
			index2 = 2*choose2 + wc2;
		}

		/* Now copying over neutral fragment, accounting for recombination */
		if(R != 0){
			/* Choosing number of recombination events on arms 1, 2 */
			recix1 = 0;
			recix2 = 0;
			nrec1 = gsl_ran_poisson(r,R);
			nrec2 = gsl_ran_poisson(r,R);
			double *recev1 = calloc(nrec1 + 1,sizeof(double));
			double *recev2 = calloc(nrec2 + 1,sizeof(double));
			for(k = 0; k < nrec1; k++){
				*(recev1 + k) = gsl_ran_flat(r,0,1);
			}
			for(k = 0; k < nrec2; k++){
				*(recev2 + k) = gsl_ran_flat(r,0,1);
			}
			*(recev1 + nrec1) = 1.01;
			*(recev2 + nrec2) = 1.01;
			rec_sort(recev1,nrec1 + 1);
			rec_sort(recev2,nrec2 + 1);

			for(j = 0; j < npoly; j++){
				if( (*(polypos + j)) >= *(recev1 + recix1)){
					wc1 = bitswitch(wc1);
					index1 = 2*choose1 + wc1;
					recix1++;
				}
				if( (*(polypos + j)) >= *(recev2 + recix2)){
					wc2 = bitswitch(wc2);
					if(isself == 1){
						index2 = 2*choose1 + wc2;
					}else if(isself == 0){
						index2 = 2*choose2 + wc2;
					}
					recix2++;
				}

				*((*(neutout + 2*i)) + j) = *((*(neutin + index1)) + j);
				*((*(neutout + 2*i + 1)) + j) = *((*(neutin + index2)) + j);

			}

			free(recev1);
			free(recev2);

		}else if (R == 0){
			for(j = 0; j < npoly; j++){
				*((*(neutout + 2*i)) + j) = *((*(neutin + index1)) + j);
				*((*(neutout + 2*i + 1)) + j) = *((*(neutin + index2)) + j);
			}
		}
	}

	free(selfev);

}	/* End of reproduction routine */

/* Adding neutral polymorphism routine */
void addpoly(unsigned int N, unsigned int **neutin, double *polyloc, unsigned int *npoly, double theta, const gsl_rng *r){
	unsigned int i, j;
	int k;
	unsigned int newpoly = 0;			/* Number of new polymorphisms added */
	unsigned int pchr = 0;				/* Chromosome of appearance */
	unsigned int pfound = 0;			/* Flag to denote if right index found */
	double ploc = 0;					/* Location of polymorphism */
/*	FILE *ofp_poly;				 Pointer for polymorphisms output */

	newpoly = gsl_ran_poisson(r,(0.5*theta));
/*	printf("Add %d new polymorphisms\n",newpoly);*/

	for(j = 0; j < newpoly; j++){
		ploc = gsl_ran_flat(r,0.0,1.0);
/*		printf("ploc is %lf\n",ploc);	*/
		pchr = gsl_rng_uniform_int(r,2*N);
/*		printf("pchr is %d\n",pchr);*/

		/* Inserting new polymorphism */
		if(*(npoly) == 0){
/*			printf("zero entry here\n");	*/
			*(polyloc + 0) = ploc;
				for(i = 0; i < 2*N; i++){
					*((*(neutin + i)) + 0) = 0;
				}
			*((*(neutin + pchr)) + 0) = 1;
		}else if(*(npoly) > 0){
/*			printf("non-zero entry here\n");		*/
			pfound = 0;
			for(k = (*(npoly) - 1); k >= 0; k--){
/*				printf("k is %d\n",k);*/
				if(*(polyloc + k) > ploc){
/*					printf("is here 1\n");*/
					*(polyloc + k + 1) = *(polyloc + k);
					for(i = 0; i < 2*N; i++){
						*((*(neutin + i)) + k + 1) = *((*(neutin + i)) + k);
					}
				}else if(*(polyloc + k) <= ploc){
/*					printf("is here 2\n");				*/
					pfound = 1;
					*(polyloc + k + 1) = ploc;
					for(i = 0; i < 2*N; i++){
						*((*(neutin + i)) + k + 1) = 0;
					}
					*((*(neutin + pchr)) + k + 1) = 1;
					break;
				}
			}

			if(pfound == 0){		/* In this case polymorphism is inserted at the start */
				*(polyloc + 0) = ploc;
				for(i = 0; i < 2*N; i++){
					*((*(neutin + i)) + 0) = 0;
				}
				*((*(neutin + pchr)) + 0) = 1;
			}

		}

		*(npoly) += 1;
/*		printf("n poly is %d\n",*(npoly));			*/
		if(*(npoly) == INITTS){
			fprintf(stderr,"Too many neutral polymorphisms.\n");
			exit(1);
		}

	}
}

/* Reassigning new population routine (after trimming) */
void reassign2(unsigned int **neutin, unsigned int **neutout, double *posin, double *posout, unsigned int npoly, unsigned int N, unsigned int aft){
	unsigned int i, j;		/* pop counter, poly counter */

	for(i = 0; i < 2*N; i++){
		for(j = 0; j < npoly; j++){
			*((*(neutout + i)) + j) = *((*(neutin + i)) + j);
			if(i == 0 && aft == 1){
				*(posout + j) = *(posin + j);
			}
		}
	}

}	/* End of reassignment routine */

/* Printing population state to file */
void popprint(unsigned int **neutin, double *posin, unsigned int npoly, unsigned int N, unsigned int ei, unsigned int suffix){
	unsigned int i, j;		/* pop counter, poly counter */
	char Pout[32];				 /* String to hold filename in (Mutations) */
	FILE *ofp_poly = NULL;		 /* Pointer for data output */

	/* Printing out polymorphism table */
	sprintf(Pout,"Pop_%d.dat",ei + suffix);
	ofp_poly = fopen(Pout,"w");
	for(j = 0; j < npoly; j++){
		fprintf(ofp_poly,"%lf ",*(posin + j));
		for(i = 0; i < 2*N; i++){
			fprintf(ofp_poly,"%d ",*((*(neutin + i)) + j));
		}
		fprintf(ofp_poly,"\n");
	}
	fclose(ofp_poly);

}	/* End of population printing routine */

/* Printing polymorphism routines */
void polyprint(unsigned int **neutin, double *posin, unsigned int npoly, unsigned int N){
	unsigned int i, j;		/* pop counter, poly counter */

	for(j = 0; j < npoly; j++){
		printf("%lf ",*(posin + j));
		for(i = 0; i < 2*N; i++){
			printf("%d ",*((*(neutin + i)) + j));
		}
		printf("\n");
	}
	Wait();

}	/* End of reassignment routine */

/* Cutting out non-polymorphic sites */
void ptrim(unsigned int size, unsigned int **neutin, unsigned int **neutout, double *posin, double *posout, unsigned int npoly, unsigned int *npolyT, double *avpi){
	unsigned int i, j;
	unsigned int newp = 0;		/* New number of polymorphisms */
	unsigned int count = 0;		/* Polymorphism count */
	double cumpi = 0;			/* Cumulative pi (diversity) */

/*	printf("Trimming routine activated\n");	*/
	for(j = 0; j < npoly; j++){
		count = 0;
		*(posout + newp) = *(posin + j);
		for(i = 0; i < size; i++){
			*((*(neutout + i)) + newp) = *((*(neutin + i)) + j);
			count += *((*(neutin + i)) + j);
		}
/*		printf("For poly %d, count is %d\n",j,count);*/
		if((count > 0) && (count < size)){
/*			printf("Keeping poly %d at location %lf\n",j,*(posin + j));	*/
			newp++;
			cumpi += ((count*(size-count))/(1.0*size*size));
/*			printf("Count is %d, freq is %lf, cumulative value is %lf\n",count,count/(1.0*size),cumpi);	*/
		}

		/*
		if(count==size){
			printf("A fixed neutral mutation here\n");
		}
		*/

	}

	*npolyT = newp;
	*avpi = (cumpi/(1.0*newp));
}

/* Count number of unique haplotypes */
unsigned int uniqueH(unsigned int **neutin, unsigned int npoly, unsigned int N){

	unsigned int i, j, k;
	unsigned int unH = 0;			/* Count of unique haplotypes */
	unsigned int nomatch = 0;		/* Marker to show if sites do not match */

	unsigned int **UH = calloc(2*N,sizeof(unsigned int *));
	for(i = 0; i < 2*N; i++){
		UH[i] = calloc(npoly,sizeof(unsigned int));
	}

	/* Assigning first unique haplotype */
	for(j = 0; j < npoly; j++){
		*((*(UH + 0)) + j) = *((*(neutin + 0)) + j);
	}
	unH++;

	/* Now testing other haplotypes */
	for(i = 0; i < 2*N; i++){
		nomatch = 0;
		for(k = 0; k < unH; k++){
			for(j = 0; j < npoly; j++){
				if(*((*(UH + k)) + j) != *((*(neutin + i)) + j)){
					nomatch++;
					break;
				}
			}
		}

		if(nomatch == unH){		/* If haplotype does not match ANY unique ones, then add as new */
			for(j = 0; j < npoly; j++){
				*((*(UH + unH)) + j) = *((*(neutin + i)) + j);
			}
			unH++;
		}
	}

	for(i = 0; i < 2*N; i++){
		free(UH[i]);
	}
	free(UH);

	return unH;

}

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

	/* Declare variables here */
	unsigned int a, i;			/* Counters */
	unsigned int N = 0;			/* Population Size */
	unsigned int npoly = 0;		/* Number of neutral alleles */
	unsigned int npolyT = 0;	/* Number of neutral alleles (after cutting out fixed sites) */
	unsigned int suffix = 0;	/* Number of times to subsample from population */
	unsigned int base = 0;		/* Number of baseline forward-in-time simulations */
/*	unsigned int n = 0;			sprintf counter */
	unsigned int ttot = 0;		/* Total time of simulation */
/*	unsigned int unH = 0;		 Number of unique haplotypes */
	unsigned int tbi = 0;		/* Burn-in time */
/*	double npr = 0;				 Number of times to print out number of stats */
	double self = 0;			/* Selfing rate */
	double Rin = 0;				/* Recombination rate (input) */
	double R = 0;				/* Recombination rate (at a time) */
	double theta = 0;			/* Rate of neutral mutation */
	double avpi = 0;
	char Sout[32];				/* String to hold filename in (Seed) */
	FILE *ofp_sd;				/* Pointer for seed output */

	/* GSL random number definitions */
	const gsl_rng_type * T;
	gsl_rng * r;

	/* Reading in data from command line */
	if(argc != 7){
		fprintf(stderr,"Not enough inputs (need: N self R theta basereps suffix).\n");
		exit(1);
	}

	N = atoi(argv[1]);
	if(N <= 0){
		fprintf(stderr,"Total Population size N is zero or negative, not allowed.\n");
		exit(1);
	}

	self = strtod(argv[2],NULL);
	if(self < 0 || self > 1){
		fprintf(stderr,"Selfing rate must lie between 0 and 1 (inclusive).\n");
		exit(1);
	}

	Rin = strtod(argv[3],NULL);
	if(Rin < 0){
		fprintf(stderr,"Recombination rate must be positive or zero.\n");
		exit(1);
	}
	Rin /= 2.0*N;

	theta = strtod(argv[4],NULL);
	if(theta <= 0){
		fprintf(stderr,"Mutation rate must be a positive value.\n");
		exit(1);
	}

	base = atoi(argv[5]);
	if(base <= 0){
		fprintf(stderr,"Number of baseline simulations must be greater than zero.\n");
		exit(1);
	}

	suffix = atoi(argv[6]);
	if(argv[6] < 0){
		fprintf(stderr,"File index must be greater than or equal to zero.\n");
		exit(1);
	}

	/* create a generator chosen by the
    environment variable GSL_RNG_TYPE */

    mkdir("SeedsPop/", 0777);
	gsl_rng_env_setup();
	if (!getenv("GSL_RNG_SEED")) gsl_rng_default_seed = time(0);
	T = gsl_rng_default;
	r = gsl_rng_alloc(T);
	sprintf(Sout,"SeedsPop/Seed_%d.dat",suffix);
	ofp_sd = fopen(Sout,"w");
	fprintf(ofp_sd,"%lu\n",gsl_rng_default_seed);
	fclose(ofp_sd);

	/* Vector of 0 to (N-1) for sampling random numbers */
	unsigned int *nums = calloc(N,sizeof(unsigned int));
	unsigned int *nums2 = calloc(2*N,sizeof(unsigned int));
	for(a = 0; a < 2*N; a++){
		*(nums2 + a) = a;
		if(a < N){
			*(nums + a) = a;
		}
	}

	/* Executing simulation */
	for(i = 0; i < base; i++){

/*		printf("Starting run %d\n",i);*/

		/* Setting up selection, neutral tables per individual */
		double *polypos = calloc(INITTS,sizeof(double));						/* Position of neutral mutations */
		double *polyposF = calloc(INITTS,sizeof(double));						/* Position of neutral mutations (final for polymorphism sampling) */
		unsigned int **neutindvP = calloc(2*N,sizeof(unsigned int *));			/* Table of neutral markers per individual (parents) */
		unsigned int **neutindvO = calloc(2*N,sizeof(unsigned int *));			/* Table of neutral markers per individual (offspring) */
		unsigned int **neutindvF = calloc(2*N,sizeof(unsigned int *));			/* Table of neutral markers per individual (final for sampling) */
		for(a = 0; a < 2*N; a++){
			neutindvP[a] = calloc(INITTS,sizeof(unsigned int));
			neutindvO[a] = calloc(INITTS,sizeof(unsigned int));
			neutindvF[a] = calloc(INITTS,sizeof(unsigned int));
		}
		npoly = 0;
		ttot = 0;
		R = Rin;
		tbi = 20*N;
/*		npr = 20.0;	*/

		while(ttot < tbi){

			/*
			int tick = (tbi)/(npr);
			if(ttot%tick == 0){
				unH = uniqueH(neutindvP, npoly, N);
				printf("Time is %d; npoly is %d; av pi is %lf, unique haps are %d\n",ttot,npoly,avpi,unH);
			}
			*/

			/* Generating new population */
			generation_BI(N,self,R,nums,neutindvP,neutindvO,npoly,polypos,r);

			/* Neutral Mutation phase */
			addpoly(N, neutindvO, polypos, &npoly, theta, r);

			/* Reassigning matrices */
			reassign2(neutindvO, neutindvP, polyposF, polypos, npoly, N, 0);

			ttot++;
			ptrim(2*N, neutindvP, neutindvF, polypos, polyposF, npoly, &npolyT, &avpi);
			npoly = npolyT;
			reassign2(neutindvF, neutindvP, polyposF, polypos, npoly, N, 1);

			/*
			printf("Number of unique haplotypes are %d\n",uniqueH(neutindvP, npoly, N));
			polyprint(neutindvP, polypos, npoly, N);
			*/
		}
/*		printf("End of burn-in\n");*/

		/* After burn-in, print out neutral mutants  */
		popprint(neutindvP, polypos, npoly, N, i, suffix);

		/*
		reassign2(neutindvP, neutindvA, polypos, polyposA, npoly, N);
		npolyA = npoly;
		*/

		for(a = 0; a < 2*N; a++){
			free(neutindvF[a]);
			free(neutindvO[a]);
			free(neutindvP[a]);
		}
		free(neutindvF);
		free(neutindvO);
		free(neutindvP);
		free(polyposF);
		free(polypos);

	}

	free(nums2);
	free(nums);
	gsl_rng_free(r);

	return 0;

}	/* End of main program */

/* End of File */

max_stars_repo_name: MattHartfield/DomSelfAdapt

max_issues_repo_name: MattHartfield/DomSelfAdapt

max_forks_repo_name: MattHartfield/DomSelfAdapt


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

#ifndef __GSL_VECTOR_SHORT_H__
#define __GSL_VECTOR_SHORT_H__

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_block_short.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 size;
  size_t stride;
  short *data;
  gsl_block_short *block;
  int owner;
}
gsl_vector_short;

typedef struct
{
  gsl_vector_short vector;
} _gsl_vector_short_view;

typedef _gsl_vector_short_view gsl_vector_short_view;

typedef struct
{
  gsl_vector_short vector;
} _gsl_vector_short_const_view;

typedef const _gsl_vector_short_const_view gsl_vector_short_const_view;

/* Allocation */

gsl_vector_short *gsl_vector_short_alloc (const size_t n);
gsl_vector_short *gsl_vector_short_calloc (const size_t n);

gsl_vector_short *gsl_vector_short_alloc_from_block (gsl_block_short * b,
                                                     const size_t offset,
                                                     const size_t n,
                                                     const size_t stride);

gsl_vector_short *gsl_vector_short_alloc_from_vector (gsl_vector_short * v,
                                                      const size_t offset,
                                                      const size_t n,
                                                      const size_t stride);

void gsl_vector_short_free (gsl_vector_short * v);

/* Views */

_gsl_vector_short_view
gsl_vector_short_view_array (short *v, size_t n);

_gsl_vector_short_view
gsl_vector_short_view_array_with_stride (short *base,
                                         size_t stride,
                                         size_t n);

_gsl_vector_short_const_view
gsl_vector_short_const_view_array (const short *v, size_t n);

_gsl_vector_short_const_view
gsl_vector_short_const_view_array_with_stride (const short *base,
                                               size_t stride,
                                               size_t n);

_gsl_vector_short_view
gsl_vector_short_subvector (gsl_vector_short *v,
                            size_t i,
                            size_t n);

_gsl_vector_short_view
gsl_vector_short_subvector_with_stride (gsl_vector_short *v,
                                        size_t i,
                                        size_t stride,
                                        size_t n);

_gsl_vector_short_const_view
gsl_vector_short_const_subvector (const gsl_vector_short *v,
                                  size_t i,
                                  size_t n);

_gsl_vector_short_const_view
gsl_vector_short_const_subvector_with_stride (const gsl_vector_short *v,
                                              size_t i,
                                              size_t stride,
                                              size_t n);

/* Operations */

short gsl_vector_short_get (const gsl_vector_short * v, const size_t i);
void gsl_vector_short_set (gsl_vector_short * v, const size_t i, short x);

short *gsl_vector_short_ptr (gsl_vector_short * v, const size_t i);
const short *gsl_vector_short_const_ptr (const gsl_vector_short * v, const size_t i);

void gsl_vector_short_set_zero (gsl_vector_short * v);
void gsl_vector_short_set_all (gsl_vector_short * v, short x);
int gsl_vector_short_set_basis (gsl_vector_short * v, size_t i);

int gsl_vector_short_fread (FILE * stream, gsl_vector_short * v);
int gsl_vector_short_fwrite (FILE * stream, const gsl_vector_short * v);
int gsl_vector_short_fscanf (FILE * stream, gsl_vector_short * v);
int gsl_vector_short_fprintf (FILE * stream, const gsl_vector_short * v,
                              const char *format);

int gsl_vector_short_memcpy (gsl_vector_short * dest, const gsl_vector_short * src);

int gsl_vector_short_reverse (gsl_vector_short * v);

int gsl_vector_short_swap (gsl_vector_short * v, gsl_vector_short * w);
int gsl_vector_short_swap_elements (gsl_vector_short * v, const size_t i, const size_t j);

short gsl_vector_short_max (const gsl_vector_short * v);
short gsl_vector_short_min (const gsl_vector_short * v);
void gsl_vector_short_minmax (const gsl_vector_short * v, short * min_out, short * max_out);

size_t gsl_vector_short_max_index (const gsl_vector_short * v);
size_t gsl_vector_short_min_index (const gsl_vector_short * v);
void gsl_vector_short_minmax_index (const gsl_vector_short * v, size_t * imin, size_t * imax);

int gsl_vector_short_add (gsl_vector_short * a, const gsl_vector_short * b);
int gsl_vector_short_sub (gsl_vector_short * a, const gsl_vector_short * b);
int gsl_vector_short_mul (gsl_vector_short * a, const gsl_vector_short * b);
int gsl_vector_short_div (gsl_vector_short * a, const gsl_vector_short * b);
int gsl_vector_short_scale (gsl_vector_short * a, const double x);
int gsl_vector_short_add_constant (gsl_vector_short * a, const double x);

int gsl_vector_short_isnull (const gsl_vector_short * v);

#ifdef HAVE_INLINE

extern inline
short
gsl_vector_short_get (const gsl_vector_short * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (i >= v->size)
    {
      GSL_ERROR_VAL ("index out of range", GSL_EINVAL, 0);
    }
#endif
  return v->data[i * v->stride];
}

extern inline
void
gsl_vector_short_set (gsl_vector_short * v, const size_t i, short x)
{
#if GSL_RANGE_CHECK
  if (i >= v->size)
    {
      GSL_ERROR_VOID ("index out of range", GSL_EINVAL);
    }
#endif
  v->data[i * v->stride] = x;
}

extern inline
short *
gsl_vector_short_ptr (gsl_vector_short * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (i >= v->size)
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return (short *) (v->data + i * v->stride);
}

extern inline
const short *
gsl_vector_short_const_ptr (const gsl_vector_short * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (i >= v->size)
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return (const short *) (v->data + i * v->stride);
}

#endif /* HAVE_INLINE */

__END_DECLS

#endif /* __GSL_VECTOR_SHORT_H__ */

max_stars_repo_name: andrewkern/segSiteHMM

max_issues_repo_name: andrewkern/segSiteHMM

max_forks_repo_name: andrewkern/segSiteHMM


text:
/*
  This file is part of the Astrometry.net suite.
  Copyright 2006-2008 Dustin Lang, Keir Mierle and Sam Roweis.

  The Astrometry.net suite 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, version 2.

  The Astrometry.net suite 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 the Astrometry.net suite ; if not, write to the Free Software
  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301 USA
*/

#ifndef AN_GSL_UTILS_H
#define AN_GSL_UTILS_H

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

void gslutils_use_error_system();

/**
 Solves a least-squares matrix equation
 A X_i = B_i

 For NB pairs of X_i, B_i.

 NOTE: THIS DESTROYS A!

 A: MxN matrix
 B: array of NB x length-M vectors
 X: must be an array big enough to hold NB vectors.
   (they will be length-N).
 resids: if non-NULL, must be an array big enough to hold NB vectors
   (they will be length-M).

 The result vectors are freshly allocated and should be freed with gsl_vector_free().
 */
int gslutils_solve_leastsquares(gsl_matrix* A, gsl_vector** B,
                                gsl_vector** X, gsl_vector** resids,
                                int NB);

/**
 Same as above, but using varargs.  There must be exactly 3 * NB additional
 arguments, in the order:

 B0, &X0, &resid0, B1, &X1, &resid1, ...

 ie, the types must be repeating triples of:

 gsl_vector* b, gsl_vector** x, gsl_vector** resid

 */
int gslutils_solve_leastsquares_v(gsl_matrix* A, int NB, ...);

// C = A B
void gslutils_matrix_multiply(gsl_matrix* C, const gsl_matrix* A, const gsl_matrix* B);

int gslutils_invert_3x3(const double* A, double* B);

#endif

max_stars_repo_name: simonct/CoreAstro

max_issues_repo_name: simonct/CoreAstro

max_forks_repo_name: simonct/CoreAstro


text:
#pragma once

#define GLM_ENABLE_EXPERIMENTAL
#include <glm/glm.hpp>
#include <gsl/gsl>

namespace Library
{
	class Mesh;

	template <typename T>
	struct VertexDeclaration
	{
		static constexpr uint32_t VertexSize() { return gsl::narrow_cast<uint32_t>(sizeof(T)); }
		static constexpr uint32_t VertexBufferByteWidth(size_t vertexCount) { return gsl::narrow_cast<uint32_t>(sizeof(T) * vertexCount); }
		static void CreateVertexBuffer(const gsl::span<const T>& vertices, GLuint& vertexBuffer);
	};

	struct VertexPosition final : public VertexDeclaration<VertexPosition>
    {
		glm::vec4 Position{ 0 };

		VertexPosition() = default;

		VertexPosition(const glm::vec4& position) :
			Position(position) { }

		static void CreateVertexBuffer(const Library::Mesh& mesh, GLuint& vertexBuffer);
		static void CreateVertexBuffer(const gsl::span<const VertexPosition>& vertices, GLuint& vertexBuffer)
		{
			VertexDeclaration::CreateVertexBuffer(vertices, vertexBuffer);
		}
    };

	struct VertexPositionColor final :public VertexDeclaration<VertexPositionColor>
    {
		glm::vec4 Position{ 0 };
		glm::vec4 Color{ 0 };

		VertexPositionColor() = default;

		VertexPositionColor(const glm::vec4& position, const glm::vec4& color) :
			Position(position), Color(color) { }

		static void CreateVertexBuffer(const Library::Mesh& mesh, GLuint& vertexBuffer, bool randomColor = false);
		static void CreateVertexBuffer(const gsl::span<const VertexPositionColor>& vertices, GLuint& vertexBuffer)
		{
			VertexDeclaration::CreateVertexBuffer(vertices, vertexBuffer);
		}
    };

	struct VertexPositionTexture final : public VertexDeclaration<VertexPositionTexture>
	{
		glm::vec4 Position{ 0 };
		glm::vec2 TextureCoordinates{ 0 };

		VertexPositionTexture() = default;

		VertexPositionTexture(const glm::vec4& position, const glm::vec2& textureCoordinates) :
			Position(position), TextureCoordinates(textureCoordinates) { }

		static void CreateVertexBuffer(const Library::Mesh& mesh, GLuint& vertexBuffer);
		static void CreateVertexBuffer(const gsl::span<const VertexPositionTexture>& vertices, GLuint& vertexBuffer)
		{
			VertexDeclaration::CreateVertexBuffer(vertices, vertexBuffer);
		}
	};

	struct VertexPositionNormal final : public VertexDeclaration<VertexPositionNormal>
	{
		glm::vec4 Position{ 0 };
		glm::vec3 Normal{ 0 };

		VertexPositionNormal() = default;

		VertexPositionNormal(const glm::vec4& position, const glm::vec3& normal) :
			Position(position), Normal(normal) { }

		static void CreateVertexBuffer(const Library::Mesh& mesh, GLuint& vertexBuffer);
		static void CreateVertexBuffer(const gsl::span<const VertexPositionNormal>& vertices, GLuint& vertexBuffer)
		{
			VertexDeclaration::CreateVertexBuffer(vertices, vertexBuffer);
		}
	};

	struct VertexPositionTextureNormal final : public VertexDeclaration<VertexPositionTextureNormal>
	{
		glm::vec4 Position{ 0 };
		glm::vec2 TextureCoordinates{ 0 };
		glm::vec3 Normal{ 0 };

		VertexPositionTextureNormal() = default;

		VertexPositionTextureNormal(const glm::vec4& position, const glm::vec2& textureCoordinates, const glm::vec3& normal) :
			Position(position), TextureCoordinates(textureCoordinates), Normal(normal) { }

		static void CreateVertexBuffer(const Library::Mesh& mesh, GLuint& vertexBuffer);
		static void CreateVertexBuffer(const gsl::span<const VertexPositionTextureNormal>& vertices, GLuint& vertexBuffer)
		{
			VertexDeclaration::CreateVertexBuffer(vertices, vertexBuffer);
		}
	};

	struct VertexPositionTextureNormalTangentBinormal final : public VertexDeclaration<VertexPositionTextureNormalTangentBinormal>
	{
		glm::vec4 Position{ 0 };
		glm::vec2 TextureCoordinates{ 0 };
		glm::vec3 Normal{ 0 };
		glm::vec3 Tangent{ 0 };
		glm::vec3 Binormal{ 0 };

		VertexPositionTextureNormalTangentBinormal() = default;

		VertexPositionTextureNormalTangentBinormal(const glm::vec4& position, const glm::vec2& textureCoordinates, const glm::vec3& normal, const glm::vec3& tangent, const glm::vec3& binormal) :
			Position(position), TextureCoordinates(textureCoordinates), Normal(normal), Tangent(tangent), Binormal(binormal) { }

		static void CreateVertexBuffer(const Library::Mesh& mesh, GLuint& vertexBuffer);
		static void CreateVertexBuffer(const gsl::span<const VertexPositionTextureNormalTangentBinormal>& vertices, GLuint& vertexBuffer)
		{
			VertexDeclaration::CreateVertexBuffer(vertices, vertexBuffer);
		}
	};
}

max_stars_repo_name: DakkyDaWolf/OpenGL

max_issues_repo_name: DakkyDaWolf/OpenGL

max_forks_repo_name: DakkyDaWolf/OpenGL


text:
#ifndef ConditionalTest_SERIAL_H
#define ConditionalTest_SERIAL_H
#include <ace/core/core.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_statistics_double.h>
#include <gsl/gsl_randist.h>
#include "conditionaltest.h"

class ConditionalTest::Serial : public EAbstractAnalyticSerial
{
    Q_OBJECT
public:
    explicit Serial(ConditionalTest* parent);
    virtual std::unique_ptr<EAbstractAnalyticBlock> execute(const EAbstractAnalyticBlock* block) override final;

    // helper functions
    int test(CCMatrix::Pair& ccmPair, qint32 clusterIndex, qint32& testIndex, qint32 featureIndex, qint32 labelIndex, QVector<QVector<double>>& pValues, QVector<QVector<double>>& r2);
    int prepAnxData(QString testLabel, int dataIndex, TESTTYPE testType);
    bool isEmpty(QVector<QVector<double>>& matrix);
    int clusterInfo(CCMatrix::Pair& ccmPair, int clusterIndex, QString label, TESTTYPE testType);

    // Binomial Tests
    double binomial();
    double testOne();
    double testTwo();

    // Hypergeometrix Test.
    double hypergeom(CCMatrix::Pair& ccmPair, int clusterIndex, QString test_label);

    // Regression Test
    void regression(QVector<QString> &amxInfo, CCMatrix::Pair& ccmPair, int clusterIndex, TESTTYPE testType, QVector<double>& results);
    double fTest(double chisq, gsl_matrix* X, gsl_vector* Y, gsl_matrix* cov, gsl_vector* C);

private:
    /*!
     * Pointer to the serials objects parent KNNAnalytic.
     */
    ConditionalTest* _base;
    /*!
     * Annotation matrix data for testing.
     */
    QVector<QString> _amxData;
    /*!
     * Category count.
     */
    qint32 _catCount {0};
    /*!
     * Category count in cluster.
     */
    qint32 _catInCluster {0};
    /*!
     * Size of the cluster.
     */
    qint32 _clusterSize {0};
};

#endif

max_stars_repo_name: mfayk/KINC

max_issues_repo_name: mfayk/KINC

max_forks_repo_name: mfayk/KINC


text:
#pragma once

#include <gsl/gsl>

namespace aas {

	struct SkelBone;

	/**
	* Represents a sphere for cloth collision purposes.
	*/
	struct CollisionSphere {
		int boneId;
		float radius;
		Matrix3x4 worldMatrix;
		Matrix3x4 worldMatrixInverse;
		std::unique_ptr<CollisionSphere> next;
	};

	/**
	* Represents a cylinder for cloth collision purposes.
	*/
	struct CollisionCylinder {
		int boneId;
		float radius;
		float height;
		Matrix3x4 worldMatrix;
		Matrix3x4 worldMatrixInverse;
		std::unique_ptr<CollisionCylinder> next;
	};

	struct CollisionGeometry {
		std::unique_ptr<CollisionSphere> firstSphere;
		std::unique_ptr<CollisionCylinder> firstCylinder;
	};

	/**
	 * Extract the cloth collision information from a list of bones.
	 * The information is parsed from the bone names.
	 */
	CollisionGeometry FindCollisionGeometry(gsl::span<SkelBone> bones);

}


text:
/* rng/ranf.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 James Theiler, Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 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 <math.h>
#include <stdlib.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_sys.h>

/* This is the CRAY RANF generator. The generator returns the
   upper 32 bits from each term of the sequence,

   x_{n+1} = (a x_n) mod m

   using 48-bit unsigned arithmetic, with a = 0x2875A2E7B175 and m =
   2^48. The seed specifies the lower 32 bits of the initial value,
   x_1, with the lowest bit set (to prevent the seed taking an even
   value), and the upper 16 bits set to 0.

   There is a subtlety in the implementation of the seed. The initial
   state is put one step back by multiplying by the modular inverse of
   a mod m. This is done for compatibility with the original CRAY
   implementation.

   Note, you can only seed the generator with integers up to 2^32,
   while the CRAY uses wide integers which can cover all 2^48 states
   of the generator.

   The theoretical value of x_{10001} is 141091827447341.

   The period of this generator is 2^{46}. */

static inline void ranf_advance (void *vstate);
static unsigned long int ranf_get (void *vstate);
static double ranf_get_double (void *vstate);
static void ranf_set (void *state, unsigned long int s);

static const unsigned short int a0 = 0xB175 ;
static const unsigned short int a1 = 0xA2E7 ;
static const unsigned short int a2 = 0x2875 ;

typedef struct
  {
    unsigned short int x0, x1, x2;
  }
ranf_state_t;

static inline void
ranf_advance (void *vstate)
{
  ranf_state_t *state = (ranf_state_t *) vstate;

  const unsigned long int x0 = (unsigned long int) state->x0 ;
  const unsigned long int x1 = (unsigned long int) state->x1 ;
  const unsigned long int x2 = (unsigned long int) state->x2 ;

  unsigned long int r ;

  r = a0 * x0 ;
  state->x0 = (r & 0xFFFF) ;

  r >>= 16 ;
  r += a0 * x1 + a1 * x0 ;
  state->x1 = (r & 0xFFFF) ;

  r >>= 16 ;
  r += a0 * x2 + a1 * x1 + a2 * x0 ;
  state->x2 = (r & 0xFFFF) ;
}

static unsigned long int
ranf_get (void *vstate)
{
  unsigned long int x1, x2;

  ranf_state_t *state = (ranf_state_t *) vstate;
  ranf_advance (state) ;

  x1 = (unsigned long int) state->x1;
  x2 = (unsigned long int) state->x2;

  return (x2 << 16) + x1;
}

static double
ranf_get_double (void * vstate)
{
  ranf_state_t *state = (ranf_state_t *) vstate;

  ranf_advance (state) ;

  return (ldexp((double) state->x2, -16)
          + ldexp((double) state->x1, -32)
          + ldexp((double) state->x0, -48)) ;
}

static void
ranf_set (void *vstate, unsigned long int s)
{
  ranf_state_t *state = (ranf_state_t *) vstate;

  unsigned short int x0, x1, x2 ;
  unsigned long int r ;

  unsigned long int b0 = 0xD6DD ;
  unsigned long int b1 = 0xB894 ;
  unsigned long int b2 = 0x5CEE ;

  if (s == 0)  /* default seed */
    {
      x0 = 0x9CD1 ;
      x1 = 0x53FC ;
      x2 = 0x9482 ;
    }
  else
    {
      x0 = (s | 1) & 0xFFFF ;
      x1 = s >> 16 & 0xFFFF ;
      x2 = 0 ;
    }

  r = b0 * x0 ;
  state->x0 = (r & 0xFFFF) ;

  r >>= 16 ;
  r += b0 * x1 + b1 * x0 ;
  state->x1 = (r & 0xFFFF) ;

  r >>= 16 ;
  r += b0 * x2 + b1 * x1 + b2 * x0 ;
  state->x2 = (r & 0xFFFF) ;

  return;
}

static const gsl_rng_type ranf_type =
{"ranf",                        /* name */
 0xffffffffUL,                  /* RAND_MAX */
 0,                             /* RAND_MIN */
 sizeof (ranf_state_t),
 &ranf_set,
 &ranf_get,
 &ranf_get_double
};

const gsl_rng_type *gsl_rng_ranf = &ranf_type;

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
/**
 * \author Sylvain Marsat, University of Maryland - NASA GSFC
 *
 * \brief C header for the instrumental noise for LISA-type detectors.
 *
 * Formulas taken from Królak&al gr-qc/0401108 (c.f. section III).
 *
 */

#ifndef _LISANOISE_H
#define _LISANOISE_H

#define _XOPEN_SOURCE 500

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

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

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

#include "constants.h"
#include "LISAgeometry.h"

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

/************************************************************************/
/****** Global variables storing min and max f for the noise PSD  *******/

/* Defines bounds in frequency beyond which we don't trust the instrument model anymore - all waveforms will be cut to this range */
/* Here extended range - allows for instance to taper the FD signal for f>1Hz */
#define __LISASimFD_Noise_fLow 1.e-6
#define __LISASimFD_Noise_fHigh 5.
/* Original, more conservative bounds */
//#define __LISASimFD_Noise_fLow 1.e-5
//#define __LISASimFD_Noise_fHigh 1.

/**************************************************************/
/****** Prototypes: functions evaluating the noise PSD  *******/

/* Function returning the relevant noise function, given a set of TDI observables and a channel */
ObjectFunction NoiseFunction(const LISAconstellation *variant, const TDItag tditag, const int nchan);

/* Noise Sn for TDI observables - factors have been scaled out both in the response and the noise */
double SnXYZ(const LISAconstellation *variant, double f);
double Snalphabetagamma(const LISAconstellation *variant, double f);
double SnAXYZ(const LISAconstellation *variant, double f);
double SnEXYZ(const LISAconstellation *variant, double f);
double SnTXYZ(const LISAconstellation *variant, double f);
double SnAalphabetagamma(const LISAconstellation *variant, double f);
double SnEalphabetagamma(const LISAconstellation *variant, double f);
double SnTalphabetagamma(const LISAconstellation *variant, double f);

/* Noise functions for AET(XYZ) without rescaling */
double SnAXYZNoRescaling(const LISAconstellation *variant, double f);
double SnEXYZNoRescaling(const LISAconstellation *variant, double f);
double SnTXYZNoRescaling(const LISAconstellation *variant, double f);

/* Function returning the relevant noise function, given a set of TDI observables and a channel */
/* double (*NoiseFunction(const TDItag tditag, const int chan))(double); */

/* The noise functions themselves */
/* double NoiseSnA(const double f); */
/* double NoiseSnE(const double f); */
/* double NoiseSnT(const double f); */

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

#endif /* _LISANOISE_H */

max_stars_repo_name: titodalcanton/flare

max_issues_repo_name: titodalcanton/flare

max_forks_repo_name: titodalcanton/flare


text:
#ifndef CONTROLLERS_IMAGEOPERATIONACTIONCONTROLLER_H
#define CONTROLLERS_IMAGEOPERATIONACTIONCONTROLLER_H

#include <QObject>

#include "s3d/cv/image_operation/image_operation.h"

#include <QAction>

#include <gsl/gsl>

class ImageOperationActionController : public QObject {
  Q_OBJECT

public:
  ImageOperationActionController(gsl::not_null<QAction*> action,
                                 gsl::not_null<s3d::image_operation::ImageOperation*> imageOperation);

  virtual void onActionTriggered() = 0;

signals:
  void imageOperationTriggered();

protected:
  QAction* m_action;
  s3d::image_operation::ImageOperation* m_imageOperation;
};

#endif //CONTROLLERS_IMAGEOPERATIONACTIONCONTROLLER_H


text:
/* movstat/gsl_movstat.h
 *
 * Copyright (C) 2018 Patrick Alken
 *
 * 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_MOVSTAT_H__
#define __GSL_MOVSTAT_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 <gsl/gsl_math.h>
#include <gsl/gsl_vector.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 enum
{
  GSL_MOVSTAT_END_PADZERO,
  GSL_MOVSTAT_END_PADVALUE,
  GSL_MOVSTAT_END_TRUNCATE
} gsl_movstat_end_t;

/* accumulator struct
 * size          - return number of bytes needed for accumulator with maximum of n elements
 * init          - initialize accumulator state
 * insert        - insert a single sample into accumulator; if there are already n
 *                 samples in accumulator, oldest sample is overwritten
 * delete_oldest - delete oldest sample from accumulator
 * get           - return accumulated value
 */
typedef struct
{
  size_t (*size) (const size_t n);
  int (*init) (const size_t n, void * vstate);
  int (*insert) (const double x, void * vstate);
  int (*delete_oldest) (void * vstate);
  int (*get) (void * params, double * result, const void * vstate);
} gsl_movstat_accum;

typedef struct
{
  double (* function) (const size_t n, double x[], void * params);
  void * params;
} gsl_movstat_function;

#define GSL_MOVSTAT_FN_EVAL(F,n,x) (*((F)->function))((n),(x),(F)->params)

/* workspace for moving window statistics */

typedef struct
{
  size_t H;          /* number of previous samples in window */
  size_t J;          /* number of after samples in window */
  size_t K;          /* window size K = H + J + 1 */
  double *work;      /* workspace, size K */
  void *state;       /* state workspace for various accumulators */
  size_t state_size; /* bytes allocated for 'state' */
} gsl_movstat_workspace;

/* alloc.c */

GSL_FUN gsl_movstat_workspace *gsl_movstat_alloc(const size_t K);
GSL_FUN gsl_movstat_workspace *gsl_movstat_alloc2(const size_t H, const size_t J);
GSL_FUN gsl_movstat_workspace *gsl_movstat_alloc_with_size(const size_t accum_state_size, const size_t H, const size_t J);
GSL_FUN void gsl_movstat_free(gsl_movstat_workspace * w);

/* apply.c */
GSL_FUN int gsl_movstat_apply_accum(const gsl_movstat_end_t endtype, const gsl_vector * x,
                            const gsl_movstat_accum * accum, void * accum_params,
                            gsl_vector * y, gsl_vector * z,
                            gsl_movstat_workspace * w);
GSL_FUN int gsl_movstat_apply(const gsl_movstat_end_t endtype, const gsl_movstat_function * F,
                      const gsl_vector * x, gsl_vector * y, gsl_movstat_workspace * w);

/* fill.c */
GSL_FUN size_t gsl_movstat_fill(const gsl_movstat_end_t endtype, const gsl_vector * x, const size_t idx,
                        const size_t H, const size_t J, double * window);

GSL_FUN int gsl_movstat_mean(const gsl_movstat_end_t endtype, const gsl_vector * x, gsl_vector * y, gsl_movstat_workspace * w);
GSL_FUN int gsl_movstat_variance(const gsl_movstat_end_t endtype, const gsl_vector * x, gsl_vector * y, gsl_movstat_workspace * w);
GSL_FUN int gsl_movstat_sd(const gsl_movstat_end_t endtype, const gsl_vector * x, gsl_vector * y, gsl_movstat_workspace * w);
GSL_FUN int gsl_movstat_median(const gsl_movstat_end_t endtype, const gsl_vector * x, gsl_vector * y, gsl_movstat_workspace * w);
GSL_FUN int gsl_movstat_min(const gsl_movstat_end_t endtype, const gsl_vector * x, gsl_vector * y, gsl_movstat_workspace * w);
GSL_FUN int gsl_movstat_max(const gsl_movstat_end_t endtype, const gsl_vector * x, gsl_vector * y, gsl_movstat_workspace * w);
GSL_FUN int gsl_movstat_minmax(const gsl_movstat_end_t endtype, const gsl_vector * x, gsl_vector * y_min, gsl_vector * y_max, gsl_movstat_workspace * w);
GSL_FUN int gsl_movstat_mad0(const gsl_movstat_end_t endtype, const gsl_vector * x, gsl_vector * xmedian,
                     gsl_vector * xmad, gsl_movstat_workspace * w);
GSL_FUN int gsl_movstat_mad(const gsl_movstat_end_t endtype, const gsl_vector * x, gsl_vector * xmedian,
                    gsl_vector * xmad, gsl_movstat_workspace * w);
GSL_FUN int gsl_movstat_qqr(const gsl_movstat_end_t endtype, const gsl_vector * x, const double q,
                    gsl_vector * xqqr, gsl_movstat_workspace * w);
GSL_FUN int gsl_movstat_Sn(const gsl_movstat_end_t endtype, const gsl_vector * x,
                   gsl_vector * xscale, gsl_movstat_workspace * w);
GSL_FUN int gsl_movstat_Qn(const gsl_movstat_end_t endtype, const gsl_vector * x,
                   gsl_vector * xscale, gsl_movstat_workspace * w);
GSL_FUN int gsl_movstat_sum(const gsl_movstat_end_t endtype, const gsl_vector * x, gsl_vector * y, gsl_movstat_workspace * w);

/* accumulator variables */

GSL_VAR const gsl_movstat_accum * gsl_movstat_accum_mad;
GSL_VAR const gsl_movstat_accum * gsl_movstat_accum_max;
GSL_VAR const gsl_movstat_accum * gsl_movstat_accum_mean;
GSL_VAR const gsl_movstat_accum * gsl_movstat_accum_median;
GSL_VAR const gsl_movstat_accum * gsl_movstat_accum_min;
GSL_VAR const gsl_movstat_accum * gsl_movstat_accum_minmax;
GSL_VAR const gsl_movstat_accum * gsl_movstat_accum_sd;
GSL_VAR const gsl_movstat_accum * gsl_movstat_accum_Sn;
GSL_VAR const gsl_movstat_accum * gsl_movstat_accum_sum;
GSL_VAR const gsl_movstat_accum * gsl_movstat_accum_Qn;
GSL_VAR const gsl_movstat_accum * gsl_movstat_accum_qqr;
GSL_VAR const gsl_movstat_accum * gsl_movstat_accum_userfunc;
GSL_VAR const gsl_movstat_accum * gsl_movstat_accum_variance;

__END_DECLS

#endif /* __GSL_MOVSTAT_H__ */

max_stars_repo_name: mgreter/astrometrylib

max_issues_repo_name: zzpwahaha/VimbaCamJILA

max_forks_repo_name: zzpwahaha/VimbaCamJILA


text:
#include <pygsl/utils.h>
#include <pygsl/error_helpers.h>
#include <pygsl/block_helpers.h>
#include <gsl/gsl_sf.h>
#include <gsl/gsl_nan.h>

#ifndef IMPORTALL
static PyObject *module=NULL;
#endif

typedef  int (array_p_evaluator_iid_ad)(int nmin, int nmax, double x, double * result_array);
static PyObject*
PyGSL_sf_array_evaluator_legendre_iid_ad(PyObject *self, PyObject *args,
					 array_p_evaluator_iid_ad * eval)
{
     PyArrayObject *result = NULL;
     int lmax=0, m=0, dimension = 0, ret;
     double x=0, *data=NULL;

     FUNC_MESS_BEGIN();
     if (!PyArg_ParseTuple(args, "iid", &lmax, &m, &x)){
	  return NULL;
     }
     if(m < 0){
	  PyErr_SetString(PyExc_ValueError,
			  "Nmin must be bigger than 0!");
	  return NULL;
     }
     if(lmax < m){
	  PyErr_SetString(PyExc_ValueError,
			  "Nmax must be bigger or equal to nmin!");
     }
     dimension = gsl_sf_legendre_array_size(lmax, m);
     result = (PyArrayObject *) PyArray_FromDims(1, &dimension, PyArray_DOUBLE);
     if(result == NULL)
	  return NULL;

     data = (double *) result->data;
     ret =  eval(lmax, m, x, data);
     if(PyGSL_ERROR_FLAG(ret) != GSL_SUCCESS)
	  goto fail;

     FUNC_MESS_END();
     return (PyObject *) result;

 fail:
     Py_XDECREF(result);
     return NULL;
}

static PyObject*
PyGSL_sf_array_evaluator_iid_ad(PyObject *self, PyObject *args,
				array_p_evaluator_iid_ad * eval)
{
     PyArrayObject *result = NULL;
     int nmin=0, nmax=0, dimension = 0, ret;
     double x=0, *data=NULL;

     FUNC_MESS_BEGIN();
     if (!PyArg_ParseTuple(args, "iid", &nmin, &nmax, &x)){
	  return NULL;
     }
     if(nmin < 0){
	  PyErr_SetString(PyExc_ImportError,
			  "Nmin must be bigger than 0!");
	  return NULL;
     }
     if(nmax < nmin){
	  PyErr_SetString(PyExc_ImportError,
			  "Nmax must be bigger or equal to nmin!");
     }
     dimension = nmax - nmin + 1; /* Goes form nmin to nmax, both included */
     result = (PyArrayObject *) PyArray_FromDims(1, &dimension, PyArray_DOUBLE);
     if(result == NULL)
	  return NULL;

     data = (double *) result->data;
     ret =  eval(nmin, nmax, x, data);
     if(PyGSL_ERROR_FLAG(ret) != GSL_SUCCESS)
	  goto fail;

     FUNC_MESS_END();
     return (PyObject *) result;

 fail:
     Py_XDECREF(result);
     return NULL;
}

typedef  int (array_p_evaluator_id_ad)(int nmax, double x, double * result_array);

static PyObject*
PyGSL_sf_array_evaluator_id_ad(PyObject *self, PyObject *args, array_p_evaluator_id_ad * eval)
{
     PyArrayObject *result = NULL;
     int nmin=0, nmax=0, dimension = 0, ret;
     double x=0, *data=NULL;

     FUNC_MESS_BEGIN();
     if (!PyArg_ParseTuple(args, "id", &nmax, &x)){
	  return NULL;
     }
     if(nmin < 0){
	  PyErr_SetString(PyExc_ImportError,
			  "Nmin must be bigger than 0!");
	  return NULL;
     }
     dimension = nmax - nmin + 1; /* Goes form nmin to nmax, both included */
     result = (PyArrayObject *) PyArray_FromDims(1, &dimension, PyArray_DOUBLE);
     if(result == NULL)
	  return NULL;

     data = (double *) result->data;
     ret =  eval(nmax, x, data);
     if(PyGSL_ERROR_FLAG(ret) != GSL_SUCCESS)
	  goto fail;

     FUNC_MESS_END();
     return (PyObject *) result;

 fail:
     Py_XDECREF(result);
     return NULL;
}

typedef  int (array_p_evaluator_idd_ad)(int nmax, double x1, double x2, double * result_array);
static PyObject*
PyGSL_sf_array_evaluator_idd_ad(PyObject *self, PyObject *args, array_p_evaluator_idd_ad * eval)
{
     PyArrayObject *result = NULL;
     int nmin=0, nmax=0, dimension = 0, ret;
     double x=0, x1=0, *data=NULL;

     FUNC_MESS_BEGIN();
     if (!PyArg_ParseTuple(args, "idd", &nmax, &x, &x1)){
	  return NULL;
     }
     if(nmin < 0){
	  PyErr_SetString(PyExc_ImportError,
			  "Nmin must be bigger than 0!");
	  return NULL;
     }
     dimension = nmax - nmin + 1; /* Goes form nmin to nmax, both included */
     result = (PyArrayObject *) PyArray_FromDims(1, &dimension, PyArray_DOUBLE);
     if(result == NULL)
	  return NULL;

     data = (double *) result->data;
     ret =  eval(nmax, x, x1, data);
     if(PyGSL_ERROR_FLAG(ret) != GSL_SUCCESS)
	  goto fail;

     FUNC_MESS_END();
     return (PyObject *) result;

 fail:
     Py_XDECREF(result);
     return NULL;
}
typedef  int (array_p_evaluator_did_ad)( double x1, int nmax, double x2, double * result_array);
static PyObject*
PyGSL_sf_array_evaluator_did_ad(PyObject *self, PyObject *args, array_p_evaluator_did_ad * eval)
{
     PyArrayObject *result = NULL;
     int nmin=0, nmax=0, dimension = 0, ret;
     double x=0, x1=0, *data=NULL;

     FUNC_MESS_BEGIN();
     if (!PyArg_ParseTuple(args, "did",&x, &nmax, &x1)){
	  return NULL;
     }
     dimension = nmax - nmin + 1; /* Goes form nmin to nmax, both included */
     result = (PyArrayObject *) PyArray_FromDims(1, &dimension, PyArray_DOUBLE);
     if(result == NULL)
	  return NULL;

     data = (double *) result->data;
     ret =  eval(nmax, x, x1, data);
     if(PyGSL_ERROR_FLAG(ret) != GSL_SUCCESS)
	  goto fail;

     FUNC_MESS_END();
     return (PyObject *) result;

 fail:
     Py_XDECREF(result);
     return NULL;
}

typedef int (array_p_evaluator_didd_add)(double , int , double , double , double * array, double*);
static PyObject*
PyGSL_sf_array_evaluator_didd_add(PyObject *self, PyObject *args, array_p_evaluator_didd_add * eval)
{
     PyArrayObject *result = NULL;
     int nmin=0, nmax=0, dimension = 0, ret;
     double x=0, x1=0, *data=NULL, l_min, exponent;

     FUNC_MESS_BEGIN();
     if (!PyArg_ParseTuple(args, "didd", &l_min, &nmax, &x, &x1)){
	  return NULL;
     }
     if(nmin < 0){
	  PyErr_SetString(PyExc_ImportError,
			  "Nmin must be bigger than 0!");
	  return NULL;
     }
     dimension = nmax - nmin + 1; /* Goes form nmin to nmax, both included */
     result = (PyArrayObject *) PyArray_FromDims(1, &dimension, PyArray_DOUBLE);
     if(result == NULL)
	  return NULL;

     data = (double *) result->data;
     ret =  eval(l_min, nmax, x, x1, data, &exponent);
     FUNC_MESS_END();
     if(PyGSL_ERROR_FLAG(ret) != GSL_SUCCESS)
	  goto fail;

     return Py_BuildValue("Od",result,exponent);
 fail:
     Py_XDECREF(result);
     return NULL;
}

typedef int (array_p_evaluator_didd_addadd)(double , int , double , double , double * array1, double*, double *array2, double*);
static PyObject*
PyGSL_sf_array_evaluator_didd_addadd(PyObject *self, PyObject *args, array_p_evaluator_didd_addadd * eval)
{
     PyArrayObject *result1 = NULL,*result2 = NULL;
     int nmin=0, nmax=0, dimension = 0, ret;
     double x=0, x1=0, *data1=NULL, *data2=NULL, l_min, exponent1,exponent2;

     FUNC_MESS_BEGIN();
     if (!PyArg_ParseTuple(args, "didd", &l_min, &nmax, &x, &x1)){
	  return NULL;
     }
     if(nmin < 0){
	  PyErr_SetString(PyExc_ImportError,
			  "Nmin must be bigger than 0!");
	  return NULL;
     }
     dimension = nmax - nmin + 1; /* Goes form nmin to nmax, both included */
     result1 = (PyArrayObject *) PyArray_FromDims(1, &dimension, PyArray_DOUBLE);
     if(result1 == NULL)
	  goto fail;
     result2 = (PyArrayObject *) PyArray_FromDims(1, &dimension, PyArray_DOUBLE);
     if(result2 == NULL)
	  goto fail;

     data1 = (double *) result1->data;
     data2 = (double *) result2->data;

     ret =  eval(l_min, nmax, x, x1, data1, &exponent1, data2, &exponent2);
     FUNC_MESS_END();
     if(PyGSL_ERROR_FLAG(ret) != GSL_SUCCESS)
	  goto fail;

     return Py_BuildValue("OdOd",result1,exponent1,result2,exponent2);
 fail:
     Py_XDECREF(result1);
     Py_XDECREF(result2);
     return NULL;
}

typedef int (array_p_evaluator_didd_adadadaddd)(double , int , double , double , double * a1, double * a2, double * a3, double * a4,
						double*, double*);
static PyObject*
PyGSL_sf_array_evaluator_didd_adadadaddd(PyObject *self, PyObject *args, array_p_evaluator_didd_adadadaddd * eval)
{
     PyArrayObject *result1 = NULL,*result2 = NULL, *result3=NULL,*result4=NULL;
     int nmin=0, nmax=0, dimension = 0, ret;
     double x=0, x1=0, l_min, exponent1, exponent2;

     FUNC_MESS_BEGIN();
     if (!PyArg_ParseTuple(args, "didd", &l_min, &nmax, &x, &x1)){
	  return NULL;
     }
     if(nmin < 0){
	  PyErr_SetString(PyExc_ImportError,
			  "Nmin must be bigger than 0!");
	  return NULL;
     }
     dimension = nmax - nmin + 1; /* Goes form nmin to nmax, both included */
     result1 = (PyArrayObject *) PyArray_FromDims(1, &dimension, PyArray_DOUBLE);
     if(result1 == NULL)
	  goto fail;
     result2 = (PyArrayObject *) PyArray_FromDims(1, &dimension, PyArray_DOUBLE);
     if(result2 == NULL)
	  goto fail;
     result3 = (PyArrayObject *) PyArray_FromDims(1, &dimension, PyArray_DOUBLE);
     if(result3 == NULL)
	  goto fail;
     result4 = (PyArrayObject *) PyArray_FromDims(1, &dimension, PyArray_DOUBLE);
     if(result4 == NULL)
	  goto fail;

     ret =  eval(l_min, nmax, x, x1, (double *) result1->data, (double *)result2->data,
		 (double *) result3->data, (double *)result4->data, &exponent1, &exponent2);
     FUNC_MESS_END();
     if(PyGSL_ERROR_FLAG(ret) != GSL_SUCCESS)
	  goto fail;

     return Py_BuildValue("OOOOdd",result1,result2,result3,result4,exponent1, exponent2);
 fail:
     Py_XDECREF(result1);
     Py_XDECREF(result2);
     Py_XDECREF(result3);
     Py_XDECREF(result4);
     return NULL;
}

#define SF_ARRAY(name, function)                                           \
static PyObject* sf_ ## name      (PyObject *self, PyObject *args)         \
{                                                                          \
     PyObject * tmp;                                                       \
     FUNC_MESS_BEGIN();                                                    \
     tmp =  PyGSL_sf_array_evaluator_ ## function (self, args, gsl_sf_ ## name); \
     if (tmp == NULL){                                                     \
	  PyGSL_add_traceback(module, __FILE__, __FUNCTION__,  __LINE__);  \
     }                                                                     \
     FUNC_MESS_END();                                                      \
     return tmp;                                                           \
}

SF_ARRAY(bessel_Jn_array, iid_ad);
SF_ARRAY(bessel_Yn_array, iid_ad);
SF_ARRAY(bessel_In_array, iid_ad);
SF_ARRAY(bessel_Kn_array, iid_ad);
SF_ARRAY(bessel_Kn_scaled_array, iid_ad);
SF_ARRAY(bessel_jl_array, id_ad);
SF_ARRAY(bessel_jl_steed_array, id_ad);
SF_ARRAY(bessel_yl_array, id_ad);
SF_ARRAY(bessel_il_scaled_array, id_ad);
SF_ARRAY(bessel_kl_scaled_array, id_ad);
SF_ARRAY(gegenpoly_array, idd_ad);
SF_ARRAY(legendre_H3d_array, idd_ad);
SF_ARRAY(coulomb_wave_F_array, didd_add);
SF_ARRAY(coulomb_wave_sphF_array, didd_add);
SF_ARRAY(coulomb_wave_FG_array, didd_addadd);
SF_ARRAY(coulomb_wave_FGp_array, didd_adadadaddd);
SF_ARRAY(coulomb_CL_array, did_ad);
SF_ARRAY(legendre_Pl_array, id_ad);
SF_ARRAY(legendre_Plm_array, legendre_iid_ad);
SF_ARRAY(legendre_sphPlm_array, legendre_iid_ad);

static PyMethodDef sf_array_functions[] = {
     {"bessel_Jn_array",         (PyCFunction) sf_bessel_Jn_array,         METH_VARARGS, NULL},
     {"bessel_Yn_array",         (PyCFunction) sf_bessel_Yn_array,         METH_VARARGS, NULL},
     {"bessel_In_array",         (PyCFunction) sf_bessel_In_array,         METH_VARARGS, NULL},
     {"bessel_Kn_array",         (PyCFunction) sf_bessel_Kn_array,         METH_VARARGS, NULL},
     {"bessel_Kn_scaled_array",  (PyCFunction) sf_bessel_Kn_scaled_array,  METH_VARARGS, NULL},
     {"bessel_jl_array",         (PyCFunction) sf_bessel_jl_array,         METH_VARARGS, NULL},
     {"bessel_jl_steed_array",   (PyCFunction) sf_bessel_jl_steed_array,   METH_VARARGS, NULL},
     {"bessel_yl_array",         (PyCFunction) sf_bessel_yl_array,         METH_VARARGS, NULL},
     {"bessel_il_scaled_array",  (PyCFunction) sf_bessel_il_scaled_array,  METH_VARARGS, NULL},
     {"bessel_kl_scaled_array",  (PyCFunction) sf_bessel_kl_scaled_array,  METH_VARARGS, NULL},
     {"gegenpoly_array",         (PyCFunction) sf_gegenpoly_array,         METH_VARARGS, NULL},
     {"legendre_H3d_array",      (PyCFunction) sf_legendre_H3d_array,      METH_VARARGS, NULL},
     {"coulomb_wave_F_array",    (PyCFunction) sf_coulomb_wave_F_array,    METH_VARARGS, NULL},
     {"coulomb_wave_sphF_array", (PyCFunction) sf_coulomb_wave_sphF_array, METH_VARARGS, NULL},
     {"coulomb_wave_FG_array",   (PyCFunction) sf_coulomb_wave_FG_array,   METH_VARARGS, NULL},
     {"coulomb_wave_FGp_array",  (PyCFunction) sf_coulomb_wave_FGp_array,  METH_VARARGS, NULL},
     {"coulomb_CL_array",        (PyCFunction) sf_coulomb_CL_array,        METH_VARARGS, NULL},
     {"legendre_Pl_array",       (PyCFunction) sf_legendre_Pl_array,       METH_VARARGS, NULL},
     {"legendre_Plm_array",      (PyCFunction) sf_legendre_Plm_array,      METH_VARARGS, NULL},
     {"legendre_sphPlm_array",   (PyCFunction) sf_legendre_sphPlm_array,   METH_VARARGS, NULL},
     {NULL, NULL, 0}
};

#ifndef IMPORTALL

DL_EXPORT(void) initsfarray(void)
{
     module = Py_InitModule("sfarray", sf_array_functions);
     import_array();
       init_pygsl();

}
#endif

max_stars_repo_name: juhnowski/FishingRod

max_issues_repo_name: juhnowski/FishingRod

max_forks_repo_name: juhnowski/FishingRod


text:
/* dht/test_dht.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 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.
 */

/* Author:  G. Jungman
 */
#include <config.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_dht.h>

/* Test exact small transform.
 */
int
test_dht_exact(void)
{
  int stat = 0;
  double f_in[3] = { 1.0, 2.0, 3.0 };
  double f_out[3];
  gsl_dht * t = gsl_dht_new(3, 1.0, 1.0);
  gsl_dht_apply(t, f_in, f_out);

  /* Check values. */
  if(fabs( f_out[0]-( 0.375254649407520))/0.375254649407520 > 1.0e-14) stat++;
  if(fabs( f_out[1]-(-0.133507872695560))/0.133507872695560 > 1.0e-14) stat++;
  if(fabs( f_out[2]-( 0.044679925143840))/0.044679925143840 > 1.0e-14) stat++;

  /* Check inverse.
   * We have to adjust the normalization
   * so we can use the same precalculated transform.
   */
  gsl_dht_apply(t, f_out, f_in);
  f_in[0] *= 13.323691936314223*13.323691936314223;  /* jzero[1,4]^2 */
  f_in[1] *= 13.323691936314223*13.323691936314223;
  f_in[2] *= 13.323691936314223*13.323691936314223;

  /* The loss of precision on the inverse
   * is a little surprising. However, this
   * thing is quite tricky since the band-limited
   * function represented by the samples {1,2,3}
   * need not be very nice. Like in any spectral
   * application, you really have to have some
   * a-priori knowledge of the underlying function.
   */
  if(fabs( f_in[0]-1.0)/1.0 > 2.0e-05) stat++;
  if(fabs( f_in[1]-2.0)/2.0 > 2.0e-05) stat++;
  if(fabs( f_in[2]-3.0)/3.0 > 2.0e-05) stat++;

  gsl_dht_free(t);

  return stat;
}

/* Test the transform
 * Integrate[x J_0(a x) / (x^2 + 1), {x,0,Inf}] = K_0(a)
 */
int
test_dht_simple(void)
{
  int stat = 0;
  int n;
  double f_in[128];
  double f_out[128];
  gsl_dht * t = gsl_dht_new(128, 0.0, 100.0);

  for(n=0; n<128; n++) {
    const double x = gsl_dht_x_sample(t, n);
    f_in[n] = 1.0/(1.0+x*x);
  }

  gsl_dht_apply(t, f_in, f_out);

  /* This is a difficult transform to calculate this way,
   * since it does not satisfy the boundary condition and
   * it dies quite slowly. So it is not meaningful to
   * compare this to high accuracy. We only check
   * that it seems to be working.
   */
  if(fabs( f_out[0]-4.00)/4.00 > 0.02) stat++;
  if(fabs( f_out[5]-1.84)/1.84 > 0.02) stat++;
  if(fabs(f_out[10]-1.27)/1.27 > 0.02) stat++;
  if(fabs(f_out[35]-0.352)/0.352 > 0.02) stat++;
  if(fabs(f_out[100]-0.0237)/0.0237 > 0.02) stat++;

  gsl_dht_free(t);

  return stat;
}

/* Test the transform
 * Integrate[ x exp(-x) J_1(a x), {x,0,Inf}] = a F(3/2, 2; 2; -a^2)
 */
int
test_dht_exp1(void)
{
  int stat = 0;
  int n;
  double f_in[128];
  double f_out[128];
  gsl_dht * t = gsl_dht_new(128, 1.0, 20.0);

  for(n=0; n<128; n++) {
    const double x = gsl_dht_x_sample(t, n);
    f_in[n] = exp(-x);
  }

  gsl_dht_apply(t, f_in, f_out);

  /* Spot check.
   * Note that the systematic errors in the calculation
   * are quite large, so it is meaningless to compare
   * to a high accuracy.
   */
  if(fabs( f_out[0]-0.181)/0.181 > 0.02) stat++;
  if(fabs( f_out[5]-0.357)/0.357 > 0.02) stat++;
  if(fabs(f_out[10]-0.211)/0.211 > 0.02) stat++;
  if(fabs(f_out[35]-0.0289)/0.0289 > 0.02) stat++;
  if(fabs(f_out[100]-0.00221)/0.00211 > 0.02) stat++;

  gsl_dht_free(t);

  return stat;
}

/* Test the transform
 * Integrate[ x^2 (1-x^2) J_1(a x), {x,0,1}] = 2/a^2 J_3(a)
 */
int
test_dht_poly1(void)
{
  int stat = 0;
  int n;
  double f_in[128];
  double f_out[128];
  gsl_dht * t = gsl_dht_new(128, 1.0, 1.0);

  for(n=0; n<128; n++) {
    const double x = gsl_dht_x_sample(t, n);
    f_in[n] = x * (1.0 - x*x);
  }

  gsl_dht_apply(t, f_in, f_out);

  /* Spot check. This function satisfies the boundary condition,
   * so the accuracy should be ok.
   */
  if(fabs( f_out[0]-0.057274214)/0.057274214    > 1.0e-07) stat++;
  if(fabs( f_out[5]-(-0.000190850))/0.000190850 > 1.0e-05) stat++;
  if(fabs(f_out[10]-0.000024342)/0.000024342    > 1.0e-04) stat++;
  if(fabs(f_out[35]-(-4.04e-07))/4.04e-07       > 1.0e-03) stat++;
  if(fabs(f_out[100]-1.0e-08)/1.0e-08           > 0.25)    stat++;

  gsl_dht_free(t);

  return stat;
}

int main()
{
  gsl_ieee_env_setup ();

  gsl_test( test_dht_exact(),   "Small Exact DHT");
  gsl_test( test_dht_simple(),  "Simple  DHT");
  gsl_test( test_dht_exp1(),    "Exp  J1 DHT");
  gsl_test( test_dht_poly1(),   "Poly J1 DHT");

  exit (gsl_test_summary());
}

max_stars_repo_name: utdsimmons/ohpc

max_issues_repo_name: utdsimmons/ohpc

max_forks_repo_name: utdsimmons/ohpc


text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <ccl.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>

//Structure containing info about a model parameter
typedef struct {
  int index;
  double value;
} LSS_param;

//Structure defining one of the tracers being correlated
//Note that each LSS tracer contains many redshift bins
//(whereas a CCL_ClTracer corresponds to a single bin)
typedef struct {
  int n_bins; //Number of redshift bins

  int *num_z; //Number of z-values for the n(z) of each bin (dims -> [nbins])
  double **z_nz_arr; //z-values for the n(z) of each bin    (dims -> [nbins,num_z])
  double **n_nz_arr; //n(z) of each bin                     (dims -> [nbins,num_z])

  int n_nodes_b; //Number of nodes defining b(z)
  double *z_bz_arr; //Redshift values for each node (dims -> [n_nodes_b])
  LSS_param *b_bz_arr; //Bias values for each node  (dims -> [n_nodes_b])
  int interp_scheme_b; //Interpolation scheme for b(z)

  //Add similar stuff for magnification bias, 2nd-order bias etc.,

  //Add similar stuff for photo-z uncertainties?

  //Flags describing theoretical calculation
  int flag_w_rsd; //Include RSDs?
  int flag_w_mag; //Include magnification?
  int non_linear_type; //Scheme for non-linearities?

  //Array of CCL_ClTracers for each redshift bin (stored here to avoid recomputation?)
  CCL_ClTracer *t; // (dims -> [n_bins])
} LSS_tracer_info;

//Structure containing info needed to compute the LSS likelihood
typedef struct {
  //Cosmological parameters
  LSS_param par_Oc;
  LSS_param par_Ob;
  LSS_param par_Ok;
  LSS_param par_h;
  LSS_param par_s8;
  LSS_param par_ns;
  LSS_param par_w0;
  LSS_param par_wa;

  //Tracers to correlate
  LSS_tracer_info *tr1;
  LSS_tracer_info *tr2;

  //Probably a lot of stuff missing here
  //...
} LSS_likelihood_workspace;

//Structure defining a set of power spectra
typedef struct {
  //Number of bins for each tracer
  int nbin_1;
  int nbin_2;

  //Number of bandpowers for each cross-power spectrum
  int **n_bpw; // (dim -> [nbin_1,nbin_2])
  //Effective ell for each bandpower
  double ***l_bpw; // (dim -> [nbin_1,nbin_2,n_bpw[nbin_1,nbin_2]])
  //Effective C_ell for each bandpower
  gsl_vector *cl_bpw; // (flattened [nbin_1,nbin_2,n_bpw[nbin_1,nbin_2]])
  //This will probably have to be more general, with a weight-per-ell
  //array for each bandpower (e.g. to account for pseudo-Cl binning).
} LSS_2pt_Cell;

//Something like the above should be defined for the precision matrix,
//or it could be read directly into a 2D matrix

//////////
// Functions defined in lss_cell_io.c
//
//Read power spectrum
LSS_2pt_Cell *lss_read_cell(char *fname);
//Write power spectrum
void lss_write_cell(char *fname,LSS_2pt_Cell *cl);
//Read precision matrix
gsl_matrix *lss_read_precision(char *fname);
//Write precision matrix
void lss_write_precision(char *fname,gsl_matrix *prec,LSS_2pt_Cell *cl);
//Destructor
void lss_free_cell(LSS_2pt_Cell *cl);

//////////
// Functions defined in lss_tracers.c
//
//Creates new tracer
LSS_tracer_info *lss_tracer_new(int n_bins,
				int *num_z,double **z_nz_arr,double **n_nz_arr,
				int n_nodes_b,double *z_bz_arr,double *b_bz_arr,
				int flag_w_rsd,int flag_w_mag,int non_linear_type);
//Destructor
void lss_tracer_free(LSS_tracer_info *tr);

//////////
// Functions defined in lss_cell_th.c
//
//Compute theory power spectrum (and return as gsl_vector)
gsl_vector *lss_cell_theory(ccl_cosmology *cosmo,LSS_likelihood_workspace *w,LSS_2pt_Cell *cl);

//////////
// Functions defined in lss_cell_like.c
//
//Computes likelihood for a set of cosmological and nuisance parameters (encoded in params)
double lss_likelihood(int n_par,double *params,LSS_2pt_Cell *data,gsl_matrix *prec,LSS_likelihood_workspace *w);
//Creates a new LSS_likelihood_workspace
LSS_likelihood_workspace *lss_like_workspace_new(ccl_cosmology *cosmo,LSS_tracer_info *tr1,LSS_tracer_info *tr2);
//Destructor
void lss_like_workspace_free(LSS_likelihood_workspace *w);

max_stars_repo_name: LSSTDESC/LSSLike

max_issues_repo_name: LSSTDESC/LSSLike

max_forks_repo_name: LSSTDESC/LSSLike


text:
//This file contains mostly untested code for computing splined
//periodic flux tube configurations

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <builtin_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv2.h>
#include "intT.h"

//extern "C"
//__device__ __host__ float interp(float rho2, float *rho2sp, float4 *coefs);

int func (double t, const double y[], double f[],
           void *params)
{
    double *p = (double *)params;
    double lambda2 = p[0];
    double a = p[1];
    f[0] = 2.0/lambda2*t*exp(-pow(sin(pi*t/a),2)/lambda2);
     return GSL_SUCCESS;
}

int jac(double t, const double y[], double *dfdy, double dfdt[], void *params)
{
    double *p = (double *)params;
    double lambda2 = p[0];
    double a = p[1];

    gsl_matrix_view dfdy_mat
        = gsl_matrix_view_array(dfdy, 1, 1);
    gsl_matrix * m = &dfdy_mat.matrix;

    double der=0.0;

    gsl_matrix_set (m, 0, 0, der);
    dfdt[0] = 2.0/lambda2*exp(-pow(sin(pi*t/a),2)/lambda2)*(1.0-
	2.0*pi/(a*lambda2)*t*sin(pi*t/a)*cos(pi/a*t));
    return GSL_SUCCESS;
}

float4* flspline(double *fli, float *rho2, int N)
//Computes the cubic spline coefficients for f_lambda(rho^2)
//uses Burden and Faires algorithm 3.5 (Clamped Cubic Spline)
//Step numbers refer to this textbook
{
  int i;
  double* h;
  double* l;
  double* mu;
  double* z;
  double* alpha;
  float4* flscoefs;
  h=(double *)malloc(N*sizeof(*h));
  l=(double *)malloc(N*sizeof(*l));
  mu=(double *)malloc(N*sizeof(*mu));
  z=(double *)malloc(N*sizeof(*z));
  alpha = (double *)malloc(N*sizeof(*alpha));
  flscoefs=(float4 *)malloc(N*sizeof(*flscoefs));

//Step 1
  for(i=0;i<N;i++)
  {
	flscoefs[i].x = fli[i];
  }
  for(i=0;i<N-1;i++) h[i]=rho2[i+1]-rho2[i];
//Step 2
  alpha[0]=3.0*(fli[1]-fli[0])/h[0];
  alpha[N-1]=-3.0*(fli[N-1]-fli[N-2])/h[N-2];
//Step 3
  for(i=1;i<N-1;i++)
  {
	//printf("%d \n",i);
	alpha[i]=3.0/h[i]*(fli[i+1]-fli[i])-3.0/h[i-1]*(fli[i]-fli[i-1]);
  }
//Step 4 -
//Steps 4-7 solve a tridiagonal system using Crout factorization (Burden and Faires alg 6.7)
  l[0]=2.0*h[0];
  mu[0]=0.5;
  z[0]=alpha[0]/l[0];
//Step 5
  for(i=1;i<N-1;i++)
  {
	l[i]=2.0*(rho2[i+1]-rho2[i-1])-h[i-1]*mu[i-1];
	mu[i]=h[i]/l[i];
	z[i]=(alpha[i]-h[i-1]*z[i-1])/l[i];
  }
//Step 6
  l[N-1]=h[N-2]*(2.0-mu[N-2]);
  z[N-1]=(alpha[N-1]-h[N-2]*z[N-2])/l[N-1];
  flscoefs[N-1].z = z[N-1];
//Step 7
  for(i=N-2;i>-1;i--)
  {
	flscoefs[i].z = z[i]-mu[i]*flscoefs[i+1].z;
	flscoefs[i].y = (fli[i+1]-fli[i])/h[i] - h[i]*(flscoefs[i+1].z+2.0*flscoefs[i].z)/3.0;
	flscoefs[i].w = (flscoefs[i+1].z-flscoefs[i].z)/(3.0*h[i]);
  }
//Free allocated memory
  free(h);free(l);free(mu);free(z);free(alpha);
//Step 8
  return flscoefs;

}

float4* getspline(float* rho2, double lambda2, double a, int Npoints)
{

    double params[2];
    const gsl_odeiv2_step_type * T = gsl_odeiv2_step_rk8pd;
    float4 *coefs;

    gsl_odeiv2_step * s  = gsl_odeiv2_step_alloc (T, 1);
    gsl_odeiv2_control * c  = gsl_odeiv2_control_y_new (1e-9,1e-7);
    gsl_odeiv2_evolve * e  = gsl_odeiv2_evolve_alloc (1);

    params[0] = lambda2;
    params[1] = a;

    gsl_odeiv2_system sys = {func, jac, 1, params};

    double rho = 0.0;
    double rhof = 20.0;
    double h = rhof/((double)Npoints);
    double y[1] = { 0.00 };
    double *flambda;
    int i=0;

    flambda=(double *)malloc(Npoints*sizeof(*flambda));
    while (i<Npoints)
    {
        int status = gsl_odeiv2_evolve_apply_fixed_step (e, c, s,
                                        &sys,
                                        &rho, h, y);
 	if(status == GSL_FAILURE)
	{
		printf("ERROR:GSL_FAILURE\n");
		break;
	}
        else if (status != GSL_SUCCESS)
            break;
        rho2[i]=rho*rho;
	flambda[i]=y[0];
	i+=1;
	//rho+=h;

    }

    gsl_odeiv2_evolve_free (e);
    gsl_odeiv2_control_free (c);
    gsl_odeiv2_step_free (s);

    coefs=flspline(flambda,rho2,Npoints);
    free(flambda);
    return coefs;
}

double getinterp(float rho2, float *rho2sp, float4 *coefs)
{
  int j;
  int upperi=Nspline-1, loweri=0;
  float rho2diff=0.0;
  double flambda;

    //Discover which interval to look in using a binary search
  if(rho2<rho2sp[Nspline-1] && rho2>rho2sp[0])
//  if(0)
  {
   	while(upperi-loweri > 1)
	{
		if(rho2 >= rho2sp[(upperi+loweri)/2]) loweri=(upperi+loweri)/2;
		else upperi=(upperi+loweri)/2;
  	}
  	//interpolate using the jth interval
  	j=loweri;
	rho2diff=rho2-rho2sp[j];
	//rho2diff=0.0;
  	flambda= coefs[j].x+rho2diff*(coefs[j].y+rho2diff*(coefs[j].z+rho2diff*coefs[j].w));
  }
  else
  {
	j=0;
	rho2diff=0.0f;
	flambda=0.0f;
  }

  //*flambda= coefs[j].x+rho2diff*(coefs[j].y+rho2diff*(coefs[j].z+rho2diff*coefs[j].w));
  //*fprime = coefs[j].y+rho2diff*(2.0f*coefs[j].z + rho2diff*(3.0f*coefs[j].w));
  //*flambda=1.0f-exp(-1.0f*rho2);
  //*fprime=1.0f*exp(-1.0f*rho2);
  return flambda;
}

double testspline(float rho2, double lambda2, double a, int Npoints, float *rho2sp, float4 *flcoefs)
{
    printf("testspline: rho2=%f\n",(double)rho2);
    double flodeiv, flspline;
    double params[2];
    const gsl_odeiv2_step_type * T = gsl_odeiv2_step_rk8pd;
    float4 *coefs;

    gsl_odeiv2_step * s  = gsl_odeiv2_step_alloc (T, 1);
    gsl_odeiv2_control * c  = gsl_odeiv2_control_y_new (1e-9,1e-7);
    gsl_odeiv2_evolve * e  = gsl_odeiv2_evolve_alloc (1);

    params[0] = lambda2;
    params[1] = a;

    gsl_odeiv2_system sys = {func, jac, 1, params};

    double rho = 0.0;
    double rhof = sqrt(rho2);
    double h = rhof/((double)Npoints);
    double y[1] = { 0.00 };
    double *flambda;
    int i=0;
    flambda=(double *)malloc(Npoints*sizeof(*flambda));

    while(rho<rhof)
    {

    	int status = gsl_odeiv2_evolve_apply (e, c, s,
                                        &sys,
                                        &rho, rhof, &h, y);
    	if(status == GSL_FAILURE)
    	{
		printf("ERROR:GSL_FAILURE\n");
    	}
    }
    flodeiv=y[0];
    //rho+=h;
    flspline=(double)interp(rho2, rho2sp, flcoefs);
    return (flodeiv-flspline)/flodeiv;
}

max_stars_repo_name: QEDan/Worldline-Flux-Tube-Lattice

max_issues_repo_name: QEDan/Worldline-Flux-Tube-Lattice

max_forks_repo_name: QEDan/Worldline-Flux-Tube-Lattice


text:
// https://www.gnu.org/software/gsl/doc/html/fft.html
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_fft_real.h>
#include <gsl/gsl_fft_halfcomplex.h>

int
main (void)
{
  int i, n = 100;
  double data[n];

  gsl_fft_real_wavetable * real;
  gsl_fft_halfcomplex_wavetable * hc;
  gsl_fft_real_workspace * work;

  for (i = 0; i < n; i++)
    {
      data[i] = 0.0;
    }

  for (i = n / 3; i < 2 * n / 3; i++)
    {
      data[i] = 1.0;
    }

  for (i = 0; i < n; i++)
    {
      printf ("%d: %e\n", i, data[i]);
    }
  printf ("\n");

  work = gsl_fft_real_workspace_alloc (n);
  real = gsl_fft_real_wavetable_alloc (n);

  gsl_fft_real_transform (data, 1, n,
                          real, work);

  gsl_fft_real_wavetable_free (real);

  for (i = 11; i < n; i++)
    {
      data[i] = 0;
    }

  hc = gsl_fft_halfcomplex_wavetable_alloc (n);

  gsl_fft_halfcomplex_inverse (data, 1, n,
                               hc, work);
  gsl_fft_halfcomplex_wavetable_free (hc);

  for (i = 0; i < n; i++)
    {
      printf ("%d: %e\n", i, data[i]);
    }

  gsl_fft_real_workspace_free (work);
  return 0;
}