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

#ifndef MCAS_HSTORE_HEAP_MC_EPHEMERAL_H
#define MCAS_HSTORE_HEAP_MC_EPHEMERAL_H

#include "heap_mm_ephemeral.h"

#include "hstore_config.h"
#include "histogram_log2.h"
#include "hop_hash_log.h"
#include <mm_plugin_itf.h>
#include "persistent.h"

#include <ccpm/cca.h> /* ctor_args */
#include <ccpm/interfaces.h> /* ownership_callback, (IHeap_expandable, region_vector_t) */
#include <common/byte_span.h>
#include <common/string_view.h>
#include <nupm/region_descriptor.h>
#include <gsl/span>

#include <algorithm> /* min, swap */
#include <cstddef> /* size_t */
#include <functional> /* function */
#include <memory> /* unique_ptr */
#include <vector>

namespace impl
{
	struct allocation_state_pin;
	struct allocation_state_emplace;
	struct allocation_state_extend;
}

struct heap_mc_ephemeral
  : public heap_mm_ephemeral
{
private:
	using byte_span = common::byte_span;
	using string_view = common::string_view;
	std::unique_ptr<
		ccpm::IHeap_expandable
	> _heap;
	impl::allocation_state_emplace *_ase;
	impl::allocation_state_pin *_aspd;
	impl::allocation_state_pin *_aspk;
	impl::allocation_state_extend *_asx;

	static constexpr unsigned log_min_alignment = 3U; /* log (sizeof(void *)) */
	static_assert(sizeof(void *) == 1U << log_min_alignment, "log_min_alignment does not match sizeof(void *)");
	/* Rca_LB seems not to allocate at or above about 2GiB. Limit reporting to 16 GiB. */
	static constexpr unsigned hist_report_upper_bound = 34U;
	explicit heap_mc_ephemeral(
		unsigned debug_level_
		, bool restore_not_clear
		, impl::allocation_state_emplace *ase
		, impl::allocation_state_pin *aspd
		, impl::allocation_state_pin *aspk
		, impl::allocation_state_extend *asx
		, std::unique_ptr<ccpm::IHeap_expandable> p
		, string_view id
		, string_view backing_file
		, const std::vector<byte_span> rv_full
		, const byte_span pool0_heap
	);

public:
	friend struct heap_mc;
	friend struct heap_mm;

	using common::log_source::debug_level;

	/* heap_mm version */
	explicit heap_mc_ephemeral(
		unsigned debug_level
		, bool restore_not_clear
		, MM_plugin_wrapper &&pw
		, impl::allocation_state_emplace *ase
		, impl::allocation_state_pin *aspd
		, impl::allocation_state_pin *aspk
		, impl::allocation_state_extend *asx
		, string_view id
		, string_view backing_file
		, const std::vector<byte_span> rv_full
		, byte_span pool0_heap_
	);

	std::size_t allocated() const override; // { return _allocated; }
	std::size_t capacity() const override; // { return _capacity; }
	void allocate(persistent_t<void *> &p, std::size_t sz, std::size_t alignment) override;
	std::size_t free(persistent_t<void *> &p_, std::size_t sz_) override;
	void free_tracked(const void *p, std::size_t sz) override;
	heap_mc_ephemeral(const heap_mc_ephemeral &) = delete;
	heap_mc_ephemeral& operator=(const heap_mc_ephemeral &) = delete;
	void add_managed_region_to_heap(byte_span r_heap) override;
	void reconstitute_managed_region_to_heap(byte_span r_heap, ccpm::ownership_callback_t f) override;
	bool is_crash_consistent() const override;
	bool can_reconstitute() const override;
};

#endif

max_stars_repo_name: omriarad/mcas

max_issues_repo_name: omriarad/mcas

max_forks_repo_name: omriarad/mcas


text:
/**
* VGMTrans (c) - 2002-2021
* Licensed under the zlib license
* See the included LICENSE for more information
*/

#pragma once

#include <vector>
#include <fluidsynth.h>
#include <gsl-lite.hpp>

class MusicPlayer {
public:
  MusicPlayer();
  ~MusicPlayer();

  /**
   * Toggles the status of the player
   * @returns true if the player is playing
   */
  bool toggle();
  /**
   * Stops the player (unloads resources)
   */
  void stop();
  /**
   * Moves the player position
   * @param position relative to song start
   */
  void seek(int position);

  /**
   * Loads SF2 and MIDI data into the player; resources are copied
   * @param soundfont_data
   * @param midi_data
   * @return true if data was loaded correctly
   */
  bool loadDataAndPlay(gsl::span<char> soundfont_data, gsl::span<char> midi_data);

  /**
   * Checks whether the player is playing
   * @return true if playing
   */
  [[nodiscard]] bool playing() const;
  /**
   * The number of MIDI ticks elapsed since the player was started
   * @return number of ticks relative to song start
   */
  [[nodiscard]] int elapsedTicks() const;
  /**
   * The total number of MIDI ticks in the song
   * @return total ticks in the song
   */
  [[nodiscard]] int totalTicks() const;

  /**
   * Gets all the audio driver the player supports
   * @return the driver list
   */
  [[nodiscard]] std::vector<const char *> getAvailableDrivers() const;
  /**
   * Changes the current audio driver and restarts playing
   * @param driver_name
   * @return true if the driver changed
   */
  bool setAudioDriver(const char* driver_name);

private:
  fluid_settings_t *m_settings = nullptr;
  fluid_synth_t *m_synth = nullptr;
  fluid_audio_driver_t *m_active_driver = nullptr;
  fluid_player_t *m_active_player = nullptr;

  void makeSettings();
  void makeSynth();
  void makePlayer();
};

ext: h

max_stars_repo_name: Oipo/vgmtrans

max_issues_repo_name: RGBA-CRT/vgmtrans

max_forks_repo_name: RGBA-CRT/vgmtrans


text:
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <time.h>

double ar1_ts (double * params, int n, gsl_rng * r)
{
    double beta = params[0];
    double alpha = params[1];
    double s = params[2];

    int i;
    double x = beta / (1 - alpha);  // Start at mean of stationary dist
    double sum = 0;
    for (i = 1; i <= n; i++) {
        sum += x;
        x = beta + alpha * x + gsl_ran_gaussian(r, s);
     }

    return sum / n;
}

int main(void)
{
    clock_t start, end;
    double cpu_time_used;

    int N = 1e7;
    double beta = 1.0;
    double alpha = 0.9;
    double s = 1;
    double params[3] = {beta, alpha, s};

    /* create a generator via GSL_RNG_TYPE */
    const gsl_rng_type * T;
    gsl_rng * r;
    gsl_rng_env_setup();
    T = gsl_rng_default;
    r = gsl_rng_alloc(T);
    gsl_rng_set(r, 1);

    start = clock();
    double sample_mean = ar1_ts(params, N, r);
    end = clock();
    cpu_time_used = ((double) (end - start)) / CLOCKS_PER_SEC;

    gsl_rng_free (r);

    printf("mean = %g\n", sample_mean);
    printf("time elapsed = %g seconds\n", cpu_time_used);
    return 0;
}

ext: c


text:
#pragma once

#include <cstddef>
#include <gsl/string_span>
#include <map>
#include <memory>
#include <string_view>
#include <vector>

namespace multihash {

using string_span = gsl::span<char>;

class algorithm {
 public:
  /** Compute the hash of an input range */
  template <typename InputIterator>
  auto digest(InputIterator first, InputIterator last) -> std::string;

  /** Write the hash of an input range into output */
  template <typename InputIterator, typename OutputIterator>
  auto digest(InputIterator first, InputIterator last,
              OutputIterator output) -> OutputIterator;

  virtual ~algorithm() = default;

 private:
  virtual void reset() = 0;
  virtual std::size_t digest_size() const = 0;
  virtual std::size_t block_size() const = 0;
  virtual void update(std::string_view data) = 0;
  virtual std::size_t digest(string_span output) = 0;
  std::string digest();

};

template <typename InputIterator>
auto algorithm::digest(InputIterator first, InputIterator last) -> std::string {
  auto result = std::string{};
  result.reserve(digest_size());
  digest(first, last, std::back_inserter(result));
  return result;
}

template <typename InputIterator, typename OutputIterator>
auto algorithm::digest(InputIterator first, InputIterator last,
                                     OutputIterator output) -> OutputIterator {
  reset();
  auto buffer = std::vector<char>{};
  auto chunk_size = block_size();
  buffer.reserve(chunk_size);
  while (first != last) {
    for (auto i = 0u; first != last && i < chunk_size; ++i, ++first) {
      buffer.emplace_back(*first);
    }
    auto view = std::string_view(buffer.data(), buffer.size());
    update(view);
    buffer.clear();
  }
  auto result = digest();
  std::copy(result.begin(), result.end(), output);
  return output;
}

}  // namespace multihash

max_stars_repo_name: lockblox/multihash

max_issues_repo_name: cpp-ipfs/cpp-multihash

max_forks_repo_name: cpp-ipfs/cpp-multihash


text:
#ifndef __SAMPLER_H__
#define __SAMPLER_H__

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

class Sampler {

  dsfmt_t dsfmt;

  void init_random(unsigned params_random_seed) {
#ifdef USE_MT_RANDOM
  dsfmt_init_gen_rand(&dsfmt, params_random_seed);
#else
  srand(params_random_seed);
#endif
  };

 protected:

  /* gsl random number generator */
  gsl_rng* _gsl_rng;

  /* default constructor */
  Sampler(unsigned int random_seed) {
    /* do i still need this ? */
    init_random(random_seed);
    /* init gsl random number generator */
    _gsl_rng = gsl_rng_alloc( gsl_rng_taus );
  };

  double sample_uniform();
  unsigned long int sample_integer_uniform(unsigned int from, unsigned int to);

};

#endif

max_stars_repo_name: ypetinot/web-summarization

max_issues_repo_name: ypetinot/web-summarization

max_forks_repo_name: ypetinot/web-summarization


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

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_multilarge_nlinear.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_blas.h>

/*
 * This module contains an implementation of the Steihaug-Toint
 * conjugate gradient algorithm for nonlinear optimization problems.
 * This implementation closely follows the following works:
 *
 * [1] T. Steihaug, The conjugate gradient method and trust regions
 *     in large scale optimization, SIAM J. Num. Anal., 20(3) 1983.
 *
 * In the below algorithm, the Jacobian and gradient are scaled
 * according to:
 *
 * J~ = J D^{-1}
 * g~ = D^{-1}
 *
 * prior to any calculations which results in better numerical
 * stability when solving for the Gauss-Newton step. The resulting
 * step vector is then backtransformed as:
 *
 * dx = D^{-1} dx~
 */

typedef struct
{
  size_t n;                  /* number of observations */
  size_t p;                  /* number of parameters */
  gsl_vector *z;             /* Gauss-Newton step, size p */
  gsl_vector *r;             /* steepest descent step, size p */
  gsl_vector *d;             /* steepest descent step, size p */
  gsl_vector *workp;         /* workspace, length p */
  gsl_vector *workn;         /* workspace, length n */
  double norm_g;             /* || g~ || */

  double cgtol;              /* tolerance for CG solution */
  size_t cgmaxit;            /* maximum CG iterations */
} cgst_state_t;

#include "common.c"

static void * cgst_alloc (const void * params, const size_t n, const size_t p);
static void cgst_free(void *vstate);
static int cgst_init(const void *vtrust_state, void *vstate);
static int cgst_preloop(const void * vtrust_state, void * vstate);
static int cgst_step(const void * vtrust_state, const double delta,
                     gsl_vector * dx, void * vstate);
static int cgst_preduction(const void * vtrust_state, const gsl_vector * dx,
                           double * pred, void * vstate);
static double cgst_calc_tau(const gsl_vector * p, const gsl_vector * d,
                            const double delta);

static void *
cgst_alloc (const void * params, const size_t n, const size_t p)
{
  const gsl_multilarge_nlinear_parameters *par = (const gsl_multilarge_nlinear_parameters *) params;
  cgst_state_t *state;

  state = calloc(1, sizeof(cgst_state_t));
  if (state == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate st state", GSL_ENOMEM);
    }

  state->z = gsl_vector_alloc(p);
  if (state->z == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate space for z", GSL_ENOMEM);
    }

  state->r = gsl_vector_alloc(p);
  if (state->r == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate space for r", GSL_ENOMEM);
    }

  state->d = gsl_vector_alloc(p);
  if (state->d == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate space for d", GSL_ENOMEM);
    }

  state->workp = gsl_vector_alloc(p);
  if (state->workp == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate space for workp", GSL_ENOMEM);
    }

  state->workn = gsl_vector_alloc(n);
  if (state->workn == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate space for workn", GSL_ENOMEM);
    }

  state->n = n;
  state->p = p;

  state->cgmaxit = par->max_iter;
  if (state->cgmaxit == 0)
    state->cgmaxit = n;

  state->cgtol = par->tol;

  return state;
}

static void
cgst_free(void *vstate)
{
  cgst_state_t *state = (cgst_state_t *) vstate;

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

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

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

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

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

  free(state);
}

/*
cgst_init()
  Initialize cgst solver

Inputs: vtrust_state - trust state
        vstate       - workspace

Return: success/error
*/

static int
cgst_init(const void *vtrust_state, void *vstate)
{
  /* nothing to do */

  (void)vtrust_state;
  (void)vstate;

  return GSL_SUCCESS;
}

static int
cgst_preloop(const void * vtrust_state, void * vstate)
{
  /* nothing to do */

  (void)vtrust_state;
  (void)vstate;

  return GSL_SUCCESS;
}

/*
cgst_step()
  Calculate a new step vector

Return:
GSL_SUCCESS if CG solution found
GSL_EMAXITER if no solution found
*/

static int
cgst_step(const void * vtrust_state, const double delta,
          gsl_vector * dx, void * vstate)
{
  int status;
  const gsl_multilarge_nlinear_trust_state *trust_state =
    (const gsl_multilarge_nlinear_trust_state *) vtrust_state;
  cgst_state_t *state = (cgst_state_t *) vstate;
  const gsl_vector * x = trust_state->x;
  const gsl_vector * f = trust_state->f;
  const gsl_vector * swts = trust_state->sqrt_wts;
  const gsl_vector * diag = trust_state->diag;
  const gsl_multilarge_nlinear_parameters * params = trust_state->params;
  gsl_multilarge_nlinear_fdf * fdf = trust_state->fdf;
  double alpha, beta, u;
  double norm_Jd;   /* || J D^{-1} d_i || */
  double norm_r;    /* || r_i || */
  double norm_rp1;  /* || r_{i+1} || */
  size_t i;

  /* Step 1 of [1], section 2; scale gradient as
   *
   * g~ = D^{-1} g
   *
   * for better numerical stability
   */

  for (i = 0; i < state->p; ++i)
    {
      double gi = gsl_vector_get(trust_state->g, i);
      double di = gsl_vector_get(trust_state->diag, i);

      gsl_vector_set(state->z, i, 0.0);
      gsl_vector_set(state->r, i, -gi / di);
      gsl_vector_set(state->d, i, -gi / di);
      gsl_vector_set(state->workp, i, gi / di);
    }

  /* compute || g~ || */
  state->norm_g = gsl_blas_dnrm2(state->workp);

  for (i = 0; i < state->cgmaxit; ++i)
    {
      /* workp := D^{-1} d_i */
      gsl_vector_memcpy(state->workp, state->d);
      gsl_vector_div(state->workp, trust_state->diag);

      /* workn := J D^{-1} d_i */
      status = gsl_multilarge_nlinear_eval_df(CblasNoTrans, x, f, state->workp,
                                              swts, params->h_df, params->fdtype,
                                              fdf, state->workn, NULL, NULL);
      if (status)
        return status;

      /* compute || J D^{-1} d_i || */
      norm_Jd = gsl_blas_dnrm2(state->workn);

      /* Step 2 of [1], section 2 */
      if (norm_Jd == 0.0)
        {
          double tau = cgst_calc_tau(state->z, state->d, delta);

          /* dx = z_i + tau*d_i */
          scaled_addition(1.0, state->z, tau, state->d, dx);
          gsl_vector_div(dx, diag);

          return GSL_SUCCESS;
        }

      /* Step 3 of [1], section 2 */

      norm_r = gsl_blas_dnrm2(state->r);
      u = norm_r / norm_Jd;
      alpha = u * u;

      /* workp <= z_{i+1} = z_i + alpha_i*d_i */
      scaled_addition(1.0, state->z, alpha, state->d, state->workp);

      u = gsl_blas_dnrm2(state->workp);
      if (u >= delta)
        {
          double tau = cgst_calc_tau(state->z, state->d, delta);

          /* dx = z_i + tau*d_i */
          scaled_addition(1.0, state->z, tau, state->d, dx);
          gsl_vector_div(dx, diag);

          return GSL_SUCCESS;
        }

      /* store z_{i+1} */
      gsl_vector_memcpy(state->z, state->workp);

      /* Step 4 of [1], section 2 */

      /* compute: workp := alpha B d_i = alpha D^{-1} J^T J D^{-1} d_i,
       * where J D^{-1} d_i is already stored in workn */
      status = gsl_multilarge_nlinear_eval_df(CblasTrans, x, f, state->workn,
                                              swts, params->h_df, params->fdtype,
                                              fdf, state->workp, NULL, NULL);
      if (status)
        return status;

      gsl_vector_div(state->workp, trust_state->diag);
      gsl_vector_scale(state->workp, alpha);

      /* r_{i+1} = r_i - alpha*B*d_i */
      gsl_vector_sub(state->r, state->workp);
      norm_rp1 = gsl_blas_dnrm2(state->r);

      u = norm_rp1 / state->norm_g;
      if (u < state->cgtol)
        {
          gsl_vector_memcpy(dx, state->z);
          gsl_vector_div(dx, diag);
          return GSL_SUCCESS;
        }

      /* Step 5 of [1], section 2 */

      /* compute u = ||r_{i+1}|| / || r_i|| */
      u = norm_rp1 / norm_r;
      beta = u * u;

      /* compute: d_{i+1} = rt_{i+1} + beta*d_i */
      scaled_addition(1.0, state->r, beta, state->d, state->d);
    }

  /* failed to converge, return current estimate */
  gsl_vector_memcpy(dx, state->z);
  gsl_vector_div(dx, diag);

  return GSL_EMAXITER;
}

static int
cgst_preduction(const void * vtrust_state, const gsl_vector * dx,
                double * pred, void * vstate)
{
  const gsl_multilarge_nlinear_trust_state *trust_state =
    (const gsl_multilarge_nlinear_trust_state *) vtrust_state;
  cgst_state_t *state = (cgst_state_t *) vstate;

  *pred = quadratic_preduction(trust_state, dx, state->workn);

  return GSL_SUCCESS;
}

/*
cgst_calc_tau()
  Compute tau > 0 such that:

|| p + tau*d || = delta
*/

static double
cgst_calc_tau(const gsl_vector * p, const gsl_vector * d,
              const double delta)
{
  double norm_p, norm_d, u;
  double t1, t2, tau;

  norm_p = gsl_blas_dnrm2(p);
  norm_d = gsl_blas_dnrm2(d);

  /* compute (p, d) */
  gsl_blas_ddot(p, d, &u);

  t1 = u / (norm_d * norm_d);
  t2 = t1*u + (delta + norm_p) * (delta - norm_p);
  tau = -t1 + sqrt(t2) / norm_d;

  return tau;
}

static const gsl_multilarge_nlinear_trs cgst_type =
{
  "steihaug-toint",
  cgst_alloc,
  cgst_init,
  cgst_preloop,
  cgst_step,
  cgst_preduction,
  cgst_free
};

const gsl_multilarge_nlinear_trs *gsl_multilarge_nlinear_trs_cgst = &cgst_type;

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#include <vector>
#include <fstream>
#include <cmath>
#include <cblas.h>
#include <fftw3.h>
#include <sys/time.h>
#include <time.h>

using namespace std;

class CochlearModel_1D{

private:
    int N;
    double Lb;
    double L;
    double W;
    double H;
    double rho;
    double dx;
    vector<double> x;

    vector<double> ch_damp;

    vector<double> k1;
    double m1;
    vector<double> c1;
    vector<double> k2;
    double m2;
    vector<double> c2;
    vector<double> k3;
    vector<double> c3;
    vector<double> k4;
    vector<double> c4;

    vector<double> k1k3;
    vector<double> c1c3;
    vector<double> k2k3;
    vector<double> c2c3;

    double g;
    double b;
    vector<double> gamma;

    double dt;
    double beta;

    vector<double> Gohc(vector<double> uc);
    vector<double> dGohc(vector<double> uc, vector<double> vc);
    void get_g(vector<double> vb, vector<double> ub, vector<double> vt, vector<double> ut, vector<double>& gb, vector<double>& gt);

public:
    CochlearModel_1D(int num_segment, vector<double> gamma_given)
    {
        N = num_segment;
        Lb = 3.5;
        L = 0.1;
        W = 0.1;
        H = 0.1;
        rho = 1.0;
        dx = Lb/N;

        m1 = 3e-3;
        m2 = 0.5e-3;
        for (int ii = 0; ii < N; ii++){
            x.push_back(dx * ii);
//            ch_damp.push_back(2.8 * exp(-0.1*x[ii]));
            ch_damp.push_back(1.0);

            k1.push_back(2.2e8 * exp(-3 * x[ii]));
            c1.push_back(6 + 670 * exp(-1.5 * x[ii]) * ch_damp[ii]);
            k2.push_back(1.4e6 * exp(-3.3 * x[ii]));
            c2.push_back(4.4 * exp(-1.65 * x[ii]) * ch_damp[ii]);
            k3.push_back(2.0e6 * exp(-3 * x[ii]));
            c3.push_back(0.8 * exp(-0.6 * x[ii]) * ch_damp[ii]);
            k4.push_back(1.15e8 * exp(-3 * x[ii]));
            c4.push_back(440.0 * exp(-1.5 * x[ii]) * ch_damp[ii]);

            k1k3.push_back(k1[ii] + k3[ii]);
            c1c3.push_back(c1[ii] + c3[ii]);
            k2k3.push_back(k2[ii] + k3[ii]);
            c2c3.push_back(c2[ii] + c3[ii]);
            gamma.push_back(gamma_given[ii]);
        }
        g = 1;
        b = 0.4;
        dt = 10e-6;
        beta = 50e-7;
    }

   void solve_time_domain(vector<double> f, vector<vector<double> >& mat_vb, vector<vector<double> >& mat_ub, vector<vector<double> >& mat_p);
};

max_stars_repo_name: y-muraka/Human-Cochlear-Model

max_issues_repo_name: y-muraka/Human-Cochlear-Model

max_forks_repo_name: y-muraka/Human-Cochlear-Model


text:
/*
FRACTAL - A program growing fractals to benchmark parallelization and drawing
libraries.

Copyright 2009-2021, Javier Burguete Tolosa.

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 Javier Burguete Tolosa ``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 Javier Burguete Tolosa 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.
*/

/**
 * \file draw.c
 * \brief Source file to define the drawing data and functions.
 * \author Javier Burguete Tolosa.
 * \copyright Copyright 2009-2021, Javier Burguete Tolosa.
 */
#define _GNU_SOURCE
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <gsl/gsl_rng.h>
#include <glib.h>
#include <png.h>
#include <ft2build.h>
#include FT_FREETYPE_H
#include <gtk/gtk.h>
#include <GL/glew.h>
#if HAVE_FREEGLUT
#include <GL/freeglut.h>
#elif HAVE_SDL
#include <SDL.h>
#elif HAVE_GLFW
#include <GLFW/glfw3.h>
#endif

#include "config.h"
#include "fractal.h"
#include "image.h"
#include "text.h"
#include "graphic.h"
#include "draw.h"
#include "simulator.h"

Graphic graphic[1];             ///< Graphic data.

/**
 * Function to resize the draw.
 */
void
draw_resize (int w,             ///< Graphic window width.
             int h)             ///< Graphic window height.
{
#if DEBUG
  printf ("draw_resize: start\n");
  fflush (stdout);
#endif
  if ((unsigned int) w < width)
    w = width;
  if ((unsigned int) h < height)
    h = height;
#if HAVE_FREEGLUT
  glutReshapeWindow (w, h);
#endif
  window_width = w;
  window_height = h;
  glViewport (0, 0, w, h);
#if DEBUG
  printf ("draw_resize: end\n");
  fflush (stdout);
#endif
}

#if HAVE_GTKGLAREA
void
resize (GtkGLArea * gl_area __attribute__((unused)), int w, int h)
{
  draw_resize (w, h);
}
#endif

/**
 * Function to render graphics.
 */
void
draw ()
{
#if HAVE_GLFW
  int graphic_width, graphic_height;
#endif

#if DEBUG
  printf ("draw: start\n");
  fflush (stdout);
#endif

#if HAVE_GLFW
  glfwGetFramebufferSize (window, &graphic_width, &graphic_height);
  draw_resize (graphic_width, graphic_height);
#endif

  graphic_render (graphic);

  // Displaying the draw
#if HAVE_GTKGLAREA
  gtk_widget_queue_draw (GTK_WIDGET (dialog_simulator->gl_area));
#elif HAVE_FREEGLUT
  glutSwapBuffers ();
#elif HAVE_SDL
  SDL_GL_SwapWindow (window);
#elif HAVE_GLFW
  glfwSwapBuffers (window);
#endif

#if DEBUG
  printf ("draw: end\n");
  fflush (stdout);
#endif
}

max_stars_repo_name: jburguete/fractal

max_issues_repo_name: jburguete/fractal

max_forks_repo_name: jburguete/fractal


text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <stdarg.h>
#include <gbpLib.h>
#include <gbpMath.h>
#include <gbpCosmo_mass_functions.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>

// Compute comoving mass function
double mass_function(double M_interp, double z, cosmo_info **cosmo, int mode, ...) {
    va_list vargs;
    va_start(vargs, mode);

    // Get passed parameters (if mode specifies they are there).
    double *P = NULL;
    if(SID_CHECK_BITFIELD_SWITCH(mode, MF_PASS_PARAMS))
        P = (double *)va_arg(vargs, double *);

    // Initialize/compute some misc. cosmology things
    double sigma_interp = sqrt(power_spectrum_variance(k_of_M(M_interp, *cosmo), z, cosmo, PSPEC_LINEAR_TF, PSPEC_ALL_MATTER));
    double Omega_M      = ((double *)ADaPS_fetch(*cosmo, "Omega_M"))[0];
    double rho_o        = Omega_M * rho_crit_z(0, *cosmo);

    // Compute mass function (eqn 41 of Lukic et al, 2007, for example)
    double dlnInvs_dlogM = dln_Inv_sigma_dlogM(cosmo, M_interp, PSPEC_LINEAR_TF, PSPEC_ALL_MATTER);
    double rval          = rho_o * scaled_mass_function(sigma_interp, mode, P) * dlnInvs_dlogM / M_interp;
    va_end(vargs);
    return (rval);
}

max_stars_repo_name: gbpoole/gbpCode

max_issues_repo_name: gbpoole/gbpCode

max_forks_repo_name: gbpoole/gbpCode