text:
#ifndef NewtonRaphsonSolver_H
#define NewtonRaphsonSolver_H

#include <stdio.h>
#include <iostream>
#include <vector>
#include <gsl/gsl_linalg.h>
#include "Node.h"
#include "ShapeBase.h"
//#include <omp.h>

/**
 *  The Newton-Raphson solver class
 *  */
class NewtonRaphsonSolver{
private:
    size_t nDim;                        ///< Dimension of the space (3D)
    size_t nNodes;                      ///< Number of nodes of the system
    double threshold;                   ///< Convergence threshold for iterations

public:
    NewtonRaphsonSolver(int nDim, int nNodes); 	///< Constructer of the N-R solver
    ~NewtonRaphsonSolver();						///< Desturctor of the N-R solver

    gsl_matrix* un;								///< The initial positions of the nodes, as calculated at the end of previous step "n"
    gsl_matrix* ge;								///< The matrix containing elastic forces on each node, size (nDim*nNodes,1). Organisation is [Node0,x ; Node0,y ; Node0,z; ... ; Noden,x ; Noden,y ; Noden,z]
    gsl_matrix* gvInternal;						///< The matrix containing internal viscous forces on each node, size (nDim*nNodes,1). Organisation is [Node0,x ; Node0,y ; Node0,z; ... ; Noden,x ; Noden,y ; Noden,z]
    gsl_matrix* gvExternal;						///< The matrix containing external viscous forces on each node, size (nDim*nNodes,1). Organisation is [Node0,x ; Node0,y ; Node0,z; ... ; Noden,x ; Noden,y ; Noden,z]
    gsl_matrix* gExt;							///< The matrix containing external forces on each node (currently includes packing forces), size (nDim*nNodes,1). Organisation is [Node0,x ; Node0,y ; Node0,z; ... ; Noden,x ; Noden,y ; Noden,z]
    gsl_vector* gSum;							///< The matrix containing sum of NewtonRaphsonSolver#ge, NewtonRaphsonSolver#gvInternal, NewtonRaphsonSolver#gvExternal, NewtonRaphsonSolver#gExt. Organisation is [Node0,x ; Node0,y ; Node0,z; ... ; Noden,x ; Noden,y ; Noden,z]
    gsl_matrix* uk;								///< The matrix storing the position of each node at iteration "k". Initiated in function NewtonRaphsonSolver#initialteUkMatrix at the beginning of each step, and updated by function NewtonRaphsonSolver#updateUkInIteration during the iteartions.
    gsl_matrix* displacementPerDt;				///< The displacement per time step of each node in current iteration "k", from its position at the end of the last time step "n"
    gsl_vector* deltaU;							///< The incremental change in positions as calculated in current iteration, resulting from the imbalance of elastic, viscous and any other external forces acting on each nodes. The solver minimises this value, convergence occurs when all incremental movements for all nodes sufficiently close to zero.
    gsl_matrix* K;								///< The Jacobian matrix, derivative of sum of forces acting on each node with respect to displacements.
    bool boundNodesWithSlaveMasterDefinition;   ///< The boolean stating if there are degrees of freedom slave to other nodes (masters).
    std::vector< std::vector<int> > slaveMasterList;    ///< The 2D integer vector storing the slave-master degrees of freedom couples, such that array [i][0] = slave to array[i][1], each i representing one dim of node position ( z of node 2 is DoF i=8);

    void setMatricesToZeroAtTheBeginningOfIteration(); 		///< The function setting the calculation matrices to zero at the beginning of each iteration.
    void setMatricesToZeroInsideIteration();												///< The function setting the relevant matrices to zero at each iteration.
    void constructUnMatrix(const std::vector<std::unique_ptr<Node> > &Nodes);               ///< This function constructs NewtonRaphsonSolver#un matrix at the beginning of the iterations.
    void initialteUkMatrix();                                                               ///< This function initiates NewtonRaphsonSolver#uk matrix at the beginning of the iterations, it is initiated to be equal to NewtonRaphsonSolver#un.
    void calculateBoundKWithSlavesMasterDoF();                                              ///< This function updates the Jacobian of the system, NewtonRaphsonSolver#K, to reflect degrees of freedom binding.
    void equateSlaveDisplacementsToMasters();                                               ///< This function moves the slaves of bound couples with a displacement equivalent to the masters'.
    void calculateDisplacementMatrix(double dt);											///< This function calculates the displacement of each node in current iteration "k", from their positions at the end of the previous step "n" (NewtonRaphsonSolver#uk - NewtonRaphsonSolver#un)
    void calcutateFixedK(const std::vector <std::unique_ptr<Node>>& Nodes);					///< This function updates the Jacobian to account for nodes  that are fixed in certain dimensions in space, as part of boundary conditions.
    void calculateForcesAndJacobianMatrixNR(const std::vector <std::unique_ptr<Node>>& Nodes, const std::vector <std::unique_ptr<ShapeBase>>&  Elements, double dt );	///< This function calculates elemental forces and Jacobians, later to be combined in NewtonRaphsonSolver#K and NewtonRaphsonSolver#gSum
    void writeForcesTogeAndgvInternal(const std::vector <std::unique_ptr<Node>>& Nodes, const std::vector <std::unique_ptr<ShapeBase>>&  Elements, std::vector<std::array<double,3>>& SystemForces);	///< This function writes the values of elemental elastic (ShapeBase#ge) and internal viscous forces (ShapeBase#gvInternal) into the system elastic and internal viscous forces, NewtonRaphsonSolver#ge, and NewtonRaphsonSolver#gvInternal, respectively.
    void writeImplicitElementalKToJacobian(const std::vector <std::unique_ptr<ShapeBase>>&  Elements);	///< This function writes the elemental values for elastic part of the Jacobian - stiffness matrix - (ShapeBase#TriPointKe) and for viscous part of Jacobian (ShapeBase#TriPointKv) into the system Jacobian NewtonRaphsonSolver#K.
    void calculateExternalViscousForcesForNR(const std::vector <std::unique_ptr<Node>>& Nodes);         ///< This function calculates the external viscous forces acting on each node, the values are sotred in NewtonRaphsonSolver#gvExternal
    void addImplicitKViscousExternalToJacobian(const std::vector <std::unique_ptr<Node>>& Nodes, double dt);	///< This function adds the external related terms of the Jacobian to the system Jacobian NewtonRaphsonSolver#K.
    void checkJacobianForAblatedNodes(std::vector <int> & AblatedNodes);                    ///< This functions checks the Jacobian to ensure the diagonal terms are non-zero for ablated nodes.
    void calculateSumOfInternalForces();                                                    ///< This function adds the ealsticity and viscosity related forces (NewtonRaphsonSolver#ge, NewtonRaphsonSolver#gvInternal, NewtonRaphsonSolver#gvExternal) to sum of forces, NewtonRaphsonSolver#gSum.
	void addExernalForces();

    void solveForDeltaU();                                                                  ///< This function solves for the displacements within the N-R step.
    //raw pointers necessary for Pardiso
	int  solveWithPardiso(double* a, double*b, int* ia, int* ja, const int n_variables);
    void constructiaForPardiso(int* ia, const int nmult, std::vector<int> &ja_vec, std::vector<double> &a_vec);
    void writeKinPardisoFormat(const int nNonzero, std::vector<int> &ja_vec, std::vector<double> &a_vec, int* ja, double* a);
	void writeginPardisoFormat(double* b, const int n);
    bool checkConvergenceViaDeltaU();                                                      ///< Check for cenvergence with the norm of displacements vector,against the threshold NewtonRaphsonSolver#threshold.
    bool checkConvergenceViaForce();                                                       ///< Check for cenvergence with the norm of forces vector,against the threshold NewtonRaphsonSolver#threshold.
    void updateUkInIteration();                                                            ///< Calulate the nodal displacemetns at the kth iteration of NR dolver.

    void displayMatrix(gsl_matrix* mat, std::string matname);
    void displayMatrix(gsl_vector* mat, std::string matname);

	bool checkIfCombinationExists(int dofSlave, int dofMaster);
    void checkMasterUpdate(int& dofMaster, int& masterId); ///< This function takes a degree of freedom number as input. This DOF is supposed to be a master. If, the dof is already a slave to another dof, then update the masetr dof. Anything that would be bound to the input dof can be bound to the already existing master of the input dof.
    bool checkIfSlaveIsAlreadyMasterOfOthers(int dofSlave, int dofMaster); ///< This function checks if the slave DOF is already master of others, if so, updates the master of said slave to the new master the current slave will be bound to.

    void updateElementPositions(const std::vector<std::unique_ptr<Node> > &Nodes, const std::vector <std::unique_ptr<ShapeBase>>&  Elements);

};

#endif


text:
/* statistics/test_robust.c
 *
 * 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.
 */

#include <config.h>
#include <stdlib.h>
#include <math.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_ieee_utils.h>

int test_robust (void);

/* random vector in [-1,1] */
static int
random_array(const size_t n, double * x, gsl_rng * r)
{
  size_t i;

  for (i = 0; i < n; ++i)
    x[i] = 2.0 * gsl_rng_uniform(r) - 1.0;

  return 0;
}

/* calculate MAD statistic for input vector using slow/naive algorithm */
static double
slow_MAD(const size_t n, const double x[])
{
  double *work = malloc(n * sizeof(double));
  double median, mad;
  size_t i;

  for (i = 0; i < n; ++i)
    work[i] = x[i];

  gsl_sort(work, 1, n);
  median = gsl_stats_median_from_sorted_data(work, 1, n);

  for (i = 0; i < n; ++i)
    work[i] = fabs(x[i] - median);

  gsl_sort(work, 1, n);
  mad = gsl_stats_median_from_sorted_data(work, 1, n);

  free(work);

  return mad;
}

/* calculate S_n statistic for input vector using slow/naive algorithm */
static double
slow_Sn0(const size_t n, const double x[])
{
  double *work1 = malloc(n * sizeof(double));
  double *work2 = malloc(n * sizeof(double));
  double Sn;
  size_t i, j;

  for (i = 0; i < n; ++i)
    {
      for (j = 0; j < n; ++j)
        work1[j] = fabs(x[i] - x[j]);

      /* find himed_j | x_i - x_j | */
      gsl_sort(work1, 1, n);
      work2[i] = work1[n / 2];
    }

  /* find lomed_i { himed_j | x_i - x_j | } */
  gsl_sort(work2, 1, n);
  Sn = work2[(n + 1) / 2 - 1];

  free(work1);
  free(work2);

  return Sn;
}

/* calculate Q_n statistic for input vector using slow/naive algorithm */
static double
slow_Qn0(const size_t n, const double x[])
{
  const size_t wsize = n * (n - 1) / 2;
  const size_t n_2 = n / 2;
  const size_t k = ((n_2 + 1) * n_2) / 2;
  double *work;
  double Qn;
  size_t idx = 0;
  size_t i, j;

  if (n < 2)
    return (0.0);

  work = malloc(wsize * sizeof(double));

  for (i = 0; i < n; ++i)
    {
      for (j = i + 1; j < n; ++j)
        work[idx++] = fabs(x[i] - x[j]);
    }

  gsl_sort(work, 1, idx);
  Qn = work[k - 1];

  free(work);

  return Qn;
}

static int
test_median(const double tol, const size_t n, gsl_rng * r)
{
  double * x = malloc(n * sizeof(double));
  double median1, median2;

  random_array(n, x, r);

  median1 = gsl_stats_median(x, 1, n);

  gsl_sort(x, 1, n);
  median2 = gsl_stats_median_from_sorted_data(x, 1, n);

  gsl_test_rel(median1, median2, tol, "test_median n=%zu", n);

  free(x);

  return 0;
}

static int
test_mad(const double tol, const size_t n, gsl_rng * r)
{
  double * x = malloc(n * sizeof(double));
  double * work = malloc(n * sizeof(double));
  double mad1, mad2;

  random_array(n, x, r);

  mad1 = slow_MAD(n, x);

  gsl_sort(x, 1, n);
  mad2 = gsl_stats_mad0(x, 1, n, work);

  gsl_test_rel(mad1, mad2, tol, "test_mad n=%zu", n);

  free(x);
  free(work);

  return 0;
}

static int
test_Sn(const double tol, const size_t n, gsl_rng * r)
{
  double * x = malloc(n * sizeof(double));
  double * work = malloc(n * sizeof(double));
  double Sn1, Sn2;

  random_array(n, x, r);

  /* compute S_n with slow/naive algorithm */
  Sn1 = slow_Sn0(n, x);

  /* compute S_n with efficient algorithm */
  gsl_sort(x, 1, n);
  Sn2 = gsl_stats_Sn0_from_sorted_data(x, 1, n, work);

  gsl_test_rel(Sn2, Sn1, tol, "test_Sn n=%zu", n);

  free(x);
  free(work);

  return 0;
}

static int
test_Qn(const double tol, const size_t n, gsl_rng * r)
{
  double * x = malloc(n * sizeof(double));
  double * work = malloc(3 * n * sizeof(double));
  int * work_int = malloc(5 * n * sizeof(int));
  double Qn1, Qn2;

  random_array(n, x, r);

  /* compute Q_n with slow/naive algorithm */
  Qn1 = slow_Qn0(n, x);

  /* compute Q_n with efficient algorithm */
  gsl_sort(x, 1, n);
  Qn2 = gsl_stats_Qn0_from_sorted_data(x, 1, n, work, work_int);

  gsl_test_rel(Qn2, Qn1, tol, "test_Qn n=%zu", n);

  free(x);
  free(work);
  free(work_int);

  return 0;
}

int
test_robust (void)
{
  const double tol = 1.0e-12;
  gsl_rng * r = gsl_rng_alloc(gsl_rng_default);

  test_median(GSL_DBL_EPSILON, 1, r);
  test_median(GSL_DBL_EPSILON, 2, r);
  test_median(GSL_DBL_EPSILON, 3, r);
  test_median(GSL_DBL_EPSILON, 100, r);
  test_median(GSL_DBL_EPSILON, 101, r);
  test_median(GSL_DBL_EPSILON, 500, r);
  test_median(GSL_DBL_EPSILON, 501, r);

  test_mad(GSL_DBL_EPSILON, 1, r);
  test_mad(GSL_DBL_EPSILON, 2, r);
  test_mad(GSL_DBL_EPSILON, 3, r);
  test_mad(GSL_DBL_EPSILON, 100, r);
  test_mad(GSL_DBL_EPSILON, 101, r);
  test_mad(GSL_DBL_EPSILON, 500, r);
  test_mad(GSL_DBL_EPSILON, 501, r);

  test_Sn(tol, 1, r);
  test_Sn(tol, 2, r);
  test_Sn(tol, 3, r);
  test_Sn(tol, 100, r);
  test_Sn(tol, 101, r);
  test_Sn(tol, 500, r);
  test_Sn(tol, 501, r);

  test_Qn(tol, 1, r);
  test_Qn(tol, 2, r);
  test_Qn(tol, 3, r);
  test_Qn(tol, 4, r);
  test_Qn(tol, 5, r);
  test_Qn(tol, 100, r);
  test_Qn(tol, 101, r);
  test_Qn(tol, 500, r);
  test_Qn(tol, 501, r);

  gsl_rng_free(r);

  return 0;
}

ext: c

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
#pragma once

#include <vector>

#include <gsl/span>

#include "Mouse.h"
#include "Joystick.h"
#include "JoystickButton.h"
#include "Keyboard.h"
#include "RawInput.h"
#include "Xbox360Controller.h"

namespace qvr
{

class JoystickAxisThreshold
{
public:
	JoystickAxisThreshold(
		const qvr::JoystickAxis axis,
		const float threshold,
		const bool trueAboveThreshold)
		: mJoystickAxis(axis)
		, mThreshold(threshold)
		, mTrueAboveThreshold(trueAboveThreshold)
	{}
	JoystickAxisThreshold(
		const Xbox360Controller::Axis axis,
		const float threshold,
		const bool trueAboveThreshold)
		: JoystickAxisThreshold(
			Xbox360Controller::ToJoystickAxis(axis),
			threshold,
			trueAboveThreshold)
	{
	}
	bool IsActive(const qvr::Joystick& joystick) const
	{
		const float v = joystick.GetPosition(mJoystickAxis);
		if (mTrueAboveThreshold)
		{
			return v > mThreshold;
		}
		return v < mThreshold;
	}
	bool JustActive(const qvr::Joystick& joystick) const
	{
		const float current  = joystick.GetPosition(mJoystickAxis);
		const float previous = joystick.GetPreviousPosition(mJoystickAxis);
		if (mTrueAboveThreshold)
		{
			return current > mThreshold && previous < mThreshold;
		}
		return current < mThreshold && previous > mThreshold;
	}
private:
	const qvr::JoystickAxis mJoystickAxis;
	const float mThreshold;
	const bool mTrueAboveThreshold;
};

}

class BinaryInput
{
public:
	enum class Type
	{
		KeyboardKey,
		MouseButton,
		JoystickButton,
		JoystickAxisThreshold
	};

	BinaryInput(qvr::KeyboardKey key)
		: mType(Type::KeyboardKey)
		, mKeyboardKey(key)
	{}

	BinaryInput(qvr::MouseButton mouseButton)
		: mType(Type::MouseButton)
		, mMouseButton(mouseButton)
	{}

	BinaryInput(qvr::JoystickButton joystickButton)
		: mType(Type::JoystickButton)
		, mJoystickButton(joystickButton)
	{}

	BinaryInput(qvr::JoystickAxisThreshold joystickAxisThreshold)
		: mType(Type::JoystickAxisThreshold)
		, mJoystickAxisThreshold(joystickAxisThreshold)
	{}

	bool IsActive(const qvr::RawInputDevices& devices) const
	{
		const auto joystickIndex = qvr::JoystickIndex(0);
		const qvr::Joystick* joystick = devices.GetJoysticks().GetJoystick(joystickIndex);

		switch (mType)
		{
		case Type::KeyboardKey:
			return devices.GetKeyboard().IsDown(mKeyboardKey);
		case Type::MouseButton:
			return devices.GetMouse().IsDown(mMouseButton);
		// TODO: Which joysticks?
		case Type::JoystickButton:
			if (joystick) {
				return joystick->IsDown(mJoystickButton);
			}
			break;
		case Type::JoystickAxisThreshold:
			if (joystick) {
				return mJoystickAxisThreshold.IsActive(*joystick);
			}
			break;
		}
		return false;
	}

	bool JustActive(const qvr::RawInputDevices& devices) const
	{
		const auto joystickIndex = qvr::JoystickIndex(0);
		const qvr::Joystick* joystick = devices.GetJoysticks().GetJoystick(joystickIndex);

		switch (mType)
		{
		case Type::KeyboardKey:
			return devices.GetKeyboard().JustDown(mKeyboardKey);
		case Type::MouseButton:
			return devices.GetMouse().JustDown(mMouseButton);
			// TODO: Which joysticks?
		case Type::JoystickButton:
			if (joystick) {
				return joystick->JustDown(mJoystickButton);
			}
			break;
		case Type::JoystickAxisThreshold:
			if (joystick) {
				return mJoystickAxisThreshold.JustActive(*joystick);
			}
			break;
		}
		return false;
	}

private:
	const Type mType;

	union
	{
		const qvr::KeyboardKey mKeyboardKey;
		const qvr::MouseButton mMouseButton;
		const qvr::JoystickButton mJoystickButton;
		const qvr::JoystickAxisThreshold mJoystickAxisThreshold;
	};
};

inline bool AnyActive(const qvr::RawInputDevices& devices, const gsl::span<const BinaryInput> inputs) {
	for (const BinaryInput input : inputs) {
		if (input.IsActive(devices)) {
			return true;
		}
	}

	return false;
}

inline bool AnyJustActive(const qvr::RawInputDevices& devices, const gsl::span<const BinaryInput> inputs) {
	for (const BinaryInput input : inputs) {
		if (input.JustActive(devices)) {
			return true;
		}
	}

	return false;
}

ext: h

max_stars_repo_name: rachelnertia/Quarrel

max_issues_repo_name: rachelnertia/Quarrel

max_forks_repo_name: rachelnertia/Quiver


text:
/**
 * APER_CHECK - create a FITS image showing the locations of
 *               apertures.
 */
#include "aper_check.h"
#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_vector.h>
#include "spce_output.h"

void mark_trace_in_aXe_mask(ap_pixel * ap_p, aXe_mask *mask)
{
  gsl_vector_int *trace_inds;
  int i, ind;

  trace_inds = (gsl_vector_int *) get_trace_inds (ap_p);

  for (i = 0; i < (int)trace_inds->size; i++)
    {
      ind = gsl_vector_int_get (trace_inds, i);
      fprintf(stderr,"%d %d\n",ap_p[ind].p_x,ap_p[ind].p_y);
    }

  gsl_vector_int_free(trace_inds);
}

/**
   Function to create and return a gsl array of the same size as the grism
   part of n observation structure.

   @param ob a pointer to an observation structure
   @return a pointer to a newly allocated aXe_mask structure

*/
aXe_mask *
aXe_mask_init (observation * ob)
{
  aXe_mask *mask;
  long i,j;
  double v;

     /* Allocated space for the mask array */
  mask = (aXe_mask *) malloc (sizeof (aXe_mask));
  if (!mask)
    {
      aXe_message (aXe_M_FATAL, __FILE__, __LINE__, "Out of memory.");
    }

  /* Create the GSL image array */
  mask->img = gsl_matrix_alloc (ob->grism->size1, ob->grism->size2);
  //  fprintf(stderr,"%ld %ld\n",ob->grism->size1,ob->grism->size2);

  /* Initialise the mask array */
  for (i=0;i<(int)ob->grism->size1;i++) {
    for (j=0;j<(int)ob->grism->size2;j++) {
      v = gsl_matrix_get(ob->grism,i,j);
      gsl_matrix_set(mask->img, i,j, v);

    }
  }

  //gsl_matrix_set_all (mask->img, 1.0);

  return mask;
}

/**
 * This function sets all the pixels listed in an ap_pixel list
 * in an aXe_mask to 0.0

    @param a pointer to a ap_pixel array
    @param mask an existing aXe_mask structure

*/
void
add_ap_p_to_aXe_mask (ap_pixel * ap_p, aXe_mask * mask)
{
  ap_pixel *cur_p;

  if (ap_p==NULL) return;

  for (cur_p = ap_p; cur_p->p_x != -1; cur_p++)
    {

      /* Loop over pixels which are valid */
      if (!isnan (cur_p->count))
        {
          //   fprintf(stderr,"%d %d\n",cur_p->p_x, cur_p->p_y);

          gsl_matrix_set (mask->img, cur_p->p_x, cur_p->p_y, 0.0);
        }
    }
}

max_stars_repo_name: sosey/pyaxe

max_issues_repo_name: sosey/pyaxe

max_forks_repo_name: sosey/pyaxe


text:
/* bspline/gsl_bspline.h
 *
 * Copyright (C) 2006 Patrick Alken
 * Copyright (C) 2008 Rhys Ulerich
 *
 * 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_BSPLINE_H__
#define __GSL_BSPLINE_H__

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

#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.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 k;      /* spline order */
    size_t km1;    /* k - 1 (polynomial order) */
    size_t l;      /* number of polynomial pieces on interval */
    size_t nbreak; /* number of breakpoints (l + 1) */
    size_t n;      /* number of bspline basis functions (l + k - 1) */

    gsl_vector *knots;  /* knots vector */
    gsl_vector *deltal; /* left delta */
    gsl_vector *deltar; /* right delta */
    gsl_vector *B;      /* temporary spline results */

    /* bspline derivative parameters */
    gsl_matrix *A;      /* work matrix */
    gsl_matrix *dB;     /* temporary derivative results */
} gsl_bspline_workspace;

GSL_FUN gsl_bspline_workspace *
gsl_bspline_alloc(const size_t k, const size_t nbreak);

GSL_FUN void gsl_bspline_free(gsl_bspline_workspace *w);

GSL_FUN size_t gsl_bspline_ncoeffs(gsl_bspline_workspace * w);
GSL_FUN size_t gsl_bspline_order(gsl_bspline_workspace * w);
GSL_FUN size_t gsl_bspline_nbreak(gsl_bspline_workspace * w);
GSL_FUN double gsl_bspline_breakpoint(size_t i, gsl_bspline_workspace * w);
GSL_FUN double gsl_bspline_greville_abscissa(size_t i, gsl_bspline_workspace *w);

GSL_FUN int
gsl_bspline_knots(const gsl_vector *breakpts, gsl_bspline_workspace *w);

GSL_FUN int gsl_bspline_knots_uniform(const double a, const double b,
                              gsl_bspline_workspace *w);

GSL_FUN int
gsl_bspline_knots_greville(const gsl_vector *abscissae,
                           gsl_bspline_workspace *w,
                           double *abserr);

GSL_FUN int
gsl_bspline_eval(const double x, gsl_vector *B,
                 gsl_bspline_workspace *w);

GSL_FUN int
gsl_bspline_eval_nonzero(const double x,
                         gsl_vector *Bk,
                         size_t *istart,
                         size_t *iend,
                         gsl_bspline_workspace *w);

GSL_FUN int
gsl_bspline_deriv_eval(const double x,
                       const size_t nderiv,
                       gsl_matrix *dB,
                       gsl_bspline_workspace *w);

GSL_FUN int
gsl_bspline_deriv_eval_nonzero(const double x,
                               const size_t nderiv,
                               gsl_matrix *dB,
                               size_t *istart,
                               size_t *iend,
                               gsl_bspline_workspace *w);

__END_DECLS

#endif /* __GSL_BSPLINE_H__ */

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:
/*
  NAME:
     proj_gauss_mixtures
  PURPOSE:
     runs the full projected gaussian mixtures algorithms, with
     regularization and split and merge
  CALLING SEQUENCE:
     proj_gauss_mixtures(struct datapoint * data, int N,
     struct gaussian * gaussians, int K,bool * fixamp, bool * fixmean,
     bool * fixcovar, double * avgloglikedata, double tol,
     long long int maxiter, bool likeonly, double w, int splitnmerge,
     bool keeplog, FILE *logfile, FILE *convlogfile, bool noproj,
     bool diagerrs,noweight)
  INPUT:
     data        - the data
     N           - number of datapoints
     gaussians   - model gaussians (initial conditions)
     K           - number of model gaussians
     fixamp      - fix the amplitude?
     fixmean     - fix the mean?
     fixcovar    - fix the covar?
     tol         - proj_EM convergence limit
     maxiter     - maximum number of iterations in each proj_EM
     likeonly    - only compute the likelihood?
     w           - regularization paramter
     splitnmerge - split 'n' merge depth (how far down the list to go)
     keeplog     - keep a log in a logfile?
     logfile     - pointer to the logfile
     convlogfile - pointer to the convlogfile
     noproj      - don't perform any projections
     diagerrs    - the data->SS errors-squared are diagonal
     noweight    - don't use data-weights
  OUTPUT:
     updated model gaussians
     avgloglikedata - average log likelihood of the data
  REVISION HISTORY:
     2008-08-21 Written Bovy
     2010-03-01 Added noproj option - Bovy
     2010-04-01 Added noweight option - Bovy
*/
#ifdef _OPENMP
#include <omp.h>
#endif
#include <stdio.h>
#include <math.h>
#include <stdbool.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include "proj_gauss_mixtures.h"

void proj_gauss_mixtures(struct datapoint * data, int N,
			 struct gaussian * gaussians, int K,
			 bool * fixamp, bool * fixmean, bool * fixcovar,
			 double * avgloglikedata, double tol,
			 long long int maxiter, bool likeonly, double w,
			 int splitnmerge, bool keeplog, FILE *logfile,
			 FILE *convlogfile, bool noproj, bool diagerrs,
			 bool noweight){
  //Allocate some memory
  struct gaussian * startgaussians;
  startgaussians = gaussians;
  int d = (gaussians->VV)->size1;//dim of mm
  //Only give copies of the fix* vectors to the EM algorithm
  bool *fixamp_tmp, *fixmean_tmp,*fixcovar_tmp;
  fixamp_tmp = (bool *) malloc(K * sizeof (bool) );
  fixmean_tmp = (bool *) malloc(K * sizeof (bool) );
  fixcovar_tmp = (bool *) malloc(K * sizeof (bool) );
  int kk;
  for (kk = 0; kk != K; ++kk){
    *(fixamp_tmp++) = *(fixamp++);
    *(fixmean_tmp++) = *(fixmean++);
    *(fixcovar_tmp++) = *(fixcovar++);
  }
  fixamp -= K;
  fixamp_tmp -= K;
  fixmean -= K;
  fixmean_tmp -= K;
  fixcovar -= K;
  fixcovar_tmp -= K;
  //allocate the newalpha, newmm and newVV matrices
#ifdef _OPENMP
  nthreads = omp_get_max_threads();
#else
  nthreads = 1;
#endif
  newgaussians = (struct gaussian *) malloc(K * nthreads * sizeof (struct gaussian) );
  startnewgaussians = newgaussians;
  int ll;
  for (kk=0; kk != K*nthreads; ++kk){
    newgaussians->alpha = 0.0;
    newgaussians->mm = gsl_vector_calloc (d);
    newgaussians->VV = gsl_matrix_calloc (d,d);
    ++newgaussians;
  }
  newgaussians= startnewgaussians;
  double oldavgloglikedata;
  //allocate the q_ij matrix
  qij = gsl_matrix_alloc(N,K);
  double * qstarij = (double *) malloc(N * sizeof (double) );
  I = gsl_matrix_alloc(d,d);
  gsl_matrix_set_identity(I);//Unit matrix
  gsl_matrix_scale(I,w);//scaled to w
  //Also take care of the bbij's and the BBij's
  bs = (struct modelbs *) malloc(nthreads * K * sizeof (struct modelbs) );
  for (kk = 0; kk != nthreads*K; ++kk){
    bs->bbij = gsl_vector_alloc (d);
    bs->BBij = gsl_matrix_alloc (d,d);
    ++bs;
  }
  bs -= nthreads*K;
  //splitnmerge
  int maxsnm = K*(K-1)*(K-2)/2;
  int * snmhierarchy = (int *) malloc(maxsnm*3* sizeof (int) );
  int j,k,l;
  struct gaussian * oldgaussians = (struct gaussian *) malloc(K * sizeof (struct gaussian) );
  gsl_matrix * oldqij = gsl_matrix_alloc(N,K);
  for (kk=0; kk != K; ++kk){
    oldgaussians->mm = gsl_vector_calloc (d);
    oldgaussians->VV = gsl_matrix_calloc (d,d);
    ++oldgaussians;
  }
  oldgaussians -= K;

  //create temporary file to hold convergence info
  FILE *tmpconvfile;
  char * tmpfilename= tmpnam(NULL);;
  if (keeplog)
    tmpconvfile= fopen(tmpfilename,"w+");

  //proj_EM
  if (keeplog)
    fprintf(logfile,"#Initial proj_EM\n");
  //printf("Where's the segmentation fault?\n");
  proj_EM(data,N,gaussians,K,fixamp_tmp,fixmean_tmp,fixcovar_tmp,
	  avgloglikedata,tol,maxiter,likeonly,w,
	  keeplog,logfile,convlogfile,noproj,diagerrs,noweight);
  if (keeplog){
    fprintf(logfile,"\n");
    fprintf(convlogfile,"\n");
  }
  //reset fix* vectors
  for (kk = 0; kk != K; ++kk){
    *(fixamp_tmp++) = *(fixamp++);
    *(fixmean_tmp++) = *(fixmean++);
    *(fixcovar_tmp++) = *(fixcovar++);
  }
  fixamp -= K;
  fixamp_tmp -= K;
  fixmean -= K;
  fixmean_tmp -= K;
  fixcovar -= K;
  fixcovar_tmp -= K;

  //Run splitnmerge
  randgen = gsl_rng_alloc(gsl_rng_mt19937);
  bool weretrying = true;
  if (likeonly || splitnmerge == 0 || K < 3)
    ;
  else {
    while (weretrying){
      weretrying = false; /* this is set back to true if an improvement is found */
      //store avgloglike from normal EM and model parameters
      oldavgloglikedata = *avgloglikedata;
      gsl_matrix_memcpy(oldqij,qij);
      for (kk=0; kk != K; ++kk){
	oldgaussians->alpha = gaussians->alpha;
	gsl_vector_memcpy(oldgaussians->mm,gaussians->mm);
	gsl_matrix_memcpy(oldgaussians->VV,gaussians->VV);
	++oldgaussians;
	++gaussians;
      }
      gaussians -= K;
      oldgaussians -= K;
      //Then calculate the splitnmerge hierarchy
      calc_splitnmerge(data,N,gaussians,K,qij,snmhierarchy);
      //Then go through this hierarchy
      kk=0;
      while (kk != splitnmerge && kk != maxsnm){
	//splitnmerge
	j = *(snmhierarchy++);
	k = *(snmhierarchy++);
	l = *(snmhierarchy++);
	splitnmergegauss(gaussians,K,oldqij,j,k,l);
	//partial EM
	//Prepare fixed vectors for partial EM
	for (ll = 0; ll != K; ++ll){
	  *(fixamp_tmp++) = (ll != j && ll != k && ll != l);
	  *(fixmean_tmp++) = (ll != j && ll != k && ll != l);
	  *(fixcovar_tmp++) = (ll != j && ll != k && ll != l);
	}
	fixamp_tmp -= K;
	fixmean_tmp -= K;
	fixcovar_tmp -= K;
	if (keeplog)
	  fprintf(logfile,"#Merging %i and %i, splitting %i\n",j,k,l);
	proj_EM(data,N,gaussians,K,fixamp_tmp,fixmean_tmp,fixcovar_tmp,
		avgloglikedata,tol,maxiter,likeonly,w,keeplog,logfile,
		tmpconvfile,noproj,diagerrs,noweight);
	//reset fix* vectors
	for (ll = 0; ll != K; ++ll){
	    *(fixamp_tmp++) = *(fixamp++);
	    *(fixmean_tmp++) = *(fixmean++);
	    *(fixcovar_tmp++) = *(fixcovar++);
	}
	fixamp -= K;
	fixamp_tmp -= K;
	fixmean -= K;
	fixmean_tmp -= K;
	fixcovar -= K;
	fixcovar_tmp -= K;
	//Full EM
	if (keeplog){
	  fprintf(logfile,"#full EM:\n");
	  fprintf(tmpconvfile,"\n");
	}
	proj_EM(data,N,gaussians,K,fixamp_tmp,fixmean_tmp,fixcovar_tmp,
		avgloglikedata,tol,maxiter,likeonly,w,keeplog,logfile,
		tmpconvfile,noproj,diagerrs,noweight);
	if (keeplog){
	  fprintf(logfile,"\n");
	  fprintf(tmpconvfile,"\n");
	}
	//reset fix* vectors
	for (ll = 0; ll != K; ++ll){
	    *(fixamp_tmp++) = *(fixamp++);
	    *(fixmean_tmp++) = *(fixmean++);
	    *(fixcovar_tmp++) = *(fixcovar++);
	  }
	  fixamp -= K;
	  fixamp_tmp -= K;
	  fixmean -= K;
	  fixmean_tmp -= K;
	  fixcovar -= K;
	  fixcovar_tmp -= K;
	//Better?
	if (*avgloglikedata > oldavgloglikedata){
	  if (keeplog){
	    fprintf(logfile,"#accepted\n");
	    //Copy tmpfile into convfile
	    fseek(tmpconvfile,0,SEEK_SET);
	    while (feof(tmpconvfile) == 0)
	      fputc(fgetc(tmpconvfile),convlogfile);
	    fseek(convlogfile,-1,SEEK_CUR);
	    fclose(tmpconvfile);
	    tmpconvfile= fopen(tmpfilename,"w+");
	  }
	  weretrying = true;
	  ++kk;
	  break;
	}
	else {
	  if (keeplog){
	    fprintf(logfile,"#didn't improve likelihood\n");
	    fclose(tmpconvfile);
	    tmpconvfile= fopen(tmpfilename,"w+");
	  }
	  //revert back to the older solution
	  *avgloglikedata = oldavgloglikedata;
	  /*  gsl_matrix_memcpy(qij,oldqij);*/
	  for (ll=0; ll != K; ++ll){
	    gaussians->alpha = oldgaussians->alpha;
	    gsl_vector_memcpy(gaussians->mm,oldgaussians->mm);
	    gsl_matrix_memcpy(gaussians->VV,oldgaussians->VV);
	    ++oldgaussians;
	    ++gaussians;
	  }
	  gaussians -= K;
	  oldgaussians -= K;
	}
	++kk;
      }
      snmhierarchy -= 3*kk;
    }
  }

  if (keeplog){
    fclose(tmpconvfile);
    fprintf(convlogfile,"\n");
  }

  //Compute some criteria to set the number of Gaussians and print these to the logfile
  int ii;
  int npc,np;
  double pc,aic,mdl;
  if (keeplog){
    //Partition coefficient
    pc=0.;
    for (ii=0; ii != N; ++ii)
      for (kk=0; kk != K; ++kk)
	pc += pow(exp(gsl_matrix_get(qij,ii,kk)),2);
    pc /= N;
    fprintf(logfile,"Partition coefficient \t=\t%f\n",pc);
    //Akaike's information criterion
    npc = 1 + d + d * (d - 1) / 2;
    np = K * npc;
    aic= -2 * ((double) N - 1. - npc - 100)* *avgloglikedata + 3 * np;
    fprintf(logfile,"AIC \t\t\t=\t%f\n",aic);
    //MDL
    mdl = - *avgloglikedata * N + 0.5 * np * log(N);
    fprintf(logfile,"MDL \t\t\t=\t%f\n",mdl);
  }

  //Free memory
  gsl_matrix_free(I);
  gsl_matrix_free(qij);
  for (kk = 0; kk != nthreads*K; ++kk){
    gsl_vector_free(bs->bbij);
    gsl_matrix_free(bs->BBij);
    ++bs;
  }
  bs -= nthreads*K;
  free(bs);
  for (kk=0; kk != K*nthreads; ++kk){
    gsl_vector_free(newgaussians->mm);
    gsl_matrix_free(newgaussians->VV);
    ++newgaussians;
  }
  newgaussians= startnewgaussians;
  free(newgaussians);
  gsl_rng_free(randgen);

  gsl_matrix_free(oldqij);
  for (kk=0; kk != K; ++kk){
    gsl_vector_free(oldgaussians->mm);
    gsl_matrix_free(oldgaussians->VV);
    ++oldgaussians;
  }
  oldgaussians -= K;
  free(oldgaussians);
  free(qstarij);
  free(snmhierarchy);
  free(fixamp_tmp);
  free(fixmean_tmp);
  free(fixcovar_tmp);

  return;
}

max_stars_repo_name: surbut/mashr

max_issues_repo_name: surbut/mashr

max_forks_repo_name: surbut/mashr


text:
#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/matrix/gsl_matrix.h>

#define BASE_DOUBLE
#include <gsl/templates_on.h>
#include <gsl/matrix/minmax_source.c>
#include <gsl/templates_off.h>
#undef BASE_DOUBLE

max_stars_repo_name: MontyThibault/centre-of-mass-awareness

max_issues_repo_name: MontyThibault/centre-of-mass-awareness

max_forks_repo_name: MontyThibault/centre-of-mass-awareness


text:
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <limits.h>
#include <sys/time.h>
#include <inttypes.h>

#include <omp.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_heapsort.h>

#include "MZC3D64.h"
#include "MZC2D64.h"
#include "sph_data_types.h"
#include "sph_linked_list.h"

int safe_free_box(linkedListBox *box){
	kh_destroy(0, box->hbegin);
	kh_destroy(1, box->hend);
	free(box);

	return 0;
}

int compare_int64_t(const void *p,const void *q){
	int64_t *data1,*data2;
	data1 = (int64_t*)p;
	data2 = (int64_t*)q;

	if(data1[0] < data2[0])
		return -1;
	else if(data1[0] == data2[0])
		return 0;
	else
		return 1;
}

#define safe_check_alloc(ptr,N,dtype) {                                          \
                                       (ptr) = (dtype*)malloc((N)*sizeof(dtype));\
                                       if((ptr)==NULL){                          \
                                         success=0;                              \
                                         goto finishlabel;                       \
                                       }                                         \
                                      }

#define safe_check__aligned_alloc(ptr,alignment,N,dtype) {                                        \
                                       (ptr) = (dtype*)aligned_alloc(alignment,(N)*sizeof(dtype));\
                                       if((ptr)==NULL){                                           \
                                         success=0;                                               \
                                         goto finishlabel;                                        \
                                       }                                                          \
                                      }

#define safe_free(ptr) {               \
                        if(ptr != NULL)\
                          free(ptr);   \
                       }

int SPHparticle_SoA_malloc(int N,SPHparticle **lsph){
  int success=1;

  (*lsph) = (SPHparticle*)malloc(1*sizeof(SPHparticle));
  if(lsph==NULL){
    success = 0;
    goto finishlabel;
  }
  (*lsph)->x  = NULL; (*lsph)->y  = NULL; (*lsph)->z  = NULL;
  (*lsph)->ux = NULL; (*lsph)->uy = NULL; (*lsph)->uz = NULL;
  (*lsph)->Fx = NULL; (*lsph)->Fy = NULL; (*lsph)->Fz = NULL;
  (*lsph)->nu = NULL; (*lsph)->rho= NULL;
  (*lsph)->id = NULL; (*lsph)->hash= NULL;

  const int alignment = 32;
  safe_check__aligned_alloc((*lsph)->x   , alignment, N ,double);
  safe_check__aligned_alloc((*lsph)->y   , alignment, N ,double);
  safe_check__aligned_alloc((*lsph)->z   , alignment, N ,double);
  safe_check__aligned_alloc((*lsph)->ux  , alignment, N ,double);
  safe_check__aligned_alloc((*lsph)->uy  , alignment, N ,double);
  safe_check__aligned_alloc((*lsph)->uz  , alignment, N ,double);
  safe_check__aligned_alloc((*lsph)->Fx  , alignment, N ,double);
  safe_check__aligned_alloc((*lsph)->Fy  , alignment, N ,double);
  safe_check__aligned_alloc((*lsph)->Fz  , alignment, N ,double);
  safe_check__aligned_alloc((*lsph)->nu  , alignment, N ,double);
  safe_check__aligned_alloc((*lsph)->rho , alignment, N ,double);
  safe_check__aligned_alloc((*lsph)->id  , alignment, N ,int64_t);
  safe_check__aligned_alloc((*lsph)->hash, alignment,2*N,int64_t);

  /*
  safe_check_alloc((*lsph)->x   , N ,double);
  safe_check_alloc((*lsph)->y   , N ,double);
  safe_check_alloc((*lsph)->z   , N ,double);
  safe_check_alloc((*lsph)->ux  , N ,double);
  safe_check_alloc((*lsph)->uy  , N ,double);
  safe_check_alloc((*lsph)->uz  , N ,double);
  safe_check_alloc((*lsph)->Fx  , N ,double);
  safe_check_alloc((*lsph)->Fy  , N ,double);
  safe_check_alloc((*lsph)->Fz  , N ,double);
  safe_check_alloc((*lsph)->nu  , N ,double);
  safe_check_alloc((*lsph)->rho , N ,double);
  safe_check_alloc((*lsph)->id  , N ,int64_t);
  safe_check_alloc((*lsph)->hash,2*N,int64_t);*/

finishlabel:

  if(success)
    return 0;
  else{
    if(*lsph==NULL)
      return 1;

    safe_free((*lsph)->x);  safe_free((*lsph)->y);  safe_free((*lsph)->z);
    safe_free((*lsph)->ux); safe_free((*lsph)->uy); safe_free((*lsph)->uz);
    safe_free((*lsph)->Fx); safe_free((*lsph)->Fy); safe_free((*lsph)->Fz);
    safe_free((*lsph)->nu); safe_free((*lsph)->rho);
    safe_free((*lsph)->id); safe_free((*lsph)->hash);

    return 1;
  }
}

int SPHparticleSOA_safe_free(int N,SPHparticle **lsph){

  if(*lsph==NULL)
    return 1;

  safe_free((*lsph)->x);  safe_free((*lsph)->y);  safe_free((*lsph)->z);
  safe_free((*lsph)->ux); safe_free((*lsph)->uy); safe_free((*lsph)->uz);
  safe_free((*lsph)->Fx); safe_free((*lsph)->Fy); safe_free((*lsph)->Fz);
  safe_free((*lsph)->nu); safe_free((*lsph)->rho);
  safe_free((*lsph)->id); safe_free((*lsph)->hash);

  free((*lsph));

  return 0;
}

int arg_parse(int argc, char **argv, int64_t *N, double *h,
              linkedListBox *box){
  bool intern_h = true;

  box->Xmin = -1.0; box->Ymin = -1.0; box->Zmin = -1.0;
  box->Xmax =  2.0; box->Ymax =  2.0; box->Zmax =  2.0;

  if(argc%2==0){
    printf("wrong number of arguments!\n");
    printf("Maybe an option is missing a value?\n");
  }

  for(int i=1;i<argc;i+=2){
    if( strcmp(argv[i],"-N") == 0 ){
      *N = (int64_t) atol(argv[i+1]);
      printf("N particles = %ld\n",*N);
    }
    else if( strcmp(argv[i],"-h") == 0 ){
      *h = atof(argv[i+1]);
      printf("h = %lf\n",*h);
      intern_h = false;
    }
    else if( strcmp(argv[i],"-Xmin") == 0 ){
      box->Xmin = atof(argv[i+1]);
      printf("Xmin = %lf\n",box->Xmin);
    }
    else if( strcmp(argv[i],"-Ymin") == 0 ){
      box->Ymin = atof(argv[i+1]);
      printf("Ymin = %lf\n",box->Ymin);
    }
    else if( strcmp(argv[i],"-Zmin") == 0 ){
      box->Zmin = atof(argv[i+1]);
      printf("Zmin = %lf\n",box->Zmin);
    }
    else if( strcmp(argv[i],"-Xmax") == 0 ){
      box->Xmax = atof(argv[i+1]);
      printf("Xmax = %lf\n",box->Xmax);
    }
    else if( strcmp(argv[i],"-Ymax") == 0 ){
      box->Ymax = atof(argv[i+1]);
      printf("Ymax = %lf\n",box->Ymax);
    }
    else if( strcmp(argv[i],"-Zmax") == 0 ){
      box->Zmax = atof(argv[i+1]);
      printf("Zmax = %lf\n",box->Zmax);
    }
    else if( strcmp(argv[i],"-Nx") == 0 ){
      box->Nx   = atol(argv[i+1]);
      printf("Nx = %d\n",box->Nx);
    }
    else if( strcmp(argv[i],"-Ny") == 0 ){
      box->Ny   = atol(argv[i+1]);
      printf("Ny = %d\n",box->Ny);
    }
    else if( strcmp(argv[i],"-Nz") == 0 ){
      box->Nz   = atol(argv[i+1]);
      printf("Nz = %d\n",box->Nz);
    }
    else{
      printf("unknown option: %s %s\n",argv[i],argv[i+1]);
    }
  }

  if(intern_h){
    box->Nx = (int)( (box->Xmax-box->Xmin)/(2*(*h)) );
    box->Ny = (int)( (box->Ymax-box->Ymin)/(2*(*h)) );
    box->Nz = (int)( (box->Zmax-box->Zmin)/(2*(*h)) );
  }

  box->N  = (box->Nx)*(box->Ny)*(box->Nz);

  double min_val = fmin((box->Xmax-box->Xmin)/box->Nx,fmin((box->Ymax-box->Ymin)/box->Ny,(box->Zmax-box->Zmin)/box->Nz));
  box->width  = (int)( 0.5 + 2*(*h)/min_val );
  box->hbegin = kh_init(0);
  box->hend   = kh_init(1);

  return 0;
}

int gen_unif_rdn_pos(int64_t N, int seed, SPHparticle *lsph){

	const gsl_rng_type *T=NULL;
	gsl_rng *r=NULL;

	if(lsph==NULL)
		return 1;

	gsl_rng_env_setup();

	T = gsl_rng_default;
	r = gsl_rng_alloc(T);
	gsl_rng_set(r,seed);

	for(int64_t i=0;i<N;i+=1){
		lsph->x[i] = gsl_rng_uniform(r);
		lsph->y[i] = gsl_rng_uniform(r);
		lsph->z[i] = gsl_rng_uniform(r);

		lsph->ux[i] = 0.0; lsph->Fx[i] = 0.0;
		lsph->uy[i] = 0.0; lsph->Fy[i] = 0.0;
    lsph->uz[i] = 0.0; lsph->Fz[i] = 0.0;

		lsph->nu[i]   = 1.0/N;
		lsph->rho[i]  = 0.0;
		lsph->id[i]   = (int64_t) i;
		lsph->hash[2*i+0] = (int64_t) 0;
		lsph->hash[2*i+1] = (int64_t) i;
	}

	gsl_rng_free(r);

	return 0;
}

int gen_unif_rdn_pos_box(int64_t N, int seed, linkedListBox *box,SPHparticle *lsph){

  const gsl_rng_type *T=NULL;
  gsl_rng *r=NULL;

  if(lsph==NULL)
    return 1;

  gsl_rng_env_setup();

  T = gsl_rng_default;
  r = gsl_rng_alloc(T);
  gsl_rng_set(r,seed);

  for(int64_t i=0;i<N;i+=1){
    lsph->x[i] = gsl_rng_uniform(r)*(box->Xmax-box->Xmin) + box->Xmin;
    lsph->y[i] = gsl_rng_uniform(r)*(box->Ymax-box->Ymin) + box->Ymin;
    lsph->z[i] = gsl_rng_uniform(r)*(box->Zmax-box->Zmin) + box->Zmin;

    lsph->ux[i] = 0.0; lsph->Fx[i] = 0.0;
    lsph->uy[i] = 0.0; lsph->Fy[i] = 0.0;
    lsph->uz[i] = 0.0; lsph->Fz[i] = 0.0;

    lsph->nu[i]   = 1.0/N;
    lsph->rho[i]  = 0.0;
    lsph->id[i]   = (int64_t) i;
    lsph->hash[2*i+0] = (int64_t) 0;
    lsph->hash[2*i+1] = (int64_t) i;
  }

  gsl_rng_free(r);

  return 0;
}

int print_time_stats(int runs, double *times){
  double t[5], dt[5], total_time, dtotal_time;

  printf("fast neighbour search / SoA / outer-openMP / symmetric load balanced\n");

  total_time = 0.;
  for(int k=0;k<5;k+=1){
    t[k]=0.; dt[k]=0.;
    for(int run=0;run<runs;run+=1)
      t[k] += times[5*run+k];
    t[k] /= runs;
    for(int run=0;run<runs;run+=1)
      dt[k] += (times[5*run+k]-t[k])*(times[5*run+k]-t[k]);
    dt[k] /= runs;
    dt[k] = sqrt(dt[k]);

    total_time += t[k];
  }

  dtotal_time = 0.;
  for(int run=0;run<runs;run+=1){
    double rgm = 0.;
    for(int k=0;k<5;k+=1)
      rgm += times[5*run+k];

    dtotal_time += (rgm-total_time)*(rgm-total_time);
  }
  dtotal_time /= runs;
  dtotal_time = sqrt(dtotal_time);

  printf("compute_hash_MC3D calc time                 : %.5lf +- %.6lf s : %.3lg%% +- %.3lg%%\n",t[0],dt[0],100*t[0]/total_time,100*dt[0]/total_time);
  printf("qsort calc time                             : %.5lf +- %.6lf s : %.3lg%% +- %.3lg%%\n",t[1],dt[1],100*t[1]/total_time,100*dt[1]/total_time);
  printf("reorder_lsph_SoA calc time                  : %.5lf +- %.6lf s : %.3lg%% +- %.3lg%%\n",t[2],dt[2],100*t[2]/total_time,100*dt[2]/total_time);
  printf("setup_interval_hashtables calc time         : %.5lf +- %.6lf s : %.3lg%% +- %.3lg%%\n",t[3],dt[3],100*t[3]/total_time,100*dt[3]/total_time);
  printf("compute_density_3d load balanced calc time  : %.5lf +- %.6lf s : %.3lg%% +- %.3lg%%\n",t[4],dt[4],100*t[4]/total_time,100*dt[4]/total_time);
  printf("compute_density_3d load balanced total time : %.5lf +- %.6lf s : %.3lf%%\n",total_time,dtotal_time,100.);

  return 0;
}

int compute_hash_MC3D(int64_t N, SPHparticle *lsph, linkedListBox *box){

	if(lsph==NULL)
		return 1;

	if(box==NULL)
		return 1;

	const double etax = box->Nx/(box->Xmax - box->Xmin);
	const double etay = box->Ny/(box->Ymax - box->Ymin);
	const double etaz = box->Nz/(box->Zmax - box->Zmin);

	for(int64_t i=0;i<N;i+=1){
		uint32_t kx,ky,kz;
		kx = (uint32_t)(etax*(lsph->x[i] - box->Xmin));
		ky = (uint32_t)(etay*(lsph->y[i] - box->Ymin));
		kz = (uint32_t)(etaz*(lsph->z[i] - box->Zmin));

		if((kx<0)||(ky<0)||(kz<0))
			return 1;
		else if((kx>=box->Nx)||(ky>=box->Nx)||(kz>=box->Nx))
			return 1;
		else{
			lsph->hash[2*i] = ullMC3Dencode(kx,ky,kz);
		}
	}

	return 0;
}

int compute_hash_MC2D(int64_t N, SPHparticle *lsph, linkedListBox *box){

	if(lsph==NULL)
		return 1;

	if(box==NULL)
		return 1;

	const double etax = box->Nx/(box->Xmax - box->Xmin);
	const double etay = box->Ny/(box->Ymax - box->Ymin);

	for(int64_t i=0;i<N;i+=1){
		uint32_t kx,ky;
		kx = (uint32_t)(etax*(lsph->x[i] - box->Xmin));
		ky = (uint32_t)(etay*(lsph->y[i] - box->Ymin));

		if((kx<0)||(ky<0))
			return 1;
		else if((kx>=box->Nx)||(ky>=box->Nx))
			return 1;
		else
			lsph->hash[2*i] = ullMC2Dencode(kx,ky);
	}

	return 0;
}

#define swap_loop(N,lsph,temp_swap,member,type) for(int64_t i=0;i<(N);i+=1)                            \
																	 	              (temp_swap)[i] = (lsph)->member[(lsph)->hash[2*i+1]];\
																		            memcpy((lsph)->member,temp_swap,(N)*sizeof(type))

int reorder_lsph_SoA(int N, SPHparticle *lsph, void *swap_arr){

	int64_t *int64_temp_swap = (int64_t *)swap_arr;
	swap_loop(N,lsph,int64_temp_swap,id ,int64_t);
	double *double_temp_swap = (double *)swap_arr;
	swap_loop(N,lsph,double_temp_swap,nu ,double);
	swap_loop(N,lsph,double_temp_swap,rho,double);
	swap_loop(N,lsph,double_temp_swap,x  ,double);
	swap_loop(N,lsph,double_temp_swap,y  ,double);
	swap_loop(N,lsph,double_temp_swap,z  ,double);
	swap_loop(N,lsph,double_temp_swap,ux ,double);
	swap_loop(N,lsph,double_temp_swap,uy ,double);
	swap_loop(N,lsph,double_temp_swap,uz ,double);
	swap_loop(N,lsph,double_temp_swap,Fx ,double);
	swap_loop(N,lsph,double_temp_swap,Fy ,double);
	swap_loop(N,lsph,double_temp_swap,Fz ,double);

	return 0;
}

int setup_interval_hashtables(int64_t N,SPHparticle *lsph,linkedListBox *box){

	int ret;
	int64_t hash0 = lsph->hash[2*0];
	khiter_t kbegin,kend;

	if(lsph==NULL)
		return 1;

	if(box==NULL)
		return 1;

	kbegin = kh_put(0, box->hbegin, lsph->hash[2*0], &ret); kh_value(box->hbegin, kbegin) = (int64_t)0;
	for(int64_t i=0;i<N;i+=1){
		lsph->hash[i] = lsph->hash[2*i];
		if(lsph->hash[i] == hash0)
			continue;
		hash0 = lsph->hash[i];
		kend   = kh_put(1, box->hend  , lsph->hash[i-1], &ret); kh_value(box->hend  , kend)   = i;
		kbegin = kh_put(0, box->hbegin, lsph->hash[i  ], &ret); kh_value(box->hbegin, kbegin) = i;
	}
	// TODO: possible bug? verify 2*(N-1) instead of N-1
	kend   = kh_put(1, box->hend  , lsph->hash[2*(N-1)], &ret); kh_value(box->hend  , kend)   = N;

	return 0;
}

int neighbour_hash_3d(int64_t hash,int64_t *nblist,int width, linkedListBox *box){
	int idx=0,kx=0,ky=0,kz=0;

	kx = ullMC3DdecodeX(hash);
	ky = ullMC3DdecodeY(hash);
	kz = ullMC3DdecodeZ(hash);

	for(int ix=-width;ix<=width;ix+=1){
		for(int iy=-width;iy<=width;iy+=1){
			for(int iz=-width;iz<=width;iz+=1){
				if((kx+ix<0)||(ky+iy<0)||(kz+iz<0))
					nblist[idx++] = -1;
				else if( (kx+ix>=box->Nx)||(ky+iy>=box->Ny)||(kz+iz>=box->Nz) )
					nblist[idx++] = -1;
				else if( kh_get(0, box->hbegin, ullMC3Dencode(kx+ix,ky+iy,kz+iz)) == kh_end(box->hbegin) )
					nblist[idx++] = -1;
				else
					nblist[idx++] = ullMC3Dencode(kx+ix,ky+iy,kz+iz);
			}
		}
	}
	return 0;
}

int neighbour_hash_2d(int64_t hash,int64_t *nblist,int width, linkedListBox *box){
	int idx=0,kx=0,ky=0;

	kx = ullMC2DdecodeX(hash);
	ky = ullMC2DdecodeY(hash);

	for(int ix=-width;ix<=width;ix+=1){
		for(int iy=-width;iy<=width;iy+=1){
				if((kx+ix<0)||(ky+iy<0))
					nblist[idx++] = -1;
				else if((kx+ix>=box->Nx)||(ky+iy>=box->Ny))
					nblist[idx++] = -1;
				else if( kh_get(0, box->hbegin, ullMC2Dencode(kx+ix,ky+iy)) == kh_end(box->hbegin) )
					nblist[idx++] = -1;
				else
					nblist[idx++] = ullMC2Dencode(kx+ix,ky+iy);
		}
	}
	return 0;
}

int print_boxes_populations(linkedListBox *box){
	khiter_t kbegin,kend;

	for (kbegin = kh_begin(box->hbegin); kbegin != kh_end(box->hbegin); kbegin++){
		if (kh_exist(box->hbegin, kbegin)){
			kend = kh_get(1, box->hend, kh_key(box->hbegin, kbegin));
			printf("hash - (begin,end) = %lu - (%lu,%lu): %ld\n",
				kh_key(box->hbegin, kbegin),
				kh_value(box->hbegin, kbegin),
				kh_value(box->hend, kend),
				((int64_t)kh_value(box->hend, kend))-((int64_t)kh_value(box->hbegin, kbegin)));
		}
	}

	return 0;
}

int print_neighbour_list_MC3D(int64_t Nmax,unsigned int width,unsigned int stride,SPHparticle *lsph,linkedListBox *box){

	int64_t nblist[(2*width+1)*(2*width+1)*(2*width+1)];

	for(int64_t i=0;i<Nmax;i+=(int64_t)stride){
		int res = neighbour_hash_3d(lsph->hash[i],nblist,width,box);
		if(res!=0)
			printf("invalid neighbour_hash_3d calculation\n");
		printf("origin hash %lu : (%d,%d,%d) \n",lsph->hash[i],ullMC3DdecodeX(lsph->hash[i]),
															  ullMC3DdecodeY(lsph->hash[i]),
															  ullMC3DdecodeZ(lsph->hash[i]));
		for(unsigned int j=0;j<(2*width+1)*(2*width+1)*(2*width+1);j+=1){
			if(nblist[j]>=0)
				printf("    neighbour hash %lu : (%d,%d,%d) \n",nblist[j],ullMC3DdecodeX(nblist[j]),
																		  ullMC3DdecodeY(nblist[j]),
																		  ullMC3DdecodeZ(nblist[j]));
			else
				printf("no neighbour here\n");
		}
		printf("\n\n");
	}

	return 0;
}

int print_neighbour_list_MC2D(int64_t Nmax,unsigned int width,unsigned int stride,SPHparticle *lsph,linkedListBox *box){

	int64_t nblist[(2*width+1)*(2*width+1)];

	for(int64_t i=0;i<Nmax;i+=(int64_t)stride){
		int res = neighbour_hash_2d(lsph->hash[i],nblist,width,box);
		if(res!=0)
			printf("invalid neighbour_hash_3d calculation\n");

		printf("origin hash %lu : (%d,%d) \n",lsph->hash[i],ullMC2DdecodeX(lsph->hash[i]),
															  											  ullMC2DdecodeY(lsph->hash[i]));
		for(unsigned int j=0;j<(2*width+1)*(2*width+1);j+=1){
			if(nblist[j]>=0)
				printf("    neighbour hash %lu : (%d,%d) \n",nblist[j],ullMC3DdecodeX(lsph->hash[i]),
																		  												 ullMC3DdecodeY(lsph->hash[i]));
			else
				printf("no neighbour here\n");
		}
		printf("\n\n");
	}

	return 0;
}

int print_neighbour_list_MC3D_lsph_file(int64_t Nmax,unsigned int width,unsigned int stride,SPHparticle *lsph,linkedListBox *box){
	FILE *fp;
	int64_t nblist[(2*width+1)*(2*width+1)*(2*width+1)];

	fp = fopen("data/nblist_MC3D.csv","w");

	const double etax = box->Nx/(box->Xmax - box->Xmin);
	const double etay = box->Ny/(box->Ymax - box->Ymin);
	const double etaz = box->Nz/(box->Zmax - box->Zmin);

	for(int64_t i=0;i<Nmax;i+=(int64_t)stride){
		uint32_t kx,ky,kz;
		int res = neighbour_hash_3d(lsph->hash[i],nblist,width,box);
		if(res!=0)
			printf("invalid neighbour_hash_3d calculation\n");

		kx = (uint32_t)(etax*(lsph->x[i] - box->Xmin));
		ky = (uint32_t)(etay*(lsph->y[i] - box->Ymin));
		kz = (uint32_t)(etaz*(lsph->z[i] - box->Zmin));

		fprintf(fp,"base   hash %lu : (%d,%d,%d) \n",lsph->hash[i],kx,ky,kz);
		fprintf(fp,"origin hash %lu : (%d,%d,%d) \n",lsph->hash[i],ullMC3DdecodeX(lsph->hash[i]),
					         									  												 ullMC3DdecodeY(lsph->hash[i]),
							    								  													 ullMC3DdecodeZ(lsph->hash[i]));

		for(unsigned int j=0;j<(2*width+1)*(2*width+1)*(2*width+1);j+=1)
			if(nblist[j]>=0){
				fprintf(fp,"  %lu:(%d,%d,%d)",nblist[j],ullMC3DdecodeX(nblist[j]),
																								ullMC3DdecodeY(nblist[j]),
																								ullMC3DdecodeZ(nblist[j]));
				if(j<(2*width+1)*(2*width+1)*(2*width+1)-1)
					fprintf(fp,",");
			}

		fprintf(fp,"\n\n");
	}

	fclose(fp);

	return 0;
}

int print_neighbour_list_MC3D_lsph_ids_file(int N, SPHparticle *lsph, linkedListBox *box){
	int res;
  khiter_t kbegin,kend;
  int64_t node_hash=-1,node_begin=0, node_end=0;
  int64_t nb_begin= 0, nb_end = 0;
  int64_t nblist[(2*box->width+1)*(2*box->width+1)*(2*box->width+1)];

  FILE *fp = fopen("data/nblist_ll.json","w");
  fprintf(fp,"{\n");
  for (kbegin = kh_begin(box->hbegin); kbegin != kh_end(box->hbegin); kbegin++){

    if (kh_exist(box->hbegin, kbegin)){

      kend = kh_get(1, box->hend, kh_key(box->hbegin, kbegin));
      node_hash = kh_key(box->hbegin, kbegin);
      node_begin = kh_value(box->hbegin, kbegin);
      node_end   = kh_value(box->hend, kend);

      res = neighbour_hash_3d(node_hash,nblist,box->width,box);
			if(res!=0)
				printf("invalid neighbour_hash_3d calculation\n");
      for(int64_t ii=node_begin;ii<node_end;ii+=1){

      	fprintf(fp,"\"%ld\":[",lsph->id[ii]);
      	for(int j=0;j<(2*box->width+1)*(2*box->width+1)*(2*box->width+1);j+=1){
      		if(nblist[j]>=0){
      			//nb_hash  = nblist[j];
          	nb_begin = kh_value(box->hbegin, kh_get(0, box->hbegin, nblist[j]) );
          	nb_end   = kh_value(box->hend  , kh_get(1, box->hend  , nblist[j]) );

          	for(int64_t jj=nb_begin;jj<nb_end;jj+=1)
          		fprintf(fp,"%ld, ",lsph->id[jj]);
      		}
      	}
      	if((kbegin==kh_end(box->hbegin)-1)&&(ii==node_end-1))
      		fprintf(fp,"-1]\n");
    		else
    			fprintf(fp,"-1],\n");
      }
    }
  }
  fprintf(fp,"}");

  fflush(fp);
  fclose(fp);

	return 0;
}

int print_neighbour_list_MC2D_lsph_file(int64_t Nmax,unsigned int width,unsigned int stride,SPHparticle *lsph,linkedListBox *box){
	FILE *fp;
	int64_t nblist[(2*width+1)*(2*width+1)];

	fp = fopen("data/nblist_MC2D.csv","w");

	const double etax = box->Nx/(box->Xmax - box->Xmin);
	const double etay = box->Ny/(box->Ymax - box->Ymin);

	for(int64_t i=0;i<Nmax;i+=(int64_t)stride){
		uint32_t kx,ky;
		int res = neighbour_hash_2d(lsph->hash[i],nblist,width,box);
		if(res!=0)
			printf("invalid neighbour_hash_3d calculation\n");

		kx = (uint32_t)(etax*(lsph->x[i] - box->Xmin));
		ky = (uint32_t)(etay*(lsph->y[i] - box->Ymin));

		fprintf(fp,"base   hash %lu : (%d,%d) \n",lsph->hash[i],kx,ky);
		fprintf(fp,"origin hash %lu : (%d,%d) \n",lsph->hash[i],ullMC2DdecodeX(lsph->hash[i]),
					         									  											ullMC2DdecodeY(lsph->hash[i]));

		for(unsigned int j=0;j<(2*width+1)*(2*width+1);j+=1)
			if(nblist[j]>=0){
				fprintf(fp,"  %lu:(%d,%d)",nblist[j], ullMC2DdecodeX(nblist[j]),
																    					ullMC2DdecodeY(nblist[j]));
				if(j<(2*width+1)*(2*width+1)-1)
					fprintf(fp,",");
			}

		fprintf(fp,"\n\n");
	}

	fclose(fp);

	return 0;
}

/*
int main(){

  int j=0,numThreads=6,err;
  int64_t N = 100000;
  double S=0.,S2=0.;
  linkedListBox *box;
  SPHparticle *lsph;

  lsph = (SPHparticle*)malloc(N*sizeof(SPHparticle));

  err = gen_unif_rdn_pos( N,123123123,lsph);

  box = (linkedListBox*)malloc(1*sizeof(linkedListBox));

  box->Nx = box->Ny = box->Nz = 10;
  box->N  = (box->Nx)*(box->Ny)*(box->Nz);
  box->Xmin.x = 0.0; box->Xmin.y = 0.0; box->Xmin.z = 0.0;
  box->Xmax.x = 1.0; box->Xmax.y = 1.0; box->Xmax.z = 1.0;
  box->hbegin = kh_init(0);
  box->hend = kh_init(1);

  err = compute_hash_MC3D(N,lsph,box);

  qsort(lsph,N,sizeof(SPHparticle),compare_SPHparticle);

  err = setup_interval_hashtables(N,lsph,box);

  //print_boxes_populations(box);
  //print_neighbour_list_MC3D(N,1,N/13,lsph,box);
  print_neighbour_list_MC3D_lsph_file(N,1,1,lsph,box);

  free(lsph);
  safe_free_box(box);

  return 0;
}*/

////////////////////////////////////////////////////////////////////
// https://stackoverflow.com/questions/16007640/openmp-parallel-quicksort

#define TASK_SIZE 24

int64_t partition(int64_t * a, int64_t p, int64_t r)
{
  int64_t lt[2*(r-p)];
  int64_t gt[2*(r-p)];
  int64_t i;
  int64_t j;
  int64_t key0 = a[2*r+0];
  int64_t key1 = a[2*r+1];
  int64_t lt_n = 0;
  int64_t gt_n = 0;

  for(i = p; i < r; i++){
    if(a[2*i] < a[2*r]){
      lt[2*lt_n+0] = a[2*i+0];
      lt[2*lt_n+1] = a[2*i+1];
      lt_n++;
    }else{
      gt[2*gt_n+0] = a[2*i+0];
      gt[2*gt_n+1] = a[2*i+1];
      gt_n++;
    }
  }

  for(i = 0; i < lt_n; i++){
    a[2*(p + i)+0] = lt[2*i+0];
    a[2*(p + i)+1] = lt[2*i+1];
  }

  a[2*(p + lt_n)+0] = key0;
  a[2*(p + lt_n)+1] = key1;

  for(j = 0; j < gt_n; j++){
    a[2*(p + lt_n + j + 1)+0] = gt[2*j+0];
    a[2*(p + lt_n + j + 1)+1] = gt[2*j+1];
  }

  return p + lt_n;
}

void quicksort_omp(int64_t * a, int64_t p, int64_t r)
{
	int64_t div;

  if(p < r){
    div = partition(a, p, r);
    #pragma omp task shared(a) if(r - p > TASK_SIZE)
    quicksort_omp(a, p, div - 1);
    #pragma omp task shared(a) if(r - p > TASK_SIZE)
    quicksort_omp(a, div + 1, r);
  }
}

max_stars_repo_name: jhelsas/sphalerite

max_issues_repo_name: jhelsas/sphalerite

max_forks_repo_name: jhelsas/sphalerite


text:
/*
 * Kabsch's algorithm: compute least-square-best transformation
 * Copyright (C) 2003, Arno Formella (formella@ei.uvigo.es)
 *
 * 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.
 *
 * This is a gsl-based implementation of Kabsch's algorithm presented in:
 *
 *   W. Kabsch, A solution for the best rotation to relate two sets of vectors,
 *   Acta Cryst. (1976), A32, 922-923
 *
 *   W. Kabsch, A discussion of the solution for the best rotation to relate
 *   two sets of vectors, Acta Cryst. (1978), A34, 827-828
 *
 * The code is C++-robust.
 *
 * More information about GSL
 * ==========================
 *
 * The project homepage is http://www.gnu.org/software/gsl/
 * The development site is http://sources.redhat.com/gsl/
 */

#ifndef KABSCH_H
#define KABSCH_H

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

#if defined(__cplusplus)
extern "C" {
#endif

/*
   preconditions:
     size > 0
     X and Y point to (size x 3)-matrices specifying 3D points x_i and y_i
     U points to a (3 x 3)-matrix
     t points to a (3)-vector
     s points to double
   postconditions if return 1:
     X and Y will be centralized at their respective origins
     U will hold the rotation
     t will hold the translation
     s will hold the scaling value
     such that:
       sum_i (U * s * x_i + t - y_i)^2  is minimal
   postconditions if return 0:
     X and Y will be centralized at their respective origins
     U will hold the identity
     t will hold the difference vector of the centroids t = cy -cx
     s will hold 1.0
   note:
     if s == NULL, no scaling is computed
*/
int kabsch(
  unsigned int size,
  gsl_matrix *X,
  gsl_matrix *Y,
  gsl_matrix *U,
  gsl_vector *t,
  double *s
);

#if defined(__cplusplus)
}
#endif

#endif

max_stars_repo_name: robpearc/FoldDesign

max_issues_repo_name: robpearc/FoldDesign

max_forks_repo_name: robpearc/FoldDesign


text:
#include <stdio.h>
#include <stdlib.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_movstat.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_statistics.h>

double
func(const size_t n, double x[], void * params)
{
  const double alpha = *(double *) params;

  gsl_sort(x, 1, n);

  return gsl_stats_trmean_from_sorted_data(alpha, x, 1, n);
}

int
main(void)
{
  const size_t N = 1000;                /* length of time series */
  const size_t K = 11;                  /* window size */
  double alpha = 0.1;                   /* trimmed mean parameter */
  gsl_vector *x = gsl_vector_alloc(N);  /* input vector */
  gsl_vector *y = gsl_vector_alloc(N);  /* filtered output vector for alpha1 */
  gsl_rng *r = gsl_rng_alloc(gsl_rng_default);
  gsl_movstat_workspace *w = gsl_movstat_alloc(K);
  gsl_movstat_function F;
  size_t i;
  double sum = 0.0;

  /* generate input signal */
  for (i = 0; i < N; ++i)
    {
      double ui = gsl_ran_gaussian(r, 1.0);
      double outlier = (gsl_rng_uniform(r) < 0.01) ? 10.0*GSL_SIGN(ui) : 0.0;
      sum += ui;
      gsl_vector_set(x, i, sum + outlier);
    }

  /* apply moving window function */
  F.function = func;
  F.params = α
  gsl_movstat_apply(GSL_MOVSTAT_END_PADVALUE, &F, x, y, w);

  /* print results */
  for (i = 0; i < N; ++i)
    {
      double xi = gsl_vector_get(x, i);
      double yi = gsl_vector_get(y, i);

      printf("%f %f\n", xi, yi);
    }

  gsl_vector_free(x);
  gsl_vector_free(y);
  gsl_rng_free(r);
  gsl_movstat_free(w);

  return 0;
}


text:
#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION

#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_multifit.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <Python.h>
#include <numpy/arrayobject.h>

void mark_debug (void)
{
    static int number;
    printf("---> Mark %d <---\n", number);
    fflush(stdout);
    number++;
}

void label_debug (const char * label)
{
    printf("---> Mark %s <---\n", label);
    fflush(stdout);
}

/* =========================== Pure-C source code =========================== */

/* Regularization workspace */
struct rglz_wksp {
    double lamb; /* Lambda value */
    gsl_vector *Ldiag; /* Pointer to vector with diagonal values of the
			      * regularization matrix L. */
    const gsl_vector *Y; /* Pointer to original Y */
    gsl_vector *Ys; /* Ponter to Y vector of the standard form. */
    const gsl_matrix *X; /* Pointer to original X matrix  */
    gsl_matrix *Xs; /* Pointer to X matrix of the standard form. */
    gsl_multifit_linear_workspace *lwksp; /* Pointer to GSL workspace where SVDs are stored. */
    /* L-curve data */
    struct lcurve {
	gsl_vector *reg_param;
	gsl_vector *rho;
	gsl_vector *eta;
    } lcurve;
};
typedef struct rglz_wksp rglz_wksp;

/* Allocate and initialize inner variables of a regularization workspace. */
rglz_wksp * rglz_init(const double *Ldiag, const int p,
		      const gsl_vector *Y, const gsl_matrix *X, const gsl_vector *W,
		      gsl_multifit_linear_workspace *lwksp)
{
    rglz_wksp *rw;

    /* Allocate space for the workspace */
    rw = (rglz_wksp*) malloc(sizeof(rglz_wksp));

    rw->lamb = 0.0; /* Initial value */
    rw->lcurve.reg_param = NULL;
    rw->lcurve.rho = NULL;
    rw->lcurve.eta = NULL;

    /* Allocate space for the standard-form matrixes */
    rw->Xs = gsl_matrix_alloc(X->size1, X->size2);
    rw->Ys = gsl_vector_alloc(Y->size);

    /* Initialize the regularization weights Ldiag */
    rw->Ldiag = gsl_vector_alloc(p);
    for (int i = 0; i < p; i++)
        gsl_vector_set(rw->Ldiag, i, Ldiag[i]);

    /* Convert the system to standard form */
    gsl_multifit_linear_wstdform1(rw->Ldiag, X, W, Y, rw->Xs, rw->Ys, lwksp);

    /* Calculate SVD of standard-form matrix Xs */
    gsl_multifit_linear_svd(rw->Xs, lwksp);

    rw->lwksp = lwksp;
    rw->X = X;
    rw->Y = Y;
    return rw;
}

void rglz_estimate_lambda(rglz_wksp *rw, const int lcsize)
{
    size_t idx; /* Temporary index to store corner of the L-curve */
    /* Allocate vectors */
    rw->lcurve.reg_param = gsl_vector_alloc(lcsize);
    rw->lcurve.rho = gsl_vector_alloc(lcsize);
    rw->lcurve.eta = gsl_vector_alloc(lcsize);
    /* Create Lcurve */
    gsl_multifit_linear_lcurve(rw->Y, rw->lcurve.reg_param, rw->lcurve.rho, rw->lcurve.eta, rw->lwksp);
#ifdef PRINT_LCURVE
    for (int i = 0; i < lcsize; i++)
    {
        printf("%12.4e%12.4e\n",
                gsl_vector_get(rw->lcurve.rho, i),
                gsl_vector_get(rw->lcurve.eta, i));
    }
#endif
    /* Estimate corner */
    gsl_multifit_linear_lcorner(rw->lcurve.rho, rw->lcurve.eta, &idx);
    /* Set lambda estimate */
    rw->lamb = gsl_vector_get(rw->lcurve.reg_param, idx);
    fprintf(stdout, "From gslpy_ridge: lambda = %.4e\n", rw->lamb);
}

void rglz_solve(gsl_vector *c, double *rnorm, double *snorm, rglz_wksp *rw)
{
    /* Solve in standard form */
    gsl_multifit_linear_solve(rw->lamb, rw->Xs, rw->Ys, c, rnorm, snorm, rw->lwksp);
    /* Convert back to original parameters */
    gsl_multifit_linear_genform1(rw->Ldiag, c, c, rw->lwksp);
    return;
}

inline void rglz_set_lambda(rglz_wksp *rw, const int lamb)
{
    rw->lamb = lamb;
}

void rglz_close(rglz_wksp *rw)
{
    gsl_vector_free(rw->Ldiag);
    gsl_matrix_free(rw->Xs);
    gsl_vector_free(rw->Ys);
    if (rw->lcurve.reg_param != NULL)
	gsl_vector_free(rw->lcurve.reg_param);
    if (rw->lcurve.rho != NULL)
	gsl_vector_free(rw->lcurve.rho);
    if (rw->lcurve.eta != NULL)
	gsl_vector_free(rw->lcurve.eta);
    free(rw);
}

int
gslpy_ridge_regression(const gsl_vector * y,
		       const gsl_matrix * X,
		       const gsl_vector * w,
		       const double lambda,
		       const double * L,
		       const unsigned long Lcurve_size,
		       gsl_vector * coefs)
{

    unsigned long n_samples = X->size1;
    unsigned long n_features = X->size2;

    double snorm;
    double chisq;

    gsl_multifit_linear_workspace
	* gmlw = gsl_multifit_linear_alloc(n_samples, n_features);

    rglz_wksp * rw = rglz_init(L, n_features, y, X, w, gmlw);

    /* If the user specified gamma, do not attempt to optimize it. */
    if (lambda >= 0)
	rglz_set_lambda(rw, lambda);
    else
	rglz_estimate_lambda(rw, Lcurve_size);

    rglz_solve(coefs, &chisq, &snorm, rw);

    rglz_close(rw);

    gsl_multifit_linear_free(gmlw);
    return 0;
}

/* ============================= Python wrapper ============================= */

static PyObject*
gslpy_ridge_fit(PyObject* self, PyObject* args, PyObject* kwargs)
{
    static char * keywords[] = {
	"X",
	"y",
	"sample_weight",
	"lamb",
	"ldiag",
	"lcurve_size",
	NULL,
    };

    PyArrayObject
	*X=NULL,
	*y=NULL,
	*w=NULL,
	*Ldiag=NULL,
	*coefs=NULL;

    /* Negative lambda means that lambda is optimized */
    double lambda = -1;
    unsigned long Lcurve_size = 10000;

    /* Flag to check if things are set */
    int wei_set = 0;
    int L_set = 0;

    /* GSL stuff */
    gsl_vector * g_y;
    gsl_matrix * g_X;
    gsl_vector * g_w;
    gsl_vector * g_coefs;

    if (!PyArg_ParseTupleAndKeywords(args,
                                     kwargs,
                                     "O!O!|O!dO!k",
                                     keywords,
                                     &PyArray_Type, &X,
				     &PyArray_Type, &y,
				     &PyArray_Type, &w,
				     &lambda,
				     &PyArray_Type, &Ldiag,
				     &Lcurve_size)) {
        return NULL;
    }

    /* If weights are not given, set them as 1 */
    if (w == NULL)
    {
	wei_set = 1;
	w = (PyArrayObject*) PyArray_ZEROS(PyArray_NDIM(y),
					   PyArray_DIMS(y),
					   NPY_DOUBLE,
					   0);

	npy_double * w_data = PyArray_DATA(w);

	for (npy_intp i = 0; i < PyArray_SIZE(w); ++i)
	    w_data[i] = 1.0;
    }

    /* If Ldiag is not given, set it to identity. */
    if (Ldiag == NULL)
    {
	npy_intp Lshape[] = {PyArray_DIM(X,1)};
	L_set = 1;
	Ldiag = (PyArrayObject*) PyArray_ZEROS(1,
					       Lshape,
					       NPY_DOUBLE,
					       0);

	for (npy_intp i = 0; i < PyArray_DIM(X, 1); i++)
	    ((double *) PyArray_DATA(Ldiag))[i] = 1.0;
    }

    /* Convert PyArray objects to gsl types */
    g_y = gsl_vector_calloc(PyArray_DIM(y, 0));
    g_w = gsl_vector_calloc(PyArray_DIM(w, 0));
    g_X = gsl_matrix_calloc(PyArray_DIMS(X)[0], PyArray_DIMS(X)[1]);
    g_coefs = gsl_vector_calloc(PyArray_DIM(X, 1));

    for (npy_intp i = 0; i < (npy_intp) g_y->size; i++)
    {
	gsl_vector_set(g_y, i, ((double *) PyArray_DATA(y))[i]);
	gsl_vector_set(g_w, i, ((double *) PyArray_DATA(w))[i]);
    }

    for (npy_intp i = 0; i < PyArray_SIZE(X); i++)
	g_X->data[i] = ((double *) PyArray_DATA(X))[i];

    /* Run regression */
    gslpy_ridge_regression(g_y, g_X,
			   g_w, lambda,
			   (double *) PyArray_DATA(Ldiag),
			   Lcurve_size, g_coefs);

    /* Convert g_coefs back to a Python object */
    /* Local scope just to delete ``coefs_shape``after. */
    {
	npy_intp coefs_shape[] = {g_coefs->size, 1};
	coefs = (PyArrayObject *) PyArray_ZEROS(2,
						coefs_shape,
						NPY_DOUBLE,
						0);
	for (npy_intp i = 0; i < (npy_intp) g_coefs->size; i++)
	    ((double *) PyArray_DATA(coefs))[i] = gsl_vector_get(g_coefs, i);

    }

    if (wei_set)
	Py_XDECREF(w);

    if (L_set)
	Py_XDECREF(Ldiag);

    gsl_matrix_free(g_X);
    gsl_vector_free(g_y);
    gsl_vector_free(g_w);
    gsl_vector_free(g_coefs);

    return (PyObject *) coefs;
}

/* =========================== Python boilerplate =========================== */

PyDoc_STRVAR(module_doc,
	     "Ridge regression model methods using the GNU Scientific Library.");

PyDoc_STRVAR(fit_doc,
	     "Fits a ridge regression model to data.\n"
	     "fit(X, y, sample_weights=None, lambda=None, ldiag=None, Lcurve_size=10000).\n");

PyMethodDef methods[] = {
    {
	"fit",
	(PyCFunction) gslpy_ridge_fit,
	METH_VARARGS | METH_KEYWORDS,
	fit_doc,
    },
    {NULL},
};

PyModuleDef gslpyridge_module = {
    PyModuleDef_HEAD_INIT,
    "gslpyridge",
    module_doc,
    -1,
    methods,
    NULL,
    NULL,
    NULL,
    NULL
};

PyMODINIT_FUNC
PyInit_gslpyridge(void)
{
    import_array();
    return PyModule_Create(&gslpyridge_module);
}

max_stars_repo_name: mssm-labmmol/profiler

max_issues_repo_name: mssm-labmmol/profiler

max_forks_repo_name: mssm-labmmol/profiler


text:
#include "constants.h"
#include "artisoptions.h"

#ifndef TYPES_H
#define TYPES_H

#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <time.h>
#ifndef __CUDA_ARCH__
#include <gsl/gsl_rng.h>
#endif
//#include <gsl/gsl_sf_expint.h>

struct time
{
  double start; // time at start of this timestep.
  double width; // Width of timestep.
  double mid; // Mid time in step - computed logarithmically.
  double gamma_dep; // cmf gamma ray energy deposition rate
  double positron_dep; // cmf positron energy deposition rate
  double electron_dep; // cmf positron energy deposition rate
  double alpha_dep; // cmf positron energy deposition rate
  double cmf_lum; // cmf luminosity light curve
  int pellet_decays; // Number of pellets that decay in this time step.
};

/// Coolinglist
///============================================================================
enum coolingtype {
  COOLINGTYPE_FF         = 880,
  COOLINGTYPE_FB         = 881,
  COOLINGTYPE_COLLEXC    = 882,
  COOLINGTYPE_COLLION    = 883,
};

typedef struct cellhistorycoolinglist_t
{
  enum coolingtype type;
  int element;
  int ion;
  int level;
  int upperlevel;
} cellhistorycoolinglist_t;

/*enum heatingtype {
  HEATINGTYPE_FF         = 884,
  HEATINGTYPE_BF         = 885,
  HEATINGTYPE_COLLDEEXC  = 886,
  HEATINGTYPE_COLLRECOMB = 887,
};*/

typedef struct heatingcoolingrates
{
  double cooling_collisional;
  double cooling_fb;
  double cooling_ff;
  double cooling_adiabatic;
  double heating_collisional;
  double heating_bf;
  double heating_ff;
  double heating_dep;
  double nt_frac_heating; // = heatingrates.gamma / total_gamma_deposition, = get_nt_frac_heating(modelgridindex) when T_e solver runs
} heatingcoolingrates_t;

/// PHIXSLIST

typedef struct fullphixslist_t
{
  double nu_edge;
  int element;
  int ion;
  int level;
  int phixstargetindex;
  int index_in_groundphixslist;
} fullphixslist_t;

typedef struct groundphixslist_t
{
  double nu_edge;
  //double photoion_contr;
  double gamma_contr;
  //double stimrecomb_contr;
  //double bfheating_contr;
  int element;
  int ion;
  int level;
  int phixstargetindex;
} groundphixslist_t;

typedef struct phixslist_t
{
  groundphixslist_t *groundcont;
  double *kappa_bf_contr;
#if (DETAILED_BF_ESTIMATORS_ON)
  double *gamma_contr;
#endif
} phixslist_t;

enum packet_type {
  TYPE_ESCAPE = 32,
  TYPE_RADIOACTIVE_PELLET = 100,
  TYPE_GAMMA = 10,
  TYPE_RPKT = 11,
  TYPE_KPKT = 12,
  TYPE_MA = 13,
  TYPE_NTLEPTON = 20,
  TYPE_PRE_KPKT = 120,
  TYPE_GAMMA_KPKT = 121,
};

enum cell_boundary {
  NEG_X = 101,
  POS_X = 102,
  NEG_Y = 103,
  POS_Y = 104,
  NEG_Z = 105,
  POS_Z = 106,
  NONE = 107,
};

typedef struct mastate_t
{
  int element;              /// macro atom of type element (this is an element index)
  int ion;                  /// in ionstage ion (this is an ion index)
  int level;                /// and level=level (this is a level index)
  int activatingline;       /// Linelistindex of the activating line for bb activated MAs, -99 else.
} mastate_t;

typedef struct packet
{
  int where;      /// The grid cell that the packet is in.
  enum packet_type type;       /// Identifies the type of packet (k-, r-, etc.)
  enum cell_boundary last_cross; /// To avoid rounding errors on cell crossing.
  int interactions;/// debug: number of interactions the packet undergone
  int nscatterings;   /// records number of electron scatterings a r-pkt undergone since it was emitted
  int last_event;  /// debug: stores information about the packets history
  double pos[3];  /// Position of the packet (x,y,z).
  double dir[3];  /// Direction of propagation. (x,y,z). Always a unit vector.
  double e_cmf;   /// The energy the packet carries in the co-moving frame.
  double e_rf;    /// The energy the packet carries in the rest frame.
  double nu_cmf;  /// The frequency in the co-moving frame.
  double nu_rf;   /// The frequency in the rest frame.
  int next_trans;  /// This keeps track of the next possible line interaction of a rpkt by storing
                   /// its linelist index (to overcome numerical problems in propagating the rpkts).
  int emissiontype;   /// records how the packet was emitted if it is a r-pkt
  double em_pos[3];   /// Position of the packet (x,y,z).
  int em_time;
  double prop_time; // internal clock to track how far in time the packet has been propagated
  int absorptiontype;     /// records linelistindex of the last absorption
                          /// negative values give ff-abs (-1), bf-abs (-2), compton scattering of gammas (-3),
                          /// photoelectric effect of gammas (-4), pair production of gammas (-5)
                          /// decaying pellets of the 52Fe chain (-6) and pellets which decayed before the
                          /// onset of the simulation (-7)
                          /// decay of a positron pellet (-10)
  int trueemissiontype;  // emission type coming from a kpkt to rpkt (last thermal emission)
  int trueem_time;
  double absorptionfreq;  /// records nu_cmf of packet at last absorption
  double absorptiondir[3]; /// Direction of propagation (x,y,z) when a packet was last absorbed in a line. Always a unit vector.
  //short timestep;
  double stokes[3]; //I, Q and U Stokes parameters
  double pol_dir[3]; //unit vector which defines the coordinate system against which Q and U are measured; should always be perpendicular to dir
  double tdecay;  /// Time at which pellet decays.
  enum packet_type escape_type; /// Flag to tell us in which form it escaped from the grid.
  int escape_time; /// Time at which is passes out of the grid.
                   /// Pos, dir, where, e_rf, nu_rf should all remain set at the exit point.
  int scat_count;  /// WHAT'S THAT???
  int number;     /// A unique number to identify which packet caused potential troubles.
  bool originated_from_particlenotgamma;    // first-non-pellet packet type was gamma
  int pellet_decaytype;           // index into decay::decaytypes
  int pellet_nucindex;           // nuclide index of the decaying species
  float trueemissionvelocity;
  mastate_t mastate;
} PKT;

enum ma_action {
  /// Radiative deexcitation rate from this level.
  MA_ACTION_RADDEEXC = 0,
  /// Collisional deexcitation rate from this level.
  MA_ACTION_COLDEEXC = 1,
  /// Radiative recombination from this level.
  MA_ACTION_RADRECOMB = 2,
  /// Collisional recombination rate from this level.
  MA_ACTION_COLRECOMB = 3,
  /// Rate for internal downward transitions to same ionisation stage.
  MA_ACTION_INTERNALDOWNSAME = 4,
  /// Rate for internal upward transitions to same ionisation stage.
  MA_ACTION_INTERNALDOWNLOWER = 5,
  /// Rate for internal downward transitions to lower ionisation stage.
  MA_ACTION_INTERNALUPSAME = 6,
  /// Rate for internal upward transitions to higher ionisation stage.
  MA_ACTION_INTERNALUPHIGHER = 7,
  /// Rate for internal upward transitions to higher ionisation stage due to non-thermal collisions.
  MA_ACTION_INTERNALUPHIGHERNT = 8,
  MA_ACTION_COUNT = 9,
};

/// GRID
///============================================================================
typedef struct compositionlist_entry
{
  float abundance;         /// Abundance of the element (by mass!).
  float *groundlevelpop;   /// Pointer to an array of floats which contains the groundlevel populations
                           /// of all included ionisation stages for the element.
  float *partfunct;        /// Pointer to an array of floats which contains the partition functions
                           /// of all included ionisation stages for the element.
  //float *ltepartfunct;     /// Pointer to an array of floats which contains the LTE partition functions
  //                         /// of all included ionisation stages for the element.
} compositionlist_entry;

typedef struct gridcell
{
  double pos_init[3]; /// Initial co-ordinates of inner most corner of cell.
  // int xyz[3];         /// Integer position of cell in grid.
  int modelgridindex;
} CELL;

typedef struct mgicooling_t
{
  double *contrib;
} mgicooling_t;

#define NSYN 1 /* number of frequency points in syn calculation */

enum ray_status {
  ACTIVE = 1,
  WAITING = 2,
  FINISHED = 3,
};

typedef struct syn_ray
{
  double tstart; /* time at which the ray enters the grid */
  double rstart[3]; /* vector position at which the ray enters the grid */
  double pos[3]; /* current position of the ray */
  double nu_rf[NSYN]; /*rest frame frequencies of rays */
  double nu_rf_init[NSYN]; /*rest frame frequencies of rays at start - for a test only */
  double nu_cmf[NSYN]; /*cmf freqencies */
  double e_rf[NSYN]; /*rest frame energy of rays */
  double e_cmf[NSYN]; /* cmf energies of rays */
  int where; /* the grid cell the ray is in */
  int lindex[NSYN]; /* array of ray positions in the line list */
  enum cell_boundary last_cross; /* last boundary crossed */
  enum ray_status status; /*WAITING, then ACTIVE then FINISHED*/
} RAY;

/// ATOMIC DATA
///============================================================================

typedef struct phixstarget_entry
{
  double *spontrecombcoeff;
  #if (!NO_LUT_PHOTOION)
    double *corrphotoioncoeff;
  #endif
  #if (!NO_LUT_BFHEATING)
  double *bfheating_coeff;
  #endif
  double *bfcooling_coeff;

  double probability;        // fraction of phixs cross section leading to this final level
  int levelindex;         // index of upper ion level after photoionisation
} phixstarget_entry;

typedef struct levellist_entry
{
  double epsilon;                            /// Excitation energy of this level relative to the neutral ground level.
  int *uptrans_lineindicies;    /// Allowed upward transitions from this level
  int *downtrans_lineindicies;  /// Allowed downward transitions from this level
  int nuptrans;
  int ndowntrans;
  double phixs_threshold;                    /// Energy of first point in the photion_xs table
  phixstarget_entry *phixstargets;  /// pointer to table of target states and probabilities
  float *photoion_xs;               /// Pointer to a lookup-table providing photoionisation cross-sections for this level.
  int nphixstargets;                         /// length of phixstargets array:
  float stat_weight;                           /// Statistical weight of this level.

  int cont_index;                            /// Index of the continuum associated to this level. Negative number.
#if (!NO_LUT_PHOTOION || !NO_LUT_BFHEATING)
  int closestgroundlevelcont;
#endif

  bool metastable;                            ///

//  double photoion_xs_nu_edge;              /// nu of the first grid point in the photoion_xs lookup-table.

//double *spontrecombcoeff_E;
//double *photoioncoeff_below;
//double *photoioncoeff_above;
//double *corrphotoioncoeff_above;
//double *bfheating_coeff_above;
//double *stimulated_bfcooling_coeff;
//double *stimulated_recomb_coeff;

//float *modified_spontrecombcoeff;
//float *stimrecombcoeff;
//float *modified_stimrecombcoeff;

///float population;
///float sahafact;
///float spontaneousrecomb_ratecoeff;
///float modifiedspontaneousrecomb_ratecoeff;
///float corrphotoion_ratecoeff;
///float modifiedcorrphotoion_ratecoeff;

/// time dependent macroatom event rates
//double rad_deexc;                          /// Radiative deexcitation rate from this level.
//double internal_down_same;                 /// Rate for internal downward transitions within the same ionisation stage.
//double internal_up_same;                   /// Rate for internal upward transitions within the same ionisation stage.

} levellist_entry;

typedef struct ionlist_entry
{
  levellist_entry *levels;                   /// Carries information for each level: 0,1,...,nlevels-1
  int ionstage;                              /// Which ionisation stage: XI=0, XII=1, XIII=2, ...
  int nlevels;                               /// Number of levels for this ionisation stage
  int nlevels_nlte;                          /// number of nlte levels for this ion
  int first_nlte;                            /// index into nlte_pops array of a grid cell
  int ionisinglevels;                        /// Number of levels which have a bf-continuum
  int maxrecombininglevel;                   /// level index of the highest level with a non-zero recombination rate
  int nlevels_groundterm;
  int coolingoffset;
  int ncoolingterms;
  int uniqueionindex;
  float *Alpha_sp;
  double ionpot;                             /// Ionisation threshold to the next ionstage
  //int nbfcontinua;
  //ionsphixslist_t *phixslist;
} ionlist_entry;

typedef struct elementlist_entry
{
  ionlist_entry *ions;              /// Carries information for each ion: 0,1,...,nions-1
  int nions;                                 /// Number of ions for the current element
  int anumber;                               /// Atomic number
//  int uppermost_ion;                       /// Highest ionisation stage which has a decent population for a given cell
                                             /// Be aware that this must not be used outside of the update_grid routine
                                             /// and their daughters. Neither it will work with OpenMP threads.
  float abundance;                           ///
  float mass;                                /// Atomic mass number in multiple of MH
} elementlist_entry;

typedef struct linelist_entry
{
  double nu;                                 /// Frequency of the line transition
  float einstein_A;
  float osc_strength;
  float coll_str;
  int elementindex;                        /// It's a transition of element (not its atomic number,
                                           /// but the (x-1)th element included in the simulation.
  int ionindex;                            /// The same for the elements ion
  int upperlevelindex;                     /// And the participating upper
  int lowerlevelindex;                     /// and lower levels
  bool forbidden;
} linelist_entry;

typedef struct bflist_t
{
  int elementindex;
  int ionindex;
  int levelindex;
  int phixstargetindex;
} bflist_t;

typedef struct nne_solution_paras
{
  int cellnumber;
} nne_solution_paras;

typedef struct gslintegration_paras
{
  double nu_edge;
  double T;
  float *photoion_xs;
} gslintegration_paras;

typedef struct rpkt_cont_opacity_struct
{
  double nu; // frequency at which opacity was calculated
  double total;
  double es;
  double ff;
  double bf;
  double ffheating;
  //double bfheating;
  int modelgridindex;
  bool recalculate_required; // e.g. when cell or timestep has changed
} rpkt_cont_opacity_struct;

/// Cell history
///============================================================================
// typedef struct coolinglist_contributions
// {
//   double contribution;
// } coolinglist_contributions;

typedef struct
{
  double corrphotoioncoeff;
#if (SEPARATE_STIMRECOMB)
  double stimrecombcoeff;
#endif
} chphixstargets_struct;

typedef struct chlevels_struct
{
  double processrates[MA_ACTION_COUNT];
  chphixstargets_struct *chphixstargets;
  double bfheatingcoeff;
} chlevels_struct;

typedef struct chions_struct
{
  chlevels_struct *chlevels;              /// Pointer to the ions levellist.
} chions_struct;

typedef struct chelements_struct
{
  chions_struct *chions;                  /// Pointer to the elements ionlist.
} chelements_struct;

typedef struct cellhistory_struct
{
  double *cooling_contrib;    /// Cooling contributions by the different processes.
  chelements_struct *chelements;            /// Pointer to a nested list which helds compositional
                                            /// information for all the elements=0,1,...,nelements-1
  int cellnumber;                           /// Identifies the cell the data is valid for.
  int bfheating_mgi;
} cellhistory_struct;

#endif //TYPES_H

max_stars_repo_name: artis-mcrt/artis

max_issues_repo_name: artis-mcrt/artis

max_forks_repo_name: artis-mcrt/artis