text:
#if 1
#   include <boost/regex.hpp>
#   define REGEX_NAMESPACE boost
#else
#   include <regex>
#   define REGEX_NAMESPACE std
#endif

#include <fstream>
#include <numeric>
#include <boost/algorithm/string.hpp>
#include <boost/math/special_functions/gamma.hpp>   // lgamma
#include "utilfunc.hpp"
#include "arbor.hpp"
#include "xarbor.hpp"

using namespace arbor;

void Arbor::processCommandLineOptions(int argc, const char * argv[])
    {
    // Get user-supplied configuration settings from command line and/or config file
    boost::program_options::variables_map       vm;
    boost::program_options::options_description desc("Allowed options");
    desc.add_options()
        ("help,h",                                                                                "produce help message")
        ("version,v",                                                                                "show program version")
        ("datafile,d",      boost::program_options::value<std::string>(&_data_file_name),                                    "name of parameter sample file")
        ("grapefraction,f", boost::program_options::value< double >(&_grape_fraction)->default_value(_def_grapefraction),    "fraction of sample to use in forming grapes")
        ("minsamplesize,m", boost::program_options::value< unsigned >(&_min_sample_size)->default_value(_def_minsamplesize), "minimum number of sampled points for tree topology to be considered")
        ("seed,s",          boost::program_options::value< unsigned >(&_rnseed)->default_value(_def_rnseed),                 "seed for pseudorandom number generation")
        ;
    boost::program_options::store(boost::program_options::parse_command_line(argc, argv, desc), vm);
    boost::program_options::notify(vm);

    // If user specified --help on command line, output usage summary and quit
    if (vm.count("help") > 0)
        {
        std::cout << desc << "\n";
        std::exit(1);
        }

    // If user specified --version on command line, output version and quit
    if (vm.count("version") > 0)
        {
        std::cout << boost::str(boost::format("This is %s version %d.%d") % _program_name % _major_version % _minor_version) << std::endl;
        std::exit(1);
        }

    // Ensure that user specified --datafile on command line
    if (vm.count("datafile") == 0)
        {
        std::cout << "Must specify --datafile on command line\n";
        std::cout << desc << std::endl;
        std::exit(1);
        }

    // If user used --seed on command line, perform sanity check to ensure that value specified was greater than 0
    if (vm.count("seed") > 0)
        {
        if (_rnseed <= 0.0)
            {
            std::cout << boost::str(boost::format("Error: specified pseudorandom number seed (%.3f) should be greater than 0.") % _rnseed) << std::endl;
            std::exit(1);
            }
        }
    // If user used --grapefraction on command line, perform sanity check to ensure that value specified was greater than or equal to 0 and less than or equal to 1
    if (vm.count("grapefraction") > 0)
        {
        if (_grape_fraction <= 0.0 || _grape_fraction >= 1.0)
            {
            std::cout << boost::str(boost::format("Error: specified grapefraction (%.3f) should be between 0 and 1.") % _grape_fraction) << std::endl;
            std::exit(1);
            }
        }

    // If user used --minsamplesize on command line, perform sanity check to ensure that value specified was greater than 0
    if (vm.count("minsamplesize") > 0)
        {
        if (_min_sample_size <= 0)
            {
            std::cout << boost::str(boost::format("Error: specified minsamplesize (%d) should be greater than 0.") % _min_sample_size) << std::endl;
            std::exit(1);
            }
        }

    }

void Arbor::showUserSpecifiedOptions() const
    {
    std::cout << "\nUser-configurable settings:" << std::endl;
    std::cout << boost::str(boost::format("  seed:               %d")   % _rnseed) << std::endl;
    std::cout << boost::str(boost::format("  datafile:           %s")   % _data_file_name) << std::endl;
    std::cout << boost::str(boost::format("  minsamplesize:      %d")   % _min_sample_size) << std::endl;
    std::cout << boost::str(boost::format("  grapefraction:      %.3f") % _grape_fraction) << std::endl;
    }

void Arbor::calcParamTypesFromColumnHeaders()
    {
    // Fill in _param_types and _param_subsets vectors using strings in _column_headers vector

    std::string subsetstr;
    REGEX_NAMESPACE::smatch what;
    REGEX_NAMESPACE::regex rex_shape("log\\((\\d+)_gamma_shape\\)");
    REGEX_NAMESPACE::regex rex_pinvar("logit\\((\\d+)_pinvar\\)");
    REGEX_NAMESPACE::regex rex_exchangeability("log\\((\\d+)_r([ACGT][ACGT])/(\\d+)_rAC\\)");
    REGEX_NAMESPACE::regex rex_frequency("log\\((\\d+)_freq([ACGT])/(\\d+)_freqA\\)");

    _nsubsets = 0;
    _param_types.clear();
    _param_subsets.clear();
    for (unsigned i = 4; i < _column_headers.size(); ++i)
        {
        std::string & h = _column_headers[i];
        if (REGEX_NAMESPACE::regex_match(h, what, rex_shape))
            {
            _param_types.push_back(ptype_gamma_shape);
            subsetstr = std::string(what[1].first, what[1].second);
            unsigned subset = boost::lexical_cast<unsigned>(subsetstr);
            _param_subsets.push_back(subset);
            if (subset > _nsubsets)
                _nsubsets = subset;
            }
        else if (REGEX_NAMESPACE::regex_match(h, what, rex_pinvar))
            {
            _param_types.push_back(ptype_pinvar);
            subsetstr = std::string(what[1].first, what[1].second);
            unsigned subset = boost::lexical_cast<unsigned>(subsetstr);
            _param_subsets.push_back(subset);
            if (subset > _nsubsets)
                _nsubsets = subset;
            }
        else if (REGEX_NAMESPACE::regex_match(h, what, rex_exchangeability))
            {
            _param_types.push_back(ptype_exchangeability);
            subsetstr = std::string(what[1].first, what[1].second);
            unsigned subset = boost::lexical_cast<unsigned>(subsetstr);
            _param_subsets.push_back(subset);
            if (subset > _nsubsets)
                _nsubsets = subset;
            }
        else if (REGEX_NAMESPACE::regex_match(h, what, rex_frequency))
            {
            _param_types.push_back(ptype_frequency);
            subsetstr = std::string(what[1].first, what[1].second);
            unsigned subset = boost::lexical_cast<unsigned>(subsetstr);
            _param_subsets.push_back(subset);
            if (subset > _nsubsets)
                _nsubsets = subset;
            }
        else
            {
            _param_types.push_back(ptype_edge_length);
            _param_subsets.push_back(0);
            }
        }

    _nparams = (unsigned)_param_types.size();
    }

void Arbor::readData()
    {
    // Create a temporary map for storing frequency of each topology
    typedef std::map< unsigned, unsigned > topology_freq_t;
    topology_freq_t topology_freq;

    // Clear containers
    _column_headers.clear();
    _log_likelihoods.clear();
    _log_priors.clear();
    _parameters.clear();

    // Open the data file
    std::ifstream inf(_data_file_name.c_str());
    if (!inf)
        throw XArbor(boost::str(boost::format("Could not open file named \"%s\"") % _data_file_name));

    // Read contents of file into string
    std::string contents((std::istreambuf_iterator<char>(inf)), std::istreambuf_iterator<char>());

    // Regular expression that splits contents into separate lines
    const char * rex = "[\\r\\n]+";

    // Set up regular expression reader to iterate through rows
    REGEX_NAMESPACE::regex row_expr(rex);
    REGEX_NAMESPACE::sregex_token_iterator m(contents.begin(), contents.end(), row_expr, -1);
    REGEX_NAMESPACE::sregex_token_iterator m2;

    for (unsigned row = 0; m != m2; ++m, ++row)
        {
        std::string rowstr = m->str();
        boost::algorithm::trim(rowstr);
        std::vector< std::string > s;
        boost::split(s, rowstr, boost::is_any_of("\t"), boost::token_compress_on);

        unsigned nelements = (unsigned)s.size();

        if (row == 0) {
            // first row should be column headers
            _column_headers.resize(nelements);
            std::copy(s.begin(), s.end(), _column_headers.begin());

            calcParamTypesFromColumnHeaders();
        }
        else if (row > _skip) {
            // s[0] is the iteration: we ignore this

            // s[1] is the tree topology: this will be used as a key for _log_likelihoods, _log_priors, and _parameters
            unsigned tree_topology = boost::lexical_cast<unsigned>(s[1]);

            // Update the frequency of this tree topology
            topology_freq_t::iterator lower_bound = topology_freq.lower_bound(tree_topology);
            if (lower_bound != topology_freq.end() && !(topology_freq.key_comp()(tree_topology, lower_bound->first)))
                {
                // this tree topology has already been seen, so increment the count
                lower_bound->second += 1;
                }
            else
                {
                // this tree topology has not yet been seen, so create a new entry
                topology_freq.insert(lower_bound, topology_freq_t::value_type(tree_topology, 1));
                }

            // s[2] is the log-likelihood: store this in _log_likelihoods map
            double logL = boost::lexical_cast<double>(s[2]);
            _log_likelihoods[tree_topology].push_back(logL);

            // s[3] is the log-prior: store this in _log_priors map
            double logP = boost::lexical_cast<double>(s[3]);
            _log_priors[tree_topology].push_back(logP);

            // s[4] is the first log-transformed parameter: store all parameters in _parameters map
            std::vector<double> tmp;
            for (unsigned i = 4; i < nelements; ++i)
                {
                double v = boost::lexical_cast<double>(s[i]);
                tmp.push_back(v);
                }
            _parameters[tree_topology].push_back(tmp);

        }   // else row > 0
    } // loop over rows

    // Copy entries of topology_freq map to the vector _tree_topologies, each element of which is a (frequency,topology) pair
    _tree_topologies.clear();
    for (auto p : topology_freq)
        {
        _tree_topologies.push_back(std::make_pair(p.second,p.first));
        }

    // Sort the entries of _tree_topologies so that topologies with highest frequency come first
    std::sort(_tree_topologies.begin(), _tree_topologies.end(), std::greater<size_topol_pair_t>());

    // _nparams should equal the length of the first parameter vector for the most frequently sampled tree topology
    assert(_nparams == _parameters[_tree_topologies[0].second][0].size());
    }

void Arbor::initRandom()
    {
    _lot = Lot::SharedPtr(new Lot());
    _lot->setSeed(_rnseed);
    }

void Arbor::showParamTable() const
    {
    std::cout << std::endl;
    std::cout << boost::str(boost::format("  %12s %12s   %s") % "parameter" % "data subset" % "parameter type") << std::endl;
    for (unsigned k = 0; k < _param_types.size(); ++k)
        {
        if (_param_types[k] == 0)
            std::cout << boost::str(boost::format("  %12d %12s   %s") % (k+1) % "NA" % "edge length") << std::endl;
        else if (_param_types[k] == 1)
            std::cout << boost::str(boost::format("  %12d %12d   %s") % (k+1) % _param_subsets[k] % "gamma shape") << std::endl;
        else if (_param_types[k] == 2)
            std::cout << boost::str(boost::format("  %12d %12d   %s") % (k+1) % _param_subsets[k] % "proportion of invariable sites") << std::endl;
        else if (_param_types[k] == 3)
            std::cout << boost::str(boost::format("  %12d %12d   %s") % (k+1) % _param_subsets[k] % "GTR exchangeability") << std::endl;
        else if (_param_types[k] == 4)
            std::cout << boost::str(boost::format("  %12d %12d   %s") % (k+1) % _param_subsets[k] % "state frequency") << std::endl;
        else
            std::cout << boost::str(boost::format("  %12d %12s   %s") % (k+1) % "NA" % "unknown") << std::endl;
        }
    }

void Arbor::chooseGrapeSeeds(sample_vect_t & parameters)
    {
    // Builds two vectors of indices into _parameters vector:
    // _reference_indices holds indices of points used as grape seeds (a random proportion _grape_fraction of all points in _parameters);
    // _estimation_indices holds indices of all remaining points.

    // Clear vectors that will be used to store indices
    _reference_indices.clear();
    _estimation_indices.clear();

    for (unsigned i = 0; i < parameters.size(); ++i)
        {
        if (_lot->uniform() < _grape_fraction)
            _reference_indices.push_back(i);
        else
            _estimation_indices.push_back(i);
        }

    // Add to tally of total number of estimation points used over all topologies
    _N += (unsigned)_estimation_indices.size();
    }

void Arbor::calcMeansAndStdevs(sample_vect_t & parameters)
    {
    // Computes means and standard deviations of each parameter for the points whose indices are in _reference_indices.
    assert(_nparams > 0);

    // Clear and zero vectors that will be used to store means and standard deviations
    _means.resize(_nparams);
    _means.assign(_nparams, 0.0);
    _stdevs.resize(_nparams);
    _stdevs.assign(_nparams, 0.0);

    // Loop through all vectors in reference sample
    unsigned reference_sample_size = (unsigned)_reference_indices.size();
    for (auto i : _reference_indices)
        {
        std::vector<double> & v = parameters[i];

        // Update sum stored in _means and sum-of-squares stored in _stdevs
        for (unsigned j = 0; j < _nparams; ++j)
            {
            _means[j]  += v[j];
            _stdevs[j] += pow(v[j],2);
            }
        }

    // Complete _means and _stdevs calculation
    for (unsigned i = 0; i < _nparams; ++i)
        {
        _means[i] /= reference_sample_size;
        double tmp = (_stdevs[i] - pow(_means[i],2)*reference_sample_size)/(reference_sample_size - 1);
        _stdevs[i] = sqrt(tmp);
        }
    }

void Arbor::standardizeSampleAllowingCorrelation(sample_vect_t & parameters)
    {
    // Computes variance-covariance matrix for the points whose indices are in _reference_indices and standardizes sample.
    assert(_nparams > 0);

    unsigned n = (unsigned)parameters.size();
    unsigned p = _nparams;

    // calculate _means
    _Means.resize(p);
    _Means.setZero();
    for (auto & sample : parameters)
        {
        for (unsigned j = 0; j < p; ++j)
            _Means(j) += sample[j];
        }
    for (unsigned j = 0; j < p; ++j)
        _Means(j) /= n;

    // Store observations in matrix _Y and means in _M
    _Y.resize(n,p);
    _M.resize(n,p);
    unsigned i = 0;
    for (auto & sample : parameters)
        {
        for (unsigned j = 0; j < p; ++j)
            {
            _Y(i,j) = sample[j];
            _M(i,j) = _Means(j);
            }
        ++i;
        }

    // Let M hold the centered observations
    // TODO: should avoid building _M as it repeats the mean vector n times
    Eigen::MatrixXd M = _Y - _M;

    // Compute variance-covariance matrix V
    Eigen::MatrixXd V = M.transpose()*M/(n-1);

    // Compute square root matrix _S = V^{0.5} and inverse square root matrix _Sinv = V^{-0.5}
    Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eigensolver(V);
    if (eigensolver.info() != Eigen::Success)
        {
        if (eigensolver.info() == Eigen::NumericalIssue)
            throw XArbor("could not compute eigendecomposition of variance-covariance matrix due to a numerical issue");
        else if (eigensolver.info() == Eigen::NoConvergence)
            throw XArbor("could not compute eigendecomposition of variance-covariance matrix due to no convergence");
        else if (eigensolver.info() == Eigen::InvalidInput)
            throw XArbor("could not compute eigendecomposition of variance-covariance matrix due to invalid input");
        else
            throw XArbor("could not compute eigendecomposition of variance-covariance matrix due to some unknown cause");
        }

    // Compute square root and inverse square root of eigenvalues
    Eigen::VectorXd eval = eigensolver.eigenvalues();
    Eigen::VectorXd eval_sqrt = eval;
    Eigen::VectorXd eval_sqrt_inv = eval;
    unsigned index_of_largest_eval = 0;
    double largest_eval = 0.0;
    for (i = 0; i < p; ++i)
        {
        if (eval[i] > largest_eval) {
            largest_eval = eval[i];
            index_of_largest_eval = i;
        }
        eval_sqrt[i]     = sqrt(eval[i]);
        eval_sqrt_inv[i] = pow(eval[i], -0.5);
        }

    // Obtain eigenvectors
    Eigen::MatrixXd evec = eigensolver.eigenvectors();

    // Save eigenvector with largest eigenvalue
    _primary_eigenvector = evec.col(index_of_largest_eval);

    // Compute _S (square root of variance-covariance matrix)
    _S = evec * eval_sqrt.asDiagonal() * evec.inverse();

    // Compute _Sinv (inverse of square root of variance-covariance matrix)
    _Sinv = _S.inverse();

    // Let matrix X hold the standardized samples
    Eigen::MatrixXd X = M * _Sinv;

    // Let _logdetS hold the log of |determinant of _S|
    // |det(S)| is the Jacobian for the standardization transformation
    _logdetS = log(fabs(_S.determinant()));

    // These vectors will be used to compute the denominator sums for jacobians of dirichlet-like parameters
    std::vector<double> freqdenom(_nsubsets, 1.0);
    std::vector<double> exchgdenom(_nsubsets, 1.0);
    std::vector<double> pinvdenom(_nsubsets, 1.0);

    // Copy values in X back to parameters and compute log jacobian terms
    i = 0;
    _log_jacobians.clear();
    for (auto & sample : parameters)
        {
        freqdenom.assign(_nsubsets, 1.0);
        exchgdenom.assign(_nsubsets, 1.0);
        pinvdenom.assign(_nsubsets, 1.0);

        double logj = 0.0;
        for (unsigned j = 0; j < p; ++j)
            {
            double value = sample[j];

            // log-jacobian
            logj += value;   // term in log-jacobian

            if (_param_types[j] == ptype_pinvar)
                pinvdenom[_param_subsets[j] - 1] += exp(value);
            else if (_param_types[j] == ptype_exchangeability)
                exchgdenom[_param_subsets[j] - 1] += exp(value);
            else if (_param_types[j] == ptype_frequency)
                freqdenom[_param_subsets[j] - 1] += exp(value);

            sample[j] = X(i,j);
            }

        logj += _logdetS;
        for (unsigned k = 0; k < _nsubsets; ++k)
            {
            logj -= 2.0*log(pinvdenom[k]);
            logj -= 4.0*log(freqdenom[k]);
            logj -= 6.0*log(exchgdenom[k]);
            }

        _log_jacobians.push_back(logj);

        ++i;
        }
    }

void Arbor::standardizeSampleAssumingIndependence(sample_vect_t & parameters)
    {
    // Standardizes vectors in parameters by subtracting mean and dividing by standard deviation for each parameter separately.

    // Sanity checks: should have stored sample in _parameters and calculated _means and _stdevs
    assert(_nparams > 0);
    assert(_means.size() == _nparams);
    assert(_stdevs.size() == _nparams);

    double sum_log_sd = std::accumulate(_stdevs.begin(), _stdevs.end(), 0.0, pluslog());
    _log_jacobians.clear();

    // These vectors will be used to compute the denominator sums for jacobians of dirichlet-like parameters
    std::vector<double> freqdenom(_nsubsets, 1.0);
    std::vector<double> exchgdenom(_nsubsets, 1.0);
    std::vector<double> pinvdenom(_nsubsets, 1.0);

    for (auto & v : parameters)
        {
        freqdenom.assign(_nsubsets, 1.0);
        exchgdenom.assign(_nsubsets, 1.0);
        pinvdenom.assign(_nsubsets, 1.0);

        double logj = 0.0;
        for (unsigned i = 0; i < _nparams; ++i)
            {
            double value = v[i];

            // log-jacobian
            logj += value;   // term in log-jacobian

            if (_param_types[i] == ptype_pinvar)
                pinvdenom[_param_subsets[i] - 1] += exp(value);
            else if (_param_types[i] == ptype_exchangeability)
                exchgdenom[_param_subsets[i] - 1] += exp(value);
            else if (_param_types[i] == ptype_frequency)
                freqdenom[_param_subsets[i] - 1] += exp(value);

            // standardize ith parameter vector element
            v[i] = (value - _means[i])/_stdevs[i];
            }

        logj += sum_log_sd;
        for (unsigned k = 0; k < _nsubsets; ++k)
            {
            logj -= 2.0*log(pinvdenom[k]);
            logj -= 4.0*log(freqdenom[k]);
            logj -= 6.0*log(exchgdenom[k]);
            }

        _log_jacobians.push_back(logj);

        }
    }

double Arbor::calcDistance(std::vector<double> & a, std::vector<double> & b) const
    {
    // Returns Euclidean distance betweeen points a and b.
    assert(_nparams == a.size());
    assert(_nparams == b.size());
    double d = 0.0;
    for (unsigned j = 0; j < _nparams; ++j)
        {
        d += pow(a[j] - b[j], 2.0);
        }
    d = sqrt(d);
    return d;
    }

double Arbor::calcPutativeRadius(unsigned index, sample_vect_t & parameters, unsigned num_to_keep, std::vector<pair_uint_dbl_t> & neighbors)
    {
    // Returns smallest radius around _parameters[index] such that num_to_keep reference points are included.
    neighbors.clear();
    for (unsigned i = 0; i < _reference_indices.size(); ++i)
        {
        unsigned k = _reference_indices[i];
        double d = 0.0;
        for (unsigned j = 0; j < _nparams; ++j)
            {
            d += pow(parameters[k][j] - parameters[index][j], 2.0);
            }
        d = sqrt(d);
        neighbors.push_back(std::make_pair(i, d));
        }
    std::sort(neighbors.begin(), neighbors.end(), UintDblPairSmallestFirst);
    neighbors.erase(neighbors.begin() + num_to_keep + 1, neighbors.end());

    return neighbors[num_to_keep].second;
    }

void Arbor::createGrapes(sample_vect_t & parameters, dlb_vect_t & log_likelihoods, dlb_vect_t & log_priors, dlb_vect_t & log_jacobians)
    {
    // Creates grape objects from points stored in _reference_indices.
    assert(_reference_indices.size() > 0);

    // Sort _reference_indices so that reference points having highest log-kernel are first
    std::vector< pair_uint_dbl_t > tmp;
    for (auto j : _reference_indices)
        {
        double log_kernel = log_likelihoods[j] + log_priors[j] + log_jacobians[j];
        tmp.push_back(std::make_pair(j, log_kernel));
        }
    std::sort(tmp.begin(), tmp.end(), UintDblPairLargestFirst);

    // Recreate _reference_indices using tmp
    _reference_indices.clear();
    for (auto x : tmp)
        {
        _reference_indices.push_back(x.first);
        }

    // Initialize _min_radius to the largest possible radius
    _min_radius = std::numeric_limits<double>().max();

    // Each grape is centered at one reference point, and the radius of the grape is set just large enough
    // to enclose the num_to_keep closest reference points
    unsigned num_to_keep = 1;

    _grapes.clear();
    std::vector<pair_uint_dbl_t> neighbors;
    for (unsigned i = 0; i < _reference_indices.size(); ++i)
        {
        unsigned k = _reference_indices[i];
        double r = calcPutativeRadius(k, parameters, num_to_keep, neighbors);
        std::vector<double> & param_vect = parameters[k];
        _grapes.push_back(Grape(k, r, param_vect, neighbors));

        if (r < _min_radius)
            _min_radius = r;
        }

    // Debug sanity check: make sure all grapes contain at least 2 reference sample points (their center plus 1 other)
    for (auto & g : _grapes)
        {
        unsigned ninside = 0;
        for (unsigned i = 0; i < _reference_indices.size(); ++i)
            {
            unsigned k = _reference_indices[i];
            double d = calcDistance(g._v, parameters[k]);
            if (d <= g._radius + 1.e-8)
                ninside++;
            }
        if (ninside != 2)
            throw XArbor(boost::str(boost::format("at least one grape does not enclose 1 point besides itself (ninside = %d)") % ninside));
        }

    // All grapes have same radius equal to the smallest radius calculated for any grape
    for (auto & g : _grapes)
        {
        g._radius = _min_radius;
        }
    }

void Arbor::calcDeltaTerms(sample_vect_t & parameters, dlb_vect_t & log_likelihoods, dlb_vect_t & log_priors, dlb_vect_t & log_jacobians)
    {
    // Calculates hyperball and hypercylinder volumes for each grape and appends terms for current topology to to _delta_terms vector.
    double log_const_term = log(M_PI)*_nparams/2 - boost::math::lgamma<double>(0.5*_nparams + 1.0);
    for (unsigned i = 0; i < _grapes.size(); ++i)
        {
        Grape & g = _grapes[i];
        g._log_hyperball_volume     = log_const_term + log(g._radius)*_nparams;
        g._log_kernel               = log_likelihoods[g._index] + log_priors[g._index] + log_jacobians[g._index];
        g._log_hypercylinder_volume = g._log_hyperball_volume + g._log_kernel;
        _delta_terms.push_back(g._log_hypercylinder_volume);
        }
    //std::sort(_grapes.begin(), _grapes.end(), std::greater<Grape>());
    }

void Arbor::calcRatioTerms(sample_vect_t & parameters, dlb_vect_t & log_likelihoods, dlb_vect_t & log_priors, dlb_vect_t & log_jacobians)
    {
    // Appends terms for current topology to to _log_ratio_terms vector.
    // The tmp_count_num_grapes vector counts the number of grapes to which each estimation sample belongs
    std::vector<unsigned> tmp_count_num_grapes(_estimation_indices.size(), 0);

    _placed.clear();

    // Go through _estimation_indices and, for each point, determine if it falls within any grape. If so, append its log kernel minus
    // the representative (grape's) log kernel to the log_ratios vector.
    unsigned k = 0;
    for (auto i : _estimation_indices)
        {
        for (unsigned j = 0; j < _grapes.size(); ++j)
            {
            Grape & g = _grapes[j];
            double d = g.distPointToGrapeBoundary(parameters[i]);
            if (d < 0.0)
                {
                double log_kernel = log_likelihoods[i] + log_priors[i] + log_jacobians[i];

                g._placed_log_kernels.push_back(log_kernel);            // record log_kernel of point placed in grape g so that variance can later be calculated
                _placed.push_back(fabs(g._log_kernel - log_kernel));
                _log_ratio_terms.push_back(g._log_kernel - log_kernel);
                tmp_count_num_grapes[k] += 1;
                }
            }
        k += 1;
        }

    // Count number of estimation sample points belonging to multiple grapes
    _in_multiple_grapes = 0;
    for (k = 0; k < _estimation_indices.size(); ++k)
        {
        if (tmp_count_num_grapes[k] > 1)
            _in_multiple_grapes++;
        }
    }

void Arbor::margLikeOneTopology(sample_vect_t & parameters, dlb_vect_t & log_likelihoods, dlb_vect_t & log_priors)
    {
    // Builds _reference_indices and _estimation_indices for this tree topology
    // and adds number of estimation points to _N)
    chooseGrapeSeeds(parameters);

    // Calculates _means and _stdevs from reference sample only
    calcMeansAndStdevs(parameters);

    // Calculates _log_jacobians (one element for each element of parameters)
    if (_correlation)
        standardizeSampleAllowingCorrelation(parameters);
    else
        standardizeSampleAssumingIndependence(parameters);

    // Builds _grapes from reference sample only; sets _radius and _v for each grape created.
    createGrapes(parameters, log_likelihoods, log_priors, _log_jacobians);

    // Calculates _log_hyperball_volume, _log_kernel, and _log_hypercylinder_volume for each grape, then
    // sorts _grapes from largest to smallest _log_hypercylinder_volume and adds terms to _delta_terms vector
    calcDeltaTerms(parameters, log_likelihoods, log_priors, _log_jacobians);

    // Processes estimation sample, adding terms to _log_ratio_terms vector and _placed vector
    calcRatioTerms(parameters, log_likelihoods, log_priors, _log_jacobians);
    }

void Arbor::calcIndivTopolMargLikes()
    {
    _N = 0;
    std::vector<double> total_placed;
    _total_sample_size = 0;
    //_min_radius = std::numeric_limits<double>().max();

    _log_ratio_terms.clear();
    _delta_terms.clear();

    for (auto p : _tree_topologies)
        {
        unsigned topology  = p.second;
        unsigned frequency = p.first;
        GrapeData::SharedPtr topol_record = _db.addTopology(topology, frequency);

        _total_sample_size += frequency;
        if (frequency > _min_sample_size)
            {
            margLikeOneTopology(_parameters[topology], _log_likelihoods[topology], _log_priors[topology]);

            topol_record->_used = true;
            topol_record->_reference_samplesize = (unsigned)_reference_indices.size();
            topol_record->_estimation_samplesize = (unsigned)_estimation_indices.size();
            topol_record->_num_placed = (unsigned)_placed.size();
            topol_record->_radius = _min_radius;
            }
        }
    }

void Arbor::showTopoFreq()
    {
    unsigned i = 1;
    std::cout << "\n  Found " << _tree_topologies.size() << " tree topologies:" << std::endl;
    std::cout << boost::str(boost::format("  %10s %10s %10s\n") % "number" % "topology" % "frequency");
    for (auto p : _tree_topologies)
        {
        std::cout << boost::str(boost::format("  %10d %10d %10d\n") % (i+1) % p.second % p.first);
        ++i;
        }
    }

void Arbor::showTopoTable()
    {
    std::cout << "\n  Found " << _tree_topologies.size() << " tree topologies:" << std::endl;
    std::cout << boost::str(boost::format("  %10s %10s %10s %10s %10s %10s %10s\n") % "topology" % "frequency" % "rsample" % "esample" % "placed" % "%placed" % "radius");
    for (auto topol_record : _db._topologies)
        {
        if (topol_record->_used) {
            std::cout << boost::str(boost::format("  %10d %10d %10d %10d %10d %10.3f %10.3f\n") % topol_record->_id % topol_record->_frequency % topol_record->_reference_samplesize % topol_record->_estimation_samplesize % topol_record->_num_placed % topol_record->calcPercentPlaced() % topol_record->_radius);
            }
        else {
            std::cout << boost::str(boost::format("  %10d %10d below cutoff (%d samples required to consider topology)\n") % topol_record->_id % topol_record->_frequency % _min_sample_size);
            }
        }
    }

double Arbor::summationWithFloatingControl(const std::vector<double> & v) const
    {
    // Computes log of the sum (linear scale) of the (log scale) elements in supplied vector v.
    // For example, let v = [log(a), log(b), log(c)] and suppose max({a,b,c}} = c and that you desire log(a + b + c)
    // log(a + b + c) = log(c[a/c + b/c + c/c])
    //                = log(c) + log([a/c + b/c + 1])
    //                = log(c) + log[exp(log(a) - log(c)) + exp(log(b) - log(c)) + exp(log(c) - log(c))]
    assert(v.size() > 0);
    double max_v = *(std::max_element(v.begin(), v.end()));
    double sum_exp_diff = std::accumulate(v.begin(), v.end(), 0.0, plusexp(-max_v));
    return max_v + log(sum_exp_diff);
    }

void Arbor::calcOverallMargLike()
    {
    double log_denominator_eq11 = summationWithFloatingControl(_delta_terms);
    double log_marginal_likelihood = log_denominator_eq11;
    if (_log_ratio_terms.size() > 0)
        {
        double log_numerator_eq11 = summationWithFloatingControl(_log_ratio_terms);
        log_marginal_likelihood = log_denominator_eq11 - (log_numerator_eq11 - log(_N));
        }

    std::cout << boost::str(boost::format("Log marginal likelihood     = %.5f") % log_marginal_likelihood) << std::endl;
    }

void Arbor::run()
    {
    initRandom();
    showUserSpecifiedOptions();
    readData();
    showParamTable();
    showTopoFreq();
    calcIndivTopolMargLikes();
    showTopoTable();
    calcOverallMargLike();
    }

max_stars_repo_name: plewis/arbor

max_issues_repo_name: plewis/arbor

max_forks_repo_name: plewis/arbor


text:
/**
 * @example ocean/test/reflect_loss_test.cc
 */
#include <boost/test/unit_test.hpp>
#include <usml/ocean/ocean.h>
#include <fstream>

BOOST_AUTO_TEST_SUITE(reflect_loss_test)

using namespace boost::unit_test;
using namespace usml::ocean;

/**
 * @ingroup ocean_test
 * @{
 */

/**
 * Test the basic features of the reflection loss model using
 * the reflect_loss_constant model.
 * Generate errors if values differ by more that 1E-6 percent.
 */
BOOST_AUTO_TEST_CASE( constant_reflect_loss_test ) {
    cout << "=== reflect_loss_test: constant_reflect_loss_test ===" << endl ;

    // simple values for points and distance

    wposition1 points ;

    // compute reflect_loss

    seq_log freq( 10.0, 10.0, 7 ) ;
    cout << "freq:  " << freq << endl ;
    vector<double> amplitude( freq.size() ) ;

    double value = 3.0 ;
    reflect_loss_constant model( value ) ;
    model.reflect_loss( points, freq, 0.1, &amplitude ) ;
    cout << "amplitude: " << amplitude << endl ;

    // check the answer

    for ( size_t f=0 ; f < freq.size() ; ++f ) {
        BOOST_CHECK_CLOSE( amplitude(f), value, 1e-6 ) ;
    }
}

/**
 * Compare values of the Rayleigh model to test case (a)
 * from of the Rayleigh model in Jensen, Kuperman, Porter,
 * Schmidt, "Computational Ocean Acoustics",
 * Table 1.4 and Figure 1.22, Chapter 1.6.1, pp 46-47.
 * Write results to CSV file for processing in Excel or Matlab.
 */
BOOST_AUTO_TEST_CASE( rayleigh_test_a ) {
    cout << "=== reflect_loss_test: rayleigh_test_a ===" << endl ;
    const char* name = USML_TEST_DIR "/ocean/test/rayleigh_test_a.csv" ;
    std::ofstream os(name) ;
    cout << "writing tables to " << name << endl ;

    // simple values for points and distance

    wposition1 points ;
    points.altitude(-1000.0) ;

    seq_log freq( 10.0, 10.0, 7 ) ;
    cout << "freq:  " << freq << endl ;
    vector<double> amplitude( freq.size() ) ;

    // test case (a) - variations with sediment sound speed

    static double speed[] = { 1550.0, 1600.0, 1800 } ;
    os << "angle,cp=1550,cp=1600,cp=1800" << endl ;
    for ( int angle = 0 ; angle <= 90 ; ++angle ) {
        os << angle ;
        for ( int n=0 ; n < 3 ; ++n ) {
            reflect_loss_rayleigh model( 2.0, speed[n]/1500.0, 0.5 ) ;
            model.reflect_loss(
                points, freq, to_radians(angle), &amplitude ) ;
            os << "," << amplitude(0) ;
        }
        os << endl ;
    }
}

/**
 * Compare values of the Rayleigh model to test case (b)
 * from of the Rayleigh model in Jensen, Kuperman, Porter,
 * Schmidt, "Computational Ocean Acoustics",
 * Table 1.4 and Figure 1.22, Chapter 1.6.1, pp 46-47.
 * Write results to CSV file for processing in Excel or Matlab.
 */
BOOST_AUTO_TEST_CASE( rayleigh_test_b ) {
    cout << "=== reflect_loss_test: rayleigh_test_b ===" << endl ;
    const char* name = USML_TEST_DIR "/ocean/test/rayleigh_test_b.csv" ;
    std::ofstream os(name) ;
    cout << "writing tables to " << name << endl ;

    // simple values for points and distance

    wposition1 points ;
    points.altitude(-1000.0) ;

    seq_log freq( 10.0, 10.0, 7 ) ;
    cout << "freq:  " << freq << endl ;
    vector<double> amplitude( freq.size() ) ;

    // test case (a) - variations with sediment attenuation

    static double atten[] = { 1.0, 0.5, 0.0 } ;
    os << "angle,a=1.0,a=0.5,a=0.0" << endl ;
    for ( int angle = 0 ; angle <= 90 ; ++angle ) {
        os << angle ;
        for ( int n=0 ; n < 3 ; ++n ) {
            reflect_loss_rayleigh model( 2.0, 1600.0/1500.0, atten[n] ) ;
            model.reflect_loss(
                points, freq, to_radians(angle), &amplitude ) ;
            os << "," << amplitude(0) ;
        }
        os << endl ;
    }
}

/**
 * Compare values of the Rayleigh model to test case (c)
 * from of the Rayleigh model in Jensen, Kuperman, Porter,
 * Schmidt, "Computational Ocean Acoustics",
 * Table 1.4 and Figure 1.22, Chapter 1.6.1, pp 46-47.
 * Write results to CSV file for processing in Excel or Matlab.
 */
BOOST_AUTO_TEST_CASE( rayleigh_test_c ) {
    cout << "=== reflect_loss_test: rayleigh_test_c ===" << endl ;
    const char* name = USML_TEST_DIR "/ocean/test/rayleigh_test_c.csv" ;
    std::ofstream os(name) ;
    cout << "writing tables to " << name << endl ;

    // simple values for points and distance

    wposition1 points ;
    points.altitude(-1000.0) ;

    seq_log freq( 10.0, 10.0, 7 ) ;
    cout << "freq:  " << freq << endl ;
    vector<double> amplitude( freq.size() ) ;

    // test case (c) - variations with density

    static double density[] = { 1.5, 2.0, 2.5 } ;
    os << "angle,d=1.5,d=2.0,d=2.5" << endl ;
    for ( int angle = 0 ; angle <= 90 ; ++angle ) {
        os << angle ;
        for ( int n=0 ; n < 3 ; ++n ) {
            reflect_loss_rayleigh model( density[n], 1600.0/1500.0, 0.5 ) ;
            model.reflect_loss(
                points, freq, to_radians(angle), &amplitude ) ;
            os << "," << amplitude(0) ;
        }
        os << endl ;
    }
}

/**
 * Compare values of the Rayleigh model to test case (b)
 * from of the Rayleigh model in Jensen, Kuperman, Porter,
 * Schmidt, "Computational Ocean Acoustics",
 * Table 1.4 and Figure 1.22, Chapter 1.6.1, pp 46-47.
 * Write results to CSV file for processing in Excel or Matlab.
 *
 * Although Table 1.4 indicates an attenuation of 0.5,
 * you can tell from the cs=0 case that an attenuation of 0.0
 * was used to make Figure 1.22.
 */
BOOST_AUTO_TEST_CASE( rayleigh_test_d ) {
    cout << "=== reflect_loss_test: rayleigh_test_d ===" << endl ;
    const char* name = USML_TEST_DIR "/ocean/test/rayleigh_test_d.csv" ;
    std::ofstream os(name) ;
    cout << "writing tables to " << name << endl ;

    // simple values for points and distance

    wposition1 points ;
    points.altitude(-1000.0) ;

    seq_log freq( 10.0, 10.0, 7 ) ;
    cout << "freq:  " << freq << endl ;
    vector<double> amplitude( freq.size() ) ;

    // test case (d) - variations with shear speed

    static double shear[] = { 0.0, 200.0, 400.0, 600.0 } ;
    os << "angle,cs=0.0,cs=200,cs=400,cs=600" << endl ;
    for ( int angle = 0 ; angle <= 90 ; ++angle ) {
        os << angle ;
        for ( int n=0 ; n < 4 ; ++n ) {
            reflect_loss_rayleigh model(
                2.0, 1600.0/1500.0, 0.0, shear[n]/1500.0 ) ;
            model.reflect_loss(
                points, freq, to_radians(angle), &amplitude ) ;
            os << "," << amplitude(0) ;
        }
        os << endl ;
    }
}

/**
 * Compute Rayleigh model values for generic sediments.
 * Write results to CSV file for processing in Excel or Matlab.
 */
BOOST_AUTO_TEST_CASE( plot_rayleigh_sediments ) {
    cout << "=== reflect_loss_test: plot_rayleigh_sediments ===" << endl ;
    const char* name = USML_TEST_DIR "/ocean/test/rayleigh_sediments.csv" ;
    std::ofstream os(name) ;
    cout << "writing tables to " << name << endl ;

    // simple values for points and distance

    wposition1 points ;
    points.altitude(-1000.0) ;

    seq_log freq( 10.0, 10.0, 7 ) ;
    cout << "freq:  " << freq << endl ;
    vector<double> amplitude( freq.size() ) ;

    // test case (a) - variations with sediment sound speed

    static reflect_loss_rayleigh::bottom_type_enum sediment[] = {
        reflect_loss_rayleigh::CLAY,
        reflect_loss_rayleigh::SILT,
        reflect_loss_rayleigh::SAND,
        reflect_loss_rayleigh::GRAVEL,
        reflect_loss_rayleigh::MORAINE,
        reflect_loss_rayleigh::CHALK,
        reflect_loss_rayleigh::LIMESTONE,
        reflect_loss_rayleigh::BASALT
    } ;

    os << "angle,clay,silt,sand,gravel,moraine,chalk,limestone,basalt" << endl ;
    for ( int angle = 0 ; angle <= 90 ; ++angle ) {
        os << angle ;
        for ( int n=0 ; n < 8 ; ++n ) {
            reflect_loss_rayleigh model( sediment[n] ) ;
            model.reflect_loss(
                points, freq, to_radians(angle), &amplitude ) ;
            os << "," << amplitude(0) ;
        }
        os << endl ;
    }
}

/**
 * Test the basic features of the reflection loss model using
 * the netCDF bottom type file.
 * Generate errors if values differ by more that 1E-5 percent.
 */
BOOST_AUTO_TEST_CASE( reflect_loss_netcdf_test ) {
	cout << " === reflection_loss_test: reflection_loss_netcdf bottom type file === " << endl;
	reflect_loss_netcdf netcdf( USML_DATA_DIR "/bottom_province/sediment_test.nc" ) ;

    seq_linear frequency(1000.0, 1000.0, 0.01) ;
    double angle = M_PI_2 ;
    vector<double> amplitude( frequency.size() ) ;

    double limestone = 3.672875 ;
    double sand = 10.166660 ;
    double tolerance = 4e-4 ;

	/** bottom type numbers in the center of the data field top left, right, bottom left, right */
	netcdf.reflect_loss(wposition1(29.5, -83.4), frequency, angle, &amplitude) ;
	BOOST_CHECK_CLOSE(amplitude(0), limestone, tolerance) ;
	netcdf.reflect_loss(wposition1(30.5, -83.4), frequency, angle, &amplitude) ;
	BOOST_CHECK_CLOSE(amplitude(0), sand, tolerance) ;
	netcdf.reflect_loss(wposition1(29.5, -84.2), frequency, angle, &amplitude) ;
	BOOST_CHECK_CLOSE(amplitude(0), sand, tolerance) ;
	netcdf.reflect_loss(wposition1(30.5, -84.2), frequency, angle, &amplitude) ;
	BOOST_CHECK_CLOSE(amplitude(0), limestone, tolerance) ;
	/** bottom type numbers at the corners of the data field top left, bottom left, top right, bottom right */
	netcdf.reflect_loss(wposition1(26, -80), frequency, angle, &amplitude) ;
	BOOST_CHECK_CLOSE(amplitude(0), sand, tolerance) ;
	netcdf.reflect_loss(wposition1(26, -89), frequency, angle, &amplitude) ;
	BOOST_CHECK_CLOSE(amplitude(0), limestone, tolerance) ;
	netcdf.reflect_loss(wposition1(35, -80), frequency, angle, &amplitude) ;
	BOOST_CHECK_CLOSE(amplitude(0), sand, tolerance) ;
	netcdf.reflect_loss(wposition1(35, -89), frequency, angle, &amplitude) ;
	BOOST_CHECK_CLOSE(amplitude(0), limestone, tolerance) ;
}

/**
 * Test the accuracy Pierson and Moskowitz model for computing wave height
 * from wind speed. Compare to significant wave height plot from
 * http://www.wikiwaves.org/Ocean-Wave_Spectra.
 */
BOOST_AUTO_TEST_CASE( wave_height_pierson_test ) {
    cout << "=== reflect_loss_test: wave_height_pierson_test ===" << endl;

    // display plotting data

    const char* csv_name = USML_TEST_DIR "/ocean/test/wave_height_pierson_test.csv" ;
    std::ofstream os(csv_name) ;
    cout << "writing tables to " << csv_name << endl ;
	os << "wind,Hsig" << endl ;
    for ( double wind=0.0 ; wind <= 25.0 ; wind += 1.0 ) {
    	os << wind << "," << 4*wave_height_pierson(wind) << endl ;
    }

    // check the answer against key points in plot

    BOOST_CHECK_CLOSE(wave_height_pierson(0.0), 0.0, 1e-6);
    BOOST_CHECK_CLOSE(wave_height_pierson(15.0), 5.0/4.0, 5.0 );
    BOOST_CHECK_CLOSE(wave_height_pierson(25.0), 14.0/4.0, 5.0 );
}

/**
 * Compare the results of the reflect_loss_eckart surface reflection loss model
 * to Figures 2-5 in Jones et. al.  Returns error if the values at 10 degrees
 * are different from those that we hand calculated.
 *

 * Figure	Wind Speed	Frequency	Loss at 10 deg
 * ------	----------	---------	--------------
 * 2		4.6 m/s		3150 Hz		1.15 dB
 * 3		4.6 m/s		5000 Hz		2.90 dB
 * 4		14.4 m/s	1000 Hz		11.15 dB
 * 5		14.4 m/s	5000 Hz		278.76 dB
 *

 *
 * @xref Adrian D. Jones, Janice Sendt, Alec J. Duncan, Paul A. Clarke and
 * Amos Maggi, "Modelling the acoustic reflection loss at the rough
 * ocean surface," Proceedings of ACOUSTICS 2009, Australian Acoustical Society,
 * 23-25 November 2009, Adelaide, Australia.
 */
BOOST_AUTO_TEST_CASE( reflect_loss_eckart_test ) {
    cout << "=== reflect_loss_test: reflect_loss_eckart_test ===" << endl;
    const char* name = USML_TEST_DIR "/ocean/test/reflect_loss_eckart_test.csv" ;
    std::ofstream os(name) ;
    cout << "writing tables to " << name << endl ;

    // simple values for points and distance

    wposition1 points ;
    seq_log freqA( 1000.0, 1.0, 1 ) ;
    seq_log freqB( 3150.0, 1.0, 1 ) ;
    seq_log freqC( 5000.0, 1.0, 1 ) ;
    vector<double> amplitude( freqA.size() ) ;

    // variations with wind speed and frequency

    os << "angle,figure2,figure3,figure4,figure5" << endl ;

    reflect_loss_eckart trackQ( 4.6 ) ;
    reflect_loss_eckart trackO( 14.4 ) ;

    for ( double angle = 0.0 ; angle <= 10.0 ; angle += 0.1 ) {
        os << angle ;

        trackQ.reflect_loss( points, freqB, to_radians(angle), &amplitude ) ;
        os << "," << amplitude(0) ;
        if ( angle >= 10.0 ) BOOST_CHECK_CLOSE(amplitude(0), 1.15, 1.0 );

        trackQ.reflect_loss( points, freqC, to_radians(angle), &amplitude ) ;
        os << "," << amplitude(0) ;
        if ( angle >= 10.0 ) BOOST_CHECK_CLOSE(amplitude(0), 2.90, 1.0 );

        trackO.reflect_loss( points, freqA, to_radians(angle), &amplitude ) ;
        os << "," << amplitude(0) ;
        if ( angle >= 10.0 ) BOOST_CHECK_CLOSE(amplitude(0), 11.15, 1.0 );

        trackO.reflect_loss( points, freqC, to_radians(angle), &amplitude ) ;
        os << "," << amplitude(0) << endl ;
        if ( angle >= 10.0 ) BOOST_CHECK_CLOSE(amplitude(0), 278.76, 1.0 );
    }
}

/**
 * Compare the results of the reflect_loss_beckmann surface reflection loss model
 * to Figures 2-5 in Jones et. al.  Returns error if the values at 10 degrees
 * are different from those that we hand calculated.
 *

 * Figure	Wind Speed	Frequency	Loss at 10 deg
 * ------	----------	---------	--------------
 * 2		4.6 m/s		3150 Hz		1.15 dB
 * 3		4.6 m/s		5000 Hz		2.90 dB
 * 4		14.4 m/s	1000 Hz		11.15 dB
 * 5		14.4 m/s	5000 Hz		278.76 dB
 *

 *
 * @xref Adrian D. Jones, Janice Sendt, Alec J. Duncan, Paul A. Clarke and
 * Amos Maggi, "Modelling the acoustic reflection loss at the rough
 * ocean surface," Proceedings of ACOUSTICS 2009, Australian Acoustical Society,
 * 23-25 November 2009, Adelaide, Australia.
 */
BOOST_AUTO_TEST_CASE( reflect_loss_beckmann_test ) {
    cout << "=== reflect_loss_test: reflect_loss_beckmann_test ===" << endl;
    const char* name = USML_TEST_DIR "/ocean/test/reflect_loss_beckmann_test.csv" ;
    std::ofstream os(name) ;
    cout << "writing tables to " << name << endl ;

    // simple values for points and distance

    wposition1 points ;
    seq_log freqA( 1000.0, 1.0, 1 ) ;
    seq_log freqB( 3150.0, 1.0, 1 ) ;
    seq_log freqC( 5000.0, 1.0, 1 ) ;
    vector<double> amplitude( freqA.size() ) ;

    // variations with wind speed and frequency

    os << "angle,figure2,figure3,figure4,figure5" << endl ;

    reflect_loss_beckmann trackQ( 4.6 ) ;
    reflect_loss_beckmann trackO( 14.4 ) ;

    for ( double angle = 0.0 ; angle <= 90.0 ; angle += 1.0 ) {
        os << angle ;

        trackQ.reflect_loss( points, freqB, to_radians(angle), &amplitude ) ;
        os << "," << amplitude(0) ;
//        if ( angle >= 10.0 ) BOOST_CHECK_CLOSE(amplitude(0), 1.15, 1.0 );

        trackQ.reflect_loss( points, freqC, to_radians(angle), &amplitude ) ;
        os << "," << amplitude(0) ;
//        if ( angle >= 10.0 ) BOOST_CHECK_CLOSE(amplitude(0), 2.90, 1.0 );

        trackO.reflect_loss( points, freqA, to_radians(angle), &amplitude ) ;
        os << "," << amplitude(0) ;
//        if ( angle >= 10.0 ) BOOST_CHECK_CLOSE(amplitude(0), 11.15, 1.0 );

        trackO.reflect_loss( points, freqC, to_radians(angle), &amplitude ) ;
        os << "," << amplitude(0) << endl ;
//        if ( angle >= 10.0 ) BOOST_CHECK_CLOSE(amplitude(0), 278.76, 1.0 );
    }
}

/// @}

BOOST_AUTO_TEST_SUITE_END()

ext: cc

max_stars_repo_name: fraclipe/UnderSeaModelingLibrary

max_issues_repo_name: Wolframy-NUDT/UnderSeaModelingLibrary

max_forks_repo_name: Wolframy-NUDT/UnderSeaModelingLibrary


text:
// Vec2.cpp
#include "Vec2.h"

#include <armadillo>

namespace ugps
{
    num_ug Vec2::length() const
    {
        return sqrt(length2());
    }

    num_ug Vec2::length2() const
    {
        return x*x + y*y;
    }

    arma::Col<num_ug> Vec2::col() const
    {
        arma::Col<num_ug> c = {x,y};
        return c;
    }

    arma::Row<num_ug> Vec2::row() const
    {
        arma::Row<num_ug> r = {x,y};
        return r;
    }

    Vec2 Vec2::rotate(const num_ug& angle) const
    {
        return Vec2(
                cos(angle)*x - sin(angle)*y,
                sin(angle)*y + cos(angle)*x
                );
    }

    num_ug Vec2::angle() const
    {
        return atan2(y,x);
    }
}

ext: cpp

max_stars_repo_name: ScrambledEggsOnToast/UniversalGPS

max_issues_repo_name: ScrambledEggsOnToast/UniversalGPS

max_forks_repo_name: ScrambledEggsOnToast/UniversalGPS


text:
#ifndef STRINGRECORD_HPP
#define STRINGRECORD_HPP

#include <boost/archive/text_iarchive.hpp>
#include <boost/archive/text_oarchive.hpp>
#include <boost/serialization/list.hpp>
#include <boost/serialization/map.hpp>
#include <boost/serialization/serialization.hpp>
#include <boost/serialization/string.hpp>

namespace onepass
{
  namespace core
  {
    class StringRecord
    {
     public:
      std::string val_;
      StringRecord();
      StringRecord(std::string val);
      friend class boost::serialization::access;
      void serialize(boost::archive::text_iarchive &iar, const unsigned int);
      void serialize(boost::archive::text_oarchive &oar, const unsigned int);
      friend std::ostream &operator<<(std::ostream &os, const StringRecord sr);
    };

  }  // namespace core
}  // namespace onepass

#endif  // STRINGRECORD_HPP

ext: hpp

max_stars_repo_name: inql/OnePass

max_issues_repo_name: inql/OnePass

max_forks_repo_name: inql/OnePass


text:
/**
 * @file
 * @brief Implementation of IntegrationElement() test for geometry objects
 * @author Anian Ruoss
 * @date   2019-02-11 17:59:17
 * @copyright MIT License
 */

#include "check_integration_element.h"

#include <gtest/gtest.h>

#include <Eigen/LU>

namespace lf::geometry::test_utils {

void checkIntegrationElement(const lf::geometry::Geometry &geom,
                             const Eigen::MatrixXd &eval_points) {
  const size_t num_points = eval_points.cols();
  const size_t dim_local = geom.DimLocal();
  const size_t dim_global = geom.DimGlobal();

  Eigen::MatrixXd jacobians = geom.Jacobian(eval_points);
  Eigen::VectorXd integrationElements = geom.IntegrationElement(eval_points);

  EXPECT_EQ(integrationElements.rows(), num_points)
      << "IntegrationElement has " << integrationElements.rows()
      << " rows instead of " << num_points;
  EXPECT_EQ(integrationElements.cols(), 1)
      << "IntegrationElement has " << integrationElements.cols()
      << " cols instead of " << 1;

  for (int j = 0; j < num_points; ++j) {
    Eigen::MatrixXd jacobian =
        jacobians.block(0, j * dim_local, dim_global, dim_local);

    const double integrationElement = integrationElements(j);
    const double approx_integrationElement =
        std::sqrt((jacobian.transpose() * jacobian).determinant());

    EXPECT_FLOAT_EQ(integrationElement, approx_integrationElement)
        << "IntegrationElement incorrect at point " << eval_points.col(j);
  }
}

}  // namespace lf::geometry::test_utils

max_stars_repo_name: Fytch/lehrfempp

max_issues_repo_name: Fytch/lehrfempp

max_forks_repo_name: Fytch/lehrfempp


text:
#include <iostream>
#include <vector>
#include <queue>
#include <limits>
#include <cmath>

#include <boost/graph/adjacency_list.hpp>
#include <boost/graph/push_relabel_max_flow.hpp>
#include <boost/graph/dijkstra_shortest_paths.hpp>

typedef boost::adjacency_list_traits<boost::vecS, boost::vecS, boost::directedS> traits;
typedef boost::adjacency_list<boost::vecS, boost::vecS, boost::directedS, boost::no_property,
    boost::property<boost::edge_capacity_t, long,
        boost::property<boost::edge_residual_capacity_t, long,
            boost::property<boost::edge_reverse_t, traits::edge_descriptor>>>> graph;

typedef traits::vertex_descriptor vertex_desc;
typedef traits::edge_descriptor edge_desc;
typedef boost::adjacency_list<boost::vecS, boost::vecS, boost::undirectedS,
  boost::no_property, boost::property<boost::edge_weight_t, int> >      weighted_graph;
typedef boost::property_map<weighted_graph, boost::edge_weight_t>::type weight_map;

// Custom edge adder class, highly recommended
class edge_adder {
  graph &G;

 public:
  explicit edge_adder(graph &G) : G(G) {}

  void add_edge(int from, int to, long capacity) {
    auto c_map = boost::get(boost::edge_capacity, G);
    auto r_map = boost::get(boost::edge_reverse, G);
    const auto e = boost::add_edge(from, to, G).first;
    const auto rev_e = boost::add_edge(to, from, G).first;
    c_map[e] = capacity;
    c_map[rev_e] = 0; // reverse edge has no capacity!
    r_map[e] = rev_e;
    r_map[rev_e] = e;
  }
};

using namespace std;

struct edge {
  int u;
  int v;
  int length;
  int width;
};

// Strategy:
// - Find all edges that are part of some shortest path
// - Run a max-flow algorithm on those edges
void solve() {
  int n, m, s, f;
  cin >> n >> m >> s >> f;

  // Dijkstra Graph
  weighted_graph G(n);

  // Read neighbors
  vector<edge> edges;
  edges.reserve(m);
  int u, v, length, width;
  for (int i = 0; i < m; ++i) {
    cin >> u >> v >> width >> length;
    boost::add_edge(u, v, length, G);
    edges.push_back({u, v, length, width});
  }

  // Run dijkstra
  std::vector<int> d(n);
  boost::dijkstra_shortest_paths(G, s, boost::distance_map(boost::make_iterator_property_map( d.begin(), boost::get(boost::vertex_index, G))));

  // Flow graph
  graph gFlow(n);
  edge_adder adder(gFlow);

  // Find edges on a shortest path
  vector<edge> shortestPathEdges;
  for (auto e : edges) {
    int diff = d[e.v] - d[e.u];
    if (e.length == diff) {
      adder.add_edge(e.u, e.v, e.width);
    }
    else if (e.length == -diff) {
      adder.add_edge(e.v, e.u, e.width);
    }
  }

  // Compute flow
  long flow = boost::push_relabel_max_flow(gFlow, s, f);
  cout << flow << endl;
}

int main() {
  ios_base::sync_with_stdio(false);
  int t; cin >> t;
  while (t--)
  {
    solve();
  }
  return 0;
}

max_stars_repo_name: dsparber/algolab

max_issues_repo_name: dsparber/algolab

max_forks_repo_name: dsparber/algolab


text:
// Copyright András Vukics 2006–2020. Distributed under the Boost Software License, Version 1.0. (See accompanying file LICENSE.txt)
/// \briefFile{Definition of cppqedutils::SliceIterator together with its helpers}
#if !BOOST_PP_IS_ITERATING

#ifndef CPPQEDCORE_UTILS_SLICEITERATOR_TCC_INCLUDED
#define CPPQEDCORE_UTILS_SLICEITERATOR_TCC_INCLUDED

#include "SliceIterator.h"

#include <boost/preprocessor/repetition/enum.hpp>

namespace cppqedutils { namespace sliceiterator {

//////////////////////////
//
// Indexer Implementations
//
//////////////////////////

namespace details {

// This class serves as private base for Indexer
template<int RANK, typename V>
class IndexerBase
{
public:
  typedef typename boost::mpl::fold<tmptools::Ordinals<RANK>,
                                    boost::fusion::vector<>,
                                    boost::mpl::push_back<boost::mpl::_1,
                                                          boost::mpl::if_<tmptools::numerical_contains<V,boost::mpl::_2>,blitz::Range,int>
                                                          >
                                    >::type Idx;

private:
  typedef typename VecIdxTiny<RANK,V>::const_iterator CI;

public:
  static const Idx&
  fillIdxValues(const VecIdxTiny<RANK,V>& idx)
  {
    boost::fusion::fold(cache_,idx.begin(),
                        [](CI iter, auto& t) -> CI {
                          using T=std::decay_t<decltype(t)>;
                          if constexpr (std::is_same_v<T,blitz::Range>) {t=blitz::Range::all(); return iter;}
                          else if constexpr (std::is_same_v<T,int>) {t=*iter++; return iter;}
                        });
    return cache_;
  }

protected:
  static Idx cache_;
  // In this way the creation of a default-constructed Helper & Idx objects in fillIdxValues can be avoided

};

// Definition of members cache_ (the above is only a declaration!):

template<int RANK, typename V>
typename IndexerBase<RANK,V>::Idx IndexerBase<RANK,V>::cache_;

////////////////////////////
//
// Transpose Implementations
//
////////////////////////////

namespace namehider {

using namespace boost::mpl;
namespace mpl=boost::mpl;
using namespace tmptools;

template<int RANK, typename V>
struct Algorithm
  : fold<Ordinals<RANK>,
         pair<vector_c<int>,typename boost::mpl::begin<V>::type>,
         pair<push_back<mpl::first<mpl::_1>,
                        if_<numerical_contains<V,mpl::_2>,
                            deref<second<mpl::_1> >,
                            mpl::_2
                            >
                        >,
              if_<numerical_contains<V,mpl::_2>,
                  mpl::next<second<mpl::_1> >,
                  second<mpl::_1>
                  >
              >
         >
{};

} // namehider

template<int RANK, typename V>
using TransposerMeta_t = typename boost::mpl::first<typename namehider::Algorithm<RANK,V>::type>::type;

} // details

} } // cppqedutils::sliceiterator

//////////////////////////////////////////
//
// Brute Force Implementations for Indexer
//
//////////////////////////////////////////

#define BOOST_PP_ITERATION_LIMITS (1,BLITZ_ARRAY_LARGEST_RANK)
#define BOOST_PP_FILENAME_1 "SliceIterator.tcc"

#include BOOST_PP_ITERATE()

#undef BOOST_PP_FILENAME_1
#undef BOOST_PP_ITERATION_LIMITS

#endif // CPPQEDCORE_UTILS_SLICEITERATOR_TCC_INCLUDED

#else  // BOOST_PP_IS_ITERATING

#define rank BOOST_PP_ITERATION()

#define INDEXER_print1(z,m,data) boost::mpl::at_c<TM,m>::type::value
#define INDEXER_print2(z,m,data) boost::fusion::at_c<m>(Base::cache_)

namespace cppqedutils { namespace sliceiterator {

template<template <int> class ARRAY, typename V> struct Transposer<ARRAY,rank,V>
{
  typedef details::TransposerMeta_t<rank,V> TM;

  static auto& _(ARRAY<rank>& array)
  {
    array.transposeSelf(BOOST_PP_ENUM(rank,INDEXER_print1,~) );
    return array;
  }

};

template<template <int> class ARRAY, typename V> struct Indexer<ARRAY,rank,V> : private details::IndexerBase<rank,V>
{
  typedef details::IndexerBase<rank,V> Base;

  static auto& _(const ARRAY<rank>& array, ResArray<ARRAY,V>& resArray, const VecIdxTiny<rank,V>& idx)
  {
    Base::fillIdxValues(idx);
    resArray.reference(subscript(array,BOOST_PP_ENUM(rank,INDEXER_print2,~) ) );
    return resArray;
  }

};

} } // cppqedutils::sliceiterator

#undef INDEXER_print1
#undef INDEXER_print2

#undef rank

#endif // BOOST_PP_IS_ITERATING

ext: tcc

max_stars_repo_name: vukics/cppqed

max_issues_repo_name: vukics/cppqed

max_forks_repo_name: vukics/cppqed


text:
//
// Copyright (c) 2019-2020 Ruben Perez Hidalgo (rubenperez038 at gmail dot com)
//
// Distributed under the Boost Software License, Version 1.0. (See accompanying
// file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
//

#ifndef BOOST_MYSQL_DETAIL_AUXILIAR_CHECK_COMPLETION_TOKEN_HPP
#define BOOST_MYSQL_DETAIL_AUXILIAR_CHECK_COMPLETION_TOKEN_HPP

#include <type_traits>
#include <boost/asio/async_result.hpp>
#include "boost/mysql/error.hpp"

namespace boost {
namespace mysql {
namespace detail {

template <typename HandlerSignature>
struct get_handler_arg;

template <typename T>
struct get_handler_arg<void(error_code, T)>
{
    using type = T;
};

template <>
struct get_handler_arg<void(error_code)>
{
    using type = void;
};

template <typename HandlerType, typename HandlerArg>
constexpr bool is_handler_signature_ok()
{
    if constexpr (std::is_same_v<HandlerArg, void>)
    {
        return std::is_invocable_v<HandlerType, error_code>;
    }
    else
    {
        return std::is_invocable_v<HandlerType, error_code, HandlerArg>;
    }
}

template <typename CompletionToken, typename HandlerSignature>
constexpr void check_completion_token()
{
    using handler_type = BOOST_ASIO_HANDLER_TYPE(CompletionToken, HandlerSignature);
    using handler_arg = typename detail::get_handler_arg<HandlerSignature>::type;
    static_assert(
        is_handler_signature_ok<handler_type, handler_arg>(),
        "Invalid CompletionToken type. Check that CompletionToken fullfills the CompletionToken "
            "requirements or that the callback signature you passed in is correct"
    );
}

} // detail
} // mysql
} // boost

#endif /* INCLUDE_BOOST_MYSQL_DETAIL_AUXILIAR_CHECK_COMPLETION_TOKEN_HPP_ */

max_stars_repo_name: AlexAndDad/mysql-asio

max_issues_repo_name: AlexAndDad/mysql-asio

max_forks_repo_name: AlexAndDad/mysql-asio


text:
#include "http/http.h"
#include <boost/foreach.hpp>
#include <boost/test/unit_test.hpp>
#include "error/error.h"
#include <fstream>
#include <iostream>

using namespace pylon ;
using namespace pylon::http ;
using namespace std;
using namespace boost ;

bool parse_file(parser& p, const char* name)
{/*{{{*/

    response resp;
    ifstream data_file(name);
    buffer_t head_buf   = resp.alloc_buffer();
    size_t   read_len   = data_file.readsome(head_buf.unused_begin(), head_buf.unused_size());
    head_buf.mark_used(read_len);

    result_t h_res    = p.parse_header(head_buf,resp);
    if(!h_res.is_end() )
    {
        std::cout <<  "header: " <<  head_buf.begin() << endl;
        return false ;
    }

    buffer_t cur_buf = head_buf.right(h_res.parsed_length());
    result_t c_res;
    while(true)
    {
        c_res =   resp.is_chunked() ?  p.parse_chunk(cur_buf,resp ) : p.parse_body(cur_buf,resp );
        if (c_res.is_lack())
        {
            buffer_t unparsed    = cur_buf.right(c_res.parsed_length());
            cur_buf              = resp.extend_buffer(unparsed,c_res.lack_length());
            uint read_len        = data_file.readsome(cur_buf.unused_begin(), cur_buf.unused_size());
            cur_buf.mark_used(read_len) ;
            continue ;
        }
        break;
    }
    if(c_res.is_bad())
    {
        std::cout << " parse  failed! " << name << endl ;
        return false;
    }
    return true ;
}/*}}}*/

BOOST_AUTO_TEST_CASE( test_by_file)
{
    parser  p;
    BOOST_CHECK(parse_file(p,"../test_data/http/resp_1.txt"));
    BOOST_CHECK(parse_file(p,"../test_data/http/resp_2.txt"));
//    BOOST_CHECK(parse_file(p,"/home/zuowenjian/devspace/rocket/test_data/http/resp_3.txt"));
}

max_stars_repo_name: xcodecraft/youbar

max_issues_repo_name: xcodecraft/youbar

max_forks_repo_name: xcodecraft/youbar


text:
/**************************************************************************
** Copyright (C) 2016-2019 Toshinobu Hondo, Ph.D.
** Copyright (C) 2016-2019 MS-Cheminformatics LLC, Toin, Mie Japan
*
** Contact: toshi.hondo@qtplatz.com
**
** Reference: http://www.analog.com/media/en/technical-documentation/data-sheets/AD5593R.pdf
**           AD5593R: 8-Channel, 12-Bit, Configurable ADC/DAC with On-Chip Reference, I2C Interface Data Sheet
**           Rev.C (Last Content Update: 05/03/2017)
**************************************************************************/

#pragma once

#include <boost/optional.hpp>
#include <array>
#include <bitset>
#include <string>
#include <ostream>
#include <system_error>
#include <memory>

namespace i2c_linux { class i2c; }

namespace i2c_device {

    namespace ad5593 {

        enum AD5593R_IO_FUNCTION {
            ADC
            , DAC
            , GPIO
            , UNUSED_HIGH
            , UNUSED_LOW
            , UNUSED_TRISTATE
            , UNUSED_PULLDOWN
            , nFunctions
        };

        enum AD5593_GPIO_DIRECTION {
            GPIO_INPUT
            , GPIO_OUTPUT
        };

        constexpr size_t number_of_pins = 8;
        constexpr size_t number_of_regs = 7;

        class AD5593 {
            std::array< std::bitset< 8 >, number_of_regs > bitmaps_;
            std::unique_ptr< i2c_linux::i2c > i2c_;
            bool dirty_;
            double Vref_;
            std::bitset< 16 > pd_;
            std::bitset< 16 > ctrl_;
            std::bitset< 16 > adc_sequence_;

        public:
            ~AD5593();
            AD5593();
            AD5593( const std::string& device = "/dev/i2c-1", int address = 0x10 );
            bool open( const std::string& device = "/dev/i2c-1", int address = 0x10 );

            const std::error_code& error_code() const;
            operator bool () const;
            double Vref() const;

            const std::bitset< 16 >& pd() const { return pd_; }
            const std::bitset< 16 >& ctrl() const { return ctrl_; }

            inline bool is_dirty() const { return dirty_; }

        private:
            bool write( uint8_t addr, uint16_t value ) const;
            bool read( uint8_t addr, uint8_t * data, size_t size ) const;
            boost::optional< uint16_t > read( uint8_t addr ) const;

            void save();
            bool restore();

        public:
            bool dac_set_value( int pin, uint16_t value );

            bool reset_function( int pin, AD5593R_IO_FUNCTION );
            bool set_function( int pin, AD5593R_IO_FUNCTION );
            bool is_function( int pin, AD5593R_IO_FUNCTION ) const;

            bool fetch( bool force = true );
            bool commit();

            bool is_adc_en( int pin ) const;
            bool is_dac_en( int pin ) const;
            uint16_t dac_readback( int pin );
            uint16_t adc_readback( int pin );

            void print_config( std::ostream& ) const;
            void print_values( std::ostream& );
            bool set_adc_sequence( uint16_t sequence );
            std::bitset< 16 > adc_sequence( bool fetch = true ) const;
            const std::bitset< 8 >& adc_enabled() const;
            std::pair< std::string, std::bitset< 8 > > reg_enabled( size_t idx ) const;

            bool read_adc_sequence ( uint16_t * data, size_t size );
            void print_adc_sequence( std::ostream&, const uint16_t * sequence, size_t size, bool has_temp ) const;
            double toVolts( uint32_t, bool is_dac ) const;
            uint16_t fromVolts( double, bool is_dac ) const;
            double temperature( uint32_t ) const;
            bool set_vref( bool );
            bool set_ctrl( uint16_t );
            bool set_pd( uint16_t );
        };
    }
}


text:
#pragma once

#include "detail/function.hpp"
#include "detail/shared_data.hpp"
#include "error.hpp"
#include "log.hpp"

#include <algorithm>
#include <boost/algorithm/string.hpp>
#include <boost/asio.hpp>
#include <cctype>
#include <cstdio>
#include <curl/curl.h>
#include <functional>
#include <map>
#include <memory>
#include <string>
#include <string_view>
#include <utility>

namespace curlio {

class Session;

class Request {
public:
	typedef boost::asio::any_io_executor executor_type;

	Request() = default;
	Request(Request&& move) noexcept;
	~Request() noexcept;

	void set_url(const char* url);
	void append_http_field(const char* field);
	void set_method(const char* method);
	void set_content_length(curl_off_t length);
	bool is_valid() const noexcept { return !_data.expired(); }
	bool is_active() const noexcept { return is_valid() && _owner == nullptr; }
	CURL* native_handle() noexcept;
	template<typename Const_buffer_sequence, typename Token>
	auto async_write_some(const Const_buffer_sequence& buffers, Token&& token);
	/// Getting the executor of an inactive request is undefined behavior.
	boost::asio::any_io_executor get_executor() noexcept { return _data.lock()->executor; }
	Request& operator=(Request&& move) noexcept;

private:
	friend Session;

	/// The owner of a started request is the session which sets this to `nullptr`.
	std::shared_ptr<detail::Shared_data> _owner;
	std::weak_ptr<detail::Shared_data> _data;
	curl_slist* _request_headers = nullptr;
	detail::Function<std::size_t(boost::system::error_code, void*, std::size_t)> _send_handler;

	void _lazy_init();
	void _hook(detail::Shared_data& data) noexcept;
	static std::size_t _send_callback(void* data, std::size_t size, std::size_t count,
	                                  void* self_pointer) noexcept;
};

inline Request::Request(Request&& move) noexcept
{
	std::swap(_owner, move._owner);
	std::swap(_data, move._data);
	std::swap(_request_headers, move._request_headers);
	std::swap(_send_handler, move._send_handler);

	if (_owner != nullptr) {
		_hook(*_owner);
	} else {
		const auto data = _data.lock();
		if (data != nullptr) {
			_hook(*data);
		}
	}
}

inline Request::~Request() noexcept { curl_slist_free_all(_request_headers); }

inline void Request::set_url(const char* url)
{
	_lazy_init();
	if (curl_easy_setopt(native_handle(), CURLOPT_URL, url) != CURLE_OK) {
		throw std::system_error{ Code::bad_url };
	}
}

inline void Request::append_http_field(const char* field)
{
	_lazy_init();
	_request_headers = curl_slist_append(_request_headers, field);
	curl_easy_setopt(native_handle(), CURLOPT_HTTPHEADER, _request_headers);
}

inline void Request::set_method(const char* method)
{
	_lazy_init();
	std::string_view tmp = method;
	if (boost::algorithm::iequals(tmp, "get")) {
		curl_easy_setopt(native_handle(), CURLOPT_HTTPGET, 1);
	} else if (boost::algorithm::iequals(tmp, "post")) {
		curl_easy_setopt(native_handle(), CURLOPT_POST, 1);
	} else {
		curl_easy_setopt(native_handle(), CURLOPT_CUSTOMREQUEST, method);
	}
}

inline void Request::set_content_length(curl_off_t length)
{
	_lazy_init();
	curl_easy_setopt(native_handle(), CURLOPT_POSTFIELDSIZE_LARGE, length);
}

inline CURL* Request::native_handle() noexcept
{
	if (_owner != nullptr) {
		return _owner->handle;
	}
	const auto ptr = _data.lock();
	return ptr == nullptr ? nullptr : ptr->handle;
}

template<typename Const_buffer_sequence, typename Token>
inline auto Request::async_write_some(const Const_buffer_sequence& buffers, Token&& token)
{
	return boost::asio::async_initiate<Token, void(boost::system::error_code, std::size_t)>(
	  [this, buffers](auto&& handler) {
		  const auto ptr = _data.lock();
		  auto executor  = boost::asio::get_associated_executor(handler, ptr->executor);

		  // can immediately finish
		  if (ptr == nullptr || ptr->status & detail::finished) {
			  boost::asio::post(executor, std::bind(std::move(handler), boost::asio::error::eof, 0));
		  } else if (_send_handler) {
			  boost::asio::post(executor, std::bind(std::move(handler), Code::multiple_writes, 0));
		  } else {
			  // set write handler when cURL calls the write callback
			  _send_handler = [this, buffers, handler = std::move(handler),
			                   executor](boost::system::error_code ec, void* data, std::size_t size) mutable {
				  // copy data and finish
				  const std::size_t copied = boost::asio::buffer_copy(boost::asio::buffer(data, size), buffers);
				  boost::asio::post(executor, std::bind(std::move(handler), ec, copied));
				  return copied;
			  };

			  // TODO check for errors
			  ptr->pause_mask &= ~CURLPAUSE_SEND;
			  curl_easy_pause(ptr->handle, ptr->pause_mask);
		  }
	  },
	  token);
}

inline Request& Request::operator=(Request&& move) noexcept
{
	this->~Request();
	new (this) Request{ std::move(move) };
	return *this;
}

inline void Request::_lazy_init()
{
	if (!_data.expired()) {
		return;
	}
	_owner = std::make_shared<detail::Shared_data>();
	_data  = _owner;
	_hook(*_owner);
}

inline void Request::_hook(detail::Shared_data& data) noexcept
{
	curl_easy_setopt(data.handle, CURLOPT_READFUNCTION, &Request::_send_callback);
	curl_easy_setopt(data.handle, CURLOPT_READDATA, this);
	// enable all encodings
	curl_easy_setopt(data.handle, CURLOPT_ACCEPT_ENCODING, "");
	data.request_finisher = [this] {
		if (_send_handler) {
			_send_handler(boost::asio::error::eof, nullptr, 0);
			_send_handler.reset();
		}
	};
}

inline std::size_t Request::_send_callback(void* data, std::size_t size, std::size_t count,
                                           void* self_pointer) noexcept
{
	CURLIO_TRACE("Sending callback with buffer size of " << size * count << " bytes");
	const auto self = static_cast<Request*>(self_pointer);

	if (self->_send_handler) {
		const std::size_t copied = self->_send_handler({}, data, size * count);
		self->_send_handler.reset();
		if (copied > 0) {
			CURLIO_DEBUG("Sending " << copied << " bytes");
			return copied;
		}
	}

	// nobody wants to send -> pause for now
	CURLIO_TRACE("Pausing sending function");
	self->_data.lock()->pause_mask |= CURLPAUSE_SEND;
	return CURL_READFUNC_PAUSE;
}

} // namespace curlio

max_stars_repo_name: terrakuh/curlio

max_issues_repo_name: terrakuh/curlio

max_forks_repo_name: terrakuh/curlio