text:
/**
 */
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include "aXe_grism.h"
#include <math.h>
#include <string.h>
#include <ctype.h>
#include "drz2pet_utils.h"

/**
 * Function: make_spc_drztable
 * The function generates and fills an ap_pixel vector for one
 * drizzled 2D grism image.
 *
 * Parameters:
 * @param  obs     - pointer to observation structure with image values
 * @param  obs     - pointer to observation structure with image values
 * @param  ob      - pointer to the description of the current beam
 * @param  lambda0 - the starting wavelength in this drizzled image
 * @param  dlambda - the dispersion in this drizzled image
 *
 * Returns:
 * @return table   - pointer to an ap_pixel vector
 */
ap_pixel *
make_spc_drztable(observation * const obs, observation * const wobs,
                  const object *ob, const double lambda0,
                  const double dlambda){

  ap_pixel *table, *cur_ap;

  int npix1, npix2;
  int nsub1=1, nsub2=1;
  int i, j;
  int beamInt = 0;

  npix1 = obs->grism->size1;
  npix2 = obs->grism->size2;

  if (!(table = malloc ((npix1 * npix2 + 1) * sizeof (ap_pixel) * nsub1 * nsub2)))
    return NULL;

  cur_ap = table;

  for (i = 0; i < npix1; i++){
    for (j = 0; j < npix2; j++){

      //* check whether the pixel is within the extraction box
      if (fabs((double)j - ob->beams[beamInt].refpoint.y) > ob->beams[beamInt].width + .5)
        continue;

      //* fill the ap_pixel
      cur_ap->p_x  = i;
      cur_ap->p_y  = j;

      cur_ap->x       = (double)i - ob->beams[beamInt].refpoint.x;
      cur_ap->y       = (double)j - ob->beams[beamInt].refpoint.y;
      cur_ap->dist    = (double)j - ob->beams[beamInt].refpoint.y;
      cur_ap->xs      = (double)i - ob->beams[beamInt].refpoint.x;
      cur_ap->ys      = 0.0;
      cur_ap->dxs     = 0.0;  // be reminded: dxs is the LOCAL TRACEANGLE, not the pixels size..
      cur_ap->xi      = (double)i - ob->beams[beamInt].refpoint.x;
      cur_ap->lambda  = lambda0 + dlambda * ((double)i - ob->beams[beamInt].refpoint.x);
      cur_ap->dlambda = dlambda;

      // transfer the cps values and the weight values
      cur_ap->count  = gsl_matrix_get(obs->grism, i, j);
      cur_ap->weight = gsl_matrix_get(wobs->grism, i, j);

      cur_ap->error  = gsl_matrix_get(obs->pixerrs, i, j);
      cur_ap->contam = gsl_matrix_get(obs->dq, i, j);

      // dq is always set to zero, bad (=unexposed) pixels
      // have exptime=weight=0.0
      cur_ap->dq     = 0;

      // go to the next
      // PET pixel
      cur_ap++;
    }
  }

  cur_ap->p_x = -1;
  cur_ap->p_y = -1;
  cur_ap->count = -1;

  return table;
}

/**
 * Function: get_ID_num
 * The function generates the extension name used in the PET
 * for a 2D drizzled grism image
 *
 * Parameters:
 * @param  grism_image - the name of the 2D drizzled grism image
 * @param  extname     - The ID number of the image
 *
 * Returns:
 * @return WorkPtr     - Pointer to the ID-string
 */
char
*get_ID_num (char *grism_image, char *extname){

  char *WorkPtr;
  char *WorkPtr2;
  char **t_err=NULL;
  char tmp[MAXCHAR];
  char ID[4];

  int len;
  int i=0;

  strcpy(ID,"_ID");
  strcpy(tmp, grism_image);

  t_err = (char **) malloc(sizeof(char *)*1);

  WorkPtr = strstr(tmp, ID)+3;
  WorkPtr2 = WorkPtr + strlen (WorkPtr);
  while (isspace ((int) *WorkPtr2))
    *WorkPtr2-- = '\0';
  WorkPtr2 = WorkPtr2 - 1;
  len = strlen(WorkPtr);
  while  (i++ < len){
    strtol(WorkPtr,t_err,10);
    if (strlen(t_err[0])==0)
      break;
    *WorkPtr2-- = '\0';
  }
  strcpy(extname,WorkPtr);
  free(t_err);
  return WorkPtr;
}

/**
 * Function: normalize_weight
 * The function computes the weights used in the extraction
 * from the exposure times given for each pixel in a column.
 *
 * Parameters:
 * @param  wobs - pointer to the structure which carries the pixel data
 * @param  ob   - pointer to the beam description
 */
void
normalize_weight(observation *wobs, object *ob,
                 const int opt_extr)
{
  gsl_matrix *weight;
  //int i, j;
 // // check for optimal extraction
  if (opt_extr)
    {
      // let the optimal weights compute
      fprintf(stdout, "aXe_DRZ2PET: Use otpimal weighting.\n");
      weight = comp_opt_weight(wobs->pixerrs, wobs->dq, ob);
    }
  else
    {
      // let the exposure time weights compute
      fprintf(stdout, "aXe_DRZ2PET: Use normal, flat  weighting.\n");
      weight = comp_equ_weight(wobs->grism, ob);

      // code lines for exposure time weighting
      //fprintf(stdout, "use exposure time  weighting\n");
      //weight = comp_exp_weight(wobs->grism, ob);
    }

  // replace the weight extension with the new weights
  gsl_matrix_free(wobs->grism);
  wobs->grism = weight;
}

/*
 * Function: comp_equ_weight
 * The function creates and fills a weight map with flat entries.
 * Pixels with non-zero exposure time which which contribute to the
 * extracted spectrum all receive the weight 1.0, all other
 * get the weight 0.0.
 *
 * Parameters:
 * @param exp_map - the matrix with the exposure times
 * @param ob      - the object structure
 *
 * Returns:
 * @return weight - the matrix with the exposure time weights
 */
gsl_matrix *
comp_equ_weight(gsl_matrix *exp_map, const object *ob)
{
  gsl_matrix *weight;
  int i, j;
  int beamInt = 0;

  // allocate the matrix and set it to the default
  weight = gsl_matrix_alloc(exp_map->size1, exp_map->size2);
  gsl_matrix_set_all(weight, 0.0);

   // go over all columns
  for (i=0; i < (int)exp_map->size1; i++)
    // go over all rows
    for (j=0; j < (int)exp_map->size2; j++)
      // check whether the pixel is inside the extraction
      // area and has a positive exposure time
      if (fabs(ob->beams[beamInt].refpoint.y-(double)j) <= ob->beams[beamInt].width+0.5 &&
          gsl_matrix_get(exp_map, i, j) > 0.0 )
        // set the weight to one
        gsl_matrix_set(weight, i,j,1.0);

  // return the resulting weight map
  return weight;
}

/*
 * Function: comp_exp_weight
 * The function uses the exposure time map and computes
 * the weights which are proportional to the relative
 * exposure time of the pixel. Pixels outside of the extraction
 * area are set to 1000.0 Inside the default is 1.0 in unexposed
 * columns and 0.0 within exposed columns.
 *
 * Parameters:
 * @param exp_map - the matrix with the exposure times
 * @param ob      - the object structure
 *
 * Returns:
 * @return weight - the matrix with the exposure time weights
 */
gsl_matrix *
comp_exp_weight(gsl_matrix *exp_map, const object *ob)
{
  gsl_matrix *weight;
  int i, j;
  double sum, contr, norm, allweight;
  int beamInt = 0;

  // allocate the matrix and set it to the default
  weight = gsl_matrix_alloc(exp_map->size1, exp_map->size2);
  gsl_matrix_set_all(weight, 1000.0);

    //* go over all columns
  for (i=0; i < (int)exp_map->size1; i++)
    {
      sum = 0.0;
      contr = 0.0;
      allweight = 0.0;

      // determine for each column the total exposure time
      // and the number of pixels with non-zero exposure time
      for (j=0; j < (int)exp_map->size2; j++)
        {
          // check whether the pixel is inside the extraction
          // area and has a positive weight
          if (fabs(ob->beams[beamInt].refpoint.y-(double)j) <= ob->beams[beamInt].width+0.5 &&
              gsl_matrix_get(exp_map, i, j) > 0.0 )
            {
              // sum up the weigth and the number of pixels
              sum = sum + gsl_matrix_get(exp_map, i, j);
              contr = contr + 1.0;
            }
        }

      // check whether the row has exposure time at all
      if (sum > 0.0)
        {
          // if the row is exposed:
          // determine the mean exposure time
          norm = sum / contr;

          // go over each pixel
          for (j=0; j < (int)exp_map->size2; j++)
            {
              // check whether the pixel is inside the extraction are
              if (fabs(ob->beams[beamInt].refpoint.y-(double)j) <= ob->beams[beamInt].width+0.5)
                {
                  // compute and store the relative exposure time (==weight)
                  // for positive exposure time
                  if (gsl_matrix_get(exp_map, i, j) > 0.0)
                    {
                      gsl_matrix_set(weight, i,j,gsl_matrix_get(exp_map, i, j)/norm);
                      allweight = allweight + gsl_matrix_get(weight, i, j);
                    }
                  // make a default to avoid values <0.0
                  else
                    {
                      gsl_matrix_set(weight, i,j,0.0);
                    }
                }
            }
        }
      // treat unexposed columns differently:
      else
        {

          // go over each row
          for (j=0; j < (int)exp_map->size2; j++)
            {
              // check whether the pixel is inside the extraction are
              if (fabs(ob->beams[beamInt].refpoint.y-(double)j) <= ob->beams[beamInt].width+0.5)
                {
                  // set inside pixels to the defaulkt 1.0
                  gsl_matrix_set(weight, i,j,1.0);
                  allweight = allweight + gsl_matrix_get(weight, i, j);
                }
            }
        }
    }

  // return the resulting weight map
  return weight;
}

/*
 * Function: comp_opt_weight
 * The function uses the model map and and the variance
 * map and computes optimale weigths according to the
 * Hoorne method.
 * Pixels outside of the extraction
 * area are set to 1000.0
 *
 * Parameters:
 * @param mod_map - the matrix with the model values
 * @param var_map - the matrix with the variance values
 * @param ob      - the object structure
 *
 * Returns:
 * @return weight - the matrix with the exposure time weights
 */
gsl_matrix *
comp_opt_weight(gsl_matrix *mod_map,
                gsl_matrix *var_map, const object *ob)
{
  gsl_matrix *weight;
  int i, j;
  double mod_sum, weight_sum, contr, norm, allweight;
  int beamInt = 0;
  double mod_val;
  //double var_val;

  // allocate the weight matrix and set the default
  weight = gsl_matrix_alloc(mod_map->size1, mod_map->size2);
  gsl_matrix_set_all(weight, 0.0);

  //* go over all columns
  for (i=0; i < (int)mod_map->size1; i++)
    {
      mod_sum = 0.0;
      contr = 0.0;
      allweight = 0.0;
      weight_sum=0.0;

      // determine for each column the total model counts
      // and the number of pixels with non-zero model_counts
      for (j=0; j < (int)mod_map->size2; j++)
        {
          if (fabs(ob->beams[beamInt].refpoint.y-(double)j) <= ob->beams[beamInt].width+0.5 &&
              gsl_matrix_get(mod_map, i, j) > 0.0 )
            {
              mod_sum = mod_sum + gsl_matrix_get(mod_map, i, j);
              contr = contr + 1.0;
            }
        }

      // check whether the column has model values.
      // normalize the model values and compute
      // optimal weights if yes.
      if (mod_sum > 0.0)
        {
          //* determine the mean model value
          //      norm = mod_sum / contr;
          norm = mod_sum;

          // go over each row
          for (j=0; j < (int)mod_map->size2; j++)
            {
              // check whether the pixel is inside the extraction area
              if (fabs(ob->beams[beamInt].refpoint.y-(double)j) <= ob->beams[beamInt].width+0.5 &&
                  gsl_matrix_get(mod_map, i, j) > 0.0 )
                {
                  // normalize the model counts
                  mod_val = gsl_matrix_get(mod_map, i, j)/norm;

                  // store the normalized model counts
                  gsl_matrix_set(mod_map, i, j, mod_val);

                  // add up the normalization value for the optimal weights
                  weight_sum = weight_sum + mod_val*mod_val/gsl_matrix_get(var_map, i, j);
                }
            }

          // finally compute and write the weights:
          // go over each pixel
          for (j=0; j < (int)mod_map->size2; j++)
            {
              // check whether the pixel is inside the extraction area
              if (fabs(ob->beams[beamInt].refpoint.y-(double)j) <= ob->beams[beamInt].width+0.5)
                {
                  if (gsl_matrix_get(mod_map, i, j) > 0.0 )
                    {
                      // compute and set the individual pixel weight
                      gsl_matrix_set(weight, i, j, gsl_matrix_get(mod_map, i, j)/gsl_matrix_get(var_map, i, j)/weight_sum);
                      allweight = allweight + gsl_matrix_get(weight, i, j);
                    }
                  else
                    {
                      // set the default value if no model value is zero
                      gsl_matrix_set(weight, i, j,0.0);
                    }
                }
            }
        }
      // if the column does not have model values at all:
      else
        {
          // go over each row
          for (j=0; j < (int)mod_map->size2; j++)
            {
              // check whether the pixel is within the extraction area
              if (fabs(ob->beams[beamInt].refpoint.y-(double)j) <= ob->beams[beamInt].width+0.5)
                {
                  // set the inside default value
                  gsl_matrix_set(weight, i,j,1.0);
                  allweight = allweight + gsl_matrix_get(mod_map, i, j);
                }
            }
        }
    }

  // return the weight matrix
  return weight;
}

max_stars_repo_name: sosey/pyaxe

max_issues_repo_name: sosey/pyaxe

max_forks_repo_name: sosey/pyaxe


text:
/// based on https://mklimenko.github.io/english/2018/06/04/constexpr-random/
#pragma once

#include <atomic>
#include <cstdint>
#include <gsl/gsl>
#include <limits>
#include <random>

#include "HyperionUtils/Concepts.h"

namespace hyperion::math {
	using gsl::narrow_cast;
#ifndef _MSC_VER
	using std::size_t;
#endif //_MSC_VER

// clang-format off
#ifndef _MSC_VER
	// NOLINTNEXTLINE
	#define IGNORE_DATETIME_START \
		_Pragma("GCC diagnostic push") \
		_Pragma("GCC diagnostic ignored \"-Wdate-time\"")
#else
	// NOLINTNEXTLINE
	#define IGNORE_DATETIME_START
#endif
// clang-format on

// clang-format off
#ifndef _MSC_VER
	// NOLINTNEXTLINE
	#define IGNORE_DATETIME_STOP \
		_Pragma("GCC diagnostic pop")
#else
	// NOLINTNEXTLINE
	#define IGNORE_DATETIME_STOP
#endif
	// clang-format on

// clang-format off
#ifndef _MSC_VER
	// NOLINTNEXTLINE
	#define IGNORE_WEAK_VTABLES_START \
		_Pragma("GCC diagnostic push") \
		_Pragma("GCC diagnostic ignored \"-Wweak-vtables\"")
#else
	// NOLINTNEXTLINE
	#define IGNORE_WEAK_VTABLES_START
#endif
// clang-format on

// clang-format off
#ifndef _MSC_VER
	// NOLINTNEXTLINE
	#define IGNORE_WEAK_VTABLES_STOP \
		_Pragma("GCC diagnostic pop")
#else
	// NOLINTNEXTLINE
	#define IGNORE_WEAK_VTABLES_STOP
#endif
	// clang-format on

	IGNORE_WEAK_VTABLES_START
	class Engine {
	  public:
		constexpr Engine() noexcept = default;
		constexpr Engine(const Engine& engine) noexcept = default;
		constexpr Engine(Engine&& engine) noexcept = default;
		virtual constexpr ~Engine() noexcept = default;

		[[nodiscard]] virtual constexpr auto get_seed() const noexcept -> size_t = 0;
		virtual constexpr auto seed(size_t seed) noexcept -> void = 0;
		virtual constexpr auto generate() noexcept -> size_t = 0;

		template<size_t size>
		inline constexpr auto generate_array() noexcept -> std::array<size_t, size> {
			auto array = std::array<size_t, size>();
			for(auto& element : array) {
				element = generate();
			}
		}

		[[nodiscard]] virtual constexpr auto max_value() const noexcept -> size_t = 0;

		constexpr auto operator=(const Engine& engine) noexcept -> Engine& = default;
		constexpr auto operator=(Engine&& engine) noexcept -> Engine& = default;

		virtual constexpr auto operator()() noexcept -> size_t = 0;
	};
	IGNORE_WEAK_VTABLES_STOP

	template<size_t max = 714025>
	class LinearCongruentialEngine final : public Engine {
	  public:
		constexpr LinearCongruentialEngine() noexcept = default;
		constexpr LinearCongruentialEngine(
			const LinearCongruentialEngine& engine) noexcept = default;
		constexpr LinearCongruentialEngine(LinearCongruentialEngine&& engine) noexcept = default;
		constexpr ~LinearCongruentialEngine() noexcept final = default;

		IGNORE_DATETIME_START
		[[nodiscard]] inline constexpr auto get_seed() const noexcept -> size_t final {
			const auto* t = __TIME__;
			return time_from_string(t, 0) * 60 * 60 + time_from_string(t, 3) * 60
				   + time_from_string(t, 6);
		}
		IGNORE_DATETIME_STOP

		inline constexpr auto seed(size_t seed) noexcept -> void final {
			m_previous = seed;
		}

		[[nodiscard]] inline constexpr auto generate() noexcept -> size_t final {
			m_previous = (m_a * m_previous + m_c) % max;
			return m_previous;
		}

		[[nodiscard]] inline constexpr auto max_value() const noexcept -> size_t final {
			return max;
		}

		constexpr auto operator=(const LinearCongruentialEngine& engine) noexcept
			-> LinearCongruentialEngine& = default;
		constexpr auto operator=(LinearCongruentialEngine&& engine) noexcept
			-> LinearCongruentialEngine& = default;

		inline constexpr auto operator()() noexcept -> size_t final {
			return generate();
		}

	  private:
		size_t m_a = 4096;
		size_t m_c = 150889;
		size_t m_previous = get_seed();

		[[nodiscard]] static inline constexpr auto
		time_from_string(const char* string, int offset) noexcept -> size_t {
			return narrow_cast<size_t>(string[offset] - '0') * 10	// NOLINT
				   + narrow_cast<size_t>(string[offset + 1] - '0'); // NOLINT
		}
	};

	IGNORE_WEAK_VTABLES_START
	template<typename EngineType, utils::concepts::Numeric T = int>
	requires utils::concepts::Derived<EngineType, Engine>
	class Distribution {
	  public:
		constexpr Distribution() noexcept = default;
		constexpr Distribution(const Distribution& distribution) noexcept = default;
		constexpr Distribution(Distribution&& distribution) noexcept = default;
		virtual constexpr ~Distribution() noexcept = default;

		virtual constexpr auto normalized_random_value() noexcept -> double = 0;
		virtual constexpr auto random_value() noexcept -> T = 0;
		template<size_t size>
		inline constexpr auto normalized_random_values() noexcept -> std::array<double, size> {
			auto array = std::array<double, size>();
			for(auto& element : array) {
				element = normalized_random_value();
			}

			return array;
		}
		template<size_t size>
		inline constexpr auto random_values() noexcept -> std::array<T, size> {
			auto array = std::array<T, size>();
			for(auto& element : array) {
				element = random_value();
			}

			return array;
		}
		virtual constexpr auto seed(size_t seed) noexcept -> void = 0;
		[[nodiscard]] virtual constexpr auto get_seed() noexcept -> size_t = 0;

		[[nodiscard]] virtual constexpr auto get_min() const noexcept -> T = 0;
		virtual constexpr auto set_min(T min) noexcept -> void = 0;
		[[nodiscard]] virtual constexpr auto get_max() const noexcept -> T = 0;
		virtual constexpr auto set_max(T max) noexcept -> void = 0;

		constexpr auto
		operator=(const Distribution& distribution) noexcept -> Distribution& = default;
		constexpr auto operator=(Distribution&& distribution) noexcept -> Distribution& = default;

		virtual constexpr auto operator()() noexcept -> T = 0;
	};
	IGNORE_WEAK_VTABLES_STOP

	template<typename EngineType, utils::concepts::Numeric T = int>
	requires utils::concepts::Derived<EngineType, Engine>
	class UniformDistribution final : public Distribution<EngineType, T> {
	  public:
		constexpr UniformDistribution() noexcept requires
			utils::concepts::DefaultConstructible<EngineType>
			: m_engine(std::make_unique<EngineType>()) {
			m_engine->seed(m_engine->get_seed());
		}
		constexpr UniformDistribution(T min, T max) noexcept requires
			utils::concepts::DefaultConstructible<EngineType>
			: m_min(min), m_max(max), m_engine(std::make_unique<EngineType>()) {
			m_engine->seed(m_engine->get_seed());
		}
		constexpr UniformDistribution(T min, T max, std::unique_ptr<EngineType>&& engine) noexcept
			: m_min(min), m_max(max), m_engine(engine) {
			m_engine->seed(m_engine->get_seed());
		}
		explicit constexpr UniformDistribution(std::unique_ptr<EngineType>&& engine) noexcept
			: m_engine(std::move(engine)) {
			m_engine->seed(m_engine->get_seed());
		}
		template<typename... Args>
		requires utils::concepts::Derived<EngineType, Engine> && utils::concepts::
			ConstructibleFrom<EngineType, Args...>
		explicit constexpr UniformDistribution(Args&&... args) noexcept
			: m_engine(std::make_unique<EngineType>(args...)) {
			m_engine->seed(m_engine->get_seed());
		}
		template<typename... Args>
		requires utils::concepts::Derived<EngineType, Engine> && utils::concepts::
			ConstructibleFrom<EngineType, Args...>
		explicit constexpr UniformDistribution(T min, T max, Args&&... args) noexcept
			: m_min(min), m_max(max), m_engine(std::make_unique<EngineType>(args...)) {
			m_engine->seed(m_engine->get_seed());
		}
		constexpr UniformDistribution(const UniformDistribution& distribution) noexcept = default;
		constexpr UniformDistribution(UniformDistribution&& distribution) noexcept = default;

		constexpr ~UniformDistribution() noexcept final = default;

		inline constexpr auto normalized_random_value() noexcept -> double final {
			return narrow_cast<double>(m_engine->generate())
				   / narrow_cast<double>(m_engine->max_value());
		}

		inline constexpr auto random_value() noexcept -> T final {
			return narrow_cast<T>(this->normalized_random_value()) * (m_max - m_min) + m_min;
		}

		inline constexpr auto seed(size_t seed) noexcept -> void final {
			m_engine->seed(seed);
		}

		[[nodiscard]] inline constexpr auto get_seed() noexcept -> size_t final {
			return m_engine->get_seed();
		}

		inline constexpr auto set_min(T min) noexcept -> void final {
			m_min = min;
		}
		[[nodiscard]] inline constexpr auto get_min() const noexcept -> T final {
			return m_min;
		}

		inline constexpr auto set_max(T max) noexcept -> void final {
			m_max = max;
		}
		[[nodiscard]] inline constexpr auto get_max() const noexcept -> T final {
			return m_max;
		}

		constexpr auto operator=(const UniformDistribution& distribution) noexcept
			-> UniformDistribution& = default;
		constexpr auto
		operator=(UniformDistribution&& distribution) noexcept -> UniformDistribution& = default;

		inline constexpr auto operator()() noexcept -> T final {
			return random_value();
		}

	  private:
		T m_min = narrow_cast<T>(0);
		T m_max = narrow_cast<T>(1);
		std::unique_ptr<EngineType> m_engine;
	};

	template<utils::concepts::Numeric T = float>
	[[clang::no_destroy]] static math::UniformDistribution<math::LinearCongruentialEngine<>, T>
		GLOBAL_UNIFORM_DISTRIBUTION;

	template<utils::concepts::Numeric T = float>
	static std::atomic_bool GLOBAL_UNIFORM_DISTRIBUTION_INITIALIZED;

	template<utils::concepts::Numeric T = float>
	inline static constexpr auto initialize_global_uniform_distribution() noexcept -> void {
		bool initialized = false;
		if(GLOBAL_UNIFORM_DISTRIBUTION_INITIALIZED<T>.compare_exchange_strong(
			   initialized,
			   true,
			   std::memory_order_seq_cst))
		{
			GLOBAL_UNIFORM_DISTRIBUTION<
				T> = math::UniformDistribution<math::LinearCongruentialEngine<>, T>();
		}
	}

	template<utils::concepts::FloatingPoint T = float>
	inline auto random_value() noexcept -> T {
		initialize_global_uniform_distribution<T>();
		GLOBAL_UNIFORM_DISTRIBUTION<T>.set_min(narrow_cast<T>(0));
		GLOBAL_UNIFORM_DISTRIBUTION<T>.set_max(narrow_cast<T>(1));

		return GLOBAL_UNIFORM_DISTRIBUTION<T>();
	}

	template<utils::concepts::Numeric T = float>
	inline auto random_value(T min, T max) noexcept -> T {
		initialize_global_uniform_distribution<T>();

		GLOBAL_UNIFORM_DISTRIBUTION<T>.set_min(min);
		GLOBAL_UNIFORM_DISTRIBUTION<T>.set_max(max);

		return GLOBAL_UNIFORM_DISTRIBUTION<T>();
	}
} // namespace hyperion::math

ext: h


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

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

/* funcs.h */

#include <gsl/gsl_math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv2.h>

#define WS_SZ	1000  			/* size of the integration workspace */

#define POW_2	gsl_pow_2  		/* square */
#define VECTOR	gsl_vector_get		/* vector_get */

struct f_params { double omega_c ; double beta ; double Omega ;
	double omega_1; double alpha; } ;
int assign_p ( void*, double*, double*, double*, double* ) ;

struct ent_pars { const gsl_vector* rhoeq ; gsl_matrix* red_m ; } ;

/* Re_gsc functions */
double k_func ( double, void* ) ;
double k_func_1 ( double , void* ) ;
double k_func_2 ( double , void* ) ;
double k_func_3 ( double , void* ) ;
double k_func_4 ( double , void* ) ;
int first_int ( double*, double*, void* ) ;
int second_int ( double* , double* , void* ) ;
int third_int ( double* , double* , void* ) ;
int fourth_int ( double* , double* , void* ) ;
int re_gsc ( double* result, double* error, void* params ) ;
int re_gcs ( double* result, double* error, void* params ) ;

int re_gcc ( void* params, double* val ) ;
int im_gcc ( void* params, double* val, double* error ) ;
int re_gss ( void* params, double* val ) ;
int im_gss ( void* params, double* val, double* error ) ;
int im_gcs ( void* params, double* val ) ;
int im_gsc ( void* params, double* val ) ;
int re_gc0 ( void* params, double* val ) ;
int im_gc0 ( void* params, double* val, double* error ) ;
int re_gs0 ( void* params, double* val, double* error ) ;
int im_gs0 ( void* params, double* val ) ;

/* Img_c0, Img_cc, Img_ss, Img_cs functions */
int expi ( double , double*, double* ) ;
int expi_plus ( double, double*, double* ) ;
double fu ( double, void* ) ;
double ex ( double, void* ) ;

/* Reg_s0 functions */
double fu_inv ( double, void* ) ;
double fu_cau ( double, void* ) ;
double fu_coth ( double, void* ) ;

/* Functions to create matrices */
int red_mat ( gsl_matrix* , void* ) ;
int cp_mat ( gsl_matrix*, void* ) ;
int ham_gen ( gsl_matrix*, const double* ) ;
int mat_write ( gsl_matrix*, char* ) ;
int mat_read ( gsl_matrix*, char* ) ;
int integration ( void*, double* ) ;

/* Functions for evolution */
int generator ( double, const double*, double*, void* ) ;
int jac ( double, const double*, double*, double*, void* ) ;
int evol ( double, gsl_vector*, double, gsl_odeiv2_evolve*,
		gsl_odeiv2_control*, gsl_odeiv2_step*, gsl_odeiv2_system* ) ;
int red_evol ( void*, const double*, double, double, const gsl_vector*,	gsl_matrix* ) ;
int cp_evol ( void*, const double*, double, double, const gsl_vector*,	gsl_matrix* ) ;

/* Stationarity */
int stationary ( const gsl_matrix* , gsl_vector*  ) ;
double J ( double, double ) ;
double polarization ( void*, double ) ;

/* Entropy */
double entropy_of_state ( const gsl_vector* ) ;
double entropy_production ( const gsl_vector*, const gsl_vector*, const gsl_matrix* ) ;
double entropy_prod_time  ( double, void* ) ;
double entropy_integ ( double, const gsl_vector*, gsl_matrix* ) ;

/* DC current */
double current ( const gsl_vector* ) ;
double current_temp ( double ) ;

/* total current */
double tot_current ( const gsl_vector*, void* ) ;

/* writing operations */
int save_integrals ( void* ) ;
int save_matrices ( void* ) ;
int write_red_curr_T ( void* ) ;
int write_cp_curr_T ( void* ) ;
int write_red_curr_O ( void* ) ;
int write_cp_curr_O ( void* ) ;

/* currents */
double current_red_om ( double, void* ) ;
double current_red_T ( double, void* ) ;
double current_cp_om ( double, void* ) ;
double current_cp_T ( double, void* ) ;

/* positivity */
double r0_dot ( const gsl_matrix* , const gsl_vector*  ) ;
int bloch_vector ( gsl_vector* , double , double , double ) ;
int polars ( double* , double* , double* , const gsl_vector*  ) ;
int sample ( const gsl_matrix* , unsigned int ) ;

/* posit */
gsl_matrix_complex* state_matrix ( double, double ) ;
gsl_matrix_complex* rotate ( gsl_matrix_complex*, double ) ;
gsl_matrix_complex* D_psi ( gsl_matrix_complex* ) ;

/* pos_study */
double delta ( double, double*, double, double, double ) ;
int read_integs ( double* ) ;

/* entropy production with different parameters */
int ent_prod_fixed_t ( double, double ) ;
int ent_prod_red_evol ( void* , const double* , double , double , const gsl_vector* , gsl_matrix* , double ) ;
char* dtoascii ( double, int ) ;


text:
#ifndef TOPO_SL_H
#define TOPO_SL_H

#include <assert.h>

#include <vector>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_spmatrix.h>
#include <gsl/gsl_splinalg.h>

#include <apf.h>
#include <apfMDS.h>
#include <apfMesh.h>
#include <apfMesh2.h>

#include <apfNumbering.h>
#include <apfShape.h>

#include "topo_disp.h"

// Container for the gsl based linear solver:
class solv_kuprat {
  private:
    solv_kuprat(const solv_kuprat &obj);
    solv_kuprat& operator=( const solv_kuprat& obj );
  public:
    apf::Mesh2* m;
    field_calc* f_calc;

    int n;
    std::map<apf::MeshEntity*, int> v2id;
    std::map<apf::MeshEntity*, bool> v_map;
    std::map<int, apf::MeshEntity*> id2v;

    std::vector<apf::MeshEntity*>* verts;

    gsl_spmatrix* A;
    gsl_spmatrix* C;
    gsl_vector *f;
    gsl_vector *u;

    apf::GlobalNumbering* gn;

    // Functions:
    void numb_vert();
    void numb_verts_in();
    void dest_numb();
    void add_forces(apf::MeshEntity* tri);

    void collect_matrix_verts();
    void collect_matrix();
    void init(apf::Mesh2* m_in);
    void init_verts(apf::Mesh2* m_in);
    void clear();
    double solve(const double tol);

    solv_kuprat(apf::Mesh2* m, field_calc* f_calc_in);
    solv_kuprat(apf::Mesh2* m, field_calc* f_calc_in, std::vector<apf::MeshEntity*>* verts);

    ~solv_kuprat();

    friend std::ostream& operator<<(std::ostream& s, const solv_kuprat& sl);

};

// Container for the gsl based linear solver:
class solv_kuprat_NBC {
  private:
    solv_kuprat_NBC(const solv_kuprat_NBC &obj);
    solv_kuprat_NBC& operator=( const solv_kuprat_NBC& obj );
  public:
    apf::Mesh2* m;
    field_calc* f_calc;

    int n;
    std::map<apf::MeshEntity*, int> v2id;
    std::map<apf::MeshEntity*, bool> v_map;
    std::map<int, apf::MeshEntity*> id2v;

    std::map<apf::MeshEntity*, bool> v_sp;
    std::vector<apf::MeshEntity*>* verts;

    gsl_spmatrix* A;
    gsl_spmatrix* C;
    gsl_vector *f;
    gsl_vector *u;

    apf::GlobalNumbering* gn;

    // Functions:
    void numb_vert();
    void numb_verts_in();
    void dest_numb();
    void add_forces(apf::MeshEntity* tri);

    void collect_matrix_verts();
    void collect_matrix();

		void process_ext_shell();
		void process_ext_shell_verts();

    void init(apf::Mesh2* m_in);
    void init_verts(apf::Mesh2* m_in);
    void clear();
    double solve(const double tol);

    solv_kuprat_NBC(apf::Mesh2* m, field_calc* f_calc_in);
    solv_kuprat_NBC(apf::Mesh2* m, field_calc* f_calc_in, std::vector<apf::MeshEntity*>* verts);

    ~solv_kuprat_NBC();

    friend std::ostream& operator<<(std::ostream& s, const solv_kuprat_NBC& sl);

};
#endif

max_stars_repo_name: erdemeren/VDlib

max_issues_repo_name: erdemeren/VDlib

max_forks_repo_name: erdemeren/VDlib


text:
/* ----------------------------------------------------------------------------
  JAM_AXI_RMS_WMMT

    Calculates weighted second moment.

    INPUTS
      xp    : projected x' [pc]
      yp    : projected y' [pc]
      nxy   : number of x' and y' values given
      incl  : inclination [radians]
      lum   : projected luminous MGE
      pot   : projected potential MGE
      beta  : velocity anisotropy (1 - vz^2 / vr^2)
      vv    : velocity integral selector (1=xx, 2=yy, 3=zz, 4=xy, 5=xz, 6=yz)

    NOTES
    * Based on janis2_weighted_second_moment_squared IDL code by Michele
      Cappellari.

  Mark den Brok
  Laura L Watkins [lauralwatkins@gmail.com]

  This code is released under a BSD 2-clause license.
  If you use this code for your research, please cite:
  Watkins et al. 2013, MNRAS, 436, 2598
  "Discrete dynamical models of omega Centauri"
  http://adsabs.harvard.edu/abs/2013MNRAS.436.2598W
---------------------------------------------------------------------------- */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_errno.h>
#include "jam.h"
#include "../mge/mge.h"
#include "../tools/tools.h"

double *jam_axi_rms_wmmt( double *xp, double *yp, int nxy, double incl, \
        struct multigaussexp *lum, struct multigaussexp *pot, \
        double *beta, int vv , int *gslFlag_rms) {

    struct params_rmsint p;
    struct multigaussexp ilum, ipot;
    double ci, si, *kani, *s2l, *q2l, *s2q2l, *s2p, *e2p;
    double result, error, *sb_mu2;
    int i;
    int status = 0;

    // convert from projected MGEs to intrinsic MGEs
    ilum = mge_deproject( lum, incl );
    ipot = mge_deproject( pot, incl );

    // angles
    ci = cos( incl );
    si = sin( incl );

    // mge component combinations

    kani = (double *) malloc ( lum->ntotal * sizeof( double ) );
    s2l = (double *) malloc( lum->ntotal * sizeof( double ) );
    q2l = (double *) malloc( lum->ntotal * sizeof( double ) );
    s2q2l = (double *) malloc( lum->ntotal * sizeof( double ) );
    s2p = (double *) malloc( pot->ntotal * sizeof( double ) );
    e2p = (double *) malloc( pot->ntotal * sizeof( double ) );

    for ( i = 0; i < ilum.ntotal; i++ ) {
        kani[i] = 1. / ( 1. - beta[i] );
        s2l[i] = pow( ilum.sigma[i], 2. );
        q2l[i] = pow( ilum.q[i], 2. );
        s2q2l[i] = s2l[i] * q2l[i];
    }

    for ( i = 0; i < ipot.ntotal; i++ ) {
        s2p[i] = pow( ipot.sigma[i], 2. );
        e2p[i] = 1. - pow( ipot.q[i], 2. );
    }

    // parameters for the integrand function
    p.ci2 = ci * ci;
    p.si2 = si * si;
    p.cisi = ci * si;
    p.lum = &ilum;
    p.pot = &ipot;
    p.kani = kani;
    p.s2l = s2l;
    p.q2l = q2l;
    p.s2q2l = s2q2l;
    p.s2p = s2p;
    p.e2p = e2p;
    p.vv = vv;

    // perform integration

    gsl_integration_workspace *w = gsl_integration_workspace_alloc( 1000 );
    gsl_function F;
    gsl_set_error_handler_off();
    F.function = &jam_axi_rms_mgeint;

    sb_mu2 = (double *) malloc( nxy * sizeof( double ) );
    for ( i = 0; i < nxy; i++ ) {
        p.x2 = xp[i] * xp[i];
        p.y2 = yp[i] * yp[i];
        p.xy = xp[i] * yp[i];
        F.params = &p;
        status += gsl_integration_qag( &F, 0., 1., 0., 1e-5, 1000, 6, w, \
            &result, &error );
        sb_mu2[i] = result;
    }

    /*// if you want to test the gsl flagging thing uncomment these two lines
    int test = 1;
    status += test;
    */
    if (status > 0) { *gslFlag_rms += 1;   }

    gsl_integration_workspace_free( w );

    free( kani );
    free( s2l );
    free( q2l );
    free( s2q2l );
    free( s2p );
    free( e2p );

    return sb_mu2;

}

ext: c

max_stars_repo_name: tadejaversic/cjam

max_issues_repo_name: tadejaversic/cjam

max_forks_repo_name: tadejaversic/cjam


text:
#ifndef MODULE_WASM_TABLE_H
#define MODULE_WASM_TABLE_H

#include "function/TableFunction.h"
#include <gsl/span>

namespace wasm {

struct WasmTable {

	using wasm_external_kind_type = parse::Table;
	using any_function_type = std::optional<

	gsl::span<any_function_type> get_segment(std::size_t offset, std::size_t length)
	{
		assert(offset <= table_.size());
		assert((table_.size() - offset) >= lenfth);
		return gsl::span<TableFunction>(table_.data() + offset, length);
	}

	friend bool matches(const WasmTable& self, const parse::Table& tp)
	{
		return self.table_.size() >= tp.size() and (
			self.maximum_ == tp.maximum.value_or(std::numeric_limits<std::size_t>::max())
		);
	}

	const any_function_type& at(wasm_uint32_t idx) const
	{ return table_.at(idx); }

	any_function_type& at(wasm_uint32_t idx)
	{ return table_.at(idx); }

private:
	std::vector<table_function_type> table_;
	const std::size_t maximum_ = std::numeric_limits<std::size_t>::max();
};

} /* namespace wasm */

#endif /* MODULE_WASM_TABLE_H */

max_stars_repo_name: tvanslyke/wasm-cpp

max_issues_repo_name: tvanslyke/wasm-cpp

max_forks_repo_name: tvanslyke/wasm-cpp


text:
/**
 * @file objective_measure.h
 * @brief objective measure for speech quality or distance
 * @author Kenichi Kumatani
 */

#ifndef OBJECTIVE_MEASURE_H
#define OBJECTIVE_MEASURE_H

#include <stdio.h>
#include <assert.h>
#include <float.h>

#include <gsl/gsl_block.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_fft_real.h>
#include <gsl/gsl_fft_complex.h>
#include <common/refcount.h>
#include "common/jexception.h"

#include "stream/stream.h"
#include "feature/feature.h"
#include "modulated/modulated.h"

// ----- definition for class `SNR' -----
//

 /**
     @class calculate the signal-to-ratio (SNR) of two speech data
     @usage
   */
class SNR{
 public:
  //SNR();
  //SNR( const String& nm = "SNR" );
  float getSNR( const String& fn1, const String& fn2, int normalizationOption, int chX=1, int samplerate=16000, int cfrom=-1, int to=-1 );
  float getSNR2( gsl_vector_float *original, gsl_vector_float *enhanced, int normalizationOption);
};

typedef refcount_ptr<SNR> SNRPtr;

class segmentalSNR{
 public:
};

class ItakuraSaitoMeasurePS {
 public:
 ItakuraSaitoMeasurePS( unsigned fftLen,  unsigned r = 1, unsigned windowType = 1, const String& nm = "ItakuraSaitoMeasurePS" ):_fftLen(fftLen),_r(r),_windowType(windowType)
  {
    _D = _fftLen / (int)pow( (float)2, (int)_r );
  }

  float getDistance( const String& fn1, const String& fn2, int chX=1, int samplerate=16000, int bframe = 0, int eframe = -1 );
  //float getDistance( VectorFloatFeatureStreamPtr& original,
  //VectorFloatFeatureStreamPtr& enhanced );

  int frameShiftLength(){
    return _D;
  }

 private:
  unsigned _fftLen;
  unsigned _r;
  unsigned _windowType;
  unsigned _D;
};

typedef refcount_ptr<ItakuraSaitoMeasurePS> ItakuraSaitoMeasurePSPtr;

#endif


text:
#include <pygsl/error_helpers.h>
#include <pygsl/utils.h>
#include <gsl/gsl_errno.h>
#include <compile.h>
#include <frameobject.h>

enum handleflag {
     HANDLE_ERROR = 0,
     HANDLE_WARNING
};

static int
PyGSL_internal_error_handler(const char *reason, /* name of function*/
			     const char *file, /*from CPP*/
			     int line,   /*from CPP*/
			     int gsl_error,
			     enum handleflag flag);

static int
PyGSL_error_flag(long flag)
{
     FUNC_MESS_BEGIN();
     if(PyGSL_DEBUG_LEVEL() > 2){
	  fprintf(stderr,"I got an Error %ld\n", flag);
     }
     if(PyErr_Occurred()){
	     DEBUG_MESS(3, "Already a python error registered for flag %ld", flag);
	     return GSL_FAILURE;
     }
     if(flag>0){
	  /*
	   * How can I end here without an Python error?
	   *
	   * 25. October 2008
	   * Well, very simply when the GSL_ERROR_HANDLER is set to off.
	   *
	   * All GSL modules are
	   * supposed to call GSL_ERROR which should call the default error
	   * handler.
	   *
	   * 25. October 2008
	   * No not when the handler is set to off, which is necessary in
	   * a threaded python. That's why it was a bad idea to call the
	   * gsl_error_handler here!
	   */
	  /*
	   * gsl_error("Unknown Reason. It was not set by GSL.",  __FILE__,
	   *	    __LINE__, flag);
	  */
	  PyGSL_internal_error_handler("Unknown Reason. It was not set by GSL",  __FILE__,
				       __LINE__, flag, HANDLE_ERROR);
	  /*
	   * So lets keep the flag to return ... who knows what it will be used for...
	   * return GSL_FAILURE;
	   */
	  return flag;
     }
     FUNC_MESS_END();
     return GSL_SUCCESS;
}

static PyObject *
PyGSL_error_flag_to_pyint(long flag)
{
     PyObject * result = NULL;
     FUNC_MESS_BEGIN();
     if(GSL_FAILURE == PyGSL_error_flag(flag)){
	  return NULL;
     }
     result = PyInt_FromLong((long) flag);
     FUNC_MESS_END();

     return result;
}

static void
PyGSL_add_traceback(PyObject *module, const char *filename, const char *funcname, int lineno)
{
     PyObject *py_srcfile = NULL, *py_funcname = NULL, *py_globals = NULL,
	  *empty_tuple = NULL,  *empty_string = NULL;
     PyCodeObject *py_code = NULL;
     PyFrameObject *py_frame = NULL;

     FUNC_MESS_BEGIN();

     if(filename == NULL)
	  filename = "file ???";
     py_srcfile = PyString_FromString(filename);
     if (py_srcfile == NULL)
	  goto fail;

     if(funcname == NULL)
	  funcname = "function ???";

     py_funcname = PyString_FromString(funcname);
     if (py_funcname == NULL)
	  goto fail;

     /* Use the module if provided */
     if(module == NULL){
	  py_globals = PyDict_New();
     } else {
	  py_globals = PyModule_GetDict(module);
     }
     if (py_globals == NULL)
	  goto fail;

     empty_tuple = PyTuple_New(0);
     if (empty_tuple == NULL)
	  goto fail;

     empty_string = PyString_FromString("");
     if (empty_string == NULL)
	  goto fail;

     py_code = PyCode_New(
	  0,            /*int argcount,*/
	  0,            /*int nlocals,*/
	  0,            /*int stacksize,*/
	  0,            /*int flags,*/
	  empty_string, /*PyObject *code,*/
	  empty_tuple,  /*PyObject *consts,*/
	  empty_tuple,  /*PyObject *names,*/
	  empty_tuple,  /*PyObject *varnames,*/
	  empty_tuple,  /*PyObject *freevars,*/
	  empty_tuple,  /*PyObject *cellvars,*/
	  py_srcfile,   /*PyObject *filename,*/
	  py_funcname,  /*PyObject *name,*/
	  lineno,       /*int firstlineno,*/
	  empty_string  /*PyObject *lnotab*/
	  );
     if (py_code == NULL)
	  goto fail;

     py_frame = PyFrame_New(
	  PyThreadState_Get(), /*PyThreadState *tstate,*/
	  py_code,             /*PyCodeObject *code,*/
	  py_globals,          /*PyObject *globals,*/
	  0                    /*PyObject *locals*/
	  );

     if (py_frame == NULL)
	  goto fail;
     py_frame->f_lineno = lineno;
     PyTraceBack_Here(py_frame);

     FUNC_MESS_END();
     return;

 fail:
     FUNC_MESS("Handling failure");
     Py_XDECREF(py_srcfile);
     Py_XDECREF(py_funcname);
     Py_XDECREF(empty_tuple);
     Py_XDECREF(empty_string);
     Py_XDECREF(py_code);
     Py_XDECREF(py_frame);
}

#define PyGSL_ERRNO_MAX 32
static PyObject * errno_accel[PyGSL_ERRNO_MAX];
static PyObject * error_dict = NULL;
static PyObject * warning_dict = NULL;
static PyObject * unknown_error = NULL;

static void
PyGSL_print_accel_object(void)
{
     int i;
     FUNC_MESS_BEGIN();
     for(i = 0; i<PyGSL_ERRNO_MAX; ++i){
	  DEBUG_MESS(4, "errno_accel[%d] = %p", i, (void*)(errno_accel[i]));
     }
     FUNC_MESS_END();
}

static int
PyGSL_register_accel_err_object(PyObject * err_ob, long test_errno)
{
     PyObject *tmp;

     FUNC_MESS_BEGIN();
     assert(err_ob);
     tmp = errno_accel[test_errno];
     if(tmp != NULL){
	  PyErr_Format(PyExc_ValueError,
		       "In errno_accel: errno %ld already occupied with object %p!\n",
		       test_errno, (void *) tmp);
	  return -2;
     }
     Py_INCREF(err_ob);
     errno_accel[test_errno] = err_ob;
     FUNC_MESS_END();
     return 0;
}

/* register an error object to the Python dictionary */
static int
_PyGSL_register_err_object(PyObject *dict, PyObject * err_ob, PyObject *the_errno)
{
     PyObject *test;

     FUNC_MESS_BEGIN();
     assert(error_dict);
     test = PyDict_GetItem(dict, the_errno);
     if(test != NULL){
	  PyErr_Format(PyExc_ValueError,
		       "In dict %p: key %p already occupied with object %p!\n",
		       dict, the_errno, (void *) test);
	  return -2;
     }
     Py_INCREF(err_ob);
     PyDict_SetItem(dict, the_errno, err_ob);
     FUNC_MESS_END();
     return 0;
}

/* register an error object */
static int
_PyGSL_register_error(PyObject *dict, int errno_max, PyObject * err_ob)
{
     PyObject *tmp, *name;
     long test_errno;
     int flag;
     char * c_name;

     FUNC_MESS_BEGIN();
     assert(err_ob);
     tmp = PyObject_GetAttrString(err_ob, "errno");
     if(tmp == NULL){
	  name = PyObject_GetAttrString(err_ob, "__name__");

	  if(name == NULL)
	       c_name = "unknown name";
	  else if (!PyString_Check(name))
	       c_name = "name not str object!";
	  else
	       c_name = PyString_AsString (name);

	  fprintf(stderr, "failed to get errno from err_ob '%s' @ %p\n",
		  c_name, (void *) err_ob);
	  PyErr_Format(PyExc_AttributeError,
		       "err_ob '%s' @ %p missed attribue 'errno!'\n", c_name,
		       err_ob);
	  return -1;
     }

     if(!PyInt_CheckExact(tmp)){
	  fprintf(stderr, "errno %p from err_ob %p was not an exact int!\n",
		  (void *) tmp, (void *) err_ob);
	  PyErr_Format(PyExc_TypeError, "errno %p from err_ob %p was not an exact int!\n",
		       (void *) tmp, (void *) err_ob);
	  return -1;
     }

     test_errno = PyInt_AsLong(tmp);
     if((dict == error_dict) && (test_errno < PyGSL_ERRNO_MAX)){
	  flag = PyGSL_register_accel_err_object(err_ob, test_errno);
     }else{
	  flag = _PyGSL_register_err_object(dict, err_ob, tmp);
     }
     if(flag != 0)
	  fprintf(stderr, "Failed to register err_ob %p with errno %ld.\n"
		  "\tAlready registered?\n", err_ob, test_errno);
     FUNC_MESS_END();
     return flag;
}

static PyObject*
PyGSL_register_error_objs(PyObject *self, PyObject *args, PyObject *dict, int errno_max)
{
     int flag, i, len;
     PyObject *tmp;

     FUNC_MESS_BEGIN();
     if(!PySequence_Check(args))
	  return NULL;

     len = PySequence_Size(args);
     DEBUG_MESS(5, "Recieved %d error objects", len);
     for(i = 0; i < len; ++i){
	  tmp = PySequence_GetItem(args, i);
	  flag = _PyGSL_register_error(dict, errno_max, tmp);
	  if(flag != 0){
	       fprintf(stderr, "Failed to register error object %d\n", i);
	       return NULL;
	  }
     }
     PyGSL_print_accel_object();

     Py_INCREF(Py_None);
     FUNC_MESS_END();
     return Py_None;
}

static PyObject*
PyGSL_register_warnings(PyObject *self, PyObject *args)
{
     PyObject *tmp;
     FUNC_MESS_BEGIN();
     tmp = PyGSL_register_error_objs(self, args, warning_dict, 0);
     FUNC_MESS_END();
     return tmp;
}

static PyObject*
PyGSL_register_exceptions(PyObject *self, PyObject *args)
{
     PyObject *tmp;
     FUNC_MESS_BEGIN();
     tmp = PyGSL_register_error_objs(self, args, error_dict, PyGSL_ERRNO_MAX);
     FUNC_MESS_END();
     return tmp;
}

static PyObject *
PyGSL_get_error_object(int the_errno, PyObject ** accel, int accel_max, PyObject *dict)
{
     PyObject *tmp;

     FUNC_MESS_BEGIN();
     assert(the_errno >= 0);
     if (the_errno < accel_max){
	  DEBUG_MESS(4, "Trying to get an error object from accel array at %p",
		     (void*) accel);
	  tmp = accel[the_errno];
     }else{
	  DEBUG_MESS(4, "Trying to get an error object from dictonary at %p",
		     (void*) dict);
	  tmp =  PyDict_GetItem(dict, PyInt_FromLong(the_errno));
     }
     if(tmp == NULL){
	  DEBUG_MESS(3, "Could not find an error object for errno %d", the_errno);
	  PyGSL_print_accel_object();
	  return unknown_error;
     }
     FUNC_MESS_END();
     return tmp;
}

static int
PyGSL_init_errno(void)
{
     int i;
     FUNC_MESS_BEGIN();

     for(i = 0; i< PyGSL_ERRNO_MAX; ++i){
	  DEBUG_MESS(3, "setting errno_accel[%d] to NULL; was %p",
		     i, (void*) (errno_accel[i]));
	  errno_accel[i] = NULL;
     }
     error_dict = PyDict_New();
     if (error_dict == NULL)
	  return -1;

     warning_dict = PyDict_New();
     if (warning_dict == NULL)
	  return -1;

     unknown_error = PyExc_ValueError;
     FUNC_MESS_END();
     return 0;
}

/*
 * Warnings return a flag, so one can see if the warning raises an exception
 *  or not.
 */
static int
PyGSL_internal_error_handler(const char *reason, /* name of function*/
			     const char *file, /*from CPP*/
			     int line,   /*from CPP*/
			     int gsl_error,
			     enum handleflag flag)
{
  const char* error_explanation;
  char error_text[255];
  PyObject* gsl_error_object;

  FUNC_MESS_BEGIN();
  /*
   * GSL_ENOMEM is special. I am out of memory. No fancy tricks here.
   */
  if (GSL_ENOMEM == gsl_error){
       PyErr_NoMemory();
       return -1;
  }

  /*
   * some functions call error handler more than once before returning
   *  report only the first (most specific) error
   */
  if (line < 0) line = 0;
  /* test, if exception is already set */
  DEBUG_MESS(5, "Checking if python error occured, gsl error %d, line %d", gsl_error, line);
  if (PyErr_Occurred()) {
       if(PyGSL_DEBUG_LEVEL() > 0)
	    fprintf(stderr, "Another error occured: %s\n",error_text);
       FUNC_MESS("Already set python error found");
       return -1;
  }

  /*
   * Find the approbriate error
   */
  error_explanation = gsl_strerror(gsl_error);
  if (reason==NULL){
       reason = "no reason given!";
  }

  if (error_explanation==NULL){
      snprintf(error_text,sizeof(error_text),
	       "unknown error %d: %s",
	       gsl_error, reason);
  }else{
      snprintf(error_text,sizeof(error_text),
	       "%s: %s",
	       error_explanation,reason);
  }

  switch(flag){
  case HANDLE_ERROR:
       assert(gsl_error > 0);
       gsl_error_object = PyGSL_get_error_object(gsl_error, errno_accel, PyGSL_ERRNO_MAX, error_dict);
       Py_INCREF(gsl_error_object);
       PyErr_SetObject(gsl_error_object, PyString_FromString(error_text));
       FUNC_MESS("Set Python error object");
       return -1;
       break;
  case HANDLE_WARNING:
       assert(gsl_error > 0);
       gsl_error_object = PyGSL_get_error_object(gsl_error, NULL, 0, warning_dict);
       Py_INCREF(gsl_error_object);
       FUNC_MESS("Returning python warning");
       return PyErr_Warn(gsl_error_object, error_text);
       break;
  default:
       fprintf(stderr, "Unknown handle %d\n", flag);
  }
  FUNC_MESS("Should not end here!");
  return -1;
}
/*
 * sets the right exception, but does not return to python!
 */
static void
PyGSL_module_error_handler(const char *reason, /* name of function*/
			   const char *file, /*from CPP*/
			   int line,   /*from CPP*/
			   int gsl_error) /* real "reason" */
{
     FUNC_MESS_BEGIN();
     PyGSL_internal_error_handler(reason, file, line,  gsl_error, HANDLE_ERROR);
     FUNC_MESS_END();
}

static int
PyGSL_warning(const char *reason, /* name of function*/
	      const char *file, /*from CPP*/
	      int line,   /*from CPP*/
	      int gsl_error) /* real "reason" */
{
     int tmp;
     FUNC_MESS_BEGIN();
     tmp =  PyGSL_internal_error_handler(reason, file, line,  gsl_error, HANDLE_WARNING);
     FUNC_MESS_END();
     return tmp;
}

max_stars_repo_name: juhnowski/FishingRod

max_issues_repo_name: juhnowski/FishingRod

max_forks_repo_name: juhnowski/FishingRod


text:
/* multfit/lmder.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 <config.h>

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

#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_multifit_nlin.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_permutation.h>

typedef struct
  {
    size_t iter;
    double xnorm;
    double fnorm;
    double delta;
    double par;
    gsl_matrix *J;             /* Jacobian matrix */
    gsl_matrix *r;             /* R matrix in J = Q R P^T */
    gsl_vector *tau;
    gsl_vector *diag;          /* scaling matrix D = diag(d1,...,dp) */
    gsl_vector *qtf;           /* Q^T f */
    gsl_vector *newton;
    gsl_vector *gradient;      /* gradient g = J^T f */
    gsl_vector *x_trial;       /* trial step x + dx */
    gsl_vector *f_trial;       /* trial function f(x + dx) */
    gsl_vector *df;
    gsl_vector *sdiag;
    gsl_vector *rptdx;
    const gsl_vector *weights; /* data weights */
    gsl_vector *w;
    gsl_vector *work1;
    gsl_permutation * perm;
  }
lmder_state_t;

static int lmder_alloc (void *vstate, size_t n, size_t p);
static int lmder_set (void *vstate, const gsl_vector * swts, gsl_multifit_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_vector * dx);
static int lmsder_set (void *vstate, const gsl_vector * swts, gsl_multifit_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_vector * dx);
static int set (void *vstate, const gsl_vector * swts, gsl_multifit_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_vector * dx, int scale);
static int lmder_iterate (void *vstate, const gsl_vector * swts, gsl_multifit_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_vector * dx);
static void lmder_free (void *vstate);
static int iterate (void *vstate, const gsl_vector * swts, gsl_multifit_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_vector * dx, int scale);

#include "lmutil.c"
#include "lmpar.c"
#include "lmset.c"
#include "lmiterate.c"

static int
lmder_alloc (void *vstate, size_t n, size_t p)
{
  lmder_state_t *state = (lmder_state_t *) vstate;
  gsl_matrix *r, *J;
  gsl_vector *tau, *diag, *qtf, *newton, *gradient, *x_trial, *f_trial,
   *df, *sdiag, *rptdx, *w, *work1;
  gsl_permutation *perm;

  J = gsl_matrix_alloc (n, p);

  if (J == 0)
    {
      GSL_ERROR ("failed to allocate space for J", GSL_ENOMEM);
    }

  state->J = J;

  r = gsl_matrix_alloc (n, p);

  if (r == 0)
    {
      GSL_ERROR ("failed to allocate space for r", GSL_ENOMEM);
    }

  state->r = r;

  tau = gsl_vector_calloc (GSL_MIN(n, p));

  if (tau == 0)
    {
      gsl_matrix_free (r);

      GSL_ERROR ("failed to allocate space for tau", GSL_ENOMEM);
    }

  state->tau = tau;

  diag = gsl_vector_calloc (p);

  if (diag == 0)
    {
      gsl_matrix_free (r);
      gsl_vector_free (tau);

      GSL_ERROR ("failed to allocate space for diag", GSL_ENOMEM);
    }

  state->diag = diag;

  qtf = gsl_vector_calloc (n);

  if (qtf == 0)
    {
      gsl_matrix_free (r);
      gsl_vector_free (tau);
      gsl_vector_free (diag);

      GSL_ERROR ("failed to allocate space for qtf", GSL_ENOMEM);
    }

  state->qtf = qtf;

  newton = gsl_vector_calloc (p);

  if (newton == 0)
    {
      gsl_matrix_free (r);
      gsl_vector_free (tau);
      gsl_vector_free (diag);
      gsl_vector_free (qtf);

      GSL_ERROR ("failed to allocate space for newton", GSL_ENOMEM);
    }

  state->newton = newton;

  gradient = gsl_vector_calloc (p);

  if (gradient == 0)
    {
      gsl_matrix_free (r);
      gsl_vector_free (tau);
      gsl_vector_free (diag);
      gsl_vector_free (qtf);
      gsl_vector_free (newton);

      GSL_ERROR ("failed to allocate space for gradient", GSL_ENOMEM);
    }

  state->gradient = gradient;

  x_trial = gsl_vector_calloc (p);

  if (x_trial == 0)
    {
      gsl_matrix_free (r);
      gsl_vector_free (tau);
      gsl_vector_free (diag);
      gsl_vector_free (qtf);
      gsl_vector_free (newton);
      gsl_vector_free (gradient);

      GSL_ERROR ("failed to allocate space for x_trial", GSL_ENOMEM);
    }

  state->x_trial = x_trial;

  f_trial = gsl_vector_calloc (n);

  if (f_trial == 0)
    {
      gsl_matrix_free (r);
      gsl_vector_free (tau);
      gsl_vector_free (diag);
      gsl_vector_free (qtf);
      gsl_vector_free (newton);
      gsl_vector_free (gradient);
      gsl_vector_free (x_trial);

      GSL_ERROR ("failed to allocate space for f_trial", GSL_ENOMEM);
    }

  state->f_trial = f_trial;

  df = gsl_vector_calloc (n);

  if (df == 0)
    {
      gsl_matrix_free (r);
      gsl_vector_free (tau);
      gsl_vector_free (diag);
      gsl_vector_free (qtf);
      gsl_vector_free (newton);
      gsl_vector_free (gradient);
      gsl_vector_free (x_trial);
      gsl_vector_free (f_trial);

      GSL_ERROR ("failed to allocate space for df", GSL_ENOMEM);
    }

  state->df = df;

  sdiag = gsl_vector_calloc (p);

  if (sdiag == 0)
    {
      gsl_matrix_free (r);
      gsl_vector_free (tau);
      gsl_vector_free (diag);
      gsl_vector_free (qtf);
      gsl_vector_free (newton);
      gsl_vector_free (gradient);
      gsl_vector_free (x_trial);
      gsl_vector_free (f_trial);
      gsl_vector_free (df);

      GSL_ERROR ("failed to allocate space for sdiag", GSL_ENOMEM);
    }

  state->sdiag = sdiag;

  rptdx = gsl_vector_calloc (n);

  if (rptdx == 0)
    {
      gsl_matrix_free (r);
      gsl_vector_free (tau);
      gsl_vector_free (diag);
      gsl_vector_free (qtf);
      gsl_vector_free (newton);
      gsl_vector_free (gradient);
      gsl_vector_free (x_trial);
      gsl_vector_free (f_trial);
      gsl_vector_free (df);
      gsl_vector_free (sdiag);

      GSL_ERROR ("failed to allocate space for rptdx", GSL_ENOMEM);
    }

  state->rptdx = rptdx;

  w = gsl_vector_calloc (n);

  if (w == 0)
    {
      gsl_matrix_free (r);
      gsl_vector_free (tau);
      gsl_vector_free (diag);
      gsl_vector_free (qtf);
      gsl_vector_free (newton);
      gsl_vector_free (gradient);
      gsl_vector_free (x_trial);
      gsl_vector_free (f_trial);
      gsl_vector_free (df);
      gsl_vector_free (sdiag);
      gsl_vector_free (rptdx);

      GSL_ERROR ("failed to allocate space for w", GSL_ENOMEM);
    }

  state->w = w;

  work1 = gsl_vector_calloc (p);

  if (work1 == 0)
    {
      gsl_matrix_free (r);
      gsl_vector_free (tau);
      gsl_vector_free (diag);
      gsl_vector_free (qtf);
      gsl_vector_free (newton);
      gsl_vector_free (gradient);
      gsl_vector_free (x_trial);
      gsl_vector_free (f_trial);
      gsl_vector_free (df);
      gsl_vector_free (sdiag);
      gsl_vector_free (rptdx);
      gsl_vector_free (w);

      GSL_ERROR ("failed to allocate space for work1", GSL_ENOMEM);
    }

  state->work1 = work1;

  perm = gsl_permutation_calloc (p);

  if (perm == 0)
    {
      gsl_matrix_free (r);
      gsl_vector_free (tau);
      gsl_vector_free (diag);
      gsl_vector_free (qtf);
      gsl_vector_free (newton);
      gsl_vector_free (gradient);
      gsl_vector_free (x_trial);
      gsl_vector_free (f_trial);
      gsl_vector_free (df);
      gsl_vector_free (sdiag);
      gsl_vector_free (rptdx);
      gsl_vector_free (w);
      gsl_vector_free (work1);

      GSL_ERROR ("failed to allocate space for perm", GSL_ENOMEM);
    }

  state->perm = perm;

  return GSL_SUCCESS;
}

static int
lmder_set (void *vstate, const gsl_vector * swts,
           gsl_multifit_function_fdf * fdf, gsl_vector * x,
           gsl_vector * f, gsl_vector * dx)
{
  int status = set (vstate, swts, fdf, x, f, dx, 0);
  return status ;
}

static int
lmsder_set (void *vstate, const gsl_vector * swts,
            gsl_multifit_function_fdf * fdf, gsl_vector * x,
            gsl_vector * f, gsl_vector * dx)
{
  int status = set (vstate, swts, fdf, x, f, dx, 1);
  return status ;
}

static int
lmder_iterate (void *vstate, const gsl_vector * swts,
               gsl_multifit_function_fdf * fdf, gsl_vector * x,
               gsl_vector * f, gsl_vector * dx)
{
  int status = iterate (vstate, swts, fdf, x, f, dx, 0);
  return status;
}

static int
lmsder_iterate (void *vstate, const gsl_vector * swts,
                gsl_multifit_function_fdf * fdf, gsl_vector * x,
                gsl_vector * f, gsl_vector * dx)
{
  int status = iterate (vstate, swts, fdf, x, f, dx, 1);
  return status;
}

static int
lmder_gradient (void *vstate, gsl_vector * g)
{
  lmder_state_t *state = (lmder_state_t *) vstate;
  compute_gradient(state->r, state->qtf, g);
  return GSL_SUCCESS;
}

static int
lmder_jac (void *vstate, gsl_matrix * J)
{
  lmder_state_t *state = (lmder_state_t *) vstate;
  int s = gsl_matrix_memcpy(J, state->J);

  return s;
}

static void
lmder_free (void *vstate)
{
  lmder_state_t *state = (lmder_state_t *) vstate;

  if (state->perm)
    gsl_permutation_free (state->perm);

  if (state->work1)
    gsl_vector_free (state->work1);

  if (state->w)
    gsl_vector_free (state->w);

  if (state->rptdx)
    gsl_vector_free (state->rptdx);

  if (state->sdiag)
    gsl_vector_free (state->sdiag);

  if (state->df)
    gsl_vector_free (state->df);

  if (state->f_trial)
    gsl_vector_free (state->f_trial);

  if (state->x_trial)
    gsl_vector_free (state->x_trial);

  if (state->gradient)
    gsl_vector_free (state->gradient);

  if (state->newton)
    gsl_vector_free (state->newton);

  if (state->qtf)
    gsl_vector_free (state->qtf);

  if (state->diag)
    gsl_vector_free (state->diag);

  if (state->tau)
    gsl_vector_free (state->tau);

  if (state->r)
    gsl_matrix_free (state->r);

  if (state->J)
    gsl_matrix_free (state->J);
}

static const gsl_multifit_fdfsolver_type lmder_type =
{
  "lmder",                      /* name */
  sizeof (lmder_state_t),
  &lmder_alloc,
  &lmder_set,
  &lmder_iterate,
  &lmder_gradient,
  &lmder_jac,
  &lmder_free
};

static const gsl_multifit_fdfsolver_type lmsder_type =
{
  "lmsder",                     /* name */
  sizeof (lmder_state_t),
  &lmder_alloc,
  &lmsder_set,
  &lmsder_iterate,
  &lmder_gradient,
  &lmder_jac,
  &lmder_free
};

const gsl_multifit_fdfsolver_type *gsl_multifit_fdfsolver_lmder = &lmder_type;
const gsl_multifit_fdfsolver_type *gsl_multifit_fdfsolver_lmsder = &lmsder_type;

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#ifndef STDAFX_H
#define STDAFX_H
/*
 * Libs
 */
#include <qwt.h>
#include <qwt_legend.h>
#include <qwt_plot.h>
#include <qwt_plot_curve.h>
#include <qwt_plot_grid.h>
#include <qwt_plot_marker.h>
#include <qwt_plot_renderer.h>
#include <qwt_symbol.h>
#include <QtConcurrent/QtConcurrent>
#include <QtCore/QtCore>
#include <QtWidgets/QtWidgets>
#include <gsl/gsl>
#include <opencv2/core/core.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/imgproc/imgproc.hpp>
/*
 * STL
 */
#include <memory>
#include <vector>
/*
 * UTILITIES
 */
#include "macros.inl"
#endif  // !STDAFX_H


text:
#ifndef _KORALI_BASEPARAMETER_H_
#define _KORALI_BASEPARAMETER_H_

#include <string>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_sf.h>

namespace Korali::Parameter
{

class Base
{
 public:
 Base(std::string name);
 gsl_rng* _range;

 std::string _name;
 double _lowerBound;
 double _upperBound;

 void setBounds(double lowerBound, double upperBound) { _lowerBound = lowerBound;	_upperBound = upperBound; }
 void setName(std::string name)    { _name = name; }
 void initializeDistribution(int seed);
 void checkBounds();

 virtual double getDensity(double x) = 0;
 virtual double getDensityLog(double x) = 0;
 virtual double getRandomNumber() = 0;
};

} // namespace Korali

#endif // _KORALI_BASEPARAMETER_H_

max_stars_repo_name: anianruoss/HPCSE-II

max_issues_repo_name: valentinjacot/backupETHZ

max_forks_repo_name: valentinjacot/backupETHZ


text:
#pragma once
#include <gsl/gsl_rng.h>
#include <utility>

// we are using GSL random number generation because i don't trust
// random number generation to be consistent across various C++ stdlib
// implementations, and for the kind of MC simulator we are writing here,
// we want to be able to run it deterministically for testing purposes.

class Sampler {
private:
    gsl_rng *internal_rng_state;
public:
    unsigned long int seed;
    double generate() {
        return gsl_rng_uniform_pos(internal_rng_state);
    };

    Sampler(unsigned long int n) :
        seed (n) {
        internal_rng_state = gsl_rng_alloc(gsl_rng_default);
        gsl_rng_set(internal_rng_state, seed);
        };

    ~Sampler() {
        gsl_rng_free(internal_rng_state);
    };

    // since we don't have access to gsl internal state, can't write
    // copy constructor
    Sampler(Sampler &other) = delete;

    // move constructor
    Sampler(Sampler &&other) :
        internal_rng_state (std::exchange(other.internal_rng_state, nullptr)),
        seed (other.seed)
        {};

    // since we don't have access to gsl internal state, can't write
    // copy assignment operator
    Sampler &operator=(Sampler &other) = delete;

    // move assignment operator
    Sampler &operator=(Sampler &&other) {
        seed = other.seed;

        // we move the existing internal_rng_state into other so
        // it gets freed when other is dropped.
        std::swap(internal_rng_state, other.internal_rng_state);
        return *this;
    };

};

max_stars_repo_name: danielbarter/SMMC

max_issues_repo_name: danielbarter/SMMC

max_forks_repo_name: danielbarter/SMMC


text:
#ifndef KEYPOINTS_H_CFNIGHAG
#define KEYPOINTS_H_CFNIGHAG

#include <boost/histogram.hpp>
#include <gsl/gsl>
#include <opencv2/core/persistence.hpp>
#include <opencv2/core/types.hpp>
#include <sens_loc/analysis/distance.h>

namespace sens_loc::analysis {

/// Analyze the distribution and characteristics of detected keypoints.
class keypoints {
  public:
    /// Use a regular grid with no transformation, no overflow and an axis
    /// title.
    using axis_t = distance::axis_t;
    /// Histogram type for the 2D distribution of keypoints in the dataset.
    using distribution_histo_t =
        decltype(boost::histogram::make_histogram(axis_t{}, axis_t{}));
    /// Histogram type for 'size' and 'response' of the keypoints.
    using histo_t = decltype(boost::histogram::make_histogram(axis_t{}));

    keypoints() = default;
    keypoints(unsigned int img_width, unsigned int img_height)
        : _img_width{img_width}
        , _img_height{img_height} {}

    void configure_image_dimension(unsigned int img_width,
                                   unsigned int img_height) noexcept {
        _img_width  = img_width;
        _img_height = img_height;
    }

    void configure_size(unsigned int bins) noexcept { _size_bins = bins; }
    void configure_size(std::string axis_title) noexcept {
        _size_title = std::move(axis_title);
    }
    void configure_size(unsigned int bins, std::string axis_title) noexcept {
        _size_bins  = bins;
        _size_title = std::move(axis_title);
    }
    void enable_size_histo(bool enabled) noexcept {
        _size_histo_enabled = enabled;
    }

    void configure_response(unsigned int bins) noexcept {
        _response_bins = bins;
    }
    void configure_response(std::string axis_title) noexcept {
        _response_title = std::move(axis_title);
    }
    void configure_response(unsigned int bins,
                            std::string  axis_title) noexcept {
        _response_bins  = bins;
        _response_title = std::move(axis_title);
    }
    void enable_response_histo(bool enabled) noexcept {
        _response_histo_enabled = enabled;
    }

    void configure_distribution(unsigned int both_bins) noexcept {
        _dist_width_bins  = both_bins;
        _dist_height_bins = both_bins;
    }
    void configure_distribution(unsigned int w_bins,
                                unsigned int h_bins) noexcept {
        _dist_width_bins  = w_bins;
        _dist_height_bins = h_bins;
    }
    void configure_distribution(std::string w_axis_title,
                                std::string h_axis_title) noexcept {
        _dist_w_title = std::move(w_axis_title);
        _dist_h_title = std::move(h_axis_title);
    }

    /// Analyze the properties of a set of keypoints. The booleans are toggles
    /// to deactivate analysis for the specified quantity to save some
    /// computations.
    void analyze(gsl::span<const cv::KeyPoint> points,
                 bool                          distribution = true,
                 bool                          size         = true,
                 bool                          response     = true) noexcept;

    [[nodiscard]] const distribution_histo_t& distribution() const noexcept {
        return _distribution;
    }
    [[nodiscard]] const histo_t& size_histo() const noexcept {
        return _size_histo;
    }
    [[nodiscard]] const histo_t& response_histo() const noexcept {
        return _response_histo;
    }

    [[nodiscard]] const statistic& size() const noexcept { return _size; }
    [[nodiscard]] const statistic& response() const noexcept {
        return _response;
    }

  private:
    unsigned int _img_width  = 0U;
    unsigned int _img_height = 0U;

    bool         _size_histo_enabled = true;
    statistic    _size;
    histo_t      _size_histo;
    unsigned int _size_bins  = 50U;
    std::string  _size_title = "size of keypoints";

    bool         _response_histo_enabled = true;
    statistic    _response;
    histo_t      _response_histo;
    unsigned int _response_bins  = 50U;
    std::string  _response_title = "response of keypoints";

    distribution_histo_t _distribution;
    unsigned int         _dist_width_bins  = 50U;
    unsigned int         _dist_height_bins = 50U;
    std::string          _dist_w_title = "width - distribution of keypoints";
    std::string          _dist_h_title = "height - distribution of keypoints";
};

/// Write the statistics for the keypoint distribution in a file with OpenCVs
/// FileStorage API.
void write(cv::FileStorage& fs, const std::string& name, const keypoints& kp);

}  // namespace sens_loc::analysis

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

max_stars_repo_name: JonasToth/depth-conversions

max_issues_repo_name: JonasToth/depth-conversions

max_forks_repo_name: JonasToth/depth-conversions


text:
#ifndef ALM_KATYUSHA_H
#define ALM_KATYUSHA_H

#include "Primal_Dual_LOOPLESS_Katyusha0.h"

#include <string>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <stdio.h>      /* printf */
#include <time.h>
#include <fstream>
#include <algorithm>
#include <iomanip>
#include <ctime>
#include <sstream>
//#include "cmd_line.h"

//This class implements the method IPALM_Katyusha

/*
The optimization problem to solve is:

min \sum_i f^i(A_ix)+ \sum_i h^i(M_ix)+ P(x)
Assumption 1: For each i, f_i is 1-smooth, P(x) is seperable
// Each subproblem solves problem of the form \sum_i f^i(A_ix)+ \sum_i h^i_{\beta_s}(M_ix;\lambda_s^i) +P(x)+ \beta_s/2\|x-x_s\|^2 by L_Katyuhsa
*/

template<typename L, typename D>
class ALM_Katyusha: public Primal_Dual_LOOPLESS_Katyusha0<L, D>
{
private:

  std::vector<D> Ax_s;
  std::vector<D> old_Ax_s;

  std::vector<D> Mx_s;
  std::vector<D> old_Mx_s;

protected:
  Matrix<L, D>	data_A;

  Matrix<L, D> data_M;

  std::vector<D> x_s;

  std::vector<D> old_x_s;

  std::vector<D> lambda_s;

  std::vector<D> old_lambda_s;

  D beta_s;

  D epsilon_s;

  D m_s;

  D m_0;

  L m_1;

  L m_2;

  L m_3;

  D eta;

  D rho;

  D tau_s;

  std::vector<D> L_phi;

  D function_value;

  L print_every_N_ALM;

  D running_time_ALM;

  L nb_outer_iters;

  ofstream samp_ALM;

  ofstream samp_x_ALM;

public:

  D mu_g;

  D lambda1;

  D lambda2;

  D L_h;

  virtual inline D value_of_f_j(D, L){return D(NULL);}
  virtual inline D value_of_h_j(D, L){return D(NULL);}
  virtual inline D gradient_of_f_j(D, L){return D(NULL);}
  virtual inline D value_of_f_star_j(D, L){return D(NULL);}
  virtual inline D value_of_h_star_j(D, L){return D(NULL);}
  virtual inline D prox_of_h_star_j(D,D, L){return D(NULL);}
  virtual inline D value_of_P_j(D, L){return D(NULL);}
  virtual inline D prox_of_P_j(D, D, L){return D(NULL);}
  virtual inline D feasible_dual(vector<D> &){return D(NULL);}
  virtual inline void compute_just_in_time_prox_grad_without_x_s(D, D, D &, L,L, D, D &, L){}
  virtual inline void rescale(){}
  virtual inline void set_matrix_M(){}
  virtual inline void set_matrix_A(){}

  ALM_Katyusha()
  : Primal_Dual_LOOPLESS_Katyusha0<L,D>(),data_A(),data_M()
  {
  	this->gamma= 1;
  }

  ALM_Katyusha(const char* matrix_file, const char* vector_file)
  : Primal_Dual_LOOPLESS_Katyusha0<L,D>(matrix_file, vector_file),data_A(), data_M()
  {
    this->gamma=1;
  }

  ALM_Katyusha(const char* matrix_file)
  : Primal_Dual_LOOPLESS_Katyusha0<L,D>(matrix_file),data_A(), data_M()
  {
    this->gamma=1;
  }
/*
  ALM_Katyusha(const char* matrix_file, const char* matrix_file2)
  : Primal_Dual_LOOPLESS_Katyusha0<L,D>(),data_A(matrix_file), data_M(matrix_file2)
  {
  	this->matrix_merge(data_A,data_M);
    this->gamma=1;
  }
 */
  inline void set_L_phi(){
  	for(int i=m_3; i<m_2+ m_3; i++){
  		L_phi[i]= this->nsamples;
	  }
  	for(int i= 0; i<m_3; i++){
    	L_phi[i]= this->nsamples/beta_s;
	}
  }

  L get_nb_features(){
    return data_A.get_d();
  }

  inline D value_of_phi_i(D x, L i) {
    if ( i>= data_M.nsamples){
      return value_of_f_j(x, i- data_M.nsamples);
    }
    else{
      D tmp= prox_of_h_star_j(1.0/beta_s*x+ lambda_s[i],beta_s,i);
      return beta_s*(x*tmp- value_of_h_star_j(tmp,i)- beta_s/2*(tmp- lambda_s[i])*(tmp- lambda_s[i]));
    }
  }

  inline D gradient_of_phi_i(D x, L i){
    if (i>= data_M.nsamples){
      return gradient_of_f_j(x,i- data_M.nsamples);
    }
    else{
      D tmp= prox_of_h_star_j(1.0/beta_s*x+ lambda_s[i],beta_s,i);
      return beta_s*tmp;
    }

  }

  inline D value_of_phistar_i(D x, L i) {
    if ( i>= data_M.nsamples){
      return value_of_f_star_j(x,i- data_M.nsamples);
    }
    else{
      return beta_s*value_of_h_star_j(x/beta_s,i)+ 0.5*(x- beta_s*lambda_s[i])*(x- beta_s*lambda_s[i]);
    }
  }

  inline D value_of_g_j(D x, L j){
  	return value_of_P_j(x,j)+ beta_s/2*(x- x_s[j])*(x- x_s[j]);
  }

  inline D value_of_gstar(vector<D> &x){
  	D res= 0;
  	L l= x.size();
  	D tmp= 0;
  	for (L i= 0; i< l; i++){
  		tmp= prox_of_P_j(1.0/beta_s*x[i]+ x_s[i],beta_s,i);
        res+= x[i]*tmp- value_of_P_j(tmp,i)- beta_s/2*(tmp- x_s[i])*(tmp- x_s[i]);
	  }
  	return res;
  }

  inline D gradient_of_gstar_j(D x, L j){
  	return prox_of_P_j(1.0/beta_s*x+ x_s[j],beta_s,j);
  }

  inline D value_of_gstar_minus(vector<D> &x, D scal){
  	L l= x.size();
  	vector<D> tmp(l,0);
  	for (L j= 0;j< l; j++){
  		tmp[j]= -scal*x[j];
	  }
  	return value_of_gstar(tmp);
  }
 /* inline D get_lambda1(){return lambda1;}
  inline D get_lambda2(){return lambda2;}
*/

  void compute_x(){
    for (L i=0; i< this->nfeatures; i++){
      x_s[i]= this->primal_x[i];
    }
    for (L j=m_3; j< m_2+ m_3; j++){
      Ax_s[j- m_3]= compute_AiTx_s(j);
    }
    for (L j=0; j< m_3; j++){
      Mx_s[j]= compute_AiTx_s(j);
	}
  }

  D compute_AiTx_s(L i){
    D res=0;
    for (L k = this->ptr[i]; k < this->ptr[i + 1];k++)
    {
      L j=this->row_idx[k];
      res+= this->A[k]*x_s[j];
    }
    return res;
  }

  inline void compute_just_in_time_prox_grad(D tau, D u, D &x, L t0,L t1, D w, D &y, L j){
    D u_s=u-tau_s*beta_s*x_s[j];
    compute_just_in_time_prox_grad_without_x_s(tau, u_s,  x,  t0,  t1,  w, y,j);
  }

  void compute_function_value(){
    D res= 0;
    for (L i= 0; i< this->nfeatures; i++){
      res+= value_of_P_j(x_s[i],i);
    }
    for (L i= 0; i< data_M.nsamples; i++){
      res+= value_of_h_j(Mx_s[i],i);
    }
    for (L i= 0; i<data_A.nsamples; i++){
      res+= value_of_f_j(Ax_s[i],i);
    }
    function_value= res;
  }

  /*void compute_infeasibility(){
  D tmp1= 0;
  D tmp2= 0;
  for (L j=0; j< this->nsamples; j++){
  if (1- Ax_s[j]> tmp1){
  tmp1= 1- Ax_s[j];
}
tmp2+= max(0.0,1- Ax_s[j]);
}
residual1= tmp1;
residual2= tmp2/this->nsamples;
}*/

inline void compute_m0(D beta0, D val_eta, D val_rho, L val_tau){
  cout<<"beta0="<<beta0<<"val_eta="<<val_eta<<"; val_rho="<<val_rho<<"; val_tau="<<val_tau<<endl;
  m_s= this->nsamples/val_tau*1e+6;

  //D tmp7=1-val_tau/this->n*sqrt(beta0/(max_Lf_s+max_M_s/beta0+beta0));
  //cout<<"tmp7: "<<tmp7<<"; "<<val_tau/(this->n+0.)*sqrt(beta0/(max_Lf_s+max_M_s/beta0+beta0))<<endl;
  //m_0=(2*log(val_rho)+log(val_eta)+log(2))/log(tmp7);

  cout<<" m_s="<<m_s<<endl;
}

void update_y(){
  for(L j=0;j<m_3;j++)
  {
    D aiTx= Mx_s[j];
    lambda_s[j]= prox_of_h_star_j(1.0/beta_s*aiTx+ lambda_s[j],beta_s,j);
  }
}

void Initialize(D beta_0, D epsilon_0,  D val_eta, D val_rho,L val_tau, vector<D> & x0,vector<D> & y0){
  cout<<"start initializing"<<endl;
  set_matrix_M();
  set_matrix_A();

  this->tau=val_tau;
  m_1=data_A.get_d();
  m_2=data_A.get_n();
  m_3=data_M.get_n();
  cout<<"m_1="<<m_1<<endl;
  cout<<"m_2="<<m_2<<endl;
  cout<<"m_3="<<m_3<<endl;

  beta_s=beta_0;
  //tau_s= this->nsamples;
  tau_s= 1;

  epsilon_s=epsilon_0;
  compute_m0(beta_s,val_eta,val_rho,val_tau);

  eta=val_eta;
  rho=val_rho;

  x_s.resize(m_1,0);
  old_x_s.resize(m_1,0);
  lambda_s.resize(m_3,0);
  old_lambda_s.resize(m_3,0);
  for(L i=0;i<m_1;i++){
    x_s[i]=x0[i];
    old_x_s[i]= x0[i];
  }
  for(L j=0;j<m_3;j++){
    lambda_s[j]=y0[j];
    old_lambda_s[j]= y0[j];
  }

  Ax_s.clear();
  Ax_s.resize(m_2,0);
  old_Ax_s.clear();
  old_Ax_s.resize(m_2,0);
  Mx_s.clear();
  Mx_s.resize(m_3,0);
  old_Mx_s.clear();
  old_Mx_s.resize(m_3,0);
  L_phi.clear();
  L_phi.resize(m_2+ m_3,0);
  for(int i=0; i< m_3; i++){
  	L_phi[i]= this->nsamples/beta_0;
  }
  for(int i=m_3; i< m_2+ m_3; i++){
  	L_phi[i]= this->nsamples;
  }
}

void reset_everything(){
  epsilon_s*=rho;
  beta_s*=eta;

}

void update_m_s(L val_tau){
  D tmp= 0;
  D tmpx= 0;
  D tmpy= 0;
  D tmpy2= 0;
  for(L j=0;j<m_3;j++)
  {
    D tmp_p= prox_of_h_star_j(Mx_s[j]/beta_s/eta+ lambda_s[j], beta_s/eta, j);
    tmpy+=(tmp_p- lambda_s[j])*(tmp_p- lambda_s[j]);
    tmp+= (lambda_s[j]- old_lambda_s[j])*(lambda_s[j]- old_lambda_s[j]);
    tmpy2+= (beta_s*eta*lambda_s[j]- beta_s*old_lambda_s[j])*(beta_s*eta*lambda_s[j]- beta_s*old_lambda_s[j]);
  }
  for (L i= 0; i< m_1;i++){
    tmpx= (x_s[i]- old_x_s[i])*(x_s[i]- old_x_s[i]);
  }
  D tmp4= 2*epsilon_s+ beta_s*tmp;
  D tmp5= (1- eta)*beta_s/2*tmpy;
  D tmp6= sqrt(tmp)*beta_s*(1+ eta)*L_h;
  D tmp7= sqrt(tmp)*sqrt(tmpy2);
  D tmp1= log(tmp4+ tmp5+ tmp6+ tmp7+ 2/(2*eta- 1)*tau_s*beta_s/2*tmpx);
  D tmp2= log(epsilon_s)+ log(rho)- log(2);
  D tmp3= sqrt(lambda2/this->sumLi);
  m_s= ceil((tmp1- tmp2)/log(2)*4*max(1.0*this->nsamples,1/tmp3)/val_tau);
}

inline void compute_and_record_res(){
  if(nb_outer_iters%print_every_N_ALM==0){
    compute_function_value();
    cout<<setprecision(9)<<"Iteration: "<<nb_outer_iters<<"; time="<<running_time_ALM<< "; function value="<<function_value<<endl;
    samp_ALM<<setprecision(9)<<nb_outer_iters<<" "<<running_time_ALM<<" "<< function_value<<" "<<endl;
  }
}

void ALM_solve_with_L_Katyusha(D beta_0, D epsilon_0,  D eta, D rho,vector<D> & x0,vector<D> & y0, L val_tau, L max_nb_outer, L p_N_1, L p_N_2, string filename1, string filename2, D time){
  Initialize(beta_0, epsilon_0, eta, rho,val_tau, x0, y0);
  nb_outer_iters=0;
  //string sampname2= ALGparam.data_dir +"/results/L_Katyusha_"+filename2;
  string sampname2= "results/L_Katyusha_"+filename2;
  this->samp.open(sampname2.c_str());
  string sampname_x="results/L_Katyusha_x_"+filename1;
  //filename1= ALGparam.data_dir +"/results/ALM_"+filename1;
  filename1= "results/ALM_"+filename1;
  samp_x_ALM.open(sampname_x.c_str());
  samp_ALM.open(filename1.c_str());
  running_time_ALM=0;
  print_every_N_ALM=p_N_1;
  this->set_print_every_N(p_N_2);
  compute_and_record_res();
  D start;
  start = std::clock();
  //cout<<"m0="<<ceil(m_s/this->n*val_tau)<<endl;
  cout<<"m_s= "<< ceil(m_s/this->nsamples*val_tau)<<"; beta_s="<<beta_s<<"; epsilon_s="<<epsilon_s<<endl;
  //this->L_Katyusga_MU(x_s, val_tau, val_mu_f, val_mu_g, 3, p_N_2, ceil(m_s/this->n*val_tau), epsilon_s, filename2,1);
  rescale();
  set_L_phi();
  this->loopless_Katyusha(y0, x_s, filename2, L_phi, mu_g+ tau_s*beta_s, epsilon_s, ceil(m_s/this->nsamples*val_tau), val_tau, 1, 1, 0, 1);
  for(L i=0;i<m_1;i++){
    old_x_s[i]=x_s[i];
  }
  for(L i=0;i<m_2;i++){
    old_Ax_s[i]=Ax_s[i];
  }
  compute_x();
  for(L i=0;i<m_3;i++){
    old_lambda_s[i]=lambda_s[i];
  }
  update_y();
  update_m_s(val_tau);
  nb_outer_iters++;
  running_time_ALM+=( std::clock() - start ) / (double) CLOCKS_PER_SEC;
  compute_and_record_res();
  start = std::clock();
  reset_everything();
  running_time_ALM+=( std::clock() - start ) / (double) CLOCKS_PER_SEC;
  while(nb_outer_iters<max_nb_outer){
    start = std::clock();
    cout<<"ms="<<ceil(m_s/this->nsamples*val_tau)<<"; beta_s="<<beta_s<<"; epsilon_s"<<epsilon_s<<endl;
    rescale();
    set_L_phi();
    this->loopless_Katyusha(y0, x_s, filename2, L_phi, mu_g+ tau_s*beta_s, epsilon_s, ceil(m_s/this->nsamples*val_tau), val_tau, 1, 1, 0, 1);
    for(L i=0;i<m_1;i++){
      old_x_s[i]=x_s[i];
    }
    for(L i=0;i<m_2;i++){
      old_Ax_s[i]=Ax_s[i];
    }
    for(L i=0;i<m_3;i++){
      old_lambda_s[i]=lambda_s[i];
    }
    compute_x();
    update_y();
    update_m_s(val_tau);
    nb_outer_iters++;
    running_time_ALM+=( std::clock() - start ) / (double) CLOCKS_PER_SEC;
    compute_and_record_res();
    start = std::clock();
    reset_everything();
    running_time_ALM+=( std::clock() - start ) / (double) CLOCKS_PER_SEC;
    if (running_time_ALM> time){
      break;
    }
  }
  for (L i= 0; i< m_1; i++){
      	samp_x_ALM<< x_s[i]<< endl;
  }
}
};

#endif /* MIN_SMOOTH_CONVEX_H */


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

/* Author:  G. Jungman */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include "gsl_sf_exp.h"
#include "gsl_sf_bessel.h"

#include "error.h"

#include "chebyshev.h"
#include "cheb_eval.c"

/*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/

/* based on SLATEC bk0(), bk0e() */

/* chebyshev expansions

 series for bk0        on the interval  0.	    to  4.00000d+00
					with weighted error   3.57e-19
					 log weighted error  18.45
			       significant figures required  17.99
				    decimal places required  18.97

 series for ak0        on the interval  1.25000d-01 to  5.00000d-01
					with weighted error   5.34e-17
					 log weighted error  16.27
			       significant figures required  14.92
				    decimal places required  16.89

 series for ak02       on the interval  0.	    to  1.25000d-01
					with weighted error   2.34e-17
					 log weighted error  16.63
			       significant figures required  14.67
				    decimal places required  17.20
*/

static double bk0_data[11] = {
  -0.03532739323390276872,
   0.3442898999246284869,
   0.03597993651536150163,
   0.00126461541144692592,
   0.00002286212103119451,
   0.00000025347910790261,
   0.00000000190451637722,
   0.00000000001034969525,
   0.00000000000004259816,
   0.00000000000000013744,
   0.00000000000000000035
};
static cheb_series bk0_cs = {
  bk0_data,
  10,
  -1, 1,
  10
};

static double ak0_data[17] = {
  -0.07643947903327941,
  -0.02235652605699819,
   0.00077341811546938,
  -0.00004281006688886,
   0.00000308170017386,
  -0.00000026393672220,
   0.00000002563713036,
  -0.00000000274270554,
   0.00000000031694296,
  -0.00000000003902353,
   0.00000000000506804,
  -0.00000000000068895,
   0.00000000000009744,
  -0.00000000000001427,
   0.00000000000000215,
  -0.00000000000000033,
   0.00000000000000005
};
static cheb_series ak0_cs = {
  ak0_data,
  16,
  -1, 1,
  10
};

static double ak02_data[14] = {
  -0.01201869826307592,
  -0.00917485269102569,
   0.00014445509317750,
  -0.00000401361417543,
   0.00000015678318108,
  -0.00000000777011043,
   0.00000000046111825,
  -0.00000000003158592,
   0.00000000000243501,
  -0.00000000000020743,
   0.00000000000001925,
  -0.00000000000000192,
   0.00000000000000020,
  -0.00000000000000002
};
static cheb_series ak02_cs = {
  ak02_data,
  13,
  -1, 1,
  8
};

/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/

int gsl_sf_bessel_K0_scaled_e(const double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(x <= 0.0) {
    DOMAIN_ERROR(result);
  }
  else if(x <= 2.0) {
    const double lx = log(x);
    const double ex = exp(x);
    int stat_I0;
    gsl_sf_result I0;
    gsl_sf_result c;
    cheb_eval_e(&bk0_cs, 0.5*x*x-1.0, &c);
    stat_I0 = gsl_sf_bessel_I0_e(x, &I0);
    result->val  = ex * ((-lx+M_LN2)*I0.val - 0.25 + c.val);
    result->err  = ex * ((M_LN2+fabs(lx))*I0.err + c.err);
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return stat_I0;
  }
  else if(x <= 8.0) {
    const double sx = sqrt(x);
    gsl_sf_result c;
    cheb_eval_e(&ak0_cs, (16.0/x-5.0)/3.0, &c);
    result->val  = (1.25 + c.val) / sx;
    result->err  = c.err / sx;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    const double sx = sqrt(x);
    gsl_sf_result c;
    cheb_eval_e(&ak02_cs, 16.0/x-1.0, &c);
    result->val  = (1.25 + c.val) / sx;
    result->err  = (c.err + GSL_DBL_EPSILON) / sx;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
}

int gsl_sf_bessel_K0_e(const double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(x <= 0.0) {
    DOMAIN_ERROR(result);
  }
  else if(x <= 2.0) {
    const double lx = log(x);
    int stat_I0;
    gsl_sf_result I0;
    gsl_sf_result c;
    cheb_eval_e(&bk0_cs, 0.5*x*x-1.0, &c);
    stat_I0 = gsl_sf_bessel_I0_e(x, &I0);
    result->val  = (-lx+M_LN2)*I0.val - 0.25 + c.val;
    result->err  = (fabs(lx) + M_LN2) * I0.err + c.err;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return stat_I0;
  }
  else {
    gsl_sf_result K0_scaled;
    int stat_K0 = gsl_sf_bessel_K0_scaled_e(x, &K0_scaled);
    int stat_e  = gsl_sf_exp_mult_err_e(-x, GSL_DBL_EPSILON*fabs(x),
                                           K0_scaled.val, K0_scaled.err,
					   result);
    return GSL_ERROR_SELECT_2(stat_e, stat_K0);
  }
}

/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/

#include "eval.h"

double gsl_sf_bessel_K0_scaled(const double x)
{
  EVAL_RESULT(gsl_sf_bessel_K0_scaled_e(x, &result));
}

double gsl_sf_bessel_K0(const double x)
{
  EVAL_RESULT(gsl_sf_bessel_K0_e(x, &result));
}


text:
/* Copyright 2016, Alistair Boyle, 3-clause BSD License */
#include "config.h"
#include <stdio.h>
#include <stdlib.h> /* malloc, free */
#include <string.h> /* memcpy, memset */
#include <math.h> /* sqrt */
#include <assert.h> /* assert */
#include <time.h> /* clock */
//#include <cblas.h>
//#include <lapacke.h>

#include "argv.h"
#include "file.h"
#include "fwd.h"
#include "inv.h"
#include "model.h"
#include "matrix.h"

static int do_fwd_solve(const char * file, double tol)
{
    int ret = 1;
    model  * mdl = malloc_model();
    printf("-- loading from file %s --\n", file);
    int nret = readfile(file, mdl);
    const int n_meas = mdl->n_stimmeas;
    const int frames = mdl->n_params[1];
    double meas[frames * n_meas];
    if(!nret) {
        fprintf(stderr, "error: read %s failed\n", file);
        goto fwd_quit;
    }
    printf("-- forward solutions --\n");
    printf("%d parameter frame%s\n", frames, frames == 1 ? "" : "s");
    if(frames < 1) {
        fprintf(stderr, "error: nothing to do\n");
        goto fwd_quit;
    }
    if(mdl->n_data[0] > 0) {
        assert(mdl->n_data[0] == mdl->n_stimmeas);
        assert(mdl->n_data[1] == mdl->n_params[1]);
    }
    nret = fwd_solve(mdl, meas);
    if(!nret) {
        fprintf(stderr, "error: bad forward solve\n");
        goto fwd_quit;
    }
    for(int idx = 0; idx < frames; idx++) {
        /* calculated measurements to expected */
        if(mdl->n_data[0] > 0) {
            printf("frame#%d\n", idx + 1);
            printf("meas#  %18s %18s\n", "calculated", "from file");
            double rmse = 0; /* root mean squared error from expected */
            for(int i = 0; i < mdl->n_stimmeas; i++) {
                const double data_calc = meas[idx * n_meas + i];
                const double data_expect = mdl->data[i + mdl->n_data[0] * idx];
                printf("%4d  %18.8g %18.8g\n", i + 1, data_calc, data_expect);
                const double dm = data_calc - data_expect;
                rmse += dm * dm;
            }
            rmse = sqrt(rmse / (double) mdl->n_data[0]);
            printf("                      RMSE = %g\n", rmse);
            printf("                      tolerance = %g\n", tol);
            if (rmse > tol) {
                fprintf(stderr, "fail: calculated measurements do not match expected\n");
                ret = 1;
                goto fwd_quit;
            }
        }
        else {
            printf("frame#%d\n", idx + 1);
            printf("meas#  %18s\n", "calculated");
            for(int i = 0; i < mdl->n_stimmeas; i++) {
                const double data_calc = meas[idx * n_meas + i];
                printf("%4d  %18.04g\n", i + 1, data_calc);
            }
        }
    }
    ret = 0;
    printf("-- completed --\n");
fwd_quit:
    free_model(mdl);
    return ret;
}

#define max(a,b) (((a) > (b)) ? (a) : (b))
static int do_inv_solve(const char * file, const double tol)
{
    int ret = 1;
    model  * mdl = malloc_model();
    printf("-- loading from file %s --\n", file);
    int nret = readfile(file, mdl);
    const int frames = mdl->n_data[1] > 1 ? mdl->n_data[1] - 1 : 0;
    if(!nret) {
        fprintf(stderr, "error: read %s failed\n", file);
        free_model(mdl);
        return 1;
    }
    printf("-- inverse solutions --\n");
    printf("%d measurement frame%s\n", frames, frames == 1 ? "" : "s");
    if(frames < 1) {
        fprintf(stderr, "error: nothing to do\n");
        free_model(mdl);
        return 1;
    }
    printf_model(mdl, 1);
    assert(mdl->n_data[0] > 0);
    assert(mdl->n_data[0] == mdl->n_stimmeas);
    if(mdl->n_params[1] > 1) {
        assert(mdl->n_data[1] == mdl->n_params[1]);
    }
    const int is_pmap = (mdl->fwd.n_pmap > 0);
    int n_params = 0;
    if( mdl->n_params[1] > 0 ) {
        n_params = mdl->n_params[0];
    }
    else if( !is_pmap ) {
        n_params = mdl->fwd.n_elems;
    }
    else {
        for(int i = 0; i < mdl->fwd.n_pmap; i++) {
            n_params = max(n_params, mdl->fwd.pmap_param[i]);
        }
        n_params++;
    }
    mdl->n_params[0] = n_params;
    double params[frames * n_params];
    memset(params, 0, sizeof(double) * frames * n_params);
    nret = inv_solve(mdl, params);
    if(!nret) {
        fprintf(stderr, "error: bad inverse solve\n");
        goto inv_quit;
    }
    for(int idx = 0; idx < frames; idx++) {
        /* calculated measurements to expected */
        if(mdl->n_params[1] > 1) {
            printf("frame#%d\n", idx + 1);
            printf("param#%18s %18s\n", "calculated", "from file");
            const int cols = mdl->n_params[1];
            const int rows = mdl->n_params[0];
            assert(cols == 3); /* TODO currently: params = [ J_background, ignore, expected ] */
            double rmse = 0; /* root mean squared error from expected */
            for(int i = 0; i < rows; i++) {
                const double params_calc = params[i + 0 * rows];
                const double params_expect = mdl->params[i + 2 * rows];
                printf("%4d  %18.8e %18.8e\n", i + 1, params_calc, params_expect);
                const double dp = params_calc - params_expect;
                rmse += dp * dp;
            }
            rmse = sqrt(rmse / (double) n_params);
            printf("                      RMSE = %g\n", rmse);
            if (rmse > tol) {
                fprintf(stderr, "fail: calculated parameters do not match expected\n");
                goto inv_quit;
            }
            break; /* TODO only handle one frame in checking mode currently */
        }
        else {
            printf("frame#%d\n", idx + 1);
            printf("param#%18s\n", "calculated");
            const int rows = mdl->n_params[0];
            for(int i = 0; i < n_params; i++) {
                const double params_calc = params[i + idx * rows];
                printf("%4d  %18.04g\n", i + 1, params_calc);
            }
        }
    }
    ret = 0;
    printf("-- completed --\n");
inv_quit:
    free_model(mdl);
    return ret;
}

int main(int argc, char ** argv)
{
    int ret = 0;
    args_t args;
    const int N = 1;
    double dt [N];
    double mean_dt = 0;
    double std_dt = 0;
    if(parse_argv(argc, argv, &args) != 0) {
        return 1;
    }
    for(int i = 0; i < N; i++) {
        clock_t start = clock();
        printf("%s %s %d\n", args.file[0], args.file[1], args.mode);
        switch(args.mode) {
        case FORWARD_SOLVER: ret = do_fwd_solve(args.file[0], args.tol); break;
        case INVERSE_SOLVER: ret = do_inv_solve(args.file[1], args.tol); break;
        default: return 0;
        }
        clock_t stop = clock();
        dt[i] = (stop - start) / CLOCKS_PER_SEC;
        mean_dt += dt[i] / N;
        if(ret) {
            break;
        }
    }
    if(ret) {
        printf("error: solve failed\n");
        return ret;
    }
    for(int i = 0; i < N; i++) {
        const double tmp = dt[i] - mean_dt;
        std_dt += tmp * tmp;
    }
    std_dt = sqrt(std_dt / N);
    printf("%f\n", mean_dt);
    printf("runtime: %0.0f:%02.0f:%02.4fs ± %2.4fs\n", 0.0, 0.0, mean_dt, std_dt);
    printf("total:   %0.0f:%02.0f:%02.4fs\n", 0.0, 0.0, mean_dt*N);
    return 0;
}

max_stars_repo_name: boyle/zedhat

max_issues_repo_name: boyle/zedhat

max_forks_repo_name: boyle/zedhat