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

/** \file utils.h
 *
 *  \brief Contains definition and partial implementation of sirius::Utils class.
 */

#ifndef __UTILS_H__
#define __UTILS_H__

#include <gsl/gsl_sf_erf.h>
#include <fstream>
#include <string>
#include <complex>
#include "sirius_internal.h"
#include "typedefs.h"
#include "constants.h"
#include "sddk.hpp"

/// Utility class.
class Utils
{
    public:

        /// Maximum number of \f$ \ell, m \f$ combinations for a given \f$ \ell_{max} \f$
        static inline int lmmax(int lmax)
        {
            return (lmax + 1) * (lmax + 1);
        }

        static inline int lm_by_l_m(int l, int m)
        {
            return (l * l + l + m);
        }

        static inline int lmax_by_lmmax(int lmmax__)
        {
            int lmax = int(std::sqrt(double(lmmax__)) + 1e-8) - 1;
            if (lmmax(lmax) != lmmax__) {
                TERMINATE("wrong lmmax");
            }
            return lmax;
        }

        static inline bool file_exists(const std::string file_name)
        {
            std::ifstream ifs(file_name.c_str());
            return ifs.is_open();
        }

        static inline double fermi_dirac_distribution(double e)
        {
            double kT = 0.001;
            if (e > 100 * kT) return 0.0;
            if (e < -100 * kT) return 1.0;
            return (1.0 / (exp(e / kT) + 1.0));
        }

        static inline double gaussian_smearing(double e, double delta)
        {
            return 0.5 * (1 - gsl_sf_erf(e / delta));
        }

        static inline double cold_smearing(double e)
        {
            double a = -0.5634;

            if (e < -10.0) return 1.0;
            if (e > 10.0) return 0.0;

            return 0.5 * (1 - gsl_sf_erf(e)) - 1 - 0.25 * exp(-e * e) * (a + 2 * e - 2 * a * e * e) / sqrt(pi);
        }

        static std::string double_to_string(double val, int precision = -1)
        {
            char buf[100];

            double abs_val = std::abs(val);

            if (precision == -1)
            {
                if (abs_val > 1.0)
                {
                    precision = 6;
                }
                else if (abs_val > 1e-14)
                {
                    precision = int(-std::log(abs_val) / std::log(10.0)) + 7;
                }
                else
                {
                    return std::string("0.0");
                }
            }

            std::stringstream fmt;
            fmt << "%." << precision << "f";

            int len = snprintf(buf, 100, fmt.str().c_str(), val);
            for (int i = len - 1; i >= 1; i--)
            {
                if (buf[i] == '0' && buf[i - 1] == '0')
                {
                    buf[i] = 0;
                }
                else
                {
                    break;
                }
            }
            return std::string(buf);
        }

        static inline double phi_by_sin_cos(double sinp, double cosp)
        {
            double phi = std::atan2(sinp, cosp);
            if (phi < 0) {
                phi += twopi;
            }
            return phi;
        }

        static inline long double factorial(int n)
        {
            assert(n >= 0);

            long double result = 1.0L;
            for (int i = 1; i <= n; i++) {
                result *= i;
            }
            return result;
        }

        /// Simple hash function.
        /** Example: printf("hash: %16llX\n", hash()); */
        static uint64_t hash(void const* buff, size_t size, uint64_t h = 5381)
        {
            unsigned char const* p = static_cast<unsigned char const*>(buff);
            for(size_t i = 0; i < size; i++) {
                h = ((h << 5) + h) + p[i];
            }
            return h;
        }

        static void write_matrix(const std::string& fname,
                                 mdarray<double_complex, 2>& matrix,
                                 int nrow,
                                 int ncol,
                                 bool write_upper_only = true,
                                 bool write_abs_only = false,
                                 std::string fmt = "%18.12f")
        {
            static int icount = 0;

            if (nrow < 0 || nrow > (int)matrix.size(0) || ncol < 0 || ncol > (int)matrix.size(1))
                TERMINATE("wrong number of rows or columns");

            icount++;
            std::stringstream s;
            s << icount;
            std::string full_name = s.str() + "_" + fname;

            FILE* fout = fopen(full_name.c_str(), "w");

            for (int icol = 0; icol < ncol; icol++)
            {
                fprintf(fout, "column : %4i\n", icol);
                for (int i = 0; i < 80; i++) fprintf(fout, "-");
                fprintf(fout, "\n");
                if (write_abs_only)
                {
                    fprintf(fout, " row, absolute value\n");
                }
                else
                {
                    fprintf(fout, " row, real part, imaginary part, absolute value\n");
                }
                for (int i = 0; i < 80; i++) fprintf(fout, "-");
                fprintf(fout, "\n");

                int max_row = (write_upper_only) ? std::min(icol, nrow - 1) : (nrow - 1);
                for (int j = 0; j <= max_row; j++)
                {
                    if (write_abs_only)
                    {
                        std::string s = "%4i  " + fmt + "\n";
                        fprintf(fout, s.c_str(), j, abs(matrix(j, icol)));
                    }
                    else
                    {
                        fprintf(fout, "%4i  %18.12f %18.12f %18.12f\n", j, real(matrix(j, icol)), imag(matrix(j, icol)),
                                                                        abs(matrix(j, icol)));
                    }
                }
                fprintf(fout,"\n");
            }

            fclose(fout);
        }

        static void write_matrix(std::string const& fname,
                                 bool write_all,
                                 mdarray<double, 2>& matrix)
        {
            static int icount = 0;

            icount++;
            std::stringstream s;
            s << icount;
            std::string full_name = s.str() + "_" + fname;

            FILE* fout = fopen(full_name.c_str(), "w");

            for (int icol = 0; icol < (int)matrix.size(1); icol++)
            {
                fprintf(fout, "column : %4i\n", icol);
                for (int i = 0; i < 80; i++) fprintf(fout, "-");
                fprintf(fout, "\n");
                fprintf(fout, " row\n");
                for (int i = 0; i < 80; i++) fprintf(fout, "-");
                fprintf(fout, "\n");

                int max_row = (write_all) ? ((int)matrix.size(0) - 1) : std::min(icol, (int)matrix.size(0) - 1);
                for (int j = 0; j <= max_row; j++)
                {
                    fprintf(fout, "%4i  %18.12f\n", j, matrix(j, icol));
                }
                fprintf(fout,"\n");
            }

            fclose(fout);
        }

        static void write_matrix(std::string const& fname,
                                 bool write_all,
                                 matrix<double_complex> const& mtrx)
        {
            static int icount = 0;

            icount++;
            std::stringstream s;
            s << icount;
            std::string full_name = s.str() + "_" + fname;

            FILE* fout = fopen(full_name.c_str(), "w");

            for (int icol = 0; icol < (int)mtrx.size(1); icol++)
            {
                fprintf(fout, "column : %4i\n", icol);
                for (int i = 0; i < 80; i++) fprintf(fout, "-");
                fprintf(fout, "\n");
                fprintf(fout, " row\n");
                for (int i = 0; i < 80; i++) fprintf(fout, "-");
                fprintf(fout, "\n");

                int max_row = (write_all) ? ((int)mtrx.size(0) - 1) : std::min(icol, (int)mtrx.size(0) - 1);
                for (int j = 0; j <= max_row; j++)
                {
                    fprintf(fout, "%4i  %18.12f %18.12f\n", j, real(mtrx(j, icol)), imag(mtrx(j, icol)));
                }
                fprintf(fout,"\n");
            }

            fclose(fout);
        }

        template <typename T>
        static void check_hermitian(const std::string& name, matrix<T> const& mtrx, int n = -1)
        {
            assert(mtrx.size(0) == mtrx.size(1));

            double maxdiff = 0.0;
            int i0 = -1;
            int j0 = -1;

            if (n == -1) {
                n = static_cast<int>(mtrx.size(0));
            }

            for (int i = 0; i < n; i++) {
                for (int j = 0; j < n; j++) {
                    double diff = std::abs(mtrx(i, j) - type_wrapper<T>::conjugate(mtrx(j, i)));
                    if (diff > maxdiff) {
                        maxdiff = diff;
                        i0 = i;
                        j0 = j;
                    }
                }
            }

            if (maxdiff > 1e-10) {
                std::stringstream s;
                s << name << " is not a symmetric or hermitian matrix" << std::endl
                  << "  maximum error: i, j : " << i0 << " " << j0 << " diff : " << maxdiff;

                WARNING(s);
            }
        }

        template <typename T>
        static inline double check_hermitian(dmatrix<T>& mtrx__, int n__)
        {
            dmatrix<T> tmp(n__, n__, mtrx__.blacs_grid(), mtrx__.bs_row(), mtrx__.bs_col());
            linalg<CPU>::tranc(n__, n__, mtrx__, 0, 0, tmp, 0, 0);

            splindex<block_cyclic> spl_r(n__, mtrx__.blacs_grid().num_ranks_row(), mtrx__.blacs_grid().rank_row(), mtrx__.bs_row());
            splindex<block_cyclic> spl_c(n__, mtrx__.blacs_grid().num_ranks_col(), mtrx__.blacs_grid().rank_col(), mtrx__.bs_col());

            double max_diff{0};
            for (int i = 0; i < spl_c.local_size(); i++) {
                for (int j = 0; j < spl_r.local_size(); j++) {
                    max_diff = std::max(max_diff, std::abs(mtrx__(j, i) - tmp(j, i)));
                }
            }
            mtrx__.blacs_grid().comm().template allreduce<double, mpi_op_t::max>(&max_diff, 1);
            return max_diff;
        }

        static double confined_polynomial(double r, double R, int p1, int p2, int dm)
        {
            double t = 1.0 - std::pow(r / R, 2);
            switch (dm)
            {
                case 0: {
                    return (std::pow(r, p1) * std::pow(t, p2));
                }
                case 2: {
                    return (-4 * p1 * p2 * std::pow(r, p1) * std::pow(t, p2 - 1) / std::pow(R, 2) +
                            p1 * (p1 - 1) * std::pow(r, p1 - 2) * std::pow(t, p2) +
                            std::pow(r, p1) * (4 * (p2 - 1) * p2 * std::pow(r, 2) * std::pow(t, p2 - 2) / std::pow(R, 4) -
                                          2 * p2 * std::pow(t, p2 - 1) / std::pow(R, 2)));
                }
                default: {
                    TERMINATE("wrong derivative order");
                    return 0.0;
                }
            }
        }

        static std::vector<int> l_by_lm(int lmax)
        {
            std::vector<int> v(lmmax(lmax));
            for (int l = 0; l <= lmax; l++) {
                for (int m = -l; m <= l; m++) {
                    v[lm_by_l_m(l, m)] = l;
                }
            }
            return std::move(v);
        }

        static std::vector<std::pair<int, int>> l_m_by_lm(int lmax)
        {
            std::vector<std::pair<int, int>> v(lmmax(lmax));
            for (int l = 0; l <= lmax; l++) {
                for (int m = -l; m <= l; m++) {
                    int lm = lm_by_l_m(l, m);
                    v[lm].first = l;
                    v[lm].second = m;
                }
            }
            return std::move(v);
        }

        inline static double round(double a__, int n__)
        {
            double a0 = std::floor(a__);
            double b = std::round((a__ - a0) * std::pow(10, n__)) / std::pow(10, n__);
            return a0 + b;
        }

        template <typename T>
        inline static int sign(T val)
        {
            return (T(0) < val) - (val < T(0));
        }

        static json serialize_timers() // TODO: here in wrong place; move somewhere else
        {
            json dict;

            /* collect local timers */
            for (auto& it: sddk::timer::timers()) {
                sddk::timer::timer_stats ts;

                ts.count = static_cast<int>(it.second.size());
                ts.total_value = 0.0;
                ts.min_value = 1e100;
                ts.max_value = 0.0;
                for (int i = 0; i < ts.count; i++) {
                    ts.total_value += it.second[i];
                    ts.min_value = std::min(ts.min_value, it.second[i]);
                    ts.max_value = std::max(ts.max_value, it.second[i]);
                }
                ts.average_value = (ts.count == 0) ? 0.0 : ts.total_value / ts.count;
                if (ts.count == 0) {
                    ts.min_value = 0.0;
                }

                dict[it.first] = {ts.total_value, ts.average_value, ts.min_value, ts.max_value};
            }
            return std::move(dict);
        }

};

#endif

max_stars_repo_name: cocteautwins/SIRIUS-develop

max_issues_repo_name: cocteautwins/SIRIUS-develop

max_forks_repo_name: cocteautwins/SIRIUS-develop


text:
#ifndef NUFFTW_PLAN_H
#define NUFFTW_PLAN_H

#include <complex>
#include <cmath>
#include <limits>               // st::numeric_limits
#include <algorithm>            // std::max
#include <gsl/gsl_sf_bessel.h>  // Bessel functions
#include <gsl/gsl_sf_lambert.h> // Lambert W functions
#include <fftw3.h>
#include "timer.h"

namespace nufftw
{
    using namespace std::complex_literals;
    typedef std::complex<double> complex;

    struct options
    {
        double tol = std::numeric_limits<double>::epsilon();
        unsigned int fftw_flags = FFTW_ESTIMATE;
    };

    namespace details
    {
        void cheb_vals(int n, int p, const double* x, double* vals);
        void bessel_coeffs(int r, double gam, complex* cfs);
    }

    class plan_base
    {
    	public:
    		plan_base(const int* ns_, complex* in_, complex* out_, const options& opts_ = options()) :
                in(in_), out(out_), opts(opts_)
            {
                n=1;
                for (int i=0; i<N; ++i) {
                    ns[i]=ns_[i];
                    n*=ns[i];
                }
            }

            virtual ~plan_base()
            {
                fftw_destroy_plan(fp);
            }

    		virtual void execute() = 0;

        protected:
            int m;               // Number of nonuniform points
            int ns[N];           // Size of the problem in each dimension
            int n;               // Total size of the problem
            complex* in;         // Input coefficients
            complex* out;        // Output values
            fftw_plan fp;        // FFTW plans
            options opts;        // Options
            double gam;          // Perturbation parameter
            int r;               // Rank parameter
    };

    class plan1 : public plan_base<N>
    {
        public:
            plan1(const int* n_, complex* in_, complex* out_, const double* omega_, const options& opts_ = options()) :
                plan_base<N>(n_, in_, out_, opts_), omega(omega_)
            {}
            virtual ~plan1();
            virtual void execute();
        protected:
            const double* omega;
            plan1<1>* plans;
    };

    class plan1_2d : public plan_base
    {

    };

    class plan1_1d : public plan_base
    {
        public:
            plan1(int n_, complex* in_, complex* out_, const double* omega_, const options& opts_ = options()) :
                plan_base<1>(&n_, in_, out_, opts_), omega(omega_), t(new int[n]),
                er(new double[n]), temp(new double[n])
            {
                timer::tic();
                gam = 0;
                for (int i=0; i<n; ++i) {
                    int s = std::round(omega[i]);
                    t[i] = s%n;
                    er[i] = omega[i]-s;
                    gam = std::max(gam, std::abs(er[i]));
                }
                timer::toc("Compute gamma");

                timer::tic();
                r = 1;
                if (gam>opts.tol) {
                    // Asymptotic approximation to Lambert-W:
                    //   double xi = std::log(std::log(10./opts.tol)/(7.*gam);
                    //   double log_xi=std::log(xi), inv_xi=1./xi;
                    //   double lw = xi-log_xi*(1+inv_xi*(1+inv_xi*(0.5*log_xi-1)));
                    double lw = gsl_sf_lambert_W0(std::log(10./opts.tol)/(7.*gam));
                    r = std::ceil(5.*gam*std::exp(lw));
                }
                timer::toc("Compute r");

                timer::tic();
                // Plan the FFTs
                fft_data = reinterpret_cast<complex*>(fftw_alloc_complex(n*r));
                // int stride = r;
                int stride = 1;
                // int dist = 1;
                int dist = n;
                fp = fftw_plan_many_dft(1, &n, r,
                    reinterpret_cast<fftw_complex*>(fft_data), NULL, stride, dist,
                    reinterpret_cast<fftw_complex*>(fft_data), NULL, stride, dist,
                    FFTW_BACKWARD, opts.fftw_flags
                );
                timer::toc("Plan the FFTs");

                // Construct a low rank approximation, using Chebyshev expansions
                // for A_K = exp(-2*pi*1im*(x[j]-j/N)*k):

                // Construct a low rank approximation to
                //     A_{jk} = exp(-2*pi*1i*(x_j-s_j/N)*k), 0<=j,k<=N-1,
                // where |x_j-j/N|<= gam <=1/2.  See [1].
                u = new complex[n*r];
                v = new double[n*r];

                // Compute u
                timer::tic();
                double* cheb = new double[n*r];
                complex* bess = new complex[r*r];
                timer::tic();
                double fac = 1./gam;
                for (int i=0; i<n; ++i) temp[i]=er[i]*fac;
                details::cheb_vals(n, r, temp, cheb);
                timer::toc("Chebyshev");
                timer::tic();
                details::bessel_coeffs(r, gam, bess);
                timer::toc("Bessel");
                timer::tic();
                for (int i=0; i<n; ++i) {
                    complex scl = std::exp(-1i*M_PI*er[i]);
                    for (int j=0; j<r; ++j) {
                        u[r*i+j]=0;
                        for (int k=0; k<r; ++k) {
                            u[r*i+j]+=cheb[r*i+k]*bess[r*j+k];
                        }
                        u[r*i+j]*=scl;
                    }
                }
                timer::toc("Matrix multiply");
                delete [] bess;
                delete [] cheb;
                timer::toc("Compute u");

                // Compute v
                timer::tic();
                fac = 2./n;
                temp[0]=-1;
                for (int i=1; i<n; ++i) temp[i]=temp[i-1]+fac;
                details::cheb_vals(n, r, temp, v);
                timer::toc("Compute v");
            }

            virtual ~plan1()
            {
                delete [] t;
                delete [] er;
                delete [] u;
                delete [] v;
                delete [] temp;
                fftw_free(reinterpret_cast<fftw_complex*>(fft_data));
            }

            virtual void execute()
            {
                for (int i=0; i<n*r; ++i) fft_data[i]=0;
                for (int i=0; i<n; ++i) {
                    for (int j=0; j<r; ++j) {
                        // fft_data[r*t[i]+j] += std::conj(in[i]*u[r*i+j]);
                        fft_data[n*j+t[i]] += std::conj(in[i]*u[r*i+j]);
                    }
                }

                // Do the FFTs
                timer::tic();
                fftw_execute(fp); // normalize?
                timer::toc("Call FFTW");

                for (int i=0; i<n; ++i) {
                    out[i]=0;
                    for (int j=0; j<r; ++j) {
                        // out[i] += v[r*i+j]*std::conj(fft_data[r*i+j]);
                        out[i] += v[r*i+j]*std::conj(fft_data[n*j+i]);
                    }
                }
            }

        protected:
            const double* omega; // FFT samples
            int* t;              // Closest equispaced FFT samples to omega
            double* er;          // Perturbation error
            complex* u;          // Low rank approximation
            double* v;           // Low rank approximation
            double* temp;        // Temporary storage
            complex* fft_data;   // Data for FFTW
    };

    // template<int N>
    // class plan2 : public plan_base<N>
    // {
    //     public:
    //         plan2();
    //         virtual ~plan2();
    // };

    // template<>
    // class plan2<1> : public plan_base<1>
    // {
    //     public:
    //         plan2();
    //         virtual ~plan2();
    // };

    // template<int N>
    // class plan3 : public plan_base<N>
    // {
    //     public:
    //         plan3();
    //         virtual ~plan3();
    // };

    // template<>
    // class plan3<1> : public plan_base<1>
    // {
    //     public:
    //         plan3();
    //         virtual ~plan3();
    // };

    typedef plan_base* plan;

    // Methods for planning
    template<int N>
    void destroy_plan(plan<N> p)
    {
        if (p) delete p;
    }

    template<int N>
    void execute(const plan<N> p)
    {
        p->execute();
    }

    namespace details
    {
        /** Evaluate Chebyshev polynomials of degree 0,...,p-1 at x.
         *  (row-major order) */
        void cheb_vals(int n, int p, const double* x, double* vals)
        {
            for (int i=0; i<n; ++i) {
                vals[p*i]=1;
                vals[p*i+1]=x[i];
            }
            for (int i=0; i<n; ++i) {
                for (int j=1; j<p-1; ++j) {
                    vals[p*i+(j+1)] = 2.*x[i]*vals[p*i+j] - vals[p*i+(j-1)];
                }
            }
        }

        /** The bivarate Chebyshev coefficients for the function
         *  f(x,y) = exp(-i*x.*y) on the domain [-gam, gam]x[0,2*pi], as given
         *  by Lemma A.2 of Townsend's DPhil thesis. (column-major order) */
        void bessel_coeffs(int r, double gam, complex* cfs)
        {
            double arg = -0.5*gam*M_PI;
            for (int p=0; p<r; ++p) {
                for (int q=p%2; q<r; q+=2) {
                    cfs[r*q+p] = 4.*std::pow(1i,q)*gsl_sf_bessel_Jn((p+q)*0.5,arg)*gsl_sf_bessel_Jn((q-p)*0.5,arg);
                }
            }
            for (int p=0; p<r; ++p) cfs[p]*=0.5;
            for (int q=0; q<r; ++q) cfs[r*q]*=0.5;
        }
    }
}

#endif

ext: h

max_stars_repo_name: danfortunato/nufftw

max_issues_repo_name: danfortunato/nufftw

max_forks_repo_name: danfortunato/nufftw


text:
#ifndef AMICI_HDF5_H
#define AMICI_HDF5_H

#include <string>
#include <memory>
#include <vector>

#include <H5Cpp.h>

#include <gsl/gsl-lite.hpp>

#undef AMI_HDF5_H_DEBUG

/* Macros for enabling/disabling  HDF5 error auto-printing
 * AMICI_H5_SAVE_ERROR_HANDLER and AMICI_H5_RESTORE_ERROR_HANDLER must be called
 * within the same context, otherwise the stack handler is lost. */
#define AMICI_H5_SAVE_ERROR_HANDLER                                            \
    herr_t (*old_func)(void *);                                                \
    void *old_client_data;                                                     \
    H5Eget_auto1(&old_func, &old_client_data);                                 \
    H5Eset_auto1(NULL, NULL)

#define AMICI_H5_RESTORE_ERROR_HANDLER H5Eset_auto1(old_func, old_client_data)

namespace amici {

class ReturnData;
class ExpData;
class Model;
class Solver;

namespace hdf5 {

/* Functions for reading and writing AMICI data to/from HDF5 files. */

/**
 * @brief Open the given file for writing. Append if exists, create if not.
 * @param hdf5filename
 * @return
 */
H5::H5File createOrOpenForWriting(std::string const& hdf5filename);

/**
 * @brief Read solver options from HDF5 file
 * @param fileId hdf5 file handle to read from
 * @param solver solver to set options on
 * @param datasetPath Path inside the HDF5 file
 */
void readSolverSettingsFromHDF5(const H5::H5File &file, Solver& solver,
                                std::string const& datasetPath);

/**
 * @brief Read solver options from HDF5 file
 * @param hdffile Name of HDF5 file
 * @param solver solver to set options on
 * @param datasetPath Path inside the HDF5 file
 */
void readSolverSettingsFromHDF5(std::string const& hdffile, Solver& solver,
                                std::string const& datasetPath);

/**
 * @brief Read model data from HDF5 file
 * @param hdffile Name of HDF5 file
 * @param model model to set data on
 * @param datasetPath Path inside the HDF5 file
 */
void readModelDataFromHDF5(std::string const& hdffile, Model& model,
                           std::string const& datasetPath);

/**
 * @brief Read model data from HDF5 file
 * @param fileId hdf5 file handle to read from
 * @param model model to set data on
 * @param datasetPath Path inside the HDF5 file
 */
void readModelDataFromHDF5(H5::H5File const&file, Model& model,
                           std::string const& datasetPath);

/**
  * @brief Write ReturnData struct to HDF5 dataset
  * @param rdata Data to write
  * @param hdffile Filename of HDF5 file
  * @param datasetPath Full dataset path inside the HDF5 file (will be created)
  */

void writeReturnData(const ReturnData &rdata,
                     H5::H5File const& file,
                     const std::string& hdf5Location);

void writeReturnData(const ReturnData &rdata,
                     std::string const& hdf5Filename,
                     const std::string& hdf5Location);

void writeReturnDataDiagnosis(const ReturnData &rdata,
                              H5::H5File const& file,
                              const std::string& hdf5Location);

/**
 * @brief Create the given group and possibly parents
 * @param file
 * @param groupPath
 * @param recursively
 */
void createGroup(const H5::H5File &file,
                 std::string const& groupPath,
                 bool recursively = true);

/**
 * @brief readSimulationExpData reads AMICI experimental data from
 * attributes in HDF5 file.
 * @param hdf5Filename Name of HDF5 file
 * @param hdf5Root Path inside the HDF5 file to object having ExpData as
 * attributes
 * @param model The model for which data is to be read
 * @return
 */

std::unique_ptr<ExpData> readSimulationExpData(const std::string &hdf5Filename,
                                               const std::string &hdf5Root,
                                               const Model &model);

/**
 * @brief writeSimulationExpData writes AMICI experimental data to
 * attributes in HDF5 file.
 * @param edata The experimental data which is to be written
 * @param hdf5Filename Name of HDF5 file
 * @param hdf5Root Path inside the HDF5 file to object having ExpData as
 * attributes
 */

void writeSimulationExpData(const ExpData &edata,
                            H5::H5File const& file,
                            const std::string &hdf5Location);

/**
 * @brief attributeExists Check whether an attribute with the given
 * name exists on the given dataset
 * @param fileId The HDF5 file object
 * @param datasetPath Dataset of which attributes should be checked
 * @param attributeName Name of the attribute of interest
 * @return
 */
bool attributeExists(H5::H5File const& file,
                     const std::string &optionsObject,
                     const std::string &attributeName);

bool attributeExists(H5::H5Object const& object,
                     const std::string &attributeName);

void createAndWriteInt1DDataset(H5::H5File const& file,
                                std::string const& datasetName,
                                gsl::span<const int> buffer);

void createAndWriteInt2DDataset(H5::H5File const& file,
                                std::string const& datasetName,
                                gsl::span<const int> buffer, hsize_t m,
                                hsize_t n);

void createAndWriteDouble1DDataset(H5::H5File const& file,
                                   std::string const& datasetName,
                                   gsl::span<const double> buffer);

void createAndWriteDouble2DDataset(H5::H5File const& file,
                                   std::string const& datasetName,
                                   gsl::span<const double> buffer, hsize_t m,
                                   hsize_t n);

void createAndWriteDouble3DDataset(H5::H5File const& file,
                                   std::string const& datasetName,
                                   gsl::span<const double> buffer, hsize_t m,
                                   hsize_t n, hsize_t o);

double getDoubleScalarAttribute(const H5::H5File& file,
                                const std::string &optionsObject,
                                const std::string &attributeName);

int getIntScalarAttribute(const H5::H5File &file,
                          const std::string &optionsObject,
                          const std::string &attributeName);

std::vector<int> getIntDataset1D(const H5::H5File &file,
                                 std::string const& name);

std::vector<double> getDoubleDataset1D(const H5::H5File &file,
                                       std::string const& name);

std::vector<double> getDoubleDataset2D(const H5::H5File &file,
                                       std::string const& name,
                                       hsize_t &m, hsize_t &n);

std::vector<double> getDoubleDataset3D(const H5::H5File &file,
                                       std::string const& name,
                                       hsize_t &m, hsize_t &n, hsize_t &o);

/**
 * @brief Check if the given location (group, link or dataset) exists in the
 * given file
 * @param filename
 * @param location
 * @return
 */
bool locationExists(std::string const& filename, std::string const& location);

bool locationExists(H5::H5File const& file, std::string const& location);

} // namespace hdf5
} // namespace amici

#endif

max_stars_repo_name: paszkow/AMICI

max_issues_repo_name: paszkow/AMICI

max_forks_repo_name: paszkow/AMICI


text:
#pragma once

#include "options.h"
#include "predict.h"
#include "progress.h"
#include <gsl/gsl-lite.hpp>

namespace angonoka::cli {
namespace detail {
    /**
        Check if an event is the final one.

        @param evt Event

        @return True if this is the final event.
    */
    bool is_final_event(ProgressEvent& evt) noexcept;

    /**
        Helper constant that checks if a given class is
        a specialization of std::future.
    */
    template <typename T> inline constexpr bool is_future = false;
    template <typename T>
    inline constexpr bool is_future<std::future<T>> = true;
} // namespace detail

/**
    Prediction events handler.

    Prints various progress messages and results as
    the prediction algorithm runs.

    @var progress Chosen progress bar implementation
    @var options CLI options
*/
struct EventHandler {
    gsl::not_null<Progress*> progress;
    gsl::not_null<const Options*> options;

    /**
        Handler events without attributes.

        @param e Event
    */
    void operator()(const SimpleProgressEvent& e) const;

    /**
        Handle schedule optimization events.

        @param e Event
    */
    void operator()(const ScheduleOptimizationEvent& e) const;

    /**
        Handle schedule optimization completion.

        @param e Event
    */
    void operator()(const ScheduleOptimizationComplete& e) const;
};

/**
    Awaitable result of the prediction.
*/
template <typename T>
concept Prediction = detail::is_future<T>;

/**
    Consumes prediction events from the queue.

    @param queue        Prediction events queue
    @param prediction   Prediction result
    @param handler      Events handler
*/
void consume_events(
    Queue<ProgressEvent>& queue,
    Prediction auto& prediction,
    EventHandler handler)
{
    Expects(prediction.valid());

    using namespace std::literals::chrono_literals;
    using boost::variant2::visit;
    constexpr auto event_timeout = 100ms;
    ProgressEvent evt;
    while (!detail::is_final_event(evt)) {
        if (!queue.try_dequeue(evt)) {
            prediction.wait_for(event_timeout);
            continue;
        }
        visit(handler, evt);
    }

    Ensures(prediction.valid());
}
} // namespace angonoka::cli

max_stars_repo_name: coffee-lord/angonoka

max_issues_repo_name: coffee-lord/angonoka

max_forks_repo_name: coffee-lord/angonoka