text:
// Boost.Range library
//
// Copyright Neil Groves 2014. Use, modification and distribution is subject
// to 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)
//
// For more information, see http://www.boost.org/libs/range
//
#ifndef BOOST_RANGE_UNIT_TEST_ADAPTOR_MOCK_RANGE_HPP_INCLUDED
#define BOOST_RANGE_UNIT_TEST_ADAPTOR_MOCK_RANGE_HPP_INCLUDED

#include "mock_iterator.hpp"
#include <boost/range/iterator_range_core.hpp>

namespace boost
{
    namespace range
    {
        namespace unit_test
        {

// Make a non-empty range that models the corresponding range concept.
// This is only useful in unit tests. It is main use is to help test concepts
// assertions are present.
template<typename TraversalTag>
iterator_range<mock_iterator<TraversalTag> >&
    mock_range()
{
    static iterator_range<mock_iterator<TraversalTag> > instance(
                mock_iterator<TraversalTag>(0),
                mock_iterator<TraversalTag>(1));
    return instance;
}

template<typename TraversalTag>
const iterator_range<mock_iterator<TraversalTag> >&
    mock_const_range()
{
    static iterator_range<mock_iterator<TraversalTag> > instance(
                mock_iterator<TraversalTag>(0),
                mock_iterator<TraversalTag>(1));
    return instance;
}

        } // namespace unit_test
    } // namespace range
} // namespace boost

#endif // include guard

ext: hpp


text:
#include "RBGL.hpp"

#include <stdlib.h>

#include <boost/graph/simple_point.hpp>

extern "C"
{

#include <Rdefines.h>
#include <R_ext/Random.h>
#include <Rmath.h>

    using namespace std;
    using namespace boost;

    static void delta_and_tau
     (const Graph_ud& g, vector<int>& v_delta, vector<int>& v_tau)
    {
        Graph_ud::vertex_iterator vi, v_end;
        Graph_ud::adjacency_iterator ui, u_end, wi, w_end;

        int dv = 0, tv = 0;

        v_delta.clear();
        v_tau.clear();

        for ( tie(vi, v_end) = vertices(g); vi != v_end; ++vi )
        {
            // delta(v)
            dv = 0;
            for ( tie(ui, u_end) = adjacent_vertices(*vi, g);
                    ui != u_end; ++ui )
            {
                wi = ui;
                for ( ++wi; wi != u_end; ++wi )
                    if ( edge(*ui, *wi, g).second ) dv++;
            }
            v_delta.push_back(dv);

            // tau(v)
            dv = degree(*vi, g);
            tv = dv * ( dv - 1 ) / 2;
            v_tau.push_back(tv);
        }
    }

    SEXP clusteringCoef(
        SEXP num_verts_in, SEXP num_edges_in, SEXP R_edges_in,
        SEXP weighted, SEXP R_v_weights_in)
    {
        int i;

        int NV = INTEGER(num_verts_in)[0];
        vector<double> v_weight(NV, 1);

        if ( INTEGER(weighted)[0] )
        {
            double* weights = REAL(R_v_weights_in);
            for ( i = 0; i < NV; i++ ) v_weight[i] = weights[i];
        }

        Graph_ud g(num_verts_in, num_edges_in, R_edges_in);
        vector<int> v_delta, v_tau;
        delta_and_tau(g, v_delta, v_tau);

        double nn = 0;		// count nodes w/ deg(v) >= 2
        double cG = 0;
	Graph_ud::vertex_descriptor v;
        for ( i = 0; i < NV; i++ )
        {
	    v = vertex(i, g);
            if ( out_degree(v, g) >= 2 && v_tau[i] > 0 )
            {
                cG += v_weight[i] * v_delta[i] / v_tau[i];
                nn += v_weight[i];
            }
        }

        if ( nn ) cG /= nn;

        SEXP ccoef;
        PROTECT(ccoef = NEW_NUMERIC(1));
        REAL(ccoef)[0] = cG;
        UNPROTECT(1);
        return(ccoef);
    }

    SEXP transitivity( SEXP num_verts_in, SEXP num_edges_in, SEXP R_edges_in)
    {
        int NV = INTEGER(num_verts_in)[0];

        Graph_ud g(num_verts_in, num_edges_in, R_edges_in);
        vector<int> v_delta, v_tau;
        delta_and_tau(g, v_delta, v_tau);

	double tG = 0;
        double sum_dv = 0, sum_tv = 0;
        for ( int i = 0; i < NV; i++ )
        {
            sum_dv += v_delta[i];
            sum_tv += v_tau[i];
        }
        if ( sum_tv ) tG = sum_dv / sum_tv;

#if DEBUG
        cout << " sum_dv = " << sum_dv
        << " sum_tv = " << sum_tv
        << " v_delta.size() = " << v_delta.size()
        << " v_tau.size() = " << v_tau.size()
	<< " tG = " << tG
        << endl;
#endif

        SEXP tcoef;
        PROTECT(tcoef = NEW_NUMERIC(1));
        REAL(tcoef)[0] = tG;
        UNPROTECT(1);
        return(tcoef);
    }

    // uniformly sample in [1, n]
    static inline int uniformRandomNumber(const int n)
    {
	int j = (int)(n * unif_rand()) + 1;   // unif_rand in [0, 1)
	return j;
    }

    static inline int findIndex(const int r, const vector<int>& W)
    {
	unsigned int i;
	for ( i = 1; i < W.size(); i++ ) if ( r <= W[i] ) break;
	return i;
    }

    // input: a node in graph g
    // output: one neighbor of node n chosen uniformly randomly
    static inline void uniformRandomAdjacentNode
	(const Graph_ud::vertex_descriptor& v, const Graph_ud& g,
	 Graph_ud::vertex_descriptor& u,
	 Graph_ud::vertex_descriptor& w)
    {
	int nc = out_degree(v, g);

	Graph_ud::adjacency_iterator vi, v_end;
	tie(vi, v_end) = adjacent_vertices(v, g);

	switch (nc)
	{
	case 0:
	case 1: u = w = *vi;
		break;
	case 2:
		u = *vi; vi++;
		w = *vi;
		break;
	default:
		{
		int r1 = uniformRandomNumber(nc);
		int r2 = uniformRandomNumber(nc);

		while ( r1 == r2 ) r2 = uniformRandomNumber(nc);

		for ( int i = 0; vi != v_end; vi++, i++ )
		{
		    if ( i == r1 ) u = *vi;
		    if ( i == r2 ) w = *vi;
		}

		break;
		}
	}
#if DEBUG
	cout << " uniformRandomAdjacentNode: " << endl;
	cout << " n = " << n << endl;
	cout << " nc = " << nc << endl;
	cout << " *vi = " << *vi << endl;
	cout << " u = " << u << endl;
	cout << " w = " << w << endl;
#endif
    }

    static inline void uniformRandomAdjacentNode_i
	(const int n, const Graph_ud& g,
	 Graph_ud::vertex_descriptor& u,
	 Graph_ud::vertex_descriptor& w)
    {
	Graph_ud::vertex_descriptor v = vertex(n, g);
	uniformRandomAdjacentNode(v, g, u, w);
    }

    // Approximating Cw
    //    Outline of the algorithm:
    //    Input: integer k;
    //           array A[1..|V'|] of nodes V' = {v in V: d(v) >= 2}
    //           node weights w: V' -> N>0;
    //           adjacentcy array for each node
    //    Output: approximation of Cw
    //    Data: node variables: u, w;
    //           integer variables: r, l, j, W[0..|V'|]
    //    Algorithm:
    //    W[0] = 0
    //    for i = (1, ..., |V'|) do
    //       W[i] = W[i-1] + w(A[i])
    //    l = 0
    //    for i in (1, ..., k) do
    //    {
    //       r = UniformRandomNumber( {1,...,W[|V'|]} )
    //       j = FindIndex( j: W[j-1] < r <= W[j] )
    //       u = UniformRandomAdjacentNode(A[j])
    //       repeat
    //         w = UniformRandomAdjacentNode(A[j])
    //       until u != w
    //       if ( EdgeExists(u, w) then
    //          l = l + 1
    //    }
    //    return l/k

    SEXP clusteringCoefAppr(SEXP k_in,
        SEXP num_verts_in, SEXP num_edges_in, SEXP R_edges_in,
        SEXP weighted, SEXP R_v_weights_in)
    {
	// prepare for later unif_rand call
	GetRNGstate();

        Graph_ud g(num_verts_in, num_edges_in, R_edges_in);

        int i, l, r, j;

        int k = INTEGER(k_in)[0];
        int NV = INTEGER(num_verts_in)[0];
        vector<int> v_weight(NV, 1);
        vector<int> W(NV+1, 0);

#if DEBUG
	cout << " inside clusteringCoefAppr "
		<< " k = " << k
		<< " NV = " << NV
		<< endl;
#endif

        if ( INTEGER(weighted)[0] )
        {
            double* weights = REAL(R_v_weights_in);
            for ( i = 0; i < NV; i++ ) v_weight[i] = (int) weights[i];
        }

        Graph_ud::vertex_descriptor u=Graph_ud::null_vertex(), w=Graph_ud::null_vertex();

	W[0] = 0;
	for ( i = 1; i < NV+1; i++ ) W[i] = W[i-1] + v_weight[i-1];

	// TODO: limit nodes to those w/ degree >= 2
	//       pick a number within range uniformaly
	//	 pick an adjacent node uniformaly randomly
	for ( l = 0, i = 0; i < k; i++ )
	{
	   r = uniformRandomNumber(W[NV]);
	   j = findIndex(r, W);
	   uniformRandomAdjacentNode_i(j-1, g, u, w);

	   if ( edge(u, w, g).second ) l++;

#if DEBUG
	   cout << " i = " << i;
	   cout << " r = " << r;
	   cout << " j = " << j;
	   cout << " l = " << l << endl;
#endif
	}

	double cG = double(l) / double(k);

        SEXP ccoef;
        PROTECT(ccoef = NEW_NUMERIC(1));
        REAL(ccoef)[0] = cG;
        UNPROTECT(1);
        return(ccoef);

    }

    inline bool prob_cmp(const simple_point<int>& p1,
			 const simple_point<int>& p2)
	{ return p1.y > p2.y; }

    //  To find a random node w/ probability d(u) / sum(d(V))
    //  The following closely mirrors the codes on
    //  Unequal probability sampling; without-replacement case
    //
    //    /* Record element identities */
    //    for (i = 0; i < n; i++)
    //        perm[i] = i + 1;
    //
    //    //* Sort probabilities into descending order */
    //    //* Order element identities in parallel */
    //    revsort(p, perm, n);
    //
    //    //* Compute the sample */
    //    totalmass = 1;
    //    for (i = 0, n1 = n-1; i < nans; i++, n1--) {
    //        rT = totalmass * unif_rand();
    //        mass = 0;
    //        for (j = 0; j < n1; j++) {
    //            mass += p[j];
    //            if (rT <= mass)
    //                break;
    //        }
    //        ans[i] = perm[j];
    //        totalmass -= p[j];
    //        for(k = j; k < n1; k++) {
    //            p[k] = p[k + 1];
    //            perm[k] = perm[k + 1];
    //        }
    //    }
    static void ProbRandomNode
	(const Graph_ud::vertex_descriptor& v, const Graph_ud& g,
	 Graph_ud::vertex_descriptor& u)
    {

	typedef graph_traits<Graph_ud>::vertex_iterator vertex_iterator;

	vertex_iterator vi, v_end;

	int NV = num_vertices(g);
	std::vector < simple_point<int> > pp(num_vertices(g));

	int i = 0, totalmass = 0;
	for ( tie(vi, v_end) = vertices(g); vi != v_end; vi++, i++ )
	{
	   pp[i].x = i+1;
	   pp[i].y = out_degree(*vi, g);
	   totalmass += pp[i].y;
	}

	std::stable_sort(pp.begin(), pp.end(), prob_cmp);

	int j, k, n1, rT, mass;
	for ( i = 0, n1 = NV-1; i < NV; i++, n1-- )
	{
            rT = (int)(totalmass * unif_rand());
            mass = 0;
            for (j = 0; j < n1; j++) {
                mass += pp[j].y;
                if (rT <= mass) break;
            }
	    u = vertex(i, g);
	    if ( !edge(v, u, g).second ) break;

	    totalmass -= pp[j].y;
	    for ( k = j; k < n1; k++ ) pp[k] = pp[k+1];
	}
    }

    // Graph Generator:
    //    Outline of algorithm:
    //    Input: initial graph G: two connected nodes
    //           integer: n >= 3, d >= 2, o
    //    Output: graph G
    //    Algorithm:
    //    for ( i = (3, ..., n) do
    //    {
    //       v = NewNode()
    //       for 1, ..., Min(i-1, d) do
    //       {
    //          repeat
    //            u = RandomNode( with prob du / sum(d(v) )
    //          until node EdgeExists(v, u)
    //          AddEdge(v, u)
    //       }
    //       for 1, ..., o do
    //       {
    //          u = RandomAdjacentNode(v)
    //          repeat
    //             w = RandomAdjacentNode(v)
    //          until w != u
    //          if ( node EdgeExists(u, w) then
    //             AddEdge(u, w)
    //       }
    //    }

    SEXP graphGenerator(SEXP n_in, SEXP d_in, SEXP o_in)
    {
	int i, j;
        int n = INTEGER(n_in)[0];
        int d = INTEGER(d_in)[0];
        int o = INTEGER(o_in)[0];

	GetRNGstate();	// get random number generator ready

	// initial graph with 2 connected nodes
	Graph_ud g(2);
	boost::add_edge(0, 1, g);

        Graph_ud::vertex_descriptor v, u, w=Graph_ud::null_vertex();

	for ( i = 3; i <= n; i++ )
	{
	   // generate a new node
	   v = boost::add_vertex(g);

	   for ( j = 1; j <= min(i-1, d); j++ )
	   {
		ProbRandomNode(v, g, u);
		boost::add_edge(v, u, g);
	   }

	   for ( j = 1; j <= o; j++ )
	   {
	        uniformRandomAdjacentNode(v, g, u, w);

		if ( !edge(u, w, g).second )
		   boost::add_edge(u, w, g);
	   }
	}

#if DEBUG
	typedef graph_traits<Graph_ud>::vertex_iterator vertex_iterator;
	vertex_iterator vi, v_end;
	cout << " no. of vertices: " << num_vertices(g)
	     << " no. of edges:    " << num_edges(g)
	     << endl;

	for ( tie(vi, v_end) = vertices(g); vi != v_end; vi++, i++ )
	{
	   cout << " vertex: " << *vi
		<< " has degree: " << out_degree(*vi, g)
		<< endl;
	}
#endif

        int NE = num_edges(g);
        SEXP anslst, ncnt, ecnt, enlst;
        PROTECT(anslst = allocVector(VECSXP, 3));
        PROTECT(ncnt = NEW_INTEGER(1));
        PROTECT(ecnt = NEW_INTEGER(1));
        PROTECT(enlst = allocMatrix(INTSXP, 2, NE));

        INTEGER(ncnt)[0] = num_vertices(g);
        INTEGER(ecnt)[0] = NE;

	typedef graph_traits<Graph_ud>::edge_iterator edge_iterator;
	edge_iterator ei, e_end;
	for ( i = 0, tie(ei, e_end) = edges(g); ei != e_end ; ei++ )
        {
            INTEGER(enlst)[i++] = source(*ei, g);
            INTEGER(enlst)[i++] = target(*ei, g);
        }

	SET_VECTOR_ELT(anslst,0,ncnt);
	SET_VECTOR_ELT(anslst,1,ecnt);
        SET_VECTOR_ELT(anslst,2,enlst);
        UNPROTECT(4);
        return(anslst);
    }

}

ext: cpp

max_stars_repo_name: cran/RBGL

max_issues_repo_name: cran/RBGL

max_forks_repo_name: cran/RBGL


text:
//
//  replicates.cpp
//  cufflinks
//
//  Created by Cole Trapnell on 3/11/11.
//  Copyright 2011 Cole Trapnell. All rights reserved.
//

#include <boost/thread.hpp>

extern "C" {
#include "locfit/local.h"
}

#include <boost/random/gamma_distribution.hpp>
#include <boost/random/poisson_distribution.hpp>

#include "replicates.h"
//#include "rounding.h"

#if ENABLE_THREADS
boost::mutex _locfit_lock;
#endif

MassDispersionModel::MassDispersionModel(const std::string& name,
                                         const std::vector<double>& scaled_compatible_mass_means,
                                         const std::vector<double>& scaled_compatible_variances,
                                         const std::vector<double>& scaled_mass_variances)
{
    _name = name;

    if (scaled_compatible_mass_means.size() != scaled_mass_variances.size())
    {
        fprintf (stderr, "Error: dispersion model table is malformed\n");
    }

    double last_val = 0;
    for (size_t i = 0; i < scaled_compatible_mass_means.size(); i++)
    {

        if (last_val > scaled_compatible_mass_means[i])
        {
            fprintf (stderr, "Error: DispersionModel input is malformed\n");
        }

        if ( i == 0 || last_val < scaled_compatible_mass_means[i])
        {
            _scaled_compatible_mass_means.push_back(scaled_compatible_mass_means[i]);
            _scaled_compatible_variances.push_back(scaled_compatible_variances[i]);
            _scaled_mass_variances.push_back(scaled_mass_variances[i]);
        }
        else
        {
            // skip this element if it's equal to what we've already seen
        }

        last_val = scaled_compatible_mass_means[i];
    }
}

double MassDispersionModel::scale_mass_variance(double scaled_mass) const
{
    if (scaled_mass <= 0)
        return 0.0;

    if (_scaled_compatible_mass_means.size() < 2 || _scaled_mass_variances.size() < 2)
    {
        return scaled_mass; // revert to poisson.
    }
    if (scaled_mass > _scaled_compatible_mass_means.back())
    {
        // extrapolate to the right
        // off the right end
        double x1_mean = _scaled_compatible_mass_means[_scaled_compatible_mass_means.size()-2];
        double x2_mean = _scaled_compatible_mass_means[_scaled_compatible_mass_means.size()-1];

        double y1_var = _scaled_mass_variances[_scaled_compatible_mass_means.size()-2];
        double y2_var = _scaled_mass_variances[_scaled_compatible_mass_means.size()-1];
        double slope = 0.0;
        if (x2_mean != x1_mean)
        {
            slope = (y2_var - y1_var) / (x2_mean-x1_mean);
        }
        else if (y1_var == y2_var)
        {
            assert (false); // should have a unique'd table
        }
        double mean_interp = _scaled_mass_variances[_scaled_compatible_mass_means.size()-1] -
        slope*(scaled_mass - _scaled_compatible_mass_means.size()-1);
        if (mean_interp < scaled_mass)
            mean_interp = scaled_mass;
        assert (!isnan(mean_interp) && !isinf(mean_interp));
        return mean_interp;
    }
    else if (scaled_mass < _scaled_compatible_mass_means.front())
    {
        // extrapolate to the left
        // off the left end?
        double x1_mean = _scaled_compatible_mass_means[0];
        double x2_mean = _scaled_compatible_mass_means[1];

        double y1_var = _scaled_mass_variances[0];
        double y2_var = _scaled_mass_variances[1];
        double slope = 0.0;
        if (x2_mean != x1_mean)
        {
            slope = (y2_var - y1_var) / (x2_mean-x1_mean);
        }
        else if (y1_var == y2_var)
        {
            assert (false); // should have a unique'd table
        }
        double mean_interp = _scaled_mass_variances[0] - slope*(_scaled_compatible_mass_means[0] - scaled_mass);
        if (mean_interp < scaled_mass)
            mean_interp = scaled_mass;

        assert (!isnan(mean_interp) && !isinf(mean_interp));
        return mean_interp;
    }

    vector<double>::const_iterator lb;
    lb = lower_bound(_scaled_compatible_mass_means.begin(),
                     _scaled_compatible_mass_means.end(),
                     scaled_mass);
    if (lb < _scaled_compatible_mass_means.end())
    {
        int d = lb - _scaled_compatible_mass_means.begin();
        if (*lb == scaled_mass || lb == _scaled_compatible_mass_means.begin())
        {
            double var = _scaled_mass_variances[d];
            if (var < scaled_mass) // revert to poisson if underdispersed
                var = scaled_mass;
            assert (!isnan(var) && !isinf(var));
            return var;
        }

        //in between two points on the scale.
        d--;

        if (d < 0)
        {
            fprintf(stderr, "ARG d < 0, d = %d \n", d);
        }

        if (d >= _scaled_compatible_mass_means.size())
        {
            fprintf(stderr, "ARG d >= _scaled_compatible_mass_means.size(), d = %d\n", d);
        }
        if (d >= _scaled_mass_variances.size())
        {
            fprintf(stderr, "ARG d >= _scaled_mass_variances.size(), d = %d\n", d);
        }

        double x1_mean = _scaled_compatible_mass_means[d];
        double x2_mean = _scaled_compatible_mass_means[d + 1];

        double y1_var = _scaled_mass_variances[d];
        double y2_var = _scaled_mass_variances[d + 1];
        double slope = 0.0;
        if (x2_mean != x1_mean)
        {
            slope = (y2_var - y1_var) / (x2_mean-x1_mean);
        }
        else if (y1_var == y2_var)
        {
            assert (false); // should have a unique'd table
        }
        double mean_interp = _scaled_mass_variances[d] + slope*(scaled_mass - _scaled_compatible_mass_means[d]);
        if (mean_interp < scaled_mass) // revert to poisson if underdispersed
            mean_interp = scaled_mass;

        assert (!isnan(mean_interp) && !isinf(mean_interp));
        return mean_interp;
    }
    else
    {
        assert (!isnan(scaled_mass) && !isinf(scaled_mass));
        return scaled_mass; // revert to poisson assumption
    }
}

MleErrorModel::MleErrorModel(const std::string& name,
                             const std::vector<double>& scaled_compatible_mass_means,
                             const std::vector<double>& scaled_mle_variances)
{
    _name = name;

    if (scaled_compatible_mass_means.size() != scaled_mle_variances.size())
    {
        fprintf (stderr, "Error: dispersion model table is malformed\n");
    }

    double last_val = 0;
    for (size_t i = 0; i < scaled_compatible_mass_means.size(); i++)
    {

        if (last_val > scaled_compatible_mass_means[i])
        {
            fprintf (stderr, "Error: MLEError input is malformed\n");
        }

        if ( i == 0 || last_val < scaled_compatible_mass_means[i])
        {
            _scaled_compatible_mass_means.push_back(scaled_compatible_mass_means[i]);
            _scaled_mle_variances.push_back(scaled_mle_variances[i]);
        }
        else
        {
            // skip this element if it's equal to what we've already seen
        }

        last_val = scaled_compatible_mass_means[i];
    }
}

double MleErrorModel::scale_mle_variance(double scaled_mass) const
{
    if (scaled_mass <= 0)
        return 0.0;

    if (_scaled_compatible_mass_means.size() < 2 || _scaled_mle_variances.size() < 2)
    {
        return 0; // revert to poisson.
    }
    if (scaled_mass > _scaled_compatible_mass_means.back())
    {
        return 0;
    }
    else if (scaled_mass < _scaled_compatible_mass_means.front())
    {
        return 0; // we won't add anything if we're out of the range of the table
    }

    vector<double>::const_iterator lb;
    lb = lower_bound(_scaled_compatible_mass_means.begin(),
                     _scaled_compatible_mass_means.end(),
                     scaled_mass);
    if (lb < _scaled_compatible_mass_means.end())
    {
        int d = lb - _scaled_compatible_mass_means.begin();
        if (*lb == scaled_mass || lb == _scaled_compatible_mass_means.begin())
        {
            double var = _scaled_mle_variances[d];
            assert (!isnan(var) && !isinf(var));
            if (var < 0)
                var = 0;
            return var;
        }

        //in between two points on the scale.
        d--;

        if (d < 0)
        {
            fprintf(stderr, "ARG d < 0, d = %d \n", d);
        }

        if (d >= _scaled_compatible_mass_means.size())
        {
            fprintf(stderr, "ARG d >= _scaled_compatible_mass_means.size(), d = %d\n", d);
        }
        if (d >= _scaled_mle_variances.size())
        {
            fprintf(stderr, "ARG d >= _scaled_mass_variances.size(), d = %d\n", d);
        }

        double x1_mean = _scaled_compatible_mass_means[d];
        double x2_mean = _scaled_compatible_mass_means[d + 1];

        double y1_var = _scaled_mle_variances[d];
        double y2_var = _scaled_mle_variances[d + 1];
        double slope = 0.0;
        if (x2_mean != x1_mean)
        {
            slope = (y2_var - y1_var) / (x2_mean-x1_mean);
        }
        else if (y1_var == y2_var)
        {
            assert (false); // should have a unique'd table
        }
        double mean_interp = _scaled_mle_variances[d] + slope*(scaled_mass - _scaled_compatible_mass_means[d]);
        if (mean_interp < 0)
            mean_interp = 0;

        assert (!isnan(mean_interp) && !isinf(mean_interp));
        return mean_interp;
    }
    else
    {
        assert (!isnan(scaled_mass) && !isinf(scaled_mass));
        return 0; // revert to poisson assumption
    }
}

void transform_counts_to_common_scale(const vector<double>& scale_factors,
                                      vector<LocusCountList>& sample_compatible_count_table)
{
    // Transform raw counts to the common scale
    for (size_t i = 0; i < sample_compatible_count_table.size(); ++i)
    {
        LocusCountList& p = sample_compatible_count_table[i];
        for (size_t j = 0; j < p.counts.size(); ++j)
        {
            assert (scale_factors.size() > j);
            p.counts[j] *= (1.0 / scale_factors[j]);
        }
    }
}

void calc_geometric_scaling_factors(const vector<LocusCountList>& sample_compatible_count_table,
                                    vector<double>& scale_factors)
{

    vector<double> log_geom_means(sample_compatible_count_table.size(), 0.0);

    for (size_t i = 0; i < sample_compatible_count_table.size(); ++i)
    {
        const LocusCountList& p = sample_compatible_count_table[i];

        for (size_t j = 0; j < p.counts.size(); ++j)
        {
            //assert (log_geom_means.size() > j);
            //if (floor(p.counts[j]) > 0)
            //{
                log_geom_means[i] += (1.0/p.counts.size()) * log(floor(p.counts[j]));
            //}

        }
        //log_geom_means[i] = pow(log_geom_means[i], 1.0/(double)p.counts.size());
    }

    for (size_t j = 0; j < scale_factors.size(); ++j)
    {
        vector<double> tmp_counts;
        for (size_t i = 0; i < sample_compatible_count_table.size(); ++i)
        {
            if (log_geom_means[i] && !isinf(log_geom_means[i]) && !isnan(log_geom_means[i]) && floor(sample_compatible_count_table[i].counts[j]))
            {
                double gm = (double)log(floor(sample_compatible_count_table[i].counts[j])) - log_geom_means[i];
                assert (!isinf(gm));
                tmp_counts.push_back(gm);
            }
        }
        sort(tmp_counts.begin(), tmp_counts.end());
        if (!tmp_counts.empty())
            scale_factors[j] = exp(tmp_counts[tmp_counts.size()/2]);
        else
            scale_factors[j] = 1.0;
    }
}

void calc_classic_fpkm_scaling_factors(const vector<LocusCountList>& sample_compatible_count_table,
                                       vector<double>& scale_factors)
{
    vector<double> total_counts(sample_compatible_count_table.size(), 0.0);

    for (size_t i = 0; i < sample_compatible_count_table.size(); ++i)
    {
        const LocusCountList& p = sample_compatible_count_table[i];

        for (size_t j = 0; j < p.counts.size(); ++j)
        {
            total_counts[j] += floor(p.counts[j]);
        }
    }

    double all_library_counts = accumulate(total_counts.begin(), total_counts.end(), 0.0);
    if (all_library_counts == 0.0)
    {
        for (size_t i = 0; i < scale_factors.size(); ++i)
        {
            scale_factors[i] = 1.0;
        }
        return;
    }
    for (size_t i = 0; i < scale_factors.size(); ++i)
    {
        scale_factors[i] = total_counts[i] / all_library_counts;
    }

    double avg_scaling_factor = accumulate(scale_factors.begin(), scale_factors.end(), 0.0);
    avg_scaling_factor /= scale_factors.size();

    for (size_t i = 0; i < scale_factors.size(); ++i)
    {
        scale_factors[i] = scale_factors[i] / avg_scaling_factor;
    }
}

void calc_estimated_absolute_scaling_factors(vector<boost::shared_ptr<ReadGroupProperties> > & all_read_groups,
                                             vector<double>& scale_factors)
{
    assert (scale_factors.size() == all_read_groups.size());
    for (size_t i = 0; i < scale_factors.size(); ++i)
    {
        double med_cov = all_read_groups[i]->mode_transcript_coverage();
        if (med_cov == 0)
        {
            scale_factors[i] = 0.0;
        }
        else if (med_cov == -1)
        {
            fprintf(stderr, "Error: estimated-absolute requires pre-calculated CXB files\n");
            exit(1);
        }

        scale_factors[i] =  med_cov * 1000000000 / (all_read_groups[i]->normalized_map_mass()) ;
    }

//
//    void adjust_fpkms_by_median_coverage(const vector<vector<double > >& median_coverage,
//                                         FPKMTrackingTable& tracking)
//    {
//        for (FPKMTrackingTable::iterator itr = tracking.begin(); itr != tracking.end(); ++itr)
//        {
//            if (itr != tracking.end())
//            {
//                FPKMTracking& track = itr->second;
//                vector<FPKMContext>& fpkms = track.fpkm_series;
//
//                for (size_t i = 0; i < fpkms.size(); ++i)
//                {
//                    for (size_t j = 0; j != fpkms[i].tracking_info_per_rep.size();
//                         ++j)
//                    {
//                        double& FPKM = fpkms[i].tracking_info_per_rep[j].fpkm;
//                        FPKM /= median_coverage[i][j];
//                        FPKM *= fpkms[i].tracking_info_per_rep[j].rg_props->normalized_map_mass();
//                        FPKM /= 1000000000;
//                    }
//                }
//            }
//        }
//    }
//
//
}

void calc_quartile_scaling_factors(const vector<LocusCountList>& sample_compatible_count_table,
                                   vector<double>& scale_factors)
{

    if (sample_compatible_count_table.empty())
        return;

    vector<double> upper_quartiles(sample_compatible_count_table.front().counts.size(), 0.0);
    vector<double> total_common_masses(sample_compatible_count_table.front().counts.size(), 0.0);

    for (size_t i = 0; i < sample_compatible_count_table.front().counts.size(); ++i)
    {
        vector<double> common_scaled_counts;
        double total_common = 0.0;

        for (size_t j = 0; j < sample_compatible_count_table.size(); ++j)
        {

            //boost::shared_ptr<ReadGroupProperties> rg = bundle_factories[fac_idx];
            //double scaled_mass = scale_factors[fac_idx] * rg->total_map_mass();

            total_common += sample_compatible_count_table[j].counts[i];
            common_scaled_counts.push_back(sample_compatible_count_table[j].counts[i]);
        }

        sort(common_scaled_counts.begin(), common_scaled_counts.end());
        if (common_scaled_counts.empty())
            continue;

        int upper_quart_index = common_scaled_counts.size() * 0.75;
        double upper_quart_count = common_scaled_counts[upper_quart_index];
        upper_quartiles[i] = upper_quart_count;
        total_common_masses[i] = total_common;
    }

    long double total_mass = accumulate(total_common_masses.begin(), total_common_masses.end(), 0.0);
    long double total_norm_mass = accumulate(upper_quartiles.begin(), upper_quartiles.end(), 0.0);

    for (size_t i = 0; i < sample_compatible_count_table.front().counts.size(); ++i)
    {
        if (total_mass > 0)
        {
            double scaling_factor = upper_quartiles[i];
            scaling_factor /= (total_norm_mass / upper_quartiles.size());
            scale_factors[i] = scaling_factor;
        }
        else
        {
            scale_factors[i] = 1.0;
        }
    }
}

void calc_tmm_scaling_factors(const vector<LocusCountList>& sample_compatible_count_table,
                              vector<double>& scale_factors)
{
    scale_factors = vector<double>(sample_compatible_count_table.size(), 1.0);
}

static const int min_loci_for_fitting = 30;

struct SCVInterpolator
{
    void add_scv_pair(double est_scv, double true_scv)
    {
        true_scvs.push_back(true_scv);
        est_scvs.push_back(est_scv);
    }

    void finalize()
    {
        vector<pair<double, double> > pv;
        for (size_t i =0; i < true_scvs.size(); ++i)
        {
            pv.push_back(make_pair(est_scvs[i], true_scvs[i]));
        }
        sort(pv.begin(), pv.end());
    }

    // This was built from the dispersion model interpolator - we should refactor these
    // into a single routine.
    double interpolate_scv(double est_scv)
    {
//        if (est_scv <= 0)
//            return 0.0;

        if (est_scvs.size() < 2 || true_scvs.size() < 2)
        {
            return est_scv; // revert to poisson.
        }
        if (est_scv > est_scvs.back())
        {
            //fprintf(stderr, "Warning: extrapolating to the right\n");
            // extrapolate to the right
            // off the right end
            double x1_mean = est_scvs[est_scvs.size()-2];
            double x2_mean = est_scvs[est_scvs.size()-1];

            double y1_var = true_scvs[est_scvs.size()-2];
            double y2_var = true_scvs[est_scvs.size()-1];
            double slope = 0.0;
            if (x2_mean != x1_mean)
            {
                slope = (y2_var - y1_var) / (x2_mean-x1_mean);
            }
            else if (y1_var == y2_var)
            {
                assert (false); // should have a unique'd table
            }
            double mean_interp = true_scvs[est_scvs.size()-1] -
                slope*(est_scv - est_scvs.size()-1);
//            if (mean_interp < est_scv)
//                mean_interp = est_scv;
            assert (!isnan(mean_interp) && !isinf(mean_interp));
            return mean_interp;
        }
        else if (est_scv < est_scvs.front())
        {
            //fprintf(stderr, "Warning: extrapolating to the left\n");

            // If we're extrapolating to the left, our fit is too coarse, but
            // that probably means we don't need SCV bias correction at all.
            return est_scv;
        }

        vector<double>::const_iterator lb;
        lb = lower_bound(est_scvs.begin(),
                         est_scvs.end(),
                         est_scv);
        if (lb < est_scvs.end())
        {
            int d = lb - est_scvs.begin();
            if (*lb == est_scv || lb == est_scvs.begin())
            {
                double var = true_scvs[d];
//                if (var < est_scv) // revert to poisson if underdispersed
//                    var = est_scv;
                assert (!isnan(var) && !isinf(var));
                return var;
            }

            //in between two points on the scale.
            d--;

            if (d < 0)
            {
                fprintf(stderr, "ARG d < 0, d = %d \n", d);
            }

            if (d >= est_scvs.size())
            {
                fprintf(stderr, "ARG d >= est_scvs.size(), d = %d\n", d);
            }
            if (d >= true_scvs.size())
            {
                fprintf(stderr, "ARG d >= true_scvs.size(), d = %d\n", d);
            }

            double x1_mean = est_scvs[d];
            double x2_mean = est_scvs[d + 1];

            double y1_var = true_scvs[d];
            double y2_var = true_scvs[d + 1];
            double slope = 0.0;
            if (x2_mean != x1_mean)
            {
                slope = (y2_var - y1_var) / (x2_mean-x1_mean);
            }
            else if (y1_var == y2_var)
            {
                fprintf(stderr, "Warning: SCV table does not have unique keys\n!");
                assert (false); // should have a unique'd table
            }
            double mean_interp = true_scvs[d] + slope*(est_scv - est_scvs[d]);
//            if (mean_interp < est_scv) // revert to poisson if underdispersed
//                mean_interp = est_scv;

            assert (!isnan(mean_interp) && !isinf(mean_interp));
            return mean_interp;
        }
        else
        {
            assert (!isnan(est_scv) && !isinf(est_scv));
            return est_scv; // revert to poisson assumption
        }

        return est_scv;
    }

private:
    vector<double> true_scvs;
    vector<double> est_scvs;
};

void build_scv_correction_fit(int nreps, int ngenes, int mean_count, SCVInterpolator& true_to_est_scv_table)
{
    setuplf();

    vector<pair<double, double> > alpha_vs_scv;

    boost::mt19937 rng;
    vector<boost::random::negative_binomial_distribution<int, double> > nb_gens;

    vector<double> alpha_range;
    for(double a = 0.0002; a < 2.0; a += 0.0002)
    {
        alpha_range.push_back(a);
    }

    for (double a = 2; a < 100.0; a += 1)
    {
        alpha_range.push_back(a);
    }

    BOOST_FOREACH (double alpha, alpha_range)
    {
        double k = 1.0/alpha;
        double p = k / (k + mean_count);
        double r = (mean_count * p) / (1-p);

        //boost::random::negative_binomial_distribution<int, double> nb(r, p);

        boost::random::gamma_distribution<double> gamma(r, (1-p)/p);

        vector<double> scvs_for_alpha;
        vector<double> draws;
        for (size_t i = 0; i < ngenes; ++i)
        {
            LocusCountList locus_count("", nreps, 1, vector<string>(), vector<string>());
            for (size_t rep_idx = 0; rep_idx < nreps; ++rep_idx)
            {
                double gamma_draw = gamma(rng);
                if (gamma_draw == 0)
                {
                    locus_count.counts[rep_idx] = 0;
                    draws.push_back(0);
                }
                else
                {
                    boost::random::poisson_distribution<long, double> poisson(gamma_draw);
                    locus_count.counts[rep_idx] = poisson(rng);
                    draws.push_back(locus_count.counts[rep_idx]);
                    //fprintf(stderr, "%lg\t", locus_count.counts[rep_idx]);
                }
            }

            double mean = accumulate(locus_count.counts.begin(), locus_count.counts.end(), 0.0);
            if (mean == 0)
                continue;
            mean /= locus_count.counts.size();
            double var = 0.0;
            BOOST_FOREACH(double c,  locus_count.counts)
            {
                var += (c-mean)*(c-mean);
            }
            var /= locus_count.counts.size();
            var *= locus_count.counts.size() / (locus_count.counts.size() - 1);

            double scv = var / (mean*mean);
            scvs_for_alpha.push_back(scv);
            //fprintf(stderr, " : mean = %lg, var = %lg, scv = %lg\n", mean, var, scv);

            //fprintf(stderr, "\n");
        }

        double mean = accumulate(draws.begin(), draws.end(), 0.0);
        mean /= draws.size();
        double var = 0.0;
        BOOST_FOREACH(int c,  draws)
        {
            var += (c - mean)*(c-mean);
        }
        var /= draws.size();
        var *= draws.size() / (draws.size() - 1);

        //fprintf(stderr, "##########\n");
        //fprintf(stderr, "mean = %lf, var = %lg\n", mean, var);
        if (scvs_for_alpha.size() > 0)
        {
            double mean_scv = accumulate(scvs_for_alpha.begin(),scvs_for_alpha.end(), 0.0);
            mean_scv /= scvs_for_alpha.size();
            //fprintf(stderr, "alpha = %lg scv = %lg\n", alpha, mean_scv);
            alpha_vs_scv.push_back(make_pair(alpha, mean_scv));
        }
    }

    //fprintf(stderr, "$$$$$$$$$\n");

    //sort (alpha_range.begin(), alpha_range.end());

    char namebuf[256];
    sprintf(namebuf, "trueSCV");
    vari* cm = createvar(namebuf,STREGULAR,alpha_vs_scv.size(),VDOUBLE);
    for (size_t i = 0; i < alpha_vs_scv.size(); ++i)
    {
        cm->dpr[i] = alpha_vs_scv[i].first;
    }

    sprintf(namebuf, "estSCV");
    vari* cv = createvar(namebuf,STREGULAR,alpha_vs_scv.size(),VDOUBLE);
    for (size_t i = 0; i < alpha_vs_scv.size(); ++i)
    {
        cv->dpr[i] = alpha_vs_scv[i].second;
    }

    char locfit_cmd[2048];
    sprintf(locfit_cmd, "locfit trueSCV~estSCV");

    locfit_dispatch(locfit_cmd);

    sprintf(namebuf, "domainSCV");
    vari* cd = createvar(namebuf,STREGULAR,alpha_vs_scv.size(),VDOUBLE);
    for (size_t i = 0; i < alpha_vs_scv.size(); ++i)
    {
        cd->dpr[i] = alpha_vs_scv[i].second;
    }

    sprintf(locfit_cmd, "fittedSCV=predict domainSCV");
    locfit_dispatch(locfit_cmd);

    int n = 0;
    sprintf(namebuf, "fittedSCV");
    vari* cp = findvar(namebuf, 1, &n);
    assert(cp != NULL);

    for (size_t i = 0; i < cp->n; ++i)
    {
        //fprintf(stderr, "%lg\t%lg\n",alpha_range[i], cp->dpr[i]);
        true_to_est_scv_table.add_scv_pair(alpha_range[i], cp->dpr[i]);
    }
    true_to_est_scv_table.finalize();
}

void calculate_count_means_and_vars(const vector<LocusCountList>& sample_compatible_count_table,
                                    vector<pair<double, double> >& means_and_vars)
{

    for (size_t i = 0; i < sample_compatible_count_table.size(); ++i)
    {
        const LocusCountList& p = sample_compatible_count_table[i];
        double mean = accumulate(p.counts.begin(), p.counts.end(), 0.0);
        if (mean > 0.0 && p.counts.size() > 0)
            mean /= p.counts.size();

        double var = 0.0;
        double num_non_zero = 0;
        BOOST_FOREACH (double d, p.counts)
        {
            if (d > 0)
                num_non_zero++;
            var += (d - mean) * (d - mean);
        }
        if (var > 0.0 && p.counts.size())
        {
            var /= p.counts.size();
            var *= p.counts.size() / (p.counts.size() - 1);
        }
        means_and_vars.push_back(make_pair(mean, var));
    }
}

boost::shared_ptr<MassDispersionModel>
fit_dispersion_model_helper(const string& condition_name,
                            const vector<double>& scale_factors,
                            const vector<LocusCountList>& sample_compatible_count_table)
{
    vector<pair<double, double> > compatible_means_and_vars;

    SCVInterpolator true_to_est_scv_table;

    int num_samples = sample_compatible_count_table.front().counts.size();
    if (no_scv_correction == false)
    {
        build_scv_correction_fit(num_samples, 10000, 100000, true_to_est_scv_table);
    }

    setuplf();

    calculate_count_means_and_vars(sample_compatible_count_table, compatible_means_and_vars);

    sort(compatible_means_and_vars.begin(), compatible_means_and_vars.end());

    vector<double> compatible_count_means;
    vector<double> raw_variances;

    for(size_t i = 0; i < compatible_means_and_vars.size(); ++i)
    {
        if (compatible_means_and_vars[i].first > 0 && compatible_means_and_vars[i].second > 0.0)
        {
            compatible_count_means.push_back(compatible_means_and_vars[i].first);
            raw_variances.push_back(compatible_means_and_vars[i].second);
        }
    }

    if (compatible_count_means.size() < min_loci_for_fitting)
    {
        boost::shared_ptr<MassDispersionModel> disperser;
        disperser = boost::shared_ptr<MassDispersionModel>(new PoissonDispersionModel(condition_name));

        return disperser;
    }

    vector<double> fitted_values;

    // WARNING: locfit doesn't like undescores - need camel case for
    // variable names

    char namebuf[256];
    sprintf(namebuf, "countMeans");
    vari* cm = createvar(namebuf,STREGULAR,compatible_count_means.size(),VDOUBLE);
    for (size_t i = 0; i < compatible_count_means.size(); ++i)
    {
        cm->dpr[i] = log(compatible_count_means[i]);
    }

    //sprintf(namebuf, "countSCV");
    sprintf(namebuf, "countVariances");
    vari* cv = createvar(namebuf,STREGULAR,raw_variances.size(),VDOUBLE);
    for (size_t i = 0; i < raw_variances.size(); ++i)
    {
        cv->dpr[i] = raw_variances[i];
        //cv->dpr[i] = raw_scvs[i];
    }

    char locfit_cmd[2048];
    //sprintf(locfit_cmd, "locfit countVariances~countMeans family=gamma");
    sprintf(locfit_cmd, "locfit countVariances~countMeans family=gamma");

    locfit_dispatch(locfit_cmd);

    sprintf(locfit_cmd, "fittedVars=predict countMeans");
    locfit_dispatch(locfit_cmd);

    //sprintf(locfit_cmd, "prfit x fhat h nlx");
    //locfit_dispatch(locfit_cmd);

    double xim = 0;
    BOOST_FOREACH(double s, scale_factors)
    {
        if (s)
            xim += 1.0 / s;
    }
    xim /= scale_factors.size();

    int n = 0;
    sprintf(namebuf, "fittedVars");
    vari* cp = findvar(namebuf, 1, &n);
    assert(cp != NULL);
    for (size_t i = 0; i < cp->n; ++i)
    {
//        if (cp->dpr[i] >= 0)
//        {
            double mean = exp(cm->dpr[i]);
            double fitted_scv = (cp->dpr[i] - mean) / (mean * mean);
            double corrected_scv = true_to_est_scv_table.interpolate_scv(fitted_scv);
            double corrected_variance = mean + (corrected_scv * (mean * mean));
            double uncorrected_variance = mean + (fitted_scv * (mean * mean));
            //fitted_values.push_back(mean + (cp->dpr[i] - xim * mean));
            if (no_scv_correction == false && corrected_variance > uncorrected_variance)
                fitted_values.push_back(corrected_variance);
            else if (uncorrected_variance > 0)
                fitted_values.push_back(uncorrected_variance);
            else
                fitted_values.push_back(compatible_count_means[i]);

//        }
//        else
//        {
//            fitted_values.push_back(compatible_count_means[i]);
//        }
    }

    boost::shared_ptr<MassDispersionModel> disperser;
    disperser = boost::shared_ptr<MassDispersionModel>(new MassDispersionModel(condition_name, compatible_count_means, raw_variances, fitted_values));
    if (dispersion_method == POISSON)
        disperser = boost::shared_ptr<MassDispersionModel>(new PoissonDispersionModel(condition_name));

//    for (map<string, pair<double, double> >::iterator itr = labeled_mv_table.begin();
//         itr != labeled_mv_table.end();
//         ++itr)
//    {
//        string label = itr->first;
//        disperser->set_compatible_mean_and_var(itr->first, itr->second);
//    }

    return disperser;
}

boost::shared_ptr<MassDispersionModel>
fit_dispersion_model(const string& condition_name,
                     const vector<double>& scale_factors,
                     const vector<LocusCountList>& sample_compatible_count_table)
{
//
//#if ENABLE_THREADS
//	boost::mutex::scoped_lock lock(_locfit_lock);
//#endif
    for (size_t i = 0; i < sample_compatible_count_table.size(); ++i)
    {
        if (sample_compatible_count_table[i].counts.size() <= 1)
        {
            // only one replicate - no point in fitting variance
            return boost::shared_ptr<MassDispersionModel>(new PoissonDispersionModel(condition_name));
        }
    }
#if ENABLE_THREADS
    _locfit_lock.lock();
#endif

    ProgressBar p_bar("Modeling fragment count overdispersion.",0);

    int max_transcripts = 0;
    BOOST_FOREACH(const LocusCountList& L, sample_compatible_count_table)
    {
        if (L.num_transcripts > max_transcripts)
        {
            max_transcripts = L.num_transcripts;
        }
    }

    boost::shared_ptr<MassDispersionModel>  model = fit_dispersion_model_helper(condition_name, scale_factors, sample_compatible_count_table);

#if ENABLE_THREADS
    _locfit_lock.unlock();
#endif
    return model;
}

void build_norm_table(const vector<LocusCountList>& full_count_table,
                      boost::shared_ptr<const map<string, LibNormStandards> > normalizing_standards,
                      vector<LocusCountList>& norm_table)
{
    // If we're using housekeeping genes or spike-in controls, select the rows we'll be using from the full count table.
    if (normalizing_standards)
    {
        for (size_t i = 0; i < full_count_table.size(); ++i)
        {
            const vector<string>& gene_ids = full_count_table[i].gene_ids;
            const vector<string>& gene_short_names = full_count_table[i].gene_short_names;

            // If the row has an ID that's in the table, take it.
            map<string, LibNormStandards>::const_iterator g_id_itr = normalizing_standards->end();
            map<string, LibNormStandards>::const_iterator g_name_itr = normalizing_standards->end();

            for (size_t j = 0; j < gene_ids.size(); ++j)
            {
                g_id_itr = normalizing_standards->find(gene_ids[j]);
                if (g_id_itr != normalizing_standards->end())
                {
                    break;
                }
            }

            if (g_id_itr != normalizing_standards->end())
            {
                norm_table.push_back(full_count_table[i]);
                continue;
            }

            for (size_t j = 0; j < gene_short_names.size(); ++j)
            {
                g_name_itr = normalizing_standards->find(gene_short_names[j]);
                if (g_name_itr != normalizing_standards->end())
                {
                    break;
                }
            }

            if (g_name_itr != normalizing_standards->end())
            {
                norm_table.push_back(full_count_table[i]);
                continue;
            }

        }
    }
    else // otherwise, just take all rows.
    {
        norm_table = full_count_table;
    }
}

void normalize_counts(vector<boost::shared_ptr<ReadGroupProperties> > & all_read_groups)
{
    vector<LocusCountList> sample_compatible_count_table;
    vector<LocusCountList> sample_total_count_table;

    for (size_t i = 0; i < all_read_groups.size(); ++i)
    {
        boost::shared_ptr<ReadGroupProperties> rg_props = all_read_groups[i];
        const vector<LocusCount>& raw_compatible_counts = rg_props->raw_compatible_counts();
        const vector<LocusCount>& raw_total_counts = rg_props->raw_total_counts();

        for (size_t j = 0; j < raw_compatible_counts.size(); ++j)
        {
            if (sample_compatible_count_table.size() == j)
            {
                const string& locus_id = raw_compatible_counts[j].locus_desc;
                int num_transcripts = raw_compatible_counts[j].num_transcripts;

                const vector<string>& gene_ids = raw_compatible_counts[j].gene_ids;
                const vector<string>& gene_short_names = raw_compatible_counts[j].gene_short_names;

                sample_compatible_count_table.push_back(LocusCountList(locus_id,all_read_groups.size(), num_transcripts, gene_ids, gene_short_names));
                sample_total_count_table.push_back(LocusCountList(locus_id,all_read_groups.size(), num_transcripts, gene_ids, gene_short_names));
            }
            double scaled = raw_compatible_counts[j].count;
            //sample_compatible_count_table[j].counts[i] = scaled * unscaling_factor;
            sample_compatible_count_table[j].counts[i] = floor(scaled);
            sample_total_count_table[j].counts[i] = floor(raw_total_counts[j].count);

            assert(sample_compatible_count_table[j].counts[i] >= 0 && !isinf(sample_compatible_count_table[j].counts[i]));
        }
    }

    vector<double> scale_factors(all_read_groups.size(), 0.0);

    vector<LocusCountList> norm_table;

    if (use_compat_mass)
    {
        build_norm_table(sample_compatible_count_table, lib_norm_standards, norm_table);
    }
    else // use_total_mass
    {
        assert(use_total_mass);
        build_norm_table(sample_total_count_table, lib_norm_standards, norm_table);
    }

    if (lib_norm_method == GEOMETRIC)
    {
        calc_geometric_scaling_factors(norm_table, scale_factors);
    }
    else if (lib_norm_method == CLASSIC_FPKM)
    {
        calc_classic_fpkm_scaling_factors(norm_table, scale_factors);
    }
    else if (lib_norm_method == QUARTILE)
    {
        calc_quartile_scaling_factors(norm_table, scale_factors);
    }
    else if (lib_norm_method == TMM)
    {
        calc_tmm_scaling_factors(norm_table, scale_factors);
    }
    else if (lib_norm_method == ESTIMATED_ABSOLUTE)
    {
        calc_estimated_absolute_scaling_factors(all_read_groups, scale_factors);
    }
    else
    {
        assert (false);
    }

    for (size_t i = 0; i < all_read_groups.size(); ++i)
    {
        boost::shared_ptr<ReadGroupProperties> rg_props = all_read_groups[i];
        rg_props->internal_scale_factor(scale_factors[i]);
    }

    assert(sample_compatible_count_table.size() == sample_total_count_table.size());

    // Transform raw counts to the common scale
    for (size_t i = 0; i < sample_compatible_count_table.size(); ++i)
    {
        LocusCountList& p = sample_compatible_count_table[i];
        for (size_t j = 0; j < p.counts.size(); ++j)
        {
            assert (scale_factors.size() > j);
            p.counts[j] *= (1.0 / scale_factors[j]);
        }

        LocusCountList& t = sample_total_count_table[i];
        for (size_t j = 0; j < t.counts.size(); ++j)
        {
            assert (scale_factors.size() > j);
            t.counts[j] *= (1.0 / scale_factors[j]);
        }
    }

    for (size_t i = 0; i < all_read_groups.size(); ++i)
    {
        boost::shared_ptr<ReadGroupProperties> rg_props = all_read_groups[i];
        vector<LocusCount> scaled_compatible_counts;
        for (size_t j = 0; j < sample_compatible_count_table.size(); ++j)
        {
            string& locus_id = sample_compatible_count_table[j].locus_desc;
            double count = sample_compatible_count_table[j].counts[i];
            int num_transcripts = sample_compatible_count_table[j].num_transcripts;

            const vector<string>& gids = sample_compatible_count_table[j].gene_ids;
            const vector<string>& gnms = sample_compatible_count_table[j].gene_short_names;

            LocusCount locus_count(locus_id, count, num_transcripts, gids, gnms);
            scaled_compatible_counts.push_back(locus_count);
        }
        rg_props->common_scale_compatible_counts(scaled_compatible_counts);
    }

    for (size_t i = 0; i < all_read_groups.size(); ++i)
    {
        boost::shared_ptr<ReadGroupProperties> rg_props = all_read_groups[i];
        vector<LocusCount> scaled_total_counts;
        for (size_t j = 0; j < sample_total_count_table.size(); ++j)
        {
            string& locus_id = sample_total_count_table[j].locus_desc;
            double count = sample_total_count_table[j].counts[i];
            int num_transcripts = sample_total_count_table[j].num_transcripts;

            const vector<string>& gids = sample_total_count_table[j].gene_ids;
            const vector<string>& gnms = sample_total_count_table[j].gene_short_names;

            LocusCount locus_count(locus_id, count, num_transcripts, gids, gnms);
            scaled_total_counts.push_back(locus_count);
        }
        rg_props->common_scale_total_counts(scaled_total_counts);
    }

    double avg_total_common_scaled_count = 0.0;

    for (size_t fac_idx = 0; fac_idx < all_read_groups.size(); ++fac_idx)
    {
        double total_common = 0.0;
        if (use_compat_mass)
        {
            for (size_t j = 0; j < sample_compatible_count_table.size(); ++j)
            {
                total_common += sample_compatible_count_table[j].counts[fac_idx];
            }
        }
        else
        {
            for (size_t j = 0; j < sample_compatible_count_table.size(); ++j)
            {
                total_common += sample_total_count_table[j].counts[fac_idx];
            }

        }

        avg_total_common_scaled_count += (1.0/all_read_groups.size()) * total_common;
    }

    BOOST_FOREACH(boost::shared_ptr<ReadGroupProperties> rg, all_read_groups)
    {
        rg->normalized_map_mass(avg_total_common_scaled_count);
    }
}

max_stars_repo_name: idiv-biodiversity/cufflinks

max_issues_repo_name: idiv-biodiversity/cufflinks

max_forks_repo_name: idiv-biodiversity/cufflinks


text:
#include <eosio/chain/authorization_manager.hpp>
#include <eosio/chain/exceptions.hpp>
#include <eosio/chain/permission_object.hpp>
#include <eosio/chain/permission_link_object.hpp>
#include <eosio/chain/authority_checker.hpp>
#include <eosio/chain/controller.hpp>
#include <eosio/chain/global_property_object.hpp>
#include <eosio/chain/contract_types.hpp>
#include <eosio/chain/generated_transaction_object.hpp>
#include <boost/tuple/tuple_io.hpp>
#include <eosio/chain/database_utils.hpp>
#include <eosio/chain/protocol_state_object.hpp>

namespace eosio { namespace chain {

   using authorization_index_set = index_set<
      permission_index,
      permission_usage_index,
      permission_link_index
   >;

   authorization_manager::authorization_manager(controller& c, database& d)
   :_control(c),_db(d){}

   void authorization_manager::add_indices() {
      authorization_index_set::add_indices(_db);
   }

   void authorization_manager::initialize_database() {
      _db.create<permission_object>([](auto&){}); /// reserve perm 0 (used else where)
   }

   namespace detail {
      template<>
      struct snapshot_row_traits<permission_object> {
         using value_type = permission_object;
         using snapshot_type = snapshot_permission_object;

         static snapshot_permission_object to_snapshot_row(const permission_object& value, const chainbase::database& db) {
            snapshot_permission_object res;
            res.name = value.name;
            res.owner = value.owner;
            res.last_updated = value.last_updated;
            res.auth = value.auth.to_authority();

            // lookup parent name
            const auto& parent = db.get(value.parent);
            res.parent = parent.name;

            // lookup the usage object
            const auto& usage = db.get<permission_usage_object>(value.usage_id);
            res.last_used = usage.last_used;

            return res;
         };

         static void from_snapshot_row(snapshot_permission_object&& row, permission_object& value, chainbase::database& db) {
            value.name = row.name;
            value.owner = row.owner;
            value.last_updated = row.last_updated;
            value.auth = row.auth;

            value.parent = 0;
            if (value.id == 0) {
               EOS_ASSERT(row.parent == permission_name(), snapshot_exception, "Unexpected parent name on reserved permission 0");
               EOS_ASSERT(row.name == permission_name(), snapshot_exception, "Unexpected permission name on reserved permission 0");
               EOS_ASSERT(row.owner == name(), snapshot_exception, "Unexpected owner name on reserved permission 0");
               EOS_ASSERT(row.auth.accounts.size() == 0,  snapshot_exception, "Unexpected auth accounts on reserved permission 0");
               EOS_ASSERT(row.auth.keys.size() == 0,  snapshot_exception, "Unexpected auth keys on reserved permission 0");
               EOS_ASSERT(row.auth.waits.size() == 0,  snapshot_exception, "Unexpected auth waits on reserved permission 0");
               EOS_ASSERT(row.auth.threshold == 0,  snapshot_exception, "Unexpected auth threshold on reserved permission 0");
               EOS_ASSERT(row.last_updated == time_point(),  snapshot_exception, "Unexpected auth last updated on reserved permission 0");
               value.parent = 0;
            } else if ( row.parent != permission_name()){
               const auto& parent = db.get<permission_object, by_owner>(boost::make_tuple(row.owner, row.parent));

               EOS_ASSERT(parent.id != 0, snapshot_exception, "Unexpected mapping to reserved permission 0");
               value.parent = parent.id;
            }

            if (value.id != 0) {
               // create the usage object
               const auto& usage = db.create<permission_usage_object>([&](auto& p) {
                  p.last_used = row.last_used;
               });
               value.usage_id = usage.id;
            } else {
               value.usage_id = 0;
            }
         }
      };
   }

   void authorization_manager::add_to_snapshot( const snapshot_writer_ptr& snapshot ) const {
      authorization_index_set::walk_indices([this, &snapshot]( auto utils ){
         using section_t = typename decltype(utils)::index_t::value_type;

         // skip the permission_usage_index as its inlined with permission_index
         if (std::is_same<section_t, permission_usage_object>::value) {
            return;
         }

         snapshot->write_section<section_t>([this]( auto& section ){
            decltype(utils)::walk(_db, [this, &section]( const auto &row ) {
               section.add_row(row, _db);
            });
         });
      });
   }

   void authorization_manager::read_from_snapshot( const snapshot_reader_ptr& snapshot ) {
      authorization_index_set::walk_indices([this, &snapshot]( auto utils ){
         using section_t = typename decltype(utils)::index_t::value_type;

         // skip the permission_usage_index as its inlined with permission_index
         if (std::is_same<section_t, permission_usage_object>::value) {
            return;
         }

         snapshot->read_section<section_t>([this]( auto& section ) {
            bool more = !section.empty();
            while(more) {
               decltype(utils)::create(_db, [this, &section, &more]( auto &row ) {
                  more = section.read_row(row, _db);
               });
            }
         });
      });
   }

   const permission_object& authorization_manager::create_permission( account_name account,
                                                                      permission_name name,
                                                                      permission_id_type parent,
                                                                      const authority& auth,
                                                                      time_point initial_creation_time
                                                                    )
   {
      for(const key_weight& k: auth.keys)
         EOS_ASSERT(k.key.which() < _db.get<protocol_state_object>().num_supported_key_types, unactivated_key_type,
           "Unactivated key type used when creating permission");

      auto creation_time = initial_creation_time;
      if( creation_time == time_point() ) {
         creation_time = _control.pending_block_time();
      }

      const auto& perm_usage = _db.create<permission_usage_object>([&](auto& p) {
         p.last_used = creation_time;
      });

      const auto& perm = _db.create<permission_object>([&](auto& p) {
         p.usage_id     = perm_usage.id;
         p.parent       = parent;
         p.owner        = account;
         p.name         = name;
         p.last_updated = creation_time;
         p.auth         = auth;
      });
      return perm;
   }

   const permission_object& authorization_manager::create_permission( account_name account,
                                                                      permission_name name,
                                                                      permission_id_type parent,
                                                                      authority&& auth,
                                                                      time_point initial_creation_time
                                                                    )
   {
      for(const key_weight& k: auth.keys)
         EOS_ASSERT(k.key.which() < _db.get<protocol_state_object>().num_supported_key_types, unactivated_key_type,
           "Unactivated key type used when creating permission");

      auto creation_time = initial_creation_time;
      if( creation_time == time_point() ) {
         creation_time = _control.pending_block_time();
      }

      const auto& perm_usage = _db.create<permission_usage_object>([&](auto& p) {
         p.last_used = creation_time;
      });

      const auto& perm = _db.create<permission_object>([&](auto& p) {
         p.usage_id     = perm_usage.id;
         p.parent       = parent;
         p.owner        = account;
         p.name         = name;
         p.last_updated = creation_time;
         p.auth         = std::move(auth);
      });
      return perm;
   }

   void authorization_manager::modify_permission( const permission_object& permission, const authority& auth ) {
      for(const key_weight& k: auth.keys)
         EOS_ASSERT(k.key.which() < _db.get<protocol_state_object>().num_supported_key_types, unactivated_key_type,
           "Unactivated key type used when modifying permission");

      _db.modify( permission, [&](permission_object& po) {
         po.auth = auth;
         po.last_updated = _control.pending_block_time();
      });
   }

   void authorization_manager::remove_permission( const permission_object& permission ) {
      const auto& index = _db.template get_index<permission_index, by_parent>();
      auto range = index.equal_range(permission.id);
      EOS_ASSERT( range.first == range.second, action_validate_exception,
                  "Cannot remove a permission which has children. Remove the children first.");

      _db.get_mutable_index<permission_usage_index>().remove_object( permission.usage_id._id );
      _db.remove( permission );
   }

   void authorization_manager::update_permission_usage( const permission_object& permission ) {
      const auto& puo = _db.get<permission_usage_object, by_id>( permission.usage_id );
      _db.modify( puo, [&](permission_usage_object& p) {
         p.last_used = _control.pending_block_time();
      });
   }

   fc::time_point authorization_manager::get_permission_last_used( const permission_object& permission )const {
      return _db.get<permission_usage_object, by_id>( permission.usage_id ).last_used;
   }

   const permission_object*  authorization_manager::find_permission( const permission_level& level )const
   { try {
      EOS_ASSERT( !level.actor.empty() && !level.permission.empty(), invalid_permission, "Invalid permission" );
      return _db.find<permission_object, by_owner>( boost::make_tuple(level.actor,level.permission) );
   } EOS_RETHROW_EXCEPTIONS( chain::permission_query_exception, "Failed to retrieve permission: ${level}", ("level", level) ) }

   const permission_object&  authorization_manager::get_permission( const permission_level& level )const
   { try {
      EOS_ASSERT( !level.actor.empty() && !level.permission.empty(), invalid_permission, "Invalid permission" );
      return _db.get<permission_object, by_owner>( boost::make_tuple(level.actor,level.permission) );
   } EOS_RETHROW_EXCEPTIONS( chain::permission_query_exception, "Failed to retrieve permission: ${level}", ("level", level) ) }

   optional<permission_name> authorization_manager::lookup_linked_permission( account_name authorizer_account,
                                                                              account_name scope,
                                                                              action_name act_name
                                                                            )const
   {
      try {
         // First look up a specific link for this message act_name
         auto key = boost::make_tuple(authorizer_account, scope, act_name);
         auto link = _db.find<permission_link_object, by_action_name>(key);
         // If no specific link found, check for a contract-wide default
         if (link == nullptr) {
            boost::get<2>(key) = {};
            link = _db.find<permission_link_object, by_action_name>(key);
         }

         // If no specific or default link found, use active permission
         if (link != nullptr) {
            return link->required_permission;
         }
         return optional<permission_name>();

       //  return optional<permission_name>();
      } FC_CAPTURE_AND_RETHROW((authorizer_account)(scope)(act_name))
   }

   optional<permission_name> authorization_manager::lookup_minimum_permission( account_name authorizer_account,
                                                                               account_name scope,
                                                                               action_name act_name
                                                                             )const
   {
      // Special case native actions cannot be linked to a minimum permission, so there is no need to check.
      if( scope == config::system_account_name ) {
          EOS_ASSERT( act_name != updateauth::get_name() &&
                     act_name != deleteauth::get_name() &&
                     act_name != linkauth::get_name() &&
                     act_name != unlinkauth::get_name() &&
                     act_name != canceldelay::get_name(),
                     unlinkable_min_permission_action,
                     "cannot call lookup_minimum_permission on native actions that are not allowed to be linked to minimum permissions" );
      }

      try {
         optional<permission_name> linked_permission = lookup_linked_permission(authorizer_account, scope, act_name);
         if( !linked_permission )
            return config::active_name;

         if( *linked_permission == config::eosio_any_name )
            return optional<permission_name>();

         return linked_permission;
      } FC_CAPTURE_AND_RETHROW((authorizer_account)(scope)(act_name))
   }

   void authorization_manager::check_updateauth_authorization( const updateauth& update,
                                                               const vector<permission_level>& auths
                                                             )const
   {
      EOS_ASSERT( auths.size() == 1, irrelevant_auth_exception,
                  "updateauth action should only have one declared authorization" );
      const auto& auth = auths[0];
      EOS_ASSERT( auth.actor == update.account, irrelevant_auth_exception,
                  "the owner of the affected permission needs to be the actor of the declared authorization" );

      const auto* min_permission = find_permission({update.account, update.permission});
      if( !min_permission ) { // creating a new permission
         min_permission = &get_permission({update.account, update.parent});
      }

      EOS_ASSERT( get_permission(auth).satisfies( *min_permission,
                                                  _db.get_index<permission_index>().indices() ),
                  irrelevant_auth_exception,
                  "updateauth action declares irrelevant authority '${auth}'; minimum authority is ${min}",
                  ("auth", auth)("min", permission_level{update.account, min_permission->name}) );
   }

   void authorization_manager::check_deleteauth_authorization( const deleteauth& del,
                                                               const vector<permission_level>& auths
                                                             )const
   {
      EOS_ASSERT( auths.size() == 1, irrelevant_auth_exception,
                  "deleteauth action should only have one declared authorization" );
      const auto& auth = auths[0];
      EOS_ASSERT( auth.actor == del.account, irrelevant_auth_exception,
                  "the owner of the permission to delete needs to be the actor of the declared authorization" );

      const auto& min_permission = get_permission({del.account, del.permission});

      EOS_ASSERT( get_permission(auth).satisfies( min_permission,
                                                  _db.get_index<permission_index>().indices() ),
                  irrelevant_auth_exception,
                  "updateauth action declares irrelevant authority '${auth}'; minimum authority is ${min}",
                  ("auth", auth)("min", permission_level{min_permission.owner, min_permission.name}) );
   }

   void authorization_manager::check_linkauth_authorization( const linkauth& link,
                                                             const vector<permission_level>& auths
                                                           )const
   {
      EOS_ASSERT( auths.size() == 1, irrelevant_auth_exception,
                  "link action should only have one declared authorization" );
      const auto& auth = auths[0];
      EOS_ASSERT( auth.actor == link.account, irrelevant_auth_exception,
                  "the owner of the linked permission needs to be the actor of the declared authorization" );

      if( link.code == config::system_account_name
            || !_control.is_builtin_activated( builtin_protocol_feature_t::fix_linkauth_restriction ) )
      {
         EOS_ASSERT( link.type != updateauth::get_name(),  action_validate_exception,
                     "Cannot link eosio::updateauth to a minimum permission" );
         EOS_ASSERT( link.type != deleteauth::get_name(),  action_validate_exception,
                     "Cannot link eosio::deleteauth to a minimum permission" );
         EOS_ASSERT( link.type != linkauth::get_name(),    action_validate_exception,
                     "Cannot link eosio::linkauth to a minimum permission" );
         EOS_ASSERT( link.type != unlinkauth::get_name(),  action_validate_exception,
                     "Cannot link eosio::unlinkauth to a minimum permission" );
         EOS_ASSERT( link.type != canceldelay::get_name(), action_validate_exception,
                     "Cannot link eosio::canceldelay to a minimum permission" );
      }

      const auto linked_permission_name = lookup_minimum_permission(link.account, link.code, link.type);

      if( !linked_permission_name ) // if action is linked to amax.any permission
         return;

      EOS_ASSERT( get_permission(auth).satisfies( get_permission({link.account, *linked_permission_name}),
                                                  _db.get_index<permission_index>().indices()              ),
                  irrelevant_auth_exception,
                  "link action declares irrelevant authority '${auth}'; minimum authority is ${min}",
                  ("auth", auth)("min", permission_level{link.account, *linked_permission_name}) );
   }

   void authorization_manager::check_unlinkauth_authorization( const unlinkauth& unlink,
                                                               const vector<permission_level>& auths
                                                             )const
   {
      EOS_ASSERT( auths.size() == 1, irrelevant_auth_exception,
                  "unlink action should only have one declared authorization" );
      const auto& auth = auths[0];
      EOS_ASSERT( auth.actor == unlink.account, irrelevant_auth_exception,
                  "the owner of the linked permission needs to be the actor of the declared authorization" );

      const auto unlinked_permission_name = lookup_linked_permission(unlink.account, unlink.code, unlink.type);
      EOS_ASSERT( unlinked_permission_name.valid(), transaction_exception,
                  "cannot unlink non-existent permission link of account '${account}' for actions matching '${code}::${action}'",
                  ("account", unlink.account)("code", unlink.code)("action", unlink.type) );

      if( *unlinked_permission_name == config::eosio_any_name )
         return;

      EOS_ASSERT( get_permission(auth).satisfies( get_permission({unlink.account, *unlinked_permission_name}),
                                                  _db.get_index<permission_index>().indices()                  ),
                  irrelevant_auth_exception,
                  "unlink action declares irrelevant authority '${auth}'; minimum authority is ${min}",
                  ("auth", auth)("min", permission_level{unlink.account, *unlinked_permission_name}) );
   }

   fc::microseconds authorization_manager::check_canceldelay_authorization( const canceldelay& cancel,
                                                                            const vector<permission_level>& auths
                                                                          )const
   {
      EOS_ASSERT( auths.size() == 1, irrelevant_auth_exception,
                  "canceldelay action should only have one declared authorization" );
      const auto& auth = auths[0];

      EOS_ASSERT( get_permission(auth).satisfies( get_permission(cancel.canceling_auth),
                                                  _db.get_index<permission_index>().indices() ),
                  irrelevant_auth_exception,
                  "canceldelay action declares irrelevant authority '${auth}'; specified authority to satisfy is ${min}",
                  ("auth", auth)("min", cancel.canceling_auth) );

      const auto& trx_id = cancel.trx_id;

      const auto& generated_transaction_idx = _control.db().get_index<generated_transaction_multi_index>();
      const auto& generated_index = generated_transaction_idx.indices().get<by_trx_id>();
      const auto& itr = generated_index.lower_bound(trx_id);
      EOS_ASSERT( itr != generated_index.end() && itr->sender == account_name() && itr->trx_id == trx_id,
                  tx_not_found,
                 "cannot cancel trx_id=${tid}, there is no deferred transaction with that transaction id",
                 ("tid", trx_id) );

      auto trx = fc::raw::unpack<transaction>(itr->packed_trx.data(), itr->packed_trx.size());
      bool found = false;
      for( const auto& act : trx.actions ) {
         for( const auto& auth : act.authorization ) {
            if( auth == cancel.canceling_auth ) {
               found = true;
               break;
            }
         }
         if( found ) break;
      }

      EOS_ASSERT( found, action_validate_exception,
                  "canceling_auth in canceldelay action was not found as authorization in the original delayed transaction" );

      return (itr->delay_until - itr->published);
   }

   void noop_checktime() {}

   std::function<void()> authorization_manager::_noop_checktime{&noop_checktime};

   void
   authorization_manager::check_authorization( const vector<action>&                actions,
                                               const flat_set<public_key_type>&     provided_keys,
                                               const flat_set<permission_level>&    provided_permissions,
                                               fc::microseconds                     provided_delay,
                                               const std::function<void()>&         _checktime,
                                               bool                                 allow_unused_keys,
                                               const flat_set<permission_level>&    satisfied_authorizations
                                             )const
   {
      const auto& checktime = ( static_cast<bool>(_checktime) ? _checktime : _noop_checktime );

      auto delay_max_limit = fc::seconds( _control.get_global_properties().configuration.max_transaction_delay );

      auto effective_provided_delay =  (provided_delay >= delay_max_limit) ? fc::microseconds::maximum() : provided_delay;

      auto checker = make_auth_checker( [&](const permission_level& p){ return get_permission(p).auth; },
                                        _control.get_global_properties().configuration.max_authority_depth,
                                        provided_keys,
                                        provided_permissions,
                                        effective_provided_delay,
                                        checktime
                                      );

      map<permission_level, fc::microseconds> permissions_to_satisfy;

      for( const auto& act : actions ) {
         bool special_case = false;
         fc::microseconds delay = effective_provided_delay;

         if( act.account == config::system_account_name ) {
            special_case = true;

            if( act.name == updateauth::get_name() ) {
               check_updateauth_authorization( act.data_as<updateauth>(), act.authorization );
            } else if( act.name == deleteauth::get_name() ) {
               check_deleteauth_authorization( act.data_as<deleteauth>(), act.authorization );
            } else if( act.name == linkauth::get_name() ) {
               check_linkauth_authorization( act.data_as<linkauth>(), act.authorization );
            } else if( act.name == unlinkauth::get_name() ) {
               check_unlinkauth_authorization( act.data_as<unlinkauth>(), act.authorization );
            } else if( act.name ==  canceldelay::get_name() ) {
               delay = std::max( delay, check_canceldelay_authorization(act.data_as<canceldelay>(), act.authorization) );
            } else {
               special_case = false;
            }
         }

         for( const auto& declared_auth : act.authorization ) {

            checktime();

            if( !special_case ) {
               auto min_permission_name = lookup_minimum_permission(declared_auth.actor, act.account, act.name);
               if( min_permission_name ) { // since special cases were already handled, it should only be false if the permission is amax.any
                  const auto& min_permission = get_permission({declared_auth.actor, *min_permission_name});
                  EOS_ASSERT( get_permission(declared_auth).satisfies( min_permission,
                                                                       _db.get_index<permission_index>().indices() ),
                              irrelevant_auth_exception,
                              "action declares irrelevant authority '${auth}'; minimum authority is ${min}",
                              ("auth", declared_auth)("min", permission_level{min_permission.owner, min_permission.name}) );
               }
            }

            if( satisfied_authorizations.find( declared_auth ) == satisfied_authorizations.end() ) {
               auto res = permissions_to_satisfy.emplace( declared_auth, delay );
               if( !res.second && res.first->second > delay) { // if the declared_auth was already in the map and with a higher delay
                  res.first->second = delay;
               }
            }
         }
      }

      // Now verify that all the declared authorizations are satisfied:

      // Although this can be made parallel (especially for input transactions) with the optimistic assumption that the
      // CPU limit is not reached, because of the CPU limit the protocol must officially specify a sequential algorithm
      // for checking the set of declared authorizations.
      // The permission_levels are traversed in ascending order, which is:
      // ascending order of the actor name with ties broken by ascending order of the permission name.
      for( const auto& p : permissions_to_satisfy ) {
         checktime(); // TODO: this should eventually move into authority_checker instead
         EOS_ASSERT( checker.satisfied( p.first, p.second ), unsatisfied_authorization,
                     "transaction declares authority '${auth}', "
                     "but does not have signatures for it under a provided delay of ${provided_delay} ms, "
                     "provided permissions ${provided_permissions}, provided keys ${provided_keys}, "
                     "and a delay max limit of ${delay_max_limit_ms} ms",
                     ("auth", p.first)
                     ("provided_delay", provided_delay.count()/1000)
                     ("provided_permissions", provided_permissions)
                     ("provided_keys", provided_keys)
                     ("delay_max_limit_ms", delay_max_limit.count()/1000)
                   );

      }

      if( !allow_unused_keys ) {
         EOS_ASSERT( checker.all_keys_used(), tx_irrelevant_sig,
                     "transaction bears irrelevant signatures from these keys: ${keys}",
                     ("keys", checker.unused_keys()) );
      }
   }

   void
   authorization_manager::check_authorization( account_name                         account,
                                               permission_name                      permission,
                                               const flat_set<public_key_type>&     provided_keys,
                                               const flat_set<permission_level>&    provided_permissions,
                                               fc::microseconds                     provided_delay,
                                               const std::function<void()>&         _checktime,
                                               bool                                 allow_unused_keys
                                             )const
   {
      const auto& checktime = ( static_cast<bool>(_checktime) ? _checktime : _noop_checktime );

      auto delay_max_limit = fc::seconds( _control.get_global_properties().configuration.max_transaction_delay );

      auto checker = make_auth_checker( [&](const permission_level& p){ return get_permission(p).auth; },
                                        _control.get_global_properties().configuration.max_authority_depth,
                                        provided_keys,
                                        provided_permissions,
                                        ( provided_delay >= delay_max_limit ) ? fc::microseconds::maximum() : provided_delay,
                                        checktime
                                      );

      EOS_ASSERT( checker.satisfied({account, permission}), unsatisfied_authorization,
                  "permission '${auth}' was not satisfied under a provided delay of ${provided_delay} ms, "
                  "provided permissions ${provided_permissions}, provided keys ${provided_keys}, "
                  "and a delay max limit of ${delay_max_limit_ms} ms",
                  ("auth", permission_level{account, permission})
                  ("provided_delay", provided_delay.count()/1000)
                  ("provided_permissions", provided_permissions)
                  ("provided_keys", provided_keys)
                  ("delay_max_limit_ms", delay_max_limit.count()/1000)
                );

      if( !allow_unused_keys ) {
         EOS_ASSERT( checker.all_keys_used(), tx_irrelevant_sig,
                     "irrelevant keys provided: ${keys}",
                     ("keys", checker.unused_keys()) );
      }
   }

   flat_set<public_key_type> authorization_manager::get_required_keys( const transaction& trx,
                                                                       const flat_set<public_key_type>& candidate_keys,
                                                                       fc::microseconds provided_delay
                                                                     )const
   {
      auto checker = make_auth_checker( [&](const permission_level& p){ return get_permission(p).auth; },
                                        _control.get_global_properties().configuration.max_authority_depth,
                                        candidate_keys,
                                        {},
                                        provided_delay,
                                        _noop_checktime
                                      );

      for (const auto& act : trx.actions ) {
         for (const auto& declared_auth : act.authorization) {
            EOS_ASSERT( checker.satisfied(declared_auth), unsatisfied_authorization,
                        "transaction declares authority '${auth}', but does not have signatures for it.",
                        ("auth", declared_auth) );
         }
      }

      return checker.used_keys();
   }

} } /// namespace eosio::chain


text:
/*********************************************************************
 * Software License Agreement (BSD License)
 *
 *  Copyright (c) 2013, PAL Robotics, S.L.
 *  All rights reserved.
 *
 *  Redistribution and use in source and binary forms, with or without
 *  modification, are permitted provided that the following conditions
 *  are met:
 *
 *   * Redistributions of source code must retain the above copyright
 *     notice, this list of conditions and the following disclaimer.
 *   * Redistributions in binary form must reproduce the above
 *     copyright notice, this list of conditions and the following
 *     disclaimer in the documentation and/or other materials provided
 *     with the distribution.
 *   * Neither the name of the PAL Robotics nor the names of its
 *     contributors may be used to endorse or promote products derived
 *     from this software without specific prior written permission.
 *
 *  THIS SOFTWARE IS PROVIDED BY THE 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 OWNER 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.
 *********************************************************************/

/*
* Author: Luca Marchionni
* Author: Bence Magyar
* Author: Enrique Fernández
* Author: Paul Mathieu
*
* Author: Braden McDorman (Semio AI, Inc.)
*/

#include <quori_holonomic_drive_controller/odometry.hpp>

#include <boost/bind.hpp>

using namespace quori_holonomic_drive_controller;

namespace bacc = boost::accumulators;

Odometry::Odometry(size_t velocity_rolling_window_size)
: timestamp_(0.0)
, x_(0.0)
, y_(0.0)
, heading_(0.0)
, linear_(0.0)
, angular_(0.0)
, wheel_separation_(0.0)
, left_wheel_radius_(0.0)
, right_wheel_radius_(0.0)
, left_wheel_old_pos_(0.0)
, right_wheel_old_pos_(0.0)
, velocity_rolling_window_size_(velocity_rolling_window_size)
, linear_acc_(RollingWindow::window_size = velocity_rolling_window_size)
, angular_acc_(RollingWindow::window_size = velocity_rolling_window_size)
, integrate_fun_(boost::bind(&Odometry::integrateExact, this, _1, _2))
{
}

void Odometry::init(const ros::Time& time)
{
  // Reset accumulators and timestamp:
  resetAccumulators();
  timestamp_ = time;
}

bool Odometry::update(double left_pos, double right_pos, const ros::Time &time)
{
  /// Get current wheel joint positions:
  const double left_wheel_cur_pos  = left_pos  * left_wheel_radius_;
  const double right_wheel_cur_pos = right_pos * right_wheel_radius_;

  /// Estimate velocity of wheels using old and current position:
  const double left_wheel_est_vel  = left_wheel_cur_pos  - left_wheel_old_pos_;
  const double right_wheel_est_vel = right_wheel_cur_pos - right_wheel_old_pos_;

  /// Update old position with current:
  left_wheel_old_pos_  = left_wheel_cur_pos;
  right_wheel_old_pos_ = right_wheel_cur_pos;

  /// Compute linear and angular diff:
  const double linear  = (right_wheel_est_vel + left_wheel_est_vel) * 0.5 ;
  const double angular = (right_wheel_est_vel - left_wheel_est_vel) / wheel_separation_;

  /// Integrate odometry:
  integrate_fun_(linear, angular);

  /// We cannot estimate the speed with very small time intervals:
  const double dt = (time - timestamp_).toSec();
  if (dt < 0.0001)
    return false; // Interval too small to integrate with

  timestamp_ = time;

  /// Estimate speeds using a rolling mean to filter them out:
  linear_acc_(linear/dt);
  angular_acc_(angular/dt);

  linear_ = bacc::rolling_mean(linear_acc_);
  angular_ = bacc::rolling_mean(angular_acc_);

  return true;
}

void Odometry::updateOpenLoop(double linear, double angular, const ros::Time &time)
{
  /// Save last linear and angular velocity:
  linear_ = linear;
  angular_ = angular;

  /// Integrate odometry:
  const double dt = (time - timestamp_).toSec();
  timestamp_ = time;
  integrate_fun_(linear * dt, angular * dt);
}

void Odometry::setWheelParams(double wheel_separation, double left_wheel_radius, double right_wheel_radius)
{
  wheel_separation_   = wheel_separation;
  left_wheel_radius_  = left_wheel_radius;
  right_wheel_radius_ = right_wheel_radius;
}

void Odometry::setVelocityRollingWindowSize(size_t velocity_rolling_window_size)
{
  velocity_rolling_window_size_ = velocity_rolling_window_size;

  resetAccumulators();
}

void Odometry::integrateRungeKutta2(double linear, double angular)
{
  const double direction = heading_ + angular * 0.5;

  /// Runge-Kutta 2nd order integration:
  x_       += linear * cos(direction);
  y_       += linear * sin(direction);
  heading_ += angular;
}

/**
 * \brief Other possible integration method provided by the class
 * \param linear
 * \param angular
 */
void Odometry::integrateExact(double linear, double angular)
{
  if (fabs(angular) < 1e-6)
    integrateRungeKutta2(linear, angular);
  else
  {
    /// Exact integration (should solve problems when angular is zero):
    const double heading_old = heading_;
    const double r = linear/angular;
    heading_ += angular;
    x_       +=  r * (sin(heading_) - sin(heading_old));
    y_       += -r * (cos(heading_) - cos(heading_old));
  }
}

void Odometry::resetAccumulators()
{
  linear_acc_ = RollingMeanAcc(RollingWindow::window_size = velocity_rolling_window_size_);
  angular_acc_ = RollingMeanAcc(RollingWindow::window_size = velocity_rolling_window_size_);
}


text: #include <boost/bind/storage.hpp>

max_stars_repo_name: miathedev/BoostForArduino

max_issues_repo_name: miathedev/BoostForArduino

max_forks_repo_name: miathedev/BoostForArduino


text:
#include <algorithm>
#include <csignal>
#include <cstdlib>
#include <fstream>
#include <optional>
#include <string>
#include <string_view>
#include <thread>

#include <boost/asio/ip/basic_endpoint.hpp>
#include <boost/asio/ip/tcp.hpp>
#include <boost/asio/signal_set.hpp>
#include <boost/lexical_cast.hpp>
#include <boost/log/utility/setup/from_stream.hpp>
#include <cxxopts.hpp>
#include <sol/sol.hpp>

#include <kr/log/logging.h>
#include <kr/utility/uri.h>

int main(int argc, char *argv[])
{
    cxxopts::Options options("Allowed options");
    options.add_options()
        ("help", "help message")
        ("config", "server config file path", cxxopts::value<std::string>()->default_value("etc/server.app.lua"));

    auto vm = options.parse(argc, argv);
    if (vm.count("help"))
    {
        std::cout << options.help() << std::endl;
        return EXIT_SUCCESS;
    }

    // 配置读取
    auto config_path = vm["config"].as<std::string>();
    sol::state config;
    config.open_libraries();
    if (auto result = config.do_file(config_path); !result.valid())
    {
        std::cout << "load lua config error:" << result.get<sol::error>().what() << std::endl;
        return EXIT_FAILURE;
    }

    // 日志配置
    if (std::optional<std::string> log_conf = config["log_conf"]; log_conf)
    {
        std::ifstream settings(*log_conf);
        boost::log::init_from_stream(settings);
    }

    // 服务配置
    auto address = config["address"].get<std::string>();
    auto thread_count = config["thread"].get<int>();

    // 地址解析
    auto to_endpoint = [](std::string_view address) {
        kr::utility::uri uri(address);
        return boost::asio::ip::tcp::endpoint(boost::asio::ip::address::from_string(uri.host().data()),
            boost::lexical_cast<std::uint16_t>(uri.port()));
    };

    // 主线程logger
    kr::log::logger logger{"Service"};
    KRLOG_INFO(logger, "Main") << "Service Start.";
    // 调度器
    boost::asio::io_context io_context;

    // 信号处理
    boost::asio::signal_set signal_set{io_context};
#if defined(SIGUSR1)
    signal_set.add(SIGUSR1);
#endif
#if defined(SIGUSR2)
    signal_set.add(SIGUSR2);
#endif

    std::function<void(const boost::system::error_code &error, int signal_number)> signal_handler;
    signal_handler = [&](const boost::system::error_code &error, int signal_number) {
#if defined(SIGUSR1)
        if (signal_number == SIGUSR1)
        {
            KRLOG_INFO(logger, "Signal") << "Siguser1";
            io_context.stop();
        }
#endif

#if defined(SIGUSR2)
        // TODO: 预留信号: 配置重读,服务存盘
        if (signal_number == SIGUSR2)
        {
            KRLOG_INFO(logger, "Signal") << "Signal2";
        }
#endif
        signal_set.async_wait(signal_handler);
    };
    signal_set.async_wait(signal_handler);

    // 开启服务
    std::vector<std::thread> threads;
    std::generate_n(std::back_insert_iterator(threads), thread_count, [&io_context]() {
        return std::thread([&io_context] { io_context.run(); });
    });
    for (auto &&thread : threads)
    {
        thread.join();
    }
    return EXIT_SUCCESS;
}

max_stars_repo_name: blackkaiserxjc/kraken

max_issues_repo_name: blackkaiserxjc/kraken

max_forks_repo_name: blackkaiserxjc/kraken


text:
//
//  Copyright (c) 2009-2011 Artyom Beilis (Tonkikh)
//
//  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)
//
#ifdef BOOST_LOCALE_NO_WINAPI_BACKEND
#include <iostream>
int main()
{
        std::cout << "WinAPI Backend is not build... Skipping" << std::endl;
}
#else

#include <boost/locale/conversion.hpp>
#include <boost/locale/localization_backend.hpp>
#include <boost/locale/generator.hpp>
#include <boost/locale/info.hpp>
#include <iomanip>
#include "test_locale.hpp"
#include "test_locale_tools.hpp"
#include <iostream>

template<typename CharType>
void test_one(std::locale const &l,std::string src,std::string tgtl,std::string tgtu)
{
    TEST(boost::locale::to_upper(to_correct_string<CharType>(src,l),l) == to_correct_string<CharType>(tgtu,l));
    TEST(boost::locale::to_lower(to_correct_string<CharType>(src,l),l) == to_correct_string<CharType>(tgtl,l));
    TEST(boost::locale::fold_case(to_correct_string<CharType>(src,l),l) == to_correct_string<CharType>(tgtl,l));
}

template<typename CharType>
void test_char()
{
    boost::locale::generator gen;

    std::cout << "- Testing at least C" << std::endl;

    std::locale l = gen("en_US.UTF-8");

    test_one<CharType>(l,"Hello World i","hello world i","HELLO WORLD I");

    std::string name = "en_US.UTF-8";

    std::cout << "- Testing " << name << std::endl;
    l=gen(name);
    test_one<CharType>(l,"Façade","façade","FAÇADE");

    name = "tr_TR.UTF-8";
    std::cout << "Testing " << name << std::endl;
    test_one<CharType>(gen(name),"i","i","İ");

}

template<typename Char>
void test_normc(std::basic_string<Char> orig,std::basic_string<Char> normal,boost::locale::norm_type type)
{
    std::locale l = boost::locale::generator().generate("en_US.UTF-8");
    TEST(boost::locale::normalize(orig,type,l)==normal);
    TEST(boost::locale::normalize(orig.c_str(),type,l)==normal);
    TEST(boost::locale::normalize(orig.c_str(),orig.c_str()+orig.size(),type,l)==normal);
}

void test_norm(std::string orig,std::string normal,boost::locale::norm_type type)
{
    test_normc<char>(orig,normal,type);
    test_normc<wchar_t>(to<wchar_t>(orig),to<wchar_t>(normal),type);
}

int main()
{
    try {
        boost::locale::localization_backend_manager mgr = boost::locale::localization_backend_manager::global();
        mgr.select("winapi");
        boost::locale::localization_backend_manager::global(mgr);

        std::cout << "Testing char" << std::endl;
        test_char<char>();
        std::cout << "Testing wchar_t" << std::endl;
        test_char<wchar_t>();

        std::cout << "Testing Unicode normalization" << std::endl;
        test_norm("\xEF\xAC\x81","\xEF\xAC\x81",boost::locale::norm_nfd); /// ligature fi
        test_norm("\xEF\xAC\x81","\xEF\xAC\x81",boost::locale::norm_nfc);
        #if defined(_WIN32_NT) && _WIN32_NT >= 0x600
        test_norm("\xEF\xAC\x81","fi",boost::locale::norm_nfkd);
        test_norm("\xEF\xAC\x81","fi",boost::locale::norm_nfkc);
        #endif
        test_norm("ä","ä",boost::locale::norm_nfd); // ä to a and accent
        test_norm("ä","ä",boost::locale::norm_nfc);
    }
    catch(std::exception const &e) {
        std::cerr << "Failed " << e.what() << std::endl;
        return EXIT_FAILURE;
    }
    FINALIZE();

}

#endif // no winapi
// vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4

// boostinspect:noascii