text:
#pragma once

#ifndef _NOGSL
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>
#else
#include "CustomSolver.h"
#endif

#include "BooleanNodes.h"
#include "BooleanDAG.h"
#include "NodeVisitor.h"
#include <map>
#include "SymbolicEvaluator.h"
#include "Util.h"

#include <iostream>

using namespace std;

// Goes through the DAG and figures out where the gradient vanishes
class GradientAnalyzer: public NodeVisitor
{
	BooleanDAG& bdag;
    SymbolicEvaluator* eval;

    double threshold = 0.01;

    double upThreshold = 1000;

    bool printVals = true;
    bool printGrads = false;
    bool printGradLoss = false;
    bool printGradBlowup = false;
    bool printSqrt = false;

public:

    GradientAnalyzer(BooleanDAG& bdag_p, SymbolicEvaluator* eval_p): bdag(bdag_p), eval(eval_p) {}
    ~GradientAnalyzer(void) {}

	virtual void visit( SRC_node& node );
	virtual void visit( DST_node& node );
	virtual void visit( CTRL_node& node );
	virtual void visit( PLUS_node& node );
	virtual void visit( TIMES_node& node );
	virtual void visit( ARRACC_node& node );
	virtual void visit( DIV_node& node );
	virtual void visit( MOD_node& node );
	virtual void visit( NEG_node& node );
	virtual void visit( CONST_node& node );
	virtual void visit( LT_node& node );
	virtual void visit( EQ_node& node );
	virtual void visit( AND_node& node );
	virtual void visit( OR_node& node );
	virtual void visit( NOT_node& node );
	virtual void visit( ARRASS_node& node );
	virtual void visit( UFUN_node& node );
	virtual void visit( TUPLE_R_node& node );
	virtual void visit( ASSERT_node& node );

	virtual void run();
    virtual void run(bool_node* n);

};

max_stars_repo_name: natebragg/sketch-backend

max_issues_repo_name: natebragg/sketch-backend

max_forks_repo_name: natebragg/sketch-backend


text:
#include <gsl/gsl>

#include <functional>
#include <mutex>
#include <condition_variable>
#include <future>

#include <thread>

namespace lom {

using executor_t = void (*)(std::function<void()>);
auto default_executor_impl = [](std::function<void()> f) {
    std::thread(std::move(f)).detach();
};

extern executor_t g_default_executor;

//
//
/**
 * used to make the library header only. You'll need to add the following in 1 file in an app to use the continuation feature
 * ```
 * #define LOM_FUTURE_IMPLEMENT
 * #include "lom/future.h"
 * ```
 */
#ifdef LOM_FUTURE_IMPLEMENT
executor_t g_default_executor = default_executor_impl;
#endif

template<class TYPE>
class state : public std::enable_shared_from_this<state<TYPE>> {
    union {
        TYPE value_;
        std::exception_ptr exception_;
    };

    enum {
        NONE,
        HAS_VALUE,
        HAS_EXCEPTION,
    } status_{NONE};

    mutable std::mutex mutex_;
    mutable std::condition_variable cv_;
    std::function<void(std::shared_ptr<state>)> continuation_;

public:
    void reset() {
        switch (status_) {
            case NONE:
                break;
            case HAS_VALUE:
                value_.TYPE::~TYPE();
                break;
            case HAS_EXCEPTION:
                exception_.std::exception_ptr::~exception_ptr();
                break;
        }

        status_ = NONE;
    }

    // for make_ready_future
    template<class ...ARGS>
    state(ARGS &&...args) {
        status_ = HAS_VALUE;
        new(&value_) TYPE(std::forward<ARGS>(args)...);
    }

    state() {}

    ~state() {
        reset();
    }

    template<class ...ARGS>
    void set_value(ARGS &&...args) {
        if (status_ != NONE)
            throw std::future_error{std::future_errc::promise_already_satisfied};

        {
            std::lock_guard<std::mutex> lock(mutex_);
            status_ = HAS_VALUE;
            new(&value_) TYPE(std::forward<ARGS>(args)...);
        }
        cv_.notify_one();

        if (continuation_) {
            continuation_(this->shared_from_this());
        }
    }

    void set_exception(std::exception_ptr p) {
        if (status_ != NONE)
            throw std::future_error{std::future_errc::promise_already_satisfied};

        {
            std::lock_guard<std::mutex> lock(mutex_);
            status_ = HAS_EXCEPTION;
            new(&exception_) std::exception_ptr(p);
        }
        cv_.notify_one();

        if (continuation_)
            continuation_(this->shared_from_this());
    }

    TYPE get() {
        std::unique_lock<decltype(mutex_)> lock{mutex_};

        cv_.wait(lock, [this]() {
            return this->status_ != NONE;
        });

        if (status_ == HAS_EXCEPTION)
            std::rethrow_exception(exception_);

        return std::move(value_);
    }

    std::future_status wait() const {
        std::unique_lock<decltype(mutex_)> lock{mutex_};

        cv_.wait(lock, [this]() {
            return this->status_ != NONE;
        });

        return std::future_status::ready;
    }

    template<class Rep, class Period>
    std::future_status wait_for(const std::chrono::duration<Rep, Period> &rel_time) const {
        std::unique_lock<decltype(mutex_)> lock{mutex_};

        bool time_out = cv_.wait_for(lock, rel_time, [this]() {
            return this->status_ != NONE;
        });

        if (time_out)
            return std::future_status::timeout;

        return std::future_status::ready;
    }

    template<class Clock, class Duration>
    std::future_status wait_until(const std::chrono::time_point<Clock, Duration> &abs_time) const {
        std::unique_lock<decltype(mutex_)> lock{mutex_};

        bool time_out = cv_.wait_until(lock, abs_time, [this]() {
            return this->status_ != NONE;
        });

        if (time_out)
            return std::future_status::timeout;

        return std::future_status::ready;
    }

    void then(std::function<void(std::shared_ptr<state>)> &&func) {
        {
            std::lock_guard<std::mutex> lock(mutex_);
            if (status_ == NONE) {
                continuation_ = std::move(func);
                return;
            }
        }

        // don't need to lock cause the value is already set and it will through an exception if you try to set it again
        //
        func(this->shared_from_this());
    }
};

template<class TYPE>
class future;

template<class TYPE, class ...ARGS>
future<TYPE> create_ready_future(ARGS &&... args);

template<class future_t>
class future {
public:
    using type = future_t;

private:
    template<class ...ARGS>
    friend future<future_t> create_ready_future(ARGS &&... args);

    std::shared_ptr<state<future_t>> state_;

private:
    template<typename> struct is_future : std::false_type {};
    template<typename T> struct is_future<future<T>> : std::true_type {};
public:

    template<class ...ARGS>
    future(ARGS &&... args): state_(std::make_shared<future_t>(std::forward<ARGS>(args)...)) {}

    future() {}

    future(const future &) = delete;

    future &operator=(const future &) = default;

    future(future &&p) : state_(std::move(p.state_)) {
        p.state_ = nullptr;
    }

    future &operator=(future &&p) {
        this->state_ = std::move(p.state_);
        p.state_ = nullptr;
    }

    future(std::shared_ptr<state<future_t>> state) : state_(std::move(state)) {}

    future_t get() {
        if (!state_)
            throw std::future_error{std::future_errc::no_state};

        auto rm_state = gsl::finally([this]() {
            this->state_->reset();
            this->state_ = nullptr;
        });

        return state_->get();
    }

    std::future_status wait() {
        if (!state_)
            throw std::future_error{std::future_errc::no_state};

        return state_->wait();
    }

    template<class Rep, class Period>
    std::future_status wait_for(const std::chrono::duration<Rep, Period> &rel_time) const {
        if (!state_)
            throw std::future_error{std::future_errc::no_state};

        return state_->wait_for(rel_time);
    }

    template<class Clock, class Duration>
    std::future_status wait_until(const std::chrono::time_point<Clock, Duration> &abs_time) const {
        if (!state_)
            throw std::future_error{std::future_errc::no_state};

        return state_->wait_until(abs_time);
    }

    bool valid() {
        return state_ != nullptr;
    }

    //function that returns void
    template <class callback_t,
            class return_t = typename std::result_of<callback_t(future_t)>::type,
            typename  = typename std::enable_if<std::is_void<return_t>::value >::type
    >
    void then(callback_t cb_function){
        if (!state_)
            throw std::future_error{std::future_errc::no_state};

        state_->then([cb_function = std::move(cb_function)](std::shared_ptr<state<future_t>> state) mutable {
            g_default_executor([cb_function = std::move(cb_function), state = std::move(state)]() mutable {
                cb_function(future<future_t> {state});
            });
        });
    }

    //function that returns a non future value
    template <class callback_t,
            class inner_return_t = typename std::result_of<callback_t(future_t)>::type,
            class return_t = future<inner_return_t>,
            typename  = typename std::enable_if<!std::is_void<inner_return_t>::value >::type,
            typename  = typename std::enable_if<!is_future<inner_return_t>::value>::type
    >
    return_t then(callback_t cb_function) {
        if (!state_)
            throw std::future_error{std::future_errc::no_state};

        auto ret_state = std::make_shared<state<inner_return_t>>();

        state_->then([cb_function = std::move(cb_function), ret_state](std::shared_ptr<state<future_t>> state) mutable {
            g_default_executor([cb_function = std::move(cb_function), state = std::move(state), ret_state = std::move(ret_state)]() mutable {
                ret_state->set_value(cb_function(future<future_t> {state}));

            });
        });

        return return_t {ret_state};
    }

    // function that returns a lom::future
    template <class callback_t,
            class return_t = typename std::result_of<callback_t(future_t)>::type,
            class inner_return_t = typename return_t::type,
            typename  = typename std::enable_if<is_future<return_t>::value>::type
    >
    return_t then(callback_t cb_function) {
        static_assert(std::is_same<return_t, lom::future<double>>::value, "qwert");
        static_assert(std::is_same<inner_return_t,double>::value, "qewr");

        if (!state_)
            throw std::future_error{std::future_errc::no_state};

        auto ret_state = std::make_shared<state<inner_return_t>>();

        state_->then([cb_function = std::move(cb_function), ret_state](std::shared_ptr<state<future_t>> state) mutable {
            g_default_executor([cb_function = std::move(cb_function), state = std::move(state), ret_state = std::move(ret_state)]() mutable {
                return_t rt = cb_function(future<future_t> {state});
                rt.then([ret_state = std::move(ret_state)](return_t r){
                    ret_state->set_value(r.get());
                });
                //ret_state->set_value

            });
        });

        return return_t {ret_state};
    }
};

template<class TYPE, class ...ARGS>
future<TYPE> create_ready_future(ARGS &&... args) {
    return {std::forward<ARGS>(args)...};
}

template<class TYPE>
class promise {
    std::shared_ptr<state<TYPE>> state_;
    bool has_future = false;

public:
    promise() : state_(std::make_shared<state<TYPE>>()) {}

    promise(const promise &) = delete;

    promise &operator=(const promise &) = delete;

    promise(promise &&p) : state_(std::move(p.state_)) {}

    promise &operator=(promise &&p) noexcept {
        this->state_ = std::move(p.state_);
    }

    future<TYPE> get_future() {
        if(has_future)
            throw std::future_error{std::future_errc::future_already_retrieved};

        has_future = true;

        return {state_};
    }

    template<class ...ARGS>
    void set_value(ARGS &&...args) {
        if (!state_)
            throw std::future_error{std::future_errc::no_state};

        state_->set_value(std::forward<ARGS>(args)...);
    }

    void set_exception(std::exception_ptr p) {
        if (!state_)
            throw std::future_error{std::future_errc::no_state};

        state_->set_exception(std::move(p));
    }

};

/*
template<class future_t>
template<class callback_t, class return_t>
auto future<future_t>::then(callback_t cb_function) -> return_t {
    static_assert(!std::is_same<return_t, void>::value, "is null");

    if (!state_)
        throw std::future_error{std::future_errc::no_state};

    state_->then([cb_function = std::move(cb_function)](std::shared_ptr<state<future_t>> state) mutable {
        g_default_executor([cb_function = std::move(cb_function), state = std::move(state)]() mutable {
            cb_function(future<future_t> {state});
        });
    });
}
*/

}

ext: h

max_stars_repo_name: liamom/future

max_issues_repo_name: liamom/future

max_forks_repo_name: liamom/future


text:
#pragma once

#include <gsl/gsl>
#include <vector>
#include <sstream>

namespace multiformats {

    //
    // Raw memory byte
    //
    using byte_t = uint8_t;

    //
    // Raw memory buffer
    //   Represents a contiguous run of byte_t elements, or a const view of it.
    //
    using buffer_t = std::vector<byte_t>;
    using bufferview_t = gsl::span<const byte_t>;

    //    Appends a buffer or a byte to a buffer
    inline buffer_t& operator += (buffer_t& _Left, bufferview_t _Right) {
        _Left.insert(_Left.end(), _Right.begin(), _Right.end());
        return _Left;
    }
    inline buffer_t& operator += (buffer_t& _Left, byte_t _Value) {
        _Left.push_back(_Value);
        return _Left;
    }

    //    Concatenates two buffers
    inline buffer_t operator + (bufferview_t _Left, bufferview_t _Right) {
        auto result = buffer_t(_Left.size() + _Right.size());
        std::copy(_Right.begin(), _Right.end(), std::copy(_Left.begin(), _Left.end(), result.begin()));
        return result;
    }

    //    Buffer comparison
    inline bool operator == (bufferview_t _Left, bufferview_t _Right) {
        return std::equal(_Left.begin(), _Left.end(), _Right.begin(), _Right.end());
    }

    inline bool operator != (bufferview_t _Left, bufferview_t _Right) {
        return !(_Left == _Right);
    }

    //
    // String and constant view of strings
    //
    using string_t = std::string;
    using stringview_t = gsl::cstring_span<>;

    inline std::vector<stringview_t> split(stringview_t s, char delim) {
        auto elems = std::vector<stringview_t>{};

        auto begin = s.begin();
        auto end = s.end();
        auto first = begin;
        while (first != end) {
            auto last = std::find(first, end, delim);
            elems.push_back(s.subspan(first - begin, last - first));
            if (last == end) break;
            first = last + 1;
        }

        return elems;
    }

    inline string_t operator+(stringview_t _Left, const string_t& _Right) { return (to_string(_Left) + _Right); }
    inline string_t operator+(const string_t& _Left, stringview_t _Right) { return (_Left + to_string(_Right)); }

    //
    // Conversions
    //   to_* methods create a new container and copy the content
    //   as_* methods return a const view of the container without copy
    //
    inline buffer_t to_buffer(stringview_t s) { return { s.begin(), s.end() }; }
    inline buffer_t to_buffer(const char* s)  { return to_buffer(gsl::ensure_z(s)); }

    inline bufferview_t as_buffer(stringview_t s) { return { reinterpret_cast<const byte_t*>(s.data()), gsl::narrow<ptrdiff_t>(s.size()) }; }
    inline stringview_t as_string(bufferview_t b) { return { reinterpret_cast<const char*>(b.data()), b.size() }; }

    inline string_t to_string(bufferview_t b) { return to_string(as_string(b)); }
}

max_stars_repo_name: cedrou/cpp-multiformats

max_issues_repo_name: cedrou/multiformats

max_forks_repo_name: cedrou/multiformats


text:
/* Copyright (c) 2011-2012, Jérémy Fix. All rights reserved. */

/* Redistribution and use in source and binary forms, with or without */
/* modification, are permitted provided that the following conditions are met: */

/* * Redistributions of source code must retain the above copyright notice, */
/* this list of conditions and the following disclaimer. */
/* * 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. */
/* * None of the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. */

/* THIS SOFTWARE IS PROVIDED BY THE AUTHOR 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. */

#ifndef UKF_NDIM_STATE_H
#define UKF_NDIM_STATE_H

#include <gsl/gsl_linalg.h> // For the Cholesky decomposition
#include <gsl/gsl_math.h>
#include <gsl/gsl_blas.h>

#include "ukf_types.h"

namespace ukf
{

  /**
   * @short UKF for state estimation, additive noise case
   * The notations follow "Sigma-Point Kalman Filters for Probabilistic Inference in Dynamic State-Space Models",p108, PhD, van Der Merwe
   */
  namespace state
  {
    /**
     * @short Allocation of the vectors/matrices and initialization
     *
     */
    void ukf_init(ukf_param &p, ukf_state &s)
    {
      // Parameters for the sigma points of the process equation
      p.nbSamples = 2 * p.n + 1;
      p.lambda = p.alpha * p.alpha * (p.n + p.kpa) - p.n;
      p.gamma = sqrt(p.n + p.lambda);

      // Parameters for the sigma points of the observation equation
      p.nbSamplesMeasure = 4 * p.n + 1;
      p.lambda_aug = p.alpha * p.alpha * (2*p.n + p.kpa) - 2*p.n;
      p.gamma_aug = sqrt(2*p.n + p.lambda_aug);

      // Init the matrices used to iterate
      s.xi = gsl_vector_alloc(p.n);
      gsl_vector_set_zero(s.xi);

      s.xi_prediction = gsl_matrix_alloc(p.n, p.nbSamples);
      gsl_matrix_set_zero(s.xi_prediction);

      s.xi_mean = gsl_vector_alloc(p.n);
      gsl_vector_set_zero(s.xi_mean);

      s.Pxxi = gsl_matrix_alloc(p.n,p.n);
      gsl_matrix_set_identity(s.Pxxi);
      gsl_matrix_scale(s.Pxxi, p.prior_x);

      s.cholPxxi = gsl_matrix_alloc(p.n,p.n);
      gsl_matrix_set_zero(s.cholPxxi);

      s.Pvvi = gsl_matrix_alloc(p.n,p.n);
      s.cholPvvi = gsl_matrix_alloc(p.n,p.n);
      p.evolution_noise->init(p,s);

      s.yi_prediction = gsl_matrix_alloc(p.no, p.nbSamplesMeasure);
      gsl_matrix_set_zero(s.yi_prediction);

      s.yi_mean = gsl_vector_alloc(p.no);
      gsl_vector_set_zero(s.yi_mean);

      s.ino_yi = gsl_vector_alloc(p.no);
      gsl_vector_set_zero(s.ino_yi);

      s.Pyyi = gsl_matrix_alloc(p.no, p.no);
      gsl_matrix_set_zero(s.Pyyi);

      s.Pnni = gsl_matrix_alloc(p.no,p.no);
      gsl_matrix_set_identity(s.Pnni);
      gsl_matrix_scale(s.Pnni, p.measurement_noise);

      s.Pxyi = gsl_matrix_alloc(p.n, p.no);
      gsl_matrix_set_zero(s.Pxyi);

      s.sigmaPoint = gsl_vector_alloc(p.n);
      gsl_vector_set_zero(s.sigmaPoint);

      s.sigmaPoints = gsl_matrix_alloc(p.n, p.nbSamples);
      gsl_matrix_set_zero(s.sigmaPoints);

      s.sigmaPointMeasure = gsl_vector_alloc(p.n);
      gsl_vector_set_zero(s.sigmaPoint);

      s.sigmaPointsMeasure = gsl_matrix_alloc(p.n, p.nbSamplesMeasure);
      gsl_matrix_set_zero(s.sigmaPointsMeasure);

      // Weights used to update the statistics
      s.wm_j  = gsl_vector_alloc(p.nbSamples); // Weights used to compute the mean of the sigma points images
      s.wc_j = gsl_vector_alloc(p.nbSamples); // Weights used to update the covariance matrices

      // Set the weights
      gsl_vector_set(s.wm_j, 0, p.lambda / (p.n + p.lambda));
      gsl_vector_set(s.wc_j, 0, p.lambda / (p.n + p.lambda) + (1.0 - p.alpha*p.alpha + p.beta));
      for(int j = 1 ; j < p.nbSamples; j ++)
	{
	  gsl_vector_set(s.wm_j, j, 1.0 / (2.0 * (p.n + p.lambda)));
	  gsl_vector_set(s.wc_j, j, 1.0 / (2.0 * (p.n + p.lambda)));
	}

      // Set the weights
      s.wm_aug_j  = gsl_vector_alloc(p.nbSamplesMeasure); // Weights used to compute the mean of the sigma points images
      s.wc_aug_j = gsl_vector_alloc(p.nbSamplesMeasure); // Weights used to update the covariance matrices
      gsl_vector_set(s.wm_aug_j, 0, p.lambda_aug / (2*p.n + p.lambda_aug));
      gsl_vector_set(s.wc_aug_j, 0, p.lambda_aug / (2*p.n + p.lambda_aug) + (1.0 - p.alpha*p.alpha + p.beta));
      for(int j = 1 ; j < p.nbSamplesMeasure; j ++)
	{
	  gsl_vector_set(s.wm_aug_j, j, 1.0 / (2.0 * (2*p.n + p.lambda_aug)));
	  gsl_vector_set(s.wc_aug_j, j, 1.0 / (2.0 * (2*p.n + p.lambda_aug)));
	}

      s.Ki = gsl_matrix_alloc(p.n, p.no);
      s.Ki_T = gsl_matrix_alloc(p.no, p.n);

      // Allocate temporary matrices
      s.temp_n = gsl_vector_alloc(p.n);

      s.temp_n_1 = gsl_matrix_alloc(p.n,1);
      s.temp_1_n = gsl_matrix_alloc(1,p.n);
      s.temp_n_n = gsl_matrix_alloc(p.n, p.n);
      s.temp_n_no = gsl_matrix_alloc(p.n, p.no);
      s.temp_no_1 = gsl_matrix_alloc(p.no,1);
      s.temp_1_no = gsl_matrix_alloc(1,p.no);
      s.temp_no_no = gsl_matrix_alloc(p.no, p.no);
    }

    /**
     * @short Free of memory allocation
     *
     */
    void ukf_free(ukf_param &p, ukf_state &s)
    {
      gsl_vector_free(s.xi);
      gsl_matrix_free(s.xi_prediction);
      gsl_vector_free(s.xi_mean);
      gsl_matrix_free(s.Pxxi);
      gsl_matrix_free(s.cholPxxi);
      gsl_matrix_free(s.Pvvi);
      gsl_matrix_free(s.cholPvvi);

      gsl_matrix_free(s.yi_prediction);
      gsl_vector_free(s.yi_mean);
      gsl_vector_free(s.ino_yi);
      gsl_matrix_free(s.Pyyi);
      gsl_matrix_free(s.Pnni);

      gsl_matrix_free(s.Pxyi);

      gsl_vector_free(s.sigmaPoint);
      gsl_matrix_free(s.sigmaPoints);

      gsl_vector_free(s.sigmaPointMeasure);
      gsl_matrix_free(s.sigmaPointsMeasure);

      gsl_vector_free(s.wm_j);
      gsl_vector_free(s.wc_j);

      gsl_vector_free(s.wm_aug_j);
      gsl_vector_free(s.wc_aug_j);

      gsl_matrix_free(s.Ki);
      gsl_matrix_free(s.Ki_T);

      gsl_vector_free(s.temp_n);
      gsl_matrix_free(s.temp_n_1);
      gsl_matrix_free(s.temp_1_n);
      gsl_matrix_free(s.temp_n_n);

      gsl_matrix_free(s.temp_n_no);
      gsl_matrix_free(s.temp_no_1);
      gsl_matrix_free(s.temp_1_no);
      gsl_matrix_free(s.temp_no_no);
    }

    /**
     * @short UKF-additive (zero-mean) noise case, "Kalman Filtering and Neural Networks", p.233
     *
     */
    template<typename FFUNC,
      typename HFUNC>
    void ukf_iterate(ukf_param &p, ukf_state &s,
		     FFUNC f,
		     HFUNC h,
		     gsl_vector* yi)
    {
      int i,j,k;

      // ************************************************** //
      // ************ Compute the sigma points ************ //
      // ************************************************** //

      // 0 - Compute the Cholesky decomposition of s.Pxxi
      gsl_matrix_memcpy(s.cholPxxi, s.Pxxi);
      gsl_linalg_cholesky_decomp(s.cholPxxi);
      // Set all the elements of cholPvvi strictly above the diagonal to zero
      for(j = 0 ; j < p.n ; j++)
	for(k = j+1 ; k < p.n ; k++)
	  gsl_matrix_set(s.cholPxxi,j,k,0.0);

      // 1- Compute the sigma points,
      // Equation (3.170)
      // sigmapoint_j = x_(i-1)
      // sigmapoint_j = x_(i-1) + gamma * sqrt(P_i-1)_j for 1 <= j <= n
      // sigmapoint_j = x_(i-1) - gamma * sqrt(P_i-1)_(j-(n+1)) for n+1 <= j <= 2n
      gsl_matrix_set_col(s.sigmaPoints, 0, s.xi);
      for(j = 1 ; j < p.n+1 ; ++j)
	for(i = 0 ; i < p.n ; ++i)
	  {
	    gsl_matrix_set(s.sigmaPoints,i,j, s.xi->data[i] + p.gamma * gsl_matrix_get(s.cholPxxi, i, j-1));
	    gsl_matrix_set(s.sigmaPoints,i,j+p.n, s.xi->data[i] - p.gamma * gsl_matrix_get(s.cholPxxi, i, j-1));
	  }

      /**********************************/
      /***** Time update equations  *****/
      /**********************************/

      // Time update equations
      // 0 - Compute the image of the sigma points and the mean of these images
      gsl_vector_set_zero(s.xi_mean);
      gsl_vector_view vec_view;
      for(j = 0 ; j < p.nbSamples ; ++j)
	{
	  gsl_matrix_get_col(s.sigmaPoint, s.sigmaPoints, j);
	  vec_view = gsl_matrix_column(s.xi_prediction,j);
	  f(s.params, s.sigmaPoint, &vec_view.vector);

	  // Update the mean, Eq (3.172)
	  for(i = 0 ; i < p.n ; ++i)
	    s.xi_mean->data[i] += s.wm_j->data[j] * gsl_matrix_get(s.xi_prediction,i,j);
	}

      // 1 - Compute the covariance of the images and add the process noise,
      // Equation (3.173)
      // Warning, s.Pxxi will now hold P_xk^-
      gsl_matrix_set_zero(s.Pxxi);
      for(j = 0 ; j < p.nbSamples ; ++j)
	{
	  for(i = 0 ; i < p.n ; ++i)
	    s.temp_n_1->data[i] = gsl_matrix_get(s.xi_prediction,i,j) - s.xi_mean->data[i];

	  gsl_blas_dgemm(CblasNoTrans, CblasTrans, s.wc_j->data[j] , s.temp_n_1, s.temp_n_1, 0, s.temp_n_n);
	  gsl_matrix_add(s.Pxxi, s.temp_n_n);
	}
      // Add the covariance of the evolution noise
      gsl_matrix_add(s.Pxxi, s.Pvvi);

      // Augment sigma points
      // Equation 3.174
      // First put the images of the initial sigma points
      gsl_matrix_view mat_view;
      mat_view = gsl_matrix_submatrix(s.sigmaPointsMeasure, 0, 0, p.n, p.nbSamples);
      gsl_matrix_memcpy(&mat_view.matrix, s.xi_prediction);
      // And add the additional sigma points eq. (7.56)
      for(j = 0 ; j < p.n ; ++j)
	{
	  for(i = 0 ; i < p.n ; ++i)
	    {
	      gsl_matrix_set(s.sigmaPointsMeasure, i, j+p.nbSamples, gsl_matrix_get(s.xi_prediction,i,0)+p.gamma_aug*gsl_matrix_get(s.cholPvvi,i,j));
	      gsl_matrix_set(s.sigmaPointsMeasure, i, j+p.nbSamples+p.n, gsl_matrix_get(s.xi_prediction,i,0)-p.gamma_aug*gsl_matrix_get(s.cholPvvi,i,j));
	    }
	}

      // Compute the image of the sigma points through the observation equation
      // eq (3.175)
      gsl_vector_set_zero(s.yi_mean);
      for(j = 0 ; j < p.nbSamplesMeasure ; ++j)
	{
	  gsl_matrix_get_col(s.sigmaPointMeasure, s.sigmaPointsMeasure, j);
	  vec_view = gsl_matrix_column(s.yi_prediction,j);
	  h(s.sigmaPointMeasure, &vec_view.vector);

	  // Update the mean , eq (3.176)
	  for(i = 0 ; i < p.no ; ++i)
	    s.yi_mean->data[i] += s.wm_aug_j->data[j] * gsl_matrix_get(s.yi_prediction,i,j);
	}

      /*****************************************/
      /***** Measurement update equations  *****/
      /*****************************************/

      // Compute the covariance of the observations
      // Eq. (3.177)
      // Initialize with the observation noise covariance
      gsl_matrix_memcpy(s.Pyyi, s.Pnni);
      for(j = 0 ; j < p.nbSamplesMeasure ; ++j)
	{
	  for(i = 0 ; i < p.no ; ++i)
	    s.temp_no_1->data[i] = gsl_matrix_get(s.yi_prediction,i,j) - s.yi_mean->data[i];

	  gsl_blas_dgemm(CblasNoTrans, CblasTrans, s.wc_aug_j->data[j] , s.temp_no_1, s.temp_no_1, 0, s.temp_no_no);
	  gsl_matrix_add(s.Pyyi, s.temp_no_no);
	}

      // Compute the state/observation covariance
      // Eq (3.178)
      gsl_matrix_set_zero(s.Pxyi);
      for(j = 0 ; j < p.nbSamplesMeasure ; ++j)
	{
	  for(i = 0 ; i < p.n ; ++i)
	    s.temp_n_1->data[i] = gsl_matrix_get(s.sigmaPointsMeasure,i,j) - s.xi_mean->data[i];

	  for(i = 0 ; i < p.no ; ++i)
	    s.temp_1_no->data[i] = gsl_matrix_get(s.yi_prediction,i,j) - s.yi_mean->data[i];

	  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, s.wc_aug_j->data[j] , s.temp_n_1, s.temp_1_no, 0, s.temp_n_no);
	  gsl_matrix_add(s.Pxyi, s.temp_n_no);
	}

      // Compute the Kalman gain, eq (3.179)
      // 0- Compute the inverse of Pyyi
      gsl_matrix_memcpy(s.temp_no_no, s.Pyyi);
      gsl_linalg_cholesky_decomp(s.temp_no_no);
      gsl_linalg_cholesky_invert(s.temp_no_no);

      // 1- Compute the Kalman gain
      gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0 , s.Pxyi, s.temp_no_no, 0, s.Ki);

      // Correction : correct the estimation of the state
      // Eq. 3.180
      // Compute the innovations
      for(i = 0 ; i < p.no ; ++i)
	s.ino_yi->data[i] = gsl_vector_get(yi, i) - gsl_vector_get(s.yi_mean, i);
      gsl_vector_memcpy(s.xi, s.xi_mean);
      gsl_blas_dgemv(CblasNoTrans, 1.0 , s.Ki, s.ino_yi, 1.0, s.xi);

      // Correction : Update the covariance matrix Pk
      // Eq. 3.181
      gsl_matrix_transpose_memcpy(s.Ki_T, s.Ki);
      gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0 , s.Ki, s.Pyyi, 0, s.temp_n_no);
      gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, -1.0 , s.temp_n_no, s.Ki_T, 1.0, s.Pxxi);

      // Update of the process noise
      p.evolution_noise->updateEvolutionNoise(p,s);
      //            switch(p.process_noise_type)
      //            {
      //            case ukf::UKF_PROCESS_FIXED:
      //                //nothing to do
      //                break;
      //            case ukf::UKF_PROCESS_RLS:
      //                gsl_matrix_memcpy(s.Pvvi, s.Pxxi);
      //                gsl_matrix_scale(s.Pvvi, 1.0/p.process_noise-1.0);
      //                gsl_matrix_memcpy(s.cholPvvi, s.Pvvi);
      //                gsl_linalg_cholesky_decomp(s.cholPvvi);
      //                for(j = 0 ; j < p.n ; j++)
      //                    for(k = j+1 ; k < p.n ; k++)
      //                        gsl_matrix_set(s.cholPvvi,j,k,0.0);
      //                break;
      //            default:
      //                printf("Warning : Unrecognized process noise type\n");
      //            }

    }

    /**
     * @short Evaluation of the output from the sigma points
     *
     */
    void ukf_evaluate(ukf_param &p, ukf_state &s,
		      void (*f)(gsl_vector*, gsl_vector *, gsl_vector *),
		      void (*h)(gsl_vector*, gsl_vector *), gsl_vector* yi)
    {

      int i,j,k;
      gsl_matrix_view mat_view;
      gsl_vector_view vec_view;
      // ************************************************** //
      // ************ Compute the sigma points ************ //
      // ************************************************** //

      // 0 - Compute the Cholesky decomposition of s.Pxxi
      gsl_matrix_memcpy(s.cholPxxi, s.Pxxi);
      gsl_linalg_cholesky_decomp(s.cholPxxi);
      // Set all the elements of cholPvvi strictly above the diagonal to zero
      for(j = 0 ; j < p.n ; j++)
	for(k = j+1 ; k < p.n ; k++)
	  gsl_matrix_set(s.cholPxxi,j,k,0.0);

      // 1- Compute the sigma points,
      // Equation (3.170)
      // sigmapoint_j = x_(i-1)
      // sigmapoint_j = x_(i-1) + gamma * sqrt(P_i-1)_j for 1 <= j <= n
      // sigmapoint_j = x_(i-1) - gamma * sqrt(P_i-1)_(j-(n+1)) for n+1 <= j <= 2n
      gsl_matrix_set_col(s.sigmaPoints, 0, s.xi);
      for(j = 1 ; j < p.n+1 ; ++j)
	for(i = 0 ; i < p.n ; ++i)
	  {
	    gsl_matrix_set(s.sigmaPoints,i,j, s.xi->data[i] + p.gamma * gsl_matrix_get(s.cholPxxi, i, j-1));
	    gsl_matrix_set(s.sigmaPoints,i,j+p.n, s.xi->data[i] - p.gamma * gsl_matrix_get(s.cholPxxi, i, j-1));
	  }

      /**********************************/
      /***** Time update equations  *****/
      /**********************************/

      // Time update equations
      // 0 - Compute the image of the sigma points and the mean of these images
      gsl_vector_set_zero(s.xi_mean);
      for(j = 0 ; j < p.nbSamples ; ++j)
	{
	  gsl_matrix_get_col(s.sigmaPoint, s.sigmaPoints, j);
	  vec_view = gsl_matrix_column(s.xi_prediction,j);
	  f(s.params, s.sigmaPoint, &vec_view.vector);

	  // Update the mean, Eq (3.172)
	  for(i = 0 ; i < p.n ; ++i)
	    s.xi_mean->data[i] += s.wm_j->data[j] * gsl_matrix_get(s.xi_prediction,i,j);
	}

      // 1 - Compute the covariance of the images and add the process noise,
      // Equation (3.173)
      // Warning, s.Pxxi will now hold P_xk^-
      gsl_matrix_set_zero(s.Pxxi);
      for(j = 0 ; j < p.nbSamples ; ++j)
	{
	  for(i = 0 ; i < p.n ; ++i)
	    s.temp_n_1->data[i] = gsl_matrix_get(s.xi_prediction,i,j) - s.xi_mean->data[i];

	  gsl_blas_dgemm(CblasNoTrans, CblasTrans, s.wc_j->data[j] , s.temp_n_1, s.temp_n_1, 0, s.temp_n_n);
	  gsl_matrix_add(s.Pxxi, s.temp_n_n);
	}
      // Add the covariance of the evolution noise
      gsl_matrix_add(s.Pxxi, s.Pvvi);

      // Augment sigma points
      // Equation 3.174
      // First put the images of the initial sigma points
      mat_view = gsl_matrix_submatrix(s.sigmaPointsMeasure, 0, 0, p.n, p.nbSamples);
      gsl_matrix_memcpy(&mat_view.matrix, s.xi_prediction);
      // And add the additional sigma points eq. (7.56)
      for(j = 0 ; j < p.n ; ++j)
	{
	  for(i = 0 ; i < p.n ; ++i)
	    {
	      gsl_matrix_set(s.sigmaPointsMeasure, i, j+p.nbSamples, gsl_matrix_get(s.xi_prediction,i,0)+p.gamma_aug*gsl_matrix_get(s.cholPvvi,i,j));
	      gsl_matrix_set(s.sigmaPointsMeasure, i, j+p.nbSamples+p.n, gsl_matrix_get(s.xi_prediction,i,0)-p.gamma_aug*gsl_matrix_get(s.cholPvvi,i,j));
	    }
	}

      // Compute the image of the sigma points through the observation equation
      // eq (3.175)
      gsl_vector_set_zero(yi);
      for(j = 0 ; j < p.nbSamplesMeasure ; ++j)
	{
	  gsl_matrix_get_col(s.sigmaPointMeasure, s.sigmaPointsMeasure, j);
	  vec_view = gsl_matrix_column(s.yi_prediction,j);
	  h(s.sigmaPointMeasure, &vec_view.vector);

	  // Update the mean , eq (3.176)
	  for(i = 0 ; i < p.no ; ++i)
	    yi->data[i] += s.wm_aug_j->data[j] * gsl_matrix_get(s.yi_prediction,i,j);
	}
    }
  } // state
} // ukf

#endif // UKF_NDIM_STATE_H


text:
/* histogram/test.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 <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_histogram.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_ieee_utils.h>

#define N 397
#define NR 10

void
test1d (void)
{
  double xr[NR + 1] =
  {0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0};

  gsl_histogram *h, *h1, *hr, *g;
  size_t i, j;

  gsl_ieee_env_setup ();

  h = gsl_histogram_calloc (N);
  h1 = gsl_histogram_calloc (N);
  g = gsl_histogram_calloc (N);

  gsl_test (h->range == 0, "gsl_histogram_alloc returns valid range pointer");
  gsl_test (h->bin == 0, "gsl_histogram_alloc returns valid bin pointer");
  gsl_test (h->n != N, "gsl_histogram_alloc returns valid size");

  hr = gsl_histogram_calloc_range (NR, xr);

  gsl_test (hr->range == 0, "gsl_histogram_calloc_range returns valid range pointer");
  gsl_test (hr->bin == 0, "gsl_histogram_calloc_range returns valid bin pointer");
  gsl_test (hr->n != NR, "gsl_histogram_calloc_range returns valid size");

  {
    int status = 0;
    for (i = 0; i <= NR; i++)
      {
        if (hr->range[i] != xr[i])
          {
            status = 1;
          }
      };

    gsl_test (status, "gsl_histogram_calloc_range creates range");
  }

  for (i = 0; i <= NR; i++)
    {
      hr->range[i] = 0.0;
    }

  {
    int status = gsl_histogram_set_ranges (hr, xr, NR+1);

    for (i = 0; i <= NR; i++)
      {
        if (hr->range[i] != xr[i])
          {
            status = 1;
          }
      };

    gsl_test (status, "gsl_histogram_set_range sets range");
  }

  for (i = 0; i < N; i++)
    {
      gsl_histogram_accumulate (h, (double) i, (double) i);
    };

  {
    int status = 0;

    for (i = 0; i < N; i++)
      {
        if (h->bin[i] != (double) i)
          {
            status = 1;
          }
      };

    gsl_test (status, "gsl_histogram_accumulate writes into array");
  }

  {
    int status = 0;

    for (i = 0; i < N; i++)
      {
        if (gsl_histogram_get (h, i) != i)
          status = 1;
      };
    gsl_test (status, "gsl_histogram_get reads from array");
  }

  for (i = 0; i <= N; i++)
    {
      h1->range[i] = 100.0 + i;
    }

  gsl_histogram_memcpy (h1, h);

  {
    int status = 0;
    for (i = 0; i <= N; i++)
      {
        if (h1->range[i] != h->range[i])
          status = 1;
      };
    gsl_test (status, "gsl_histogram_memcpy copies bin ranges");
  }

  {
    int status = 0;
    for (i = 0; i < N; i++)
      {
        if (gsl_histogram_get (h1, i) != gsl_histogram_get (h, i))
          status = 1;
      };
    gsl_test (status, "gsl_histogram_memcpy copies bin values");
  }

  gsl_histogram_free (h1);

  h1 = gsl_histogram_clone (h);

  {
    int status = 0;
    for (i = 0; i <= N; i++)
      {
        if (h1->range[i] != h->range[i])
          status = 1;
      };
    gsl_test (status, "gsl_histogram_clone copies bin ranges");
  }

  {
    int status = 0;
    for (i = 0; i < N; i++)
      {
        if (gsl_histogram_get (h1, i) != gsl_histogram_get (h, i))
          status = 1;
      };
    gsl_test (status, "gsl_histogram_clone copies bin values");
  }

  gsl_histogram_reset (h);

  {
    int status = 0;

    for (i = 0; i < N; i++)
      {
        if (h->bin[i] != 0)
          status = 1;
      }
    gsl_test (status, "gsl_histogram_reset zeros array");
  }

  {
    int status = 0;

    for (i = 0; i < N; i++)
      {
        gsl_histogram_increment (h, (double) i);

        for (j = 0; j <= i; j++)
          {
            if (h->bin[j] != 1)
              {
                status = 1;
              }
          }

        for (j = i + 1; j < N; j++)
          {
            if (h->bin[j] != 0)
              {
                status = 1;
              }
          }
      }

    gsl_test (status, "gsl_histogram_increment increases bin value");
  }

  {
    int status = 0;
    for (i = 0; i < N; i++)
      {
        double x0 = 0, x1 = 0;

        gsl_histogram_get_range (h, i, &x0, &x1);

        if (x0 != i || x1 != i + 1)
          {
            status = 1;
          }
      }
    gsl_test (status, "gsl_histogram_getbinrange returns bin range");
  }

  {
    int status = 0;
    if (gsl_histogram_max (h) != N)
      status = 1;
    gsl_test (status, "gsl_histogram_max returns maximum");
  }

  {
    int status = 0;
    if (gsl_histogram_min (h) != 0)
      status = 1;
    gsl_test (status, "gsl_histogram_min returns minimum");
  }

  {
    int status = 0;
    if (gsl_histogram_bins (h) != N)
      status = 1;
    gsl_test (status, "gsl_histogram_bins returns number of bins");
  }

  h->bin[2] = 123456.0;
  h->bin[4] = -654321;

  {
    double max = gsl_histogram_max_val (h);
    gsl_test (max != 123456.0, "gsl_histogram_max_val finds maximum value");
  }

  {
    double min = gsl_histogram_min_val (h);
    gsl_test (min != -654321.0, "gsl_histogram_min_val finds minimum value");
  }

  {
    size_t imax = gsl_histogram_max_bin (h);
    gsl_test (imax != 2, "gsl_histogram_max_bin finds maximum value bin");
  }

  {
    size_t imin = gsl_histogram_min_bin (h);
    gsl_test (imin != 4, "gsl_histogram_min_bin find minimum value bin");
  }

  for (i = 0; i < N; i++)
    {
      h->bin[i] = i + 27;
      g->bin[i] = (i + 27) * (i + 1);
    }

  {
    double sum=gsl_histogram_sum (h);
    gsl_test(sum != N*27+((N-1)*N)/2, "gsl_histogram_sum sums all bin values");
  }

  gsl_histogram_memcpy (h1, g);
  gsl_histogram_add (h1, h);

  {
    int status = 0;
    for (i = 0; i < N; i++)
      {
        if (h1->bin[i] != g->bin[i] + h->bin[i])
          status = 1;
      }
    gsl_test (status, "gsl_histogram_add histogram addition");
  }

  gsl_histogram_memcpy (h1, g);
  gsl_histogram_sub (h1, h);

  {
    int status = 0;
    for (i = 0; i < N; i++)
      {
        if (h1->bin[i] != g->bin[i] - h->bin[i])
          status = 1;
      }
    gsl_test (status, "gsl_histogram_sub histogram subtraction");
  }

  gsl_histogram_memcpy (h1, g);
  gsl_histogram_mul (h1, h);

  {
    int status = 0;
    for (i = 0; i < N; i++)
      {
        if (h1->bin[i] != g->bin[i] * h->bin[i])
          status = 1;
      }
    gsl_test (status, "gsl_histogram_mul histogram multiplication");
  }

  gsl_histogram_memcpy (h1, g);
  gsl_histogram_div (h1, h);

  {
    int status = 0;
    for (i = 0; i < N; i++)
      {
        if (h1->bin[i] != g->bin[i] / h->bin[i])
          status = 1;
      }
    gsl_test (status, "gsl_histogram_div histogram division");
  }

  gsl_histogram_memcpy (h1, g);
  gsl_histogram_scale (h1, 0.5);

  {
    int status = 0;
    for (i = 0; i < N; i++)
      {
        if (h1->bin[i] != 0.5 * g->bin[i])
          status = 1;
      }
    gsl_test (status, "gsl_histogram_scale histogram scaling");
  }

  gsl_histogram_memcpy (h1, g);
  gsl_histogram_shift (h1, 0.25);

  {
    int status = 0;
    for (i = 0; i < N; i++)
      {
        if (h1->bin[i] != 0.25 + g->bin[i])
          status = 1;
      }
    gsl_test (status, "gsl_histogram_shift histogram shift");
  }

  gsl_histogram_free (h);       /* free whatever is in h */

  h = gsl_histogram_calloc_uniform (N, 0.0, 1.0);

  gsl_test (h->range == 0,
            "gsl_histogram_calloc_uniform returns valid range pointer");
  gsl_test (h->bin == 0,
            "gsl_histogram_calloc_uniform returns valid bin pointer");
  gsl_test (h->n != N,
            "gsl_histogram_calloc_uniform returns valid size");

  gsl_histogram_accumulate (h, 0.0, 1.0);
  gsl_histogram_accumulate (h, 0.1, 2.0);
  gsl_histogram_accumulate (h, 0.2, 3.0);
  gsl_histogram_accumulate (h, 0.3, 4.0);

  {
    size_t i1, i2, i3, i4;
    double expected;
    int status = gsl_histogram_find (h, 0.0, &i1);
    status = gsl_histogram_find (h, 0.1, &i2);
    status = gsl_histogram_find (h, 0.2, &i3);
    status = gsl_histogram_find (h, 0.3, &i4);

    for (i = 0; i < N; i++)
      {
        if (i == i1)
          {
            expected = 1.0;
          }
        else if (i == i2)
          {
            expected = 2.0;
          }
        else if (i == i3)
          {
            expected = 3.0;
          }
        else if (i == i4)
          {
            expected = 4.0;
          }
        else
          {
            expected = 0.0;
          }

        if (h->bin[i] != expected)
          {
            status = 1;
          }

      }
    gsl_test (status, "gsl_histogram_find returns index");
  }

  {
    FILE *f = fopen ("test.txt", "w");
    gsl_histogram_fprintf (f, h, "%.19e", "%.19e");
    fclose (f);
  }

  {
    FILE *f = fopen ("test.txt", "r");
    gsl_histogram *hh = gsl_histogram_calloc (N);
    int status = 0;

    gsl_histogram_fscanf (f, hh);

    for (i = 0; i < N; i++)
      {
        if (h->range[i] != hh->range[i])
          status = 1;
        if (h->bin[i] != hh->bin[i])
          status = 1;
      }
    if (h->range[N] != hh->range[N])
      status = 1;

    gsl_test (status, "gsl_histogram_fprintf and fscanf");

    gsl_histogram_free (hh);
    fclose (f);
  }

  {
    FILE *f = fopen ("test.dat", "wb");
    gsl_histogram_fwrite (f, h);
    fclose (f);
  }

  {
    FILE *f = fopen ("test.dat", "rb");
    gsl_histogram *hh = gsl_histogram_calloc (N);
    int status = 0;

    gsl_histogram_fread (f, hh);

    for (i = 0; i < N; i++)
      {
        if (h->range[i] != hh->range[i])
          status = 1;
        if (h->bin[i] != hh->bin[i])
          status = 1;
      }
    if (h->range[N] != hh->range[N])
      status = 1;

    gsl_test (status, "gsl_histogram_fwrite and fread");

    gsl_histogram_free (hh);
    fclose (f);
  }

  gsl_histogram_free (h);
  gsl_histogram_free (g);
  gsl_histogram_free (h1);
  gsl_histogram_free (hr);
}

ext: c

max_stars_repo_name: utdsimmons/ohpc

max_issues_repo_name: utdsimmons/ohpc

max_forks_repo_name: utdsimmons/ohpc


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

#pragma once

#include "arcana/errors.h"
#include "arcana/expected.h"
#include "arcana/functional/inplace_function.h"
#include "arcana/iterators.h"
#include "arcana/type_traits.h"

#include "cancellation.h"

#include <gsl/gsl>
#include <memory>
#include <stdexcept>

#include <atomic>

namespace mira
{
    template<typename ResultT>
    class task_completion_source;

    //
    // A scheduler that will invoke the continuation inline
    // right after the previous task.
    //
    namespace
    {
        struct inline_scheduler_t
        {
            template<typename CallableT>
            constexpr void queue(CallableT&& callable) const
            {
                callable();
            }
        };

        constexpr inline_scheduler_t inline_scheduler{};
    }
}

#include "internal/internal_task.h"

namespace mira
{
    //
    //  Generic task system to run work with continuations on a generic scheduler.
    //
    //  The scheduler on which tasks are queued must satisfy this contract:
    //
    //      struct scheduler
    //      {
    //          template<CallableT>
    //          void queue(CallableT&& callable)
    //          {
    //              callable must be usable like this: callable();
    //          }
    //      };
    //
    template<typename ResultT>
    class task
    {
        using payload_t = internal::base_task_payload_with_return<ResultT>;
        using payload_ptr = std::shared_ptr<payload_t>;

        static_assert(std::is_same<typename as_expected<ResultT>::result, ResultT>::value,
            "task can't be of expected<T>");

    public:
        using result_t = ResultT;

        task() = default;

        task(const task& other) = default;
        task(task&& other) = default;

        task(const task_completion_source<ResultT>& source)
            : m_payload{ source.m_payload }
        {}

        task(task_completion_source<ResultT>&& source)
            : m_payload{ std::move(source.m_payload) }
        {}

        task& operator=(task&& other) = default;
        task& operator=(const task& other) = default;

        bool operator==(const task& other)
        {
            return m_payload == other.m_payload;
        }

        //
        // Executes a callable on this scheduler once this task is finished and
        // returns a task that represents the callable.
        //
        // Calling .then() on the returned task will queue a task to run after the
        // callable is run.
        //
        template<typename SchedulerT, typename CallableT>
        auto then(SchedulerT& scheduler, cancellation& token, CallableT&& callable)
        {
            return internal::continuation_factory<ResultT>::create_continuation_task(
                internal::input_output_wrapper<ResultT>::wrap_callable(std::forward<CallableT>(callable), token),
                scheduler,
                m_payload);
        }

    private:
        explicit task(payload_ptr payload)
            : m_payload{ std::move(payload) }
        {}

        template<typename CallableT, typename InputT>
        explicit task(CallableT&& callable, type_of<InputT> input)
        {
            auto payload = std::make_shared<internal::task_payload<ResultT, sizeof(CallableT)>>(
                std::forward<CallableT>(callable), input);

            m_payload = std::move(payload);
        }

        template<typename OtherResultT, size_t WorkSize>
        friend struct internal::task_payload;

        template<typename OtherResultT>
        friend class task;

        friend class task_completion_source<ResultT>;

        template<typename OtherResultT, typename InputT>
        friend struct internal::task_factory;

        template<typename OtherResultT>
        friend struct internal::continuation_factory;

        template<typename SchedulerT, typename CallableT>
        friend auto make_task(SchedulerT& scheduler, cancellation& token, CallableT&& callable) -> typename internal::task_factory<typename as_expected<decltype(callable())>::result, void>::task_t;

        template<typename OtherResultT>
        friend task<typename as_expected<OtherResultT>::result> task_from_result(OtherResultT&& value);

        friend task<void> task_from_result();

        template<typename OtherResultT, typename ErrorT>
        friend task<OtherResultT> task_from_error(const ErrorT& error);

        payload_ptr m_payload;
    };
}

namespace mira
{
    template<typename ResultT>
    class task_completion_source
    {
        using payload_t = internal::noop_task_payload<ResultT>;
        using payload_ptr = std::shared_ptr<payload_t>;
        using uninitialized = std::integral_constant<int, 0>;

    public:
        using result_t = ResultT;

        task_completion_source()
            : m_payload{ std::make_shared<payload_t>() }
        {}

        static task_completion_source<ResultT> make_uninitialized()
        {
            return task_completion_source<ResultT>{ uninitialized{} };
        }

        //
        // Completes the task this source represents.
        //
        void complete()
        {
            static_assert(std::is_same<ResultT, void>::value,
                "complete with no arguments can only be used with a void completion source");
            m_payload->complete(expected<void>::make_valid());
        }

        //
        // Completes the task this source represents.
        //
        template<typename ValueT>
        void complete(ValueT&& value)
        {
            m_payload->complete(std::forward<ValueT>(value));
        }

        //
        // Returns whether or not the current source has already been completed.
        //
        bool completed() const
        {
            return m_payload->completed();
        }

        //
        // Converts this task_completion_source to a task object for consumers to use.
        //
        task<ResultT> as_task() const &
        {
            return task<ResultT>{ m_payload };
        }

        task<ResultT> as_task() &&
        {
            return task<ResultT>{ std::move(m_payload) };
        }

    private:
        explicit task_completion_source(uninitialized)
        {}

        explicit task_completion_source(std::shared_ptr<payload_t> payload)
            : m_payload{ std::move(payload) }
        {}

        friend class task<ResultT>;
        friend class abstract_task_completion_source;

        template<typename R, typename I>
        friend struct internal::task_factory;

        payload_ptr m_payload;
    };

    //
    // a type erased version of task_completion_source.
    //
    class abstract_task_completion_source
    {
        using payload_t = internal::base_task_payload;
        using payload_ptr = std::shared_ptr<payload_t>;

    public:
        abstract_task_completion_source()
            : m_payload{}
        {}

        template<typename T>
        explicit abstract_task_completion_source(const task_completion_source<T>& other)
            : m_payload{ other.m_payload }
        {}

        template<typename T>
        explicit abstract_task_completion_source(task_completion_source<T>&& other)
            : m_payload{ std::move(other.m_payload) }
        {}

        //
        // Returns whether or not the current source has already been completed.
        //
        bool completed() const
        {
            return m_payload->completed();
        }

        template<typename T>
        bool operator==(const task_completion_source<T>& other)
        {
            return m_payload == other.m_payload;
        }

        template<typename T>
        task_completion_source<T> unsafe_cast()
        {
            return task_completion_source<T>{ std::static_pointer_cast<internal::noop_task_payload<T>>(m_payload) };
        }

    private:
        payload_ptr m_payload;
    };

    //
    // creates a task and queues it to run on the given scheduler
    //
    template<typename SchedulerT, typename CallableT>
    inline auto make_task(SchedulerT& scheduler, cancellation& token, CallableT&& callable) -> typename internal::task_factory<typename as_expected<decltype(callable())>::result, void>::task_t
    {
        using callable_return_t = typename as_expected<decltype(callable())>::result;
        using wrapper = internal::input_output_wrapper<void>;

        internal::task_factory<callable_return_t, void> factory(
            wrapper::wrap_callable(std::forward<CallableT>(callable), token));

        scheduler.queue([to_run = std::move(factory.to_run)]{ to_run.m_payload->run(nullptr); });

        return factory.to_return;
    }

    //
    // creates a completed task from the given result
    //
    template<typename ResultT>
    inline task<typename as_expected<ResultT>::result> task_from_result(ResultT&& value)
    {
        task_completion_source<typename as_expected<ResultT>::result> result;
        result.complete(std::forward<ResultT>(value));
        return std::move(result);
    }

    inline task<void> task_from_result()
    {
        task_completion_source<void> result;
        result.complete();
        return std::move(result);
    }

    template<typename ResultT, typename ErrorT>
    inline task<ResultT> task_from_error(const ErrorT& error)
    {
        task_completion_source<ResultT> result;
        result.complete(std::make_error_code(error));
        return std::move(result);
    }

    template<typename ResultT>
    inline task<ResultT> task_from_error(const std::error_code& error)
    {
        task_completion_source<ResultT> result;
        result.complete(error);
        return std::move(result);
    }

    inline task<void> when_all(gsl::span<task<void>> tasks)
    {
        if (tasks.empty())
        {
            return task_from_result();
        }

        struct when_all_data
        {
            std::mutex mutex;
            size_t pendingCount;
            std::error_code error;
        };

        task_completion_source<void> result;
        auto data = std::make_shared<when_all_data>();
        data->pendingCount = tasks.size();

        for (task<void>& task : tasks)
        {
            task.then(inline_scheduler, cancellation::none(), [data, result](const expected<void>& exp) mutable {
                bool last = false;
                {
                    std::lock_guard<std::mutex> guard{ data->mutex };

                    data->pendingCount -= 1;
                    last = data->pendingCount == 0;

                    if (exp.has_error() && !data->error)
                    {
                        // set the first error, as it might have cascaded
                        data->error = exp.error();
                    }
                }

                if (last) // we were the last task to complete
                {
                    if (data->error)
                    {
                        result.complete(data->error);
                    }
                    else
                    {
                        result.complete();
                    }
                }
            });
        }

        return std::move(result);
    }

    template<typename T>
    inline task<std::vector<T>> when_all(gsl::span<task<T>> tasks)
    {
        if (tasks.empty())
        {
            return task_from_result<std::vector<T>>(std::vector<T>());
        }

        struct when_all_data
        {
            std::mutex mutex;
            size_t pendingCount;
            std::error_code error;
            std::vector<T> results;
        };

        task_completion_source<std::vector<T>> result;
        auto data = std::make_shared<when_all_data>();
        data->pendingCount = tasks.size();
        data->results.resize(tasks.size());

        //using forloop with index to be able to keep proper order of results
        for (auto idx = 0U; idx < data->results.size(); idx++)
        {
            tasks[idx].then(mira::inline_scheduler, cancellation::none(), [data, result, idx](const expected<T>& exp) mutable {
                bool last = false;
                {
                    std::lock_guard<std::mutex> guard{ data->mutex };

                    data->pendingCount -= 1;
                    last = data->pendingCount == 0;

                    if (exp.has_error() && !data->error)
                    {
                        // set the first error, as it might have cascaded
                        data->error = exp.error();
                    }
                    if (exp.has_value())
                    {
                        data->results[idx] = exp.value();
                    }
                }

                if (last) // we were the last task to complete
                {
                    if (data->error)
                    {
                        result.complete(data->error);
                    }
                    else
                    {
                        result.complete(data->results);
                    }
                }
            });
        }

        return std::move(result);
    }

    template<typename... ArgTs>
    inline task<std::tuple<typename mira::void_passthrough<ArgTs>::type...>> when_all(task<ArgTs>... tasks)
    {
        using void_passthrough_tuple = std::tuple<typename mira::void_passthrough<ArgTs>::type...>;

        struct when_all_data
        {
            std::mutex mutex;
            int pending;
            std::error_code error;
            void_passthrough_tuple results;
        };

        task_completion_source<void_passthrough_tuple> result;

        auto data = std::make_shared<when_all_data>();
        data->pending = std::tuple_size<void_passthrough_tuple>::value;

        std::tuple<task<ArgTs>&...> taskrefs = std::make_tuple(std::ref(tasks)...);

        iterate_tuple(taskrefs, [&](auto& task, auto idx) {
            using task_t = std::remove_reference_t<decltype(task)>;

            task.then(inline_scheduler,
                cancellation::none(),
                [data, result](const expected<typename task_t::result_t>& exp) mutable {
                bool last = false;
                {
                    std::lock_guard<std::mutex> guard{ data->mutex };

                    data->pending -= 1;
                    last = data->pending == 0;

                    internal::write_expected_to_tuple<decltype(idx)::value>(data->results, exp);

                    if (exp.has_error() && !data->error)
                    {
                        // set the first error, as it might have cascaded
                        data->error = exp.error();
                    }
                }

                if (last) // we were the last task to complete
                {
                    if (data->error)
                    {
                        result.complete(data->error);
                    }
                    else
                    {
                        result.complete(std::move(data->results));
                    }
                }
            });
        });
        return std::move(result);
    }
}

max_stars_repo_name: syntheticmagus/mageslam

max_issues_repo_name: syntheticmagus/mageslam

max_forks_repo_name: syntheticmagus/mageslam