text:
#pragma once

#include <cstddef>
#include <cstdint>
#include <functional>
#include <initializer_list>
#include <memory>
#include <numeric>
#include <set>
#include <utility>
#include <vector>

#include <absl/container/flat_hash_map.h>
#include <gsl/gsl>

#include "chainerx/array.h"
#include "chainerx/array_body.h"
#include "chainerx/array_node.h"
#include "chainerx/constant.h"
#include "chainerx/device.h"
#include "chainerx/dtype.h"
#include "chainerx/graph.h"
#include "chainerx/macro.h"
#include "chainerx/op_node.h"
#include "chainerx/shape.h"

namespace chainerx {
namespace backward_builder_detail {

// This class is used by the BackwardBuilder to record retained inputs and outputs.
// The records are used to create outer graph edges (between op nodes and previous array nodes) when the builder is finalized.
class RetentionRecord {
public:
    explicit RetentionRecord(size_t size) : size_{size} { CHAINERX_ASSERT(size_ > 0); }

    size_t size() const { return size_; }

    void Record(size_t index) {
        if (flags_.empty()) {
            flags_.resize(size_);
        }
        gsl::at(flags_, index) = static_cast<int8_t>(true);
    }

    bool IsAnyRecorded() const { return !flags_.empty(); }

    bool IsRecorded(size_t index) const { return static_cast<bool>(flags_[index]); }

private:
    size_t size_{};
    std::vector<int8_t> flags_{};  // binary flags
};

template <typename Tag>
class RetainedArrayToken {
public:
    RetainedArrayToken(internal::ArrayBody::Params array_params, size_t index) : array_params_{std::move(array_params)}, index_{index} {}

    ~RetainedArrayToken() = default;

    RetainedArrayToken(const RetainedArrayToken&) = default;
    RetainedArrayToken(RetainedArrayToken&&) noexcept = default;
    RetainedArrayToken& operator=(const RetainedArrayToken&) = default;
    // TODO(hvy): Make the move assignment operator noexcept.
    RetainedArrayToken& operator=(RetainedArrayToken&&) = default;  // NOLINT(performance-noexcept-move-constructor)

private:
    friend class chainerx::BackwardContext;

    // Returns the array index.
    size_t index() const { return index_; }

    const internal::ArrayBody::Params& array_params() const { return array_params_; }

    internal::ArrayBody::Params array_params_;

    size_t index_;
};

}  // namespace backward_builder_detail

// An object used by op implementations to bridge between BackwardBuilder::RetainInput() and BackwardContext::GetRetainedInput().
//
// See BackwardBuilder::RetainInput() for details.
using RetainedInputToken = backward_builder_detail::RetainedArrayToken<struct InputTag>;

// An object used by op implementations to bridge between BackwardBuilder::RetainOutput() and BackwardContext::GetRetainedOutput().
//
// See BackwardBuilder::RetainOutput() for details.
using RetainedOutputToken = backward_builder_detail::RetainedArrayToken<struct OutputTag>;

// A class that is used to define backward operations and connect the graph.
//
// This class is not thread safe.
class BackwardBuilder {
public:
    // Target is responsible to define edges from OpNode to input ArrayNodes with given BackwardFunction.
    // Note that Targets built from the same BackwardBuilder share some properties not to compute again.
    class Target {
    public:
        explicit operator bool() const { return is_definition_required(); }

        // Defines a backward function with respect to specified input arrays (target).
        void Define(const BackwardFunction& backward_func);

        bool is_definition_required() const { return !graph_to_input_array_nodes_.empty(); }

    private:
        friend class BackwardBuilder;  // Only BackwardBuilder can create Target

        using InputArrayNodes = std::vector<const std::shared_ptr<internal::ArrayNode>*>;

        Target(BackwardBuilder& builder, std::vector<size_t> input_indices);

        // Collect input ArrayNodes, grouped by graph considering IsBackpropRequired.
        // This functions is only called once in the constructor.
        absl::flat_hash_map<BackpropId, InputArrayNodes> CreateInputArrayNodesMap() const;

        BackwardBuilder& builder_;
        std::vector<size_t> input_indices_;

        // TODO(hvy): Consider using linear search since elements are usually few.
        absl::flat_hash_map<BackpropId, InputArrayNodes> graph_to_input_array_nodes_;
    };

    // TODO(niboshi): Add an overload to accept `const std::vector<Array>&` as `inputs` and `outputs`
    // Note that simply overloading with the above type will results in ambiguous calls.
    // One solution is to define a type that accepts all of the expected types of inputs.
    BackwardBuilder(const char* op_name, std::vector<ConstArrayRef> inputs, std::vector<ConstArrayRef> outputs);
    BackwardBuilder(const char* op_name, const Array& input, std::vector<ConstArrayRef> outputs)
        : BackwardBuilder{op_name, std::vector<ConstArrayRef>{input}, std::move(outputs)} {}
    BackwardBuilder(const char* op_name, std::vector<ConstArrayRef> inputs, const Array& output)
        : BackwardBuilder{op_name, std::move(inputs), std::vector<ConstArrayRef>{output}} {}
    BackwardBuilder(const char* op_name, const Array& input, const Array& output)
        : BackwardBuilder{op_name, std::vector<ConstArrayRef>{input}, std::vector<ConstArrayRef>{output}} {}
    ~BackwardBuilder() { CHAINERX_ASSERT(is_finalized_); }

    BackwardBuilder(const BackwardBuilder&) = delete;
    BackwardBuilder(BackwardBuilder&&) noexcept = default;
    BackwardBuilder& operator=(const BackwardBuilder&) = delete;
    BackwardBuilder& operator=(BackwardBuilder&&) = delete;

    // Creates a backward target for the specified inputs.
    Target CreateTarget(std::vector<size_t> input_indices) {
        // input_indices shouldn't have duplicates.
        CHAINERX_ASSERT((std::set<size_t>{input_indices.begin(), input_indices.end()}.size() == input_indices.size()));

        for (size_t input_index : input_indices) {
            CHAINERX_ASSERT(input_index < inputs_target_created_.size());
            CHAINERX_ASSERT(!inputs_target_created_[input_index]);
            inputs_target_created_[input_index] = true;
        }
        return Target{*this, std::move(input_indices)};
    }

    // Creates a backward target for the specified input.
    Target CreateTarget(size_t input_index) { return CreateTarget(std::vector<size_t>{input_index}); }

    // Creates a backward target for all the inputs.
    Target CreateTarget() {
        std::vector<size_t> input_indices;
        input_indices.resize(inputs_.size());
        std::iota(input_indices.begin(), input_indices.end(), size_t{0});

        return CreateTarget(std::move(input_indices));
    }

    // TODO(hvy): Write comment.
    RetainedInputToken RetainInput(size_t input_index);

    std::vector<RetainedInputToken> RetainInput(std::vector<size_t> indices);

    // Flags an output array to be retained for use in the backward pass.
    // Op implementations can use this function in combination with BackwardContext::GetRetainedOutput() to retrieve output arrays in the
    // backward pass.
    //
    // If an op implementation requires the output array of the forward pass in the backward pass, it should call
    // BackwardBuilder::RetainOutput() in the forward pass and keep its return value (either assign a variable or capture by
    // value in a lambda expression). In the backward pass, it should call BackwardContext::GetRetainedOutput() with this token to retrieve
    // the output array.
    //
    // Capturing the output array directly with lambda expression would cause cyclic reference and therefore would lead to memory leak.
    //
    // Reusing the token for higher-order backward functions results in undefined behavior.
    //
    // `output` must be one of the arrays specified in the constructor of BackwardBuilder as output arrays.
    // If invalid array is specified, ChainerxError will be thrown.
    RetainedOutputToken RetainOutput(size_t output_index);
    std::vector<RetainedOutputToken> RetainOutput(std::vector<size_t> indices);

    // Finalizes the builder.
    //
    // This functions must be called when targets have been created for all inputs.
    void Finalize();

private:
    // Create an op node for a specific graph.
    // Edges from output nodes to the op node are connected.
    std::shared_ptr<internal::OpNode>& FindOrCreateOpNode(const BackpropId& backprop_id);

    // Add shared ptrs between op nodes and array nodes belonging to outer graphs.
    // This functions is called once when the builder is finalized.
    // These references are required to restore retained inputs/outputs.
    void AddEdgesFromOpNodeToArrayNodeOfOuterGraphsForRetention();

    void ConnectBackpropIds();

    const char* op_name_;

    Context& context_;

    // Input arrays of the op.
    std::vector<ConstArrayRef> inputs_;

    // Flags indicating whether CreateTarget has been called for each of the input arrays.
    // All of these flags must be true after all the backwards have been defined for a BackwardBuilder.
    // This can be checked by calling is_complete();
    std::vector<bool> inputs_target_created_;

    // Output arrays of the op.
    std::vector<ConstArrayRef> outputs_;

    // A collection of op nodes, each of which corresponds to a graph.
    // This record is increasingly populated as new graphs are encountered in multiple Define() calls.
    absl::flat_hash_map<BackpropId, std::shared_ptr<internal::OpNode>> op_node_map_;

    backward_builder_detail::RetentionRecord input_retention_record_;
    backward_builder_detail::RetentionRecord output_retention_record_;

    bool has_any_applicable_outputs_;
    bool is_finalized_{false};
};

}  // namespace chainerx

max_stars_repo_name: zjzh/chainer

max_issues_repo_name: zjzh/chainer

max_forks_repo_name: zjzh/chainer


text:
// Copyright (c) Microsoft Corporation. All rights reserved.
// Licensed under the MIT License.

#pragma once

#include <gsl/gsl>

#include "core/optimizer/graph_transformer.h"
#include "core/optimizer/rule_based_graph_transformer.h"
#include "core/optimizer/rewrite_rule.h"

namespace onnxruntime {
struct FreeDimensionOverride;

namespace transformer_utils {

/** Generates all predefined rules for this level.
   If rules_to_enable is not empty, it returns the intersection of predefined rules and rules_to_enable.
   TODO: This is visible for testing at the moment, but we should rather make it private. */
std::vector<std::unique_ptr<RewriteRule>> GenerateRewriteRules(TransformerLevel level,
                                                               const std::vector<std::string>& rules_to_enable = {});

/** Generates all predefined (both rule-based and non-rule-based) transformers for this level.
    If transformers_and_rules_to_enable is not empty, it returns the intersection between the predefined transformers/rules
	and the transformers_and_rules_to_enable. */
std::vector<std::unique_ptr<GraphTransformer>> GenerateTransformers(TransformerLevel level,
                                                                    gsl::span<const FreeDimensionOverride> free_dimension_overrides,
                                                                    const std::vector<std::string>& rules_and_transformers_to_enable = {});

/** Given a TransformerLevel, this method generates a name for the rule-based graph transformer of that level. */
std::string GenerateRuleBasedTransformerName(TransformerLevel level);

}  // namespace transformer_utils
}  // namespace onnxruntime

ext: h

max_stars_repo_name: csteegz/onnxruntime


text:
#ifdef HAVE_LIBGSL
/* interface_gsl.h
 * Easel's interfaces to the GNU Scientific Library
 *
 * SRE, Tue Jul 13 15:36:48 2004
 * SVN $Id: interface_gsl.h 664 2011-02-27 17:08:36Z eddys $
 * SVN $URL: https://svn.janelia.org/eddylab/eddys/easel/branches/hmmer/3.1/interface_gsl.h $
 */
#ifndef eslINTERFACE_GSL_INCLUDED
#define eslINTERFACE_GSL_INCLUDED

#include <stdlib.h>
#include <easel/easel.h>
#include <easel/dmatrix.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_eigen.h>

extern int esl_GSL_MatrixInversion(ESL_DMATRIX *A, ESL_DMATRIX **ret_Ai);

#endif /*eslINTERFACE_GSL_INCLUDED*/
#endif /*HAVE_LIBGSL*/


text:
#ifndef _EM_CTX_H_
#define _EM_CTX_H_ 1

#include "em_defs.h"

#include <petsc.h>

#include <set>
#include <map>
#include <memory>
#include <vector>

class Mesh;

struct EMContext {
  MPI_Comm world_comm, group_comm;
  int world_size, world_rank, group_size, group_rank, group_id;

  std::vector<Point> rx;
  std::vector<double> tx;
  std::vector<double> freqs;

  Eigen::VectorXcd rsp, obs, obserr;
  Eigen::VectorXi otype, fidx, tidx, ridx;

  Point top_corners[4];
  Eigen::VectorXd ztop[4], lsig[4];

  std::shared_ptr<Mesh> original_mesh, mesh;

  std::vector<int> relevant_edges;
  std::pair<int, int> local_vertices, local_edges;

  std::set<int> bdr_cells;

  int aniso_form;
  Eigen::MatrixXd rho;
  Eigen::VectorXd lb, ub;
  std::vector<std::string> rho_name;

  Vec w;
  Mat C, M, A, B;
  KSP A_ksp, B_ksp;
  PETScBlockVector s, csem_e, dual_e;
  std::map<int, PETScBlockVector> mt_e;

  PetscViewer LS_log;

  Eigen::VectorXd csem_error;
  std::map<int, Eigen::VectorXd> mt_error;

  Mat G;
  std::vector<PetscReal> v_coords;

  PetscBool use_ams;
  char iprefix[256], oprefix[256];
  PetscReal max_rx_edge_length, refine_fraction, e_rtol, dual_rtol;
  PetscInt max_adaptive_refinements, max_dofs, refine_strategy, n_groups, K_max_it, pc_threshold, inner_pc_type, direct_solver_type;

  PetscClassId EMCTX_ID;
  PetscLogEvent CreateLS, AssembleMat, AssembleRHS, SetupAMS, CreatePC, SolveLS, EstimateError,
      RefineMesh, CalculateRSP;
};

PetscErrorCode create_context(EMContext *);
PetscErrorCode destroy_context(EMContext *);
PetscErrorCode process_options(EMContext *);

#endif

max_stars_repo_name: emfem/emfem

max_issues_repo_name: emfem/emfem

max_forks_repo_name: emfem/emfem


text:
//          Copyright Jean Pierre Cimalando 2019.
// Distributed under the Boost Software License, Version 1.0.
//    (See accompanying file LICENSE.md or copy at
//          http://www.boost.org/LICENSE_1_0.txt)

#pragma once
#include <gsl/gsl>
#include <string>

std::string get_home_directory();
std::string get_current_directory();
bool is_path_absolute(gsl::cstring_span path);
bool is_path_separator(char32_t character);
void append_path_separator(std::string &path);
std::string normalize_path_separators(gsl::cstring_span path);
std::string make_path_canonical(gsl::cstring_span path);
std::string expand_path_tilde(gsl::cstring_span path);
gsl::cstring_span path_file_name(gsl::cstring_span path);
gsl::cstring_span path_directory(gsl::cstring_span path);

std::string get_display_path(gsl::cstring_span path);

#if defined(_WIN32)
std::string known_folder_path(int csidl, std::error_code &ec);
std::string known_folder_path(int csidl);
#endif

max_stars_repo_name: jpcima/smf-dsp

max_issues_repo_name: jpcima/smf-dsp

max_forks_repo_name: jpcima/smf-dsp


text:
#pragma once

#include "BoundingVolume.h"
#include "Library.h"
#include "Tile.h"
#include "TileContext.h"
#include "TileID.h"
#include "TileRefine.h"
#include "TilesetOptions.h"

#include <CesiumAsync/AsyncSystem.h>
#include <CesiumAsync/IAssetAccessor.h>
#include <CesiumAsync/IAssetRequest.h>

#include <gsl/span>
#include <spdlog/fwd.h>

#include <cstddef>
#include <memory>

namespace Cesium3DTilesSelection {
/**
 * @brief The information that is passed to a {@link TileContentLoader} to
 * create a {@link TileContentLoadResult}.
 *
 * For many types of tile content, only the `pRequest` field is required. The
 * other members are used for content that can generate child tiles, like
 * external tilesets or composite tiles. These members are usually initialized
 * from
 * the corresponding members of the {@link Tile} that the content belongs to.
 */
struct CESIUM3DTILESSELECTION_API TileContentLoadInput {

  /**
   * @brief Creates a new, uninitialized instance for the given tile.
   *
   * The `asyncSystem`, `pLogger`, `pAssetAccessor` and `pRequest` will have
   * default values, and have to be initialized before this instance is passed
   * to one of the loader functions.
   *
   * @param tile The {@link Tile} that the content belongs to
   */
  TileContentLoadInput(const Tile& tile);

  /**
   * @brief Creates a new instance
   *
   * @param asyncSystem The async system to use for tile content loading.
   * @param pLogger The logger that will be used
   * @param pAssetAccessor The asset accessor to make further requests with.
   * @param pRequest The original tile request and its response.
   * @param tile The {@link Tile} that the content belongs to.
   */
  TileContentLoadInput(
      const CesiumAsync::AsyncSystem& asyncSystem,
      const std::shared_ptr<spdlog::logger>& pLogger,
      const std::shared_ptr<CesiumAsync::IAssetAccessor>& pAssetAccessor,
      const std::shared_ptr<CesiumAsync::IAssetRequest>& pRequest,
      const Tile& tile);

  /**
   * @brief Creates a new instance.
   *
   * @param asyncSystem The async system to use for tile content loading.
   * @param pLogger The logger that will be used
   * @param pAssetAccessor The asset accessor to make further requests with.
   * @param pRequest The original tile request and its response.
   * @param tileID The {@link TileID}
   * @param tileBoundingVolume The tile {@link BoundingVolume}
   * @param tileContentBoundingVolume The tile content {@link BoundingVolume}
   * @param tileRefine The {@link TileRefine} strategy
   * @param tileGeometricError The geometric error of the tile
   * @param tileTransform The tile transform
   */
  TileContentLoadInput(
      const CesiumAsync::AsyncSystem& asyncSystem,
      const std::shared_ptr<spdlog::logger>& pLogger,
      const std::shared_ptr<CesiumAsync::IAssetAccessor>& pAssetAccessor,
      const std::shared_ptr<CesiumAsync::IAssetRequest>& pRequest,
      const TileID& tileID,
      const BoundingVolume& tileBoundingVolume,
      const std::optional<BoundingVolume>& tileContentBoundingVolume,
      TileRefine tileRefine,
      double tileGeometricError,
      const glm::dmat4& tileTransform,
      const TilesetContentOptions& contentOptions);

  /**
   * @brief The async system to use for tile content loading.
   */
  CesiumAsync::AsyncSystem asyncSystem;

  /**
   * @brief The logger that receives details of loading errors and warnings.
   */
  std::shared_ptr<spdlog::logger> pLogger;

  /**
   * @brief The asset accessor to make further requests with.
   */
  std::shared_ptr<CesiumAsync::IAssetAccessor> pAssetAccessor;

  /**
   * @brief The asset request and response data for the tile.
   */
  std::shared_ptr<CesiumAsync::IAssetRequest> pRequest;

  /**
   * @brief The {@link TileID}.
   */
  TileID tileID;

  /**
   * @brief The tile {@link BoundingVolume}.
   */
  BoundingVolume tileBoundingVolume;

  /**
   * @brief Tile content {@link BoundingVolume}.
   */
  std::optional<BoundingVolume> tileContentBoundingVolume;

  /**
   * @brief The {@link TileRefine}.
   */
  TileRefine tileRefine;

  /**
   * @brief The geometric error.
   */
  double tileGeometricError;

  /**
   * @brief The tile transform
   */
  glm::dmat4 tileTransform;

  /**
   * @brief Options for parsing content and creating Gltf models.
   */
  TilesetContentOptions contentOptions;
};
} // namespace Cesium3DTilesSelection

max_stars_repo_name: JiangMuWen/cesium-native

max_issues_repo_name: JiangMuWen/cesium-native

max_forks_repo_name: JiangMuWen/cesium-native


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

int main(int argc, char **argv)
{
	if(argc != 3) {
		fprintf(stderr, "%s <database name> <rows>\n", argv[0]);
		exit(1);
	}

	unlink(argv[1]);
	virginian virg;
	virginian *v = &virg;

	virg_init(v);
	virg_db_create(v, argv[1]);
	virg_table_create(v, "test", VIRG_INT);
	virg_table_addcolumn(v, 0, "uniformi", VIRG_INT);
	virg_table_addcolumn(v, 0, "normali5", VIRG_INT);
	virg_table_addcolumn(v, 0, "normali20", VIRG_INT);
	virg_table_addcolumn(v, 0, "uniformf", VIRG_FLOAT);
	virg_table_addcolumn(v, 0, "normalf5", VIRG_FLOAT);
	virg_table_addcolumn(v, 0, "normalf20", VIRG_FLOAT);

	int i;
	void *buff = malloc(4 * 6);
	int *buff_i = (int*)buff;
	float *buff_f = (float*)buff;

	const gsl_rng_type *type;
	gsl_rng *ran;

	type = gsl_rng_ranlxd2;
	ran = gsl_rng_alloc(type);
	gsl_rng_set(ran, time(0));

	for(i = 0; i < atoi(argv[2]); i++) {
		buff_i[0] = (int)gsl_ran_flat(ran, -100, 100);
		buff_i[1] = (int)gsl_ran_gaussian(ran, 5);
		buff_i[2] = (int)gsl_ran_gaussian(ran, 20);
		buff_f[3] = (float)gsl_ran_flat(ran, -100, 100);
		buff_f[4] = (float)gsl_ran_gaussian(ran, 5);
		buff_f[5] = (float)gsl_ran_gaussian(ran, 20);
		virg_table_insert(v, 0, (char*)&i, buff, NULL);

		if(i % 10000 == 0) {
			printf("%i,", i);
			fflush(stdout);
		}
	}
	printf("\n");

	free(buff);

	virg_db_close(v);
	virg_close(v);
	gsl_rng_free(ran);

	return 0;
}

ext: c

max_stars_repo_name: bakks/virginian

max_issues_repo_name: krunalpuri/virginian

max_forks_repo_name: krunalpuri/virginian


text:
#ifndef __INTEGRATOR_H
#define __INTEGRATOR_H

#include "autoparams.h"
#include <gsl/gsl_integration.h>
#include <gsl/gsl_errno.h>
#include <stack>
namespace integrator
{

double infty = 1000;

class WorkspaceHandle;

class Workspaces
{
public:
    static const size_t WORKSPACE_SIZE = 1000000;
    static const size_t NUM_INITIAL_WORKSPACES = 4;

    static WorkspaceHandle get_workspace_handle();

    virtual ~Workspaces()
    {
        while (!workspace_stack_.empty())
        {
            gsl_integration_workspace_free(workspace_stack_.top());
            workspace_stack_.pop();
        }
    }
private:
    Workspaces(): // private constructor
        workspace_stack_()
    {
        for (int i = 0; i < NUM_INITIAL_WORKSPACES; i++)
            workspace_stack_.push(
                    gsl_integration_workspace_alloc(WORKSPACE_SIZE));
    }
    Workspaces(Workspaces const& copy);            // not impl
    Workspaces& operator=(Workspaces const& copy); // not impl
    std::stack<gsl_integration_workspace*> workspace_stack_;
};

class WorkspaceHandle
{
public:
    WorkspaceHandle(std::stack<gsl_integration_workspace*>* workspaces):
        workspaces_(workspaces), ws_(0)
    {
        if (!workspaces->empty())
        {
            ws_ = workspaces_->top();
            workspaces_->pop();
        }
        else ws_ = gsl_integration_workspace_alloc(Workspaces::WORKSPACE_SIZE);

    }

    virtual ~WorkspaceHandle()
    {
        workspaces_->push(ws_);
    }

    gsl_integration_workspace* get()
    {
        return ws_;
    }
private:
    std::stack<gsl_integration_workspace*>* workspaces_;
    gsl_integration_workspace* ws_;
};

WorkspaceHandle Workspaces::get_workspace_handle()
{
        static Workspaces workspaces;
        return WorkspaceHandle(&workspaces.workspace_stack_);
}

struct WrapperInfo
{
    double l;
    double r;
    double ts;
    gsl_function fu;
};

double wrapper_func(const double x, const void* p)
{
    const WrapperInfo* wi = static_cast<const WrapperInfo*>(p);
    const double l = wi->l;
    const double r = wi->r;
    const double ts = wi->ts;
    return (r-l)/2 * 1./atan(ts) * 1./(1.+x*x) *
        wi->fu.function((r-l)/2*atan(x)/atan(ts)+(l+r)/2, wi->fu.params);
}

typedef double (*func_with_nonconst_args)(double, void*);

double integrate(double (*func)(const double, const void*), const double a, const double b,
        const double eps, const void* const p=0)
{
    WorkspaceHandle w = Workspaces::get_workspace_handle();
    double result = 0;
    double error = 0;
    gsl_function fu = {(func_with_nonconst_args)func, const_cast<void*>(p)};
    gsl_integration_qags(&fu, a, b, params::integ_epsabs, eps, Workspaces::WORKSPACE_SIZE,
        w.get(), &result, &error);
    return result;
}

double integrate_inf(double (*func)(double, const void*),
        const double eps, const void* const p=0)
{
    WorkspaceHandle w = Workspaces::get_workspace_handle();
    double result = 0;
    double error = 0;
    gsl_function fu = {(func_with_nonconst_args)func, const_cast<void*>(p)};
    gsl_integration_qagi(&fu, params::integ_epsabs, eps, Workspaces::WORKSPACE_SIZE,
            w.get(), &result, &error);
    return result;
}

double integrate_infu(double (*func)(double, const void*), const double a,
        const double eps, const void* const p=0)
{
    WorkspaceHandle w = Workspaces::get_workspace_handle();
    double result = 0;
    double error = 0;
    gsl_function fu = {(func_with_nonconst_args)func, const_cast<void*>(p)};
    gsl_integration_qagiu(&fu, a, params::integ_epsabs, eps, Workspaces::WORKSPACE_SIZE,
            w.get(), &result, &error);
    return result;
}

double integrate_infl(double (*func)(double, const void*), const double b,
        const double eps, const void* const p=0)
{
    WorkspaceHandle w = Workspaces::get_workspace_handle();
    double result = 0;
    double error = 0;
    gsl_function fu = {(func_with_nonconst_args)func, const_cast<void*>(p)};
    gsl_integration_qagil(&fu, b, params::integ_epsabs, eps, Workspaces::WORKSPACE_SIZE,
            w.get(), &result, &error);
    return result;
}

double edge_emph_integrate(double (*func)(double, const void*), const double l,
        const double r, const double ts, const double a, const double b,
        const double eps, const void* const p=0)
{
    WorkspaceHandle w = Workspaces::get_workspace_handle();
    double result = 0;
    double error = 0;
    gsl_function fu = {(func_with_nonconst_args)func, const_cast<void*>(p)};
    WrapperInfo wi = {l, r, ts, fu};
    gsl_function wra_fu = {(func_with_nonconst_args)wrapper_func, &wi};
    gsl_integration_qags(&wra_fu,
            tan(2./(r-l)*(a-(l+r)/2)*atan(ts)),
            tan(2./(r-l)*(b-(l+r)/2)*atan(ts)),
            params::integ_epsabs, eps, Workspaces::WORKSPACE_SIZE,
            //GSL_INTEG_GAUSS21,
            w.get(), &result, &error);

    return result;
}

};

#endif // __INTEGRATOR_H