text:
/////////1/////////2/////////3/////////4/////////5/////////6/////////7/////////8
// Name        :
// Author      : Avi
// Revision    : $Revision: #4 $
//
// Copyright 2009-2020 ECMWF.
// This software is licensed under the terms of the Apache Licence version 2.0
// which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
// In applying this licence, ECMWF does not waive the privileges and immunities
// granted to it by virtue of its status as an intergovernmental organisation
// nor does it submit to any jurisdiction.
//
// Description :
/////////1/////////2/////////3/////////4/////////5/////////6/////////7/////////8

#include <iostream>
#include <boost/lexical_cast.hpp>
#include "SClientHandleCmd.hpp"

using namespace ecf;
using namespace std;

bool SClientHandleCmd::equals(ServerToClientCmd* rhs) const
{
	return (dynamic_cast<SClientHandleCmd*>(rhs)) ? ServerToClientCmd::equals(rhs) : false;
}

std::string SClientHandleCmd::print() const
{
	string os;
	os += "cmd:SClientHandleCmd [ ";
	os += boost::lexical_cast<string>(handle_);
	os += " ]";
	return os;
}

bool SClientHandleCmd::handle_server_response( ServerReply& server_reply, Cmd_ptr cts_cmd, bool debug ) const
{
	if (debug) cout << "  SClientHandleCmd::handle_server_response handle_ = " << handle_ << "\n";
	server_reply.set_client_handle( handle_ );
	return true;
}

std::ostream& operator<<(std::ostream& os, const SClientHandleCmd& c){ os << c.print(); return os; }

ext: cpp

max_stars_repo_name: ecmwf/ecflow

max_issues_repo_name: CoollRock/ecflow

max_forks_repo_name: CoollRock/ecflow


text:
//
// Copyright (c) 2013-2017 Vinnie Falco (vinnie dot falco at gmail dot com)
//
// Distributed under the Boost Software License, Version 1.0. (See accompanying
// file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
//

#ifndef BEAST_VERSION_HPP
#define BEAST_VERSION_HPP

#include <beast/config.hpp>
#include <boost/config.hpp>

/** @def BEAST_API_VERSION

    Identifies the API version of Beast.

    This is a simple integer that is incremented by one every time
    a set of code changes is merged to the master or develop branch.
*/
#define BEAST_VERSION 55

#define BEAST_VERSION_STRING "Beast/" BOOST_STRINGIZE(BEAST_VERSION)

#endif

ext: hpp

max_stars_repo_name: vinniefalco/rippled

max_issues_repo_name: vinniefalco/rippled

max_forks_repo_name: vinniefalco/rippled


text:
#ifndef __SMT_HPP_DEFINED__
#define __SMT_HPP_DEFINED__

#include <boost/shared_ptr.hpp>
#include <ilasynth/ast.hpp>
#include <ilasynth/ast/hash.hpp>
#include <ilasynth/exception.hpp>
#include <ilasynth/util.hpp>
#include <unordered_map>
#include <z3++.h>

namespace ilasynth {
// A function object that converts nodes into Z3 expressions.
class Z3ExprAdapter {
public:
  // Define types.
  typedef std::unordered_map<const Node*, z3::expr, decltype(&nodeHash),
                             decltype(&nodeEqual)>
      expr_map_t;
  // simplify flag.
  bool simplify;

protected:
  // These are used during memoization.
  expr_map_t exprmap;
  expr_map_t cnstmap;

  // The underlying Z3 context.
  z3::context& c;
  // The "suffix" (used only for synthesis operators).
  std::string suffix;

  /// Added to every name
  std::string name_suffix;

  // Utility function extracts the Z3 expression corresponding
  // to the i-th argument of the Node n.
  z3::expr _getArg(const expr_map_t& m, const Node* n, int i);
  // This extracts expressions from exprmap.
  z3::expr getArgExpr(const Node* n, int i) { return _getArg(exprmap, n, i); }
  // This extracts expressions from cnstmap.
  z3::expr getArgCnst(const Node* n, int i) { return _getArg(cnstmap, n, i); }

  // These functions populate exprmap and cnstmap.
  void _populateExprMap(const Node* n);
  void _populateCnstMap(const Node* n);

  // Convert a boolean variable into a Z3 expression.
  virtual z3::expr getBoolVarExpr(const BoolVar* bv);
  // Convert a boolean op into a Z3 expression.
  virtual z3::expr getBoolOpExpr(const BoolOp* op);
  // for booleans
  virtual z3::expr getChoiceExpr(const BoolChoice* op);

  // Convert a bitvector variable into a Z3 expression.
  virtual z3::expr getBitvectorVarExpr(const BitvectorVar* bvv);
  // Convert a bitvector op into a Z3 expression.
  virtual z3::expr getBvOpExpr(const BitvectorOp* op);
  // Convert a choice op into a Z3 expression. (bitvectors)
  virtual z3::expr getChoiceExpr(const BitvectorChoice* op);
  // Convert a in-range expression to Z3.
  virtual z3::expr getBVInRangeExpr(const BVInRange* op);
  // Constraints for the in-range operator.
  virtual z3::expr getBVInRangeCnst(const BVInRange* op);

  // Convert a memory var into a Z3 expression.
  virtual z3::expr getMemVarExpr(const MemVar* mv);
  // Convert a mem write expression into Z3.
  virtual z3::expr getMemOpExpr(const MemOp* mw);
  // for memories.
  virtual z3::expr getChoiceExpr(const MemChoice* op);

  // Convert a function variable into a Z3 expression.
  virtual z3::expr getFuncVarExpr(const FuncVar* fv);

public:
  // Constructors.
  Z3ExprAdapter(z3::context& c, const std::string& suffix);
  Z3ExprAdapter(z3::context& c, const char* suffix);
  // Destructor.
  virtual ~Z3ExprAdapter();

  // This is used by depthFirstVisit.
  virtual void operator()(const Node* n);

  // This will call node->depthFirstVisit.
  z3::expr getExpr(const Node* n);
  // As will this.
  z3::expr getCnst(const Node* n);

  // return the context.
  z3::context& ctx() const { return c; }

  // clear the memo.
  void clear() {
    exprmap.clear();
    cnstmap.clear();
  }

  // Extract a string representation of the numeric value
  // of node r in the model m.
  std::string extractNumeralString(z3::model& m, const Node* r);
  // Extract the integer value of node r in model m.
  int getNumeralInt(z3::model& m, const Node* r);
  // Extract the arbitrary precision int value of node r in model m.
  mp_int_t getNumeralCppInt(z3::model& m, const Node* r);
  // Extract the boolean value of node r in model m.
  bool getBoolValue(z3::model& m, const Node* r);
  // Return the value of the ith choice boolean
  template <typename T>
  bool getChoiceBool(z3::model& m, const ChoiceExpr<T>* op, int i) {
    using namespace z3;
    std::string name = op->getChoiceVarName(i) + suffix;
    expr ci = c.bool_const(name.c_str());
    z3::expr mi = m.eval(ci, true);
    Z3_lbool bi = Z3_get_bool_value(c, mi);
    // FIXME avoid macro in template
    // ILA_ASSERT(bi != Z3_L_UNDEF, "Unable to extract bool from model.");
    return bi == Z3_L_TRUE;
  }

  void setNameSuffix(const std::string& ns) { name_suffix = ns; }
  const std::string& getNameSuffix() const { return name_suffix; }

private:
  // Helper to convert choices into z3 expressions.
  template <typename T> z3::expr _getChoiceExpr(const ChoiceExpr<T>* op) {
    using namespace z3;

    expr vi_ = getArgExpr(op, 0);
    unsigned nargs = op->nArgs();
    for (unsigned i = 1; i != nargs; i++) {
      std::string name = op->getChoiceVarName(i - 1) + suffix;
      expr ci = c.bool_const(name.c_str());
      expr vi = getArgExpr(op, i);
      expr vi_next = ite(ci, vi, vi_);
      vi_ = vi_next;
    }
    return vi_;
  }
};

class DistInput;

// The function object we use during synthesis to rewrite expressions.
class Z3ExprRewritingAdapter : public Z3ExprAdapter {
protected:
  // The distinguishing input this is based on.
  const DistInput* distInp;
  // Convert a boolean variable into a Z3 expression.
  virtual z3::expr getBoolVarExpr(const BoolVar* bv);
  // Convert a bitvector variable into a Z3 expression.
  virtual z3::expr getBitvectorVarExpr(const BitvectorVar* bvv);
  // Convert a mem var into Z3.
  virtual z3::expr getMemVarExpr(const MemVar* mv);

public:
  // Constructors.
  Z3ExprRewritingAdapter(z3::context& c, const DistInput* di,
                         const std::string& suffix);
  Z3ExprRewritingAdapter(z3::context& c, const DistInput* di,
                         const char* suffix);
  // Destructor.
  ~Z3ExprRewritingAdapter();
  // wrapper over getExpr and getCnst
  z3::expr getIOCnst(const Node* n, const py::object& result);
};

class Z3FixedpointAdapter : public Z3ExprAdapter {
public:
  Z3FixedpointAdapter(z3::context& c, const std::string& suffix)
      : Z3ExprAdapter(c, suffix) {}

  ~Z3FixedpointAdapter() {}
};
} // namespace ilasynth

#endif

max_stars_repo_name: Bo-Yuan-Huang/ItSy


text:
#ifndef BOOST_NETWORK_PROTOCOL_HTTP_CLIENT_MACROS_HPP_20110430
#define BOOST_NETWORK_PROTOCOL_HTTP_CLIENT_MACROS_HPP_20110430

// Copyright 2011 Dean Michael Berris <mikhailberis@gmail.com>.
// Distributed under the Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)

#include <boost/range/iterator_range.hpp>
#include <system_error>

#ifndef BOOST_NETWORK_HTTP_BODY_CALLBACK
#define BOOST_NETWORK_HTTP_BODY_CALLBACK(function_name, range_name,        \
                                         error_name)                       \
  void function_name(boost::iterator_range<const char*> const& (range_name), \
                     std::error_code const& (error_name))
#endif

#endif /* BOOST_NETWORK_PROTOCOL_HTTP_CLIENT_MACROS_HPP_20110430 */

max_stars_repo_name: ccnyou/cpp-netlib-stable

max_issues_repo_name: ccnyou/cpp-netlib-stable

max_forks_repo_name: ccnyou/cpp-netlib-stable


text:
// Boost.Geometry (aka GGL, Generic Geometry Library)

// Copyright (c) 2007-2015 Barend Gehrels, Amsterdam, the Netherlands.
// Copyright (c) 2008-2015 Bruno Lalande, Paris, France.
// Copyright (c) 2009-2015 Mateusz Loskot, London, UK.
// Copyright (c) 2014-2015 Samuel Debionne, Grenoble, France.

// This file was modified by Oracle on 2015-2020.
// Modifications copyright (c) 2015-2020, Oracle and/or its affiliates.

// Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
// Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle

// Distributed under the Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)

#ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_EXPAND_SEGMENT_HPP
#define BOOST_GEOMETRY_ALGORITHMS_DETAIL_EXPAND_SEGMENT_HPP

#include <boost/core/ignore_unused.hpp>

#include <boost/geometry/algorithms/dispatch/expand.hpp>

#include <boost/geometry/core/tags.hpp>

namespace boost { namespace geometry
{

#ifndef DOXYGEN_NO_DISPATCH
namespace dispatch
{

template
<
    typename Box, typename Segment
>
struct expand
    <
        Box, Segment,
        box_tag, segment_tag
    >
{
    template <typename Strategy>
    static inline void apply(Box& box,
                             Segment const& segment,
                             Strategy const& strategy)
    {
        strategy.expand(box, segment).apply(box, segment);
    }
};

} // namespace dispatch
#endif // DOXYGEN_NO_DISPATCH

}} // namespace boost::geometry

#endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_EXPAND_SEGMENT_HPP


text:
#include <cstdio>
#include <iostream>
#include <boost/scoped_ptr.hpp>
#include <boost/filesystem.hpp>
#include <opencv2/core/core.hpp>
#include <opencv2/highgui/highgui.hpp>

#include "MatlabIOModel.hpp"
#include "PoseEstimator.hpp"

using namespace std;
using namespace cv;

PoseEstimator::PoseEstimator(const string& model_filename) {
  boost::scoped_ptr<Model> model;
  model.reset(new MatlabIOModel);

	bool ok = model->deserialize(model_filename);
	if (!ok) {
		printf("Error deserializing file\n");
		exit(-3);
	}

  pbd.distributeModel(*model);
}

//cv::Mat copyPixelArrayToMat(int[][][] pixelArray, int rows, int cols) {
//  const cv::Mat& image(rows, cols, CV_8UC3);
//  for(int y = 0; y < image.rows; y++){
//    for(int x = 0; x < image.cols; x++){
//      image.data[image.step[0]*y + image.step[1]* x + 0] = pixelArray[y][x][0];
//      image.data[image.step[0]*y + image.step[1]* x + 1] = pixelArray[y][x][1];
//      image.data[image.step[0]*y + image.step[1]* x + 2] = pixelArray[y][x][2];
//    }
//  }
//  return image;
//}

vector<Candidate> PoseEstimator::estimate(const string& image_filename) {
	Mat_<float> depth;
	Mat im = imread(image_filename);
  if (im.empty()) {
      printf("Image not found or invalid image format\n");
      exit(-4);
  }

	vector<Candidate> candidates;
	pbd.detect(im, depth, candidates);

  return candidates;
}

//vector<Candidate> PoseEstimator::estimate(int[][][] pixelArray,
//                                          int rows, int cols) {
//  ///3D pixelArray to cv::Mat
//  cv::Mat& image = copyPixelArrayToMat(pixelArray, rows, cols);
//  //TODO show image once here
//
//  vector<Candidate> candidates;
//  pbd.detect(image, candidates);
//
//  return candidates;
//}

max_stars_repo_name: IshitaTakeshi/VirtualFitting

max_issues_repo_name: IshitaTakeshi/VirtualFitting

max_forks_repo_name: IshitaTakeshi/VirtualFitting


text:
//
// Copyright (c) 2015 CNRS
//

#ifndef __pinocchio_force_set_hpp__
#define __pinocchio_force_set_hpp__

#include "pinocchio/spatial/fwd.hpp"
#include <Eigen/Geometry>

namespace pinocchio
{
  template<typename _Scalar, int _Options>
  class ForceSetTpl
  {
  public:
    typedef _Scalar Scalar;
    enum { Options = _Options };
    typedef Eigen::Matrix<Scalar,3,1,Options> Vector3;
    typedef Eigen::Matrix<Scalar,3,3,Options> Matrix3;
    typedef Eigen::Matrix<Scalar,6,1,Options> Vector6;
    typedef Eigen::Matrix<Scalar,6,6,Options> Matrix6;
    typedef SE3Tpl<Scalar,Options> SE3;

    typedef Eigen::Matrix<Scalar,3,Eigen::Dynamic,Options> Matrix3x;
    typedef Eigen::Matrix<Scalar,6,Eigen::Dynamic,Options> Matrix6x;

  public:
    // Constructors
    ForceSetTpl(const int & ncols ) : size(ncols),m_f(3,ncols), m_n(3,ncols)
    { m_f.fill(NAN); m_n.fill(NAN); }
    ForceSetTpl(const Matrix3x & linear, const Matrix3x & angular)
      : size((int)linear.cols()),m_f(linear), m_n(angular)
    {  assert( linear.cols() == angular.cols() ); }

    Matrix6x matrix() const
    {
      Matrix6x F(6,size); F << m_f, m_n;
      // F.template topRows<3>() = m_f;
      // F.template bottomRows<3>()  = m_n;
      return F;
    }
    operator Matrix6x () const { return matrix(); }

    // Getters
    const Matrix3x & linear() const { return m_f; }
    const Matrix3x & angular() const { return m_n; }

    /// af = aXb.act(bf)
    ForceSetTpl se3Action(const SE3 & m) const
    {
      Matrix3x Rf (m.rotation()*linear());
      return ForceSetTpl(Rf,skew(m.translation())*Rf+m.rotation()*angular());
      // TODO check if nothing better than explicitely calling skew
    }
    /// bf = aXb.actInv(af)
    ForceSetTpl se3ActionInverse(const SE3 & m) const
    {
      return ForceSetTpl(m.rotation().transpose()*linear(),
		      m.rotation().transpose()*(angular() - skew(m.translation())*linear()) );
      // TODO check if nothing better than explicitely calling skew
    }

    friend std::ostream & operator << (std::ostream & os, const ForceSetTpl & phi)
    {
      os
	<< "F =\n" << phi.linear() << std::endl
	<< "Tau =\n" << phi.angular() << std::endl;
      return os;
    }

    /* --- BLOCK ------------------------------------------------------------ */
    struct Block
    {
      ForceSetTpl & ref;
      int idx,len;
      Block( ForceSetTpl & ref, const int & idx, const int & len )
	: ref(ref), idx(idx), len(len) {}

      Eigen::Block<ForceSetTpl::Matrix3x> linear() { return ref.m_f.block(0,idx,3,len); }
      Eigen::Block<ForceSetTpl::Matrix3x> angular() { return ref.m_n.block(0,idx,3,len); }
      Eigen::Block<const ForceSetTpl::Matrix3x> linear()  const
      { return ((const ForceSetTpl::Matrix3x &)(ref.m_f)).block(0,idx,3,len); }
      Eigen::Block<const ForceSetTpl::Matrix3x> angular() const
      { return ((const ForceSetTpl::Matrix3x &)(ref.m_n)).block(0,idx,3,len); }

      ForceSetTpl::Matrix6x matrix() const
      {
	ForceSetTpl::Matrix6x res(6,len); res << linear(),angular();
	return res;
      }

      Block& operator= (const ForceSetTpl & copy)
      {
	assert(copy.size == len);
	linear() = copy.linear(); //ref.m_f.block(0,idx,3,len) = copy.m_f;
	angular() = copy.angular(); //ref.m_n.block(0,idx,3,len) = copy.m_n;
	return *this;
      }

      Block& operator= (const ForceSetTpl::Block & copy)
      {
	assert(copy.len == len);
	linear() = copy.linear(); //ref.m_f.block(0,idx,3,len) = copy.ref.m_f.block(0,copy.idx,3,copy.len);
	angular() = copy.angular(); //ref.m_n.block(0,idx,3,len) = copy.ref.m_n.block(0,copy.idx,3,copy.len);
	return *this;
      }

      template <typename D>
      Block& operator= (const Eigen::MatrixBase<D> & m)
      {
        eigen_assert(D::RowsAtCompileTime == 6);
        assert(m.cols() == len);
        linear() = m.template topRows<3>();
        angular() = m.template bottomRows<3>();
        return *this;
      }

      /// af = aXb.act(bf)
      ForceSetTpl se3Action(const SE3 & m) const
      {
	// const Eigen::Block<const Matrix3x> linear = ref.linear().block(0,idx,3,len);
	// const Eigen::Block<const Matrix3x> angular = ref.angular().block(0,idx,3,len);
	Matrix3x Rf ((m.rotation()*linear()));
	return ForceSetTpl(Rf,skew(m.translation())*Rf+m.rotation()*angular());
	// TODO check if nothing better than explicitely calling skew
      }
      /// bf = aXb.actInv(af)
      ForceSetTpl se3ActionInverse(const SE3 & m) const
      {
	// const Eigen::Block<const Matrix3x> linear = ref.linear().block(0,idx,3,len);
	// const Eigen::Block<const Matrix3x> angular = ref.angular().block(0,idx,3,len);
	return ForceSetTpl(m.rotation().transpose()*linear(),
			   m.rotation().transpose()*(angular() - skew(m.translation())*linear()) );
	// TODO check if nothing better than explicitely calling skew
      }

    };

    Block block(const int & idx, const int & len) { return Block(*this,idx,len); }

    /* CRBA joint operators
     *   - ForceSet::Block = ForceSet
     *   - ForceSet operator* (Inertia Y,Constraint S)
     *   - MatrixBase operator* (Constraint::Transpose S, ForceSet::Block)
     *   - SE3::act(ForceSet::Block)
     */

  public: //
  private:
    int size;
    Matrix3x m_f,m_n;
  };

  typedef ForceSetTpl<double,0> ForceSet;

  namespace internal
  {
    template<>
    struct SE3GroupAction<ForceSet::Block>    { typedef ForceSet ReturnType; };
  }

} // namespace pinocchio

#endif // ifndef __pinocchio_force_set_hpp__

max_stars_repo_name: matthieuvigne/pinocchio

max_issues_repo_name: matthieuvigne/pinocchio

max_forks_repo_name: matthieuvigne/pinocchio


text:
//
// Copyright (c) 2020 Richard Hodges (hodges.r@gmail.com)
// Copyright (c) 2021 Sergey Ilinykh (rion4ik@gmail.com)
//
// Distributed under the Boost Software License, Version 1.0. (See accompanying
// file LICENSE or copy at http://www.boost.org/LICENSE_1_0.txt)
//

#include "http/connection_cache.hpp"

#include <boost/program_options.hpp>

#include <cstdio>
#include <iostream>

net::awaitable< void >
visit_site(http::connection_cache &cache, std::string const &url, bool verbose = false)
{
    auto then = std::chrono::steady_clock::now();

    try
    {
        auto result = co_await cache.call(http::verb::get, url);
        auto time   = std::chrono::duration_cast< std::chrono::milliseconds >(
            std::chrono::steady_clock::now() - then);
        if (verbose)
        {
            std::cout << "GET " << url << " -> " << result->base().result_int()
                      << " " << result->base().reason() << " in "
                      << (time).count() << "ms\n";
        }
        std::cout << result->body();
    }
    catch (std::exception &e)
    {
        auto time = std::chrono::duration_cast< std::chrono::milliseconds >(
            std::chrono::steady_clock::now() - then);
        if (verbose)
        {
            std::cerr << "GET " << url << " -> "
                      << " exception: " << e.what() << " in " << (time).count()
                      << "ms\n";
        }
    }
}

int
main(int argc, char **argv)
{
    namespace po = boost::program_options;
    po::options_description desc("Allowed options");
    // clang-format off
    desc.add_options()
        ("help,h", "produce a help message")
        ("verbose,v", "verbose output")
        ("url", po::value<std::vector<std::string>>(), "input file");

    po::positional_options_description p;
    p.add("url", -1);

    po::variables_map vm;
    po::store(po::command_line_parser(argc, argv).options(desc).positional(p).run(), vm);
    po::notify(vm);

    if (vm.count("help")) {
        std::cout << desc << "\n";
        return 0;
    }

    if (!vm.count("url")) {
        std::cerr << "No input file";
        return -1;
    }
    auto urls = vm["url"].as<std::vector<std::string>>();
    bool verbose = vm.count("verbose");

    net::io_context ioctx;

    http::connection_cache cache(ioctx.get_executor());
    for (auto const &url: urls) {
        net::co_spawn(
            ioctx.get_executor(), visit_site(cache, url, verbose), net::detached);
    }

    ioctx.run();
}


text:
//
//  Eigenvalues of the radial part of a 3D exponential well
//
// compile with:
//     g++ -std=c++11 exp_well_3d.cpp -o exp_well_3d -I$HOME/miniconda3/include/ -O3
//
// using the a Miniconda install of boost for the headers

#include <cstdio>
#include <iostream>
#include <cmath>
#include <boost/array.hpp>
#include <boost/numeric/odeint.hpp>

using namespace std;
using namespace boost::numeric::odeint;

typedef boost::array< long double , 2 > vec_type;
typedef runge_kutta_cash_karp54< vec_type > error_stepper_type;
typedef controlled_runge_kutta< error_stepper_type > controlled_stepper_type;
typedef runge_kutta4< vec_type > stepper;

struct ode{

    long double m;
    long double energy;
    long double a;
    long double k;
    long double l;     // angluar quantum number

    void operator()( const vec_type &vec , vec_type &dvecdx , double r)
    {
        dvecdx[0] = vec[1];
        dvecdx[1] =  -2.0*vec[1]/r - (2*m*(energy - k*(exp(a*r) - 1.0)) - (l*(l+1)/(r*r)))*vec[0];
    }
};

struct nothing{

    void operator()( const vec_type &vec , const double r){

    }
};

void print_exp_well_psi_final(double k, double a, double m){
    controlled_stepper_type controlled_stepper;

    // Energy sweep range
    double min_energy = 0.000001;
    double max_energy = 10 * 0.00094;    // kBT = 0.00094 Ha

    int l_max = 50;                       // Maximum value of the l qunatum number

    double r_init = 0.001;               // Where to start the propogation from
    double r_final = 1.7;                // Where to end   '   '  as there is a divergence at r=0

    int n_steps = 20000;

    for (int l = 0; l <= l_max; l++ ){

        double prev_psi_final = 0.0;         // Final value of Psi i.e. Psi(r_final) for each eigenvalue

        for (int i = 0; i < n_steps; i++){

            double energy = min_energy + ((max_energy - min_energy) * i) / n_steps;

            vec_type vec = { 1.0 , 1E-20 }; // initial conditions
            integrate_adaptive( make_controlled< error_stepper_type >( 1.0e-10 , 1.0e-10 ),
                               ode {m , energy, a, k, static_cast<long double>(l)},
                               vec,
                               r_init,
                               r_final,
                               1E-10,
                               nothing {}
                               );

            // If the product of the previous Psi(r_final) and the current one is negative
            // there must be a zero between them, i.e. a normalaisable WF
            if (prev_psi_final * vec[0] < 0){
                // cout << energy << endl;
                cout << energy << '\t' << l << endl;

            }

            // Reset the previous Psi(r_final)
            prev_psi_final = vec[0];

        } // Energy
    }  // l
}

int main(int argc, char **argv){

    freopen("eigenvals_co2.txt","w", stdout);
    cout << "Eigenvalue / Ha    l number" << endl;

    double k_co2 = 0.00220873;                   // prefactor for the exponential well

    double k = 0.0001;
    double a = 2.832 * 0.529177;             // exponent in the expoential well
    double m = 48 * 1822.888486;             // mass in atomic units from amu

    print_exp_well_psi_final(k_co2, a, m);

    return 0;
}


text: #include <boost/preprocessor/list/at.hpp>

max_stars_repo_name: miathedev/BoostForArduino

max_issues_repo_name: miathedev/BoostForArduino

max_forks_repo_name: miathedev/BoostForArduino


text:
//
// Copyright 2016 Pixar
//
// Licensed under the Apache License, Version 2.0 (the "Apache License")
// with the following modification; you may not use this file except in
// compliance with the Apache License and the following modification to it:
// Section 6. Trademarks. is deleted and replaced with:
//
// 6. Trademarks. This License does not grant permission to use the trade
//    names, trademarks, service marks, or product names of the Licensor
//    and its affiliates, except as required to comply with Section 4(c) of
//    the License and to reproduce the content of the NOTICE file.
//
// You may obtain a copy of the Apache License at
//
//     http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the Apache License with the above modification is
// distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
// KIND, either express or implied. See the Apache License for the specific
// language governing permissions and limitations under the Apache License.
//
#include "pxr/base/tf/weakPtr.h"

/*
 * Compile-time testing of the Tf_SupportsWeakPtr mechanism.  Change the 0
 * to 1 to enable.
 */
#if 0

#include <boost/static_assert.hpp>

struct _Tf_TestHasGetWeakBase {
    TfWeakBase const &__GetTfWeakBase__() const;
};

struct _Tf_TestHasGetWeakBaseDerived : public _Tf_TestHasGetWeakBase
{
};

struct _Tf_TestHasGetWeakBaseNot
{
};

struct _Tf_TestIsWeakBase : public TfWeakBase
{
};

BOOST_STATIC_ASSERT(TF_SUPPORTS_WEAKPTR(_Tf_TestHasGetWeakBase));
BOOST_STATIC_ASSERT(TF_SUPPORTS_WEAKPTR(_Tf_TestHasGetWeakBaseDerived));
BOOST_STATIC_ASSERT(not TF_SUPPORTS_WEAKPTR(_Tf_TestHasGetWeakBaseNot));
BOOST_STATIC_ASSERT(not TF_SUPPORTS_WEAKPTR(TfWeakPtr<_Tf_TestIsWeakBase>));
#endif // testing Tf_SupportsWeakPtr.


text:
/*
 * Model.cpp
 *
 *  Created on: Jun 2, 2009
 *      Author: David Berrios
 */

#include "Model.h"
#include <iostream>
#include <algorithm>
#include <boost/foreach.hpp>
#include "Interpolator.h"
#include "Constants.h"

namespace ccmc
{
	/**
	 * Default constructor
	 */
	Model::Model()
	{
		missingValue = ccmc::defaults::missingValue;
		units_ = "units";
		variablesAdded = 0;
	}

	/**
	 * @brief Sets the model name to modelName.
	 *
	 * @param modelName
	 */
	void Model::setModelName(std::string modelName)
	{
		this->modelName = modelName;
	}

	/**
	 * @brief Returns the modelName of the currently selected file.
	 *
	 * @return
	 */
	std::string Model::getModelName()
	{
		return modelName;
	}

	/**
	 * @brief Closes the currently selected file
	 *
	 * @return
	 */
	long Model::close()
	{
		//TODO: clean up memory

		if (variableData.size() > 0)
		{
			boost::unordered_map<std::string, std::vector<float>*>::iterator iter;

			for ( iter=variableData.begin() ; iter != variableData.end(); iter++ )
			{
				delete iter->second;
			}
		}

		if (variableDataInt.size() > 0)
		{
			boost::unordered_map<std::string, std::vector<int>*>::iterator iter;

			for (iter=variableDataInt.begin(); iter != variableDataInt.end(); iter++)
			{
				delete iter->second;
			}
		}

		//the other maps have the same pointers, so they are already deleted.
		variableData.clear();
		variableDataInt.clear();
		variableDataByID.clear();
		variableDataIntByID.clear();

		return GeneralFileReader::close();
	}

	/**
	 * @brief Sets the missing value to use when no valid data exists.
	 *
	 * @param missingValue
	 */
	void Model::setMissingValue(float missingValue)
	{
		this->missingValue = missingValue;
	}

	/**
	 * @brief Returns the missing value that will be returned when no valid data exists.
	 *
	 * This value can be used to check whether an interpolated value is a missing value or not.
	 *
	 * @return Missing value
	 */
	float Model::getMissingValue()
	{
		return missingValue;
	}

	/**
	 * @brief Load a variable into memory.
	 *
	 * Use this method when the variable to load is of type float
	 *
	 * @param variable Variable to load into memory.
	 * @return status of the operation
	 */
	long Model::loadVariable(const std::string& variable)
	{

		//first, check to determine whether variable is already loaded
		if (variableData.find(variable) != variableData.end())
			return FileReader::OK;

		//first check if the variable exists in the file!!
		if (!this->doesVariableExist(variable))
		{
			std::cerr <<"Problem: "<< variable << " does not exist" << std::endl;
			return FileReader::VARIABLE_DOES_NOT_EXIST;
		}
		std::vector<float> * data = getVariable(variable);
		long id = getVariableID(variable);
//std::cout << BOOST_CURRENT_FUNCTION << " id: " << id << std::endl;
		if (data->size() > 0)
		{

			//std::cout << "adding " << variable << " to maps" << std::endl;
			variableData[variable] = data;
			variableDataByID[id] = data;
		} //else return false;
		else
		{
			//std::cout << "not adding " << variable << " to maps" << std::endl;
			delete data;
		}

		return FileReader::OK;
	}

	/**
	 * @brief Unload a variable from memory.
	 * This will return FileReader::OK if the variable is removed from memory, and FileReader::VARIABLE_NOT_IN_MEMORY
	 * if the variable was not already in memory.
	 *
	 * @param variable Variable to unload from memory.
	 * @return status of the operation
	 */
	long Model::unloadVariable(const std::string& variable)
	{
		//bool success = false;

		//first, check to determine whether variable is already loaded
		if (variableData.find(variable) != variableData.end())
		{
			std::vector<float> * data = variableData[variable];
			long id = getVariableID(variable);
			delete data;
			variableData.erase(variable);
			variableDataByID.erase(id);
			return FileReader::OK;
		}
		return FileReader::VARIABLE_NOT_IN_MEMORY;
	}

	/**
	 * @brief Loads a variable into memory.
	 *
	 * Use this method when the variable to load is of type int
	 *
	 * @param variable
	 * @return status of the operation
	 */
	long Model::loadVariableInt(const std::string& variable)
	{

		//first, check to determine whether variable is already loaded
		if (variableDataInt.find(variable) != variableDataInt.end())
			return true;

		//first check if the variable exists in the file!!
		if (!this->doesVariableExist(variable))
		{
			return FileReader::VARIABLE_DOES_NOT_EXIST;
		}
		std::vector<int> * data = getVariableInt(variable);
		long id = getVariableID(variable);
		if (data->size() > 0)
		{
			variableDataInt[variable] = data;
			variableDataIntByID[id] = data;
		} //else return false;

		return FileReader::OK;
	}

	/**
	 * @brief Returns a pointer to the entry in the variableData map containing the variable data.
	 *
	 * This pointer can be modified.
	 *
	 * @param variable Variable to fetch from memory.  This assumes the variable has already been loaded into memory.
	 * If the variable has not been loaded, the pointer will be NULL.
	 *
	 * @return std::vector<float>* of the requested variable.  Note that the pointer points to an entry in
	 * a map. The memory pointed to by the pointer will automatically be freed when
	 * the file is closed, or the Model object is deleted.
	 */
	std::vector<float>* Model::getVariableFromMapRW(const std::string& variable)
	{
		boost::unordered_map<std::string, std::vector<float>*>::iterator iter = variableData.find(variable);

		std::vector<float>* vData = NULL;
		if (iter != variableData.end())
			vData = (*iter).second;
		else
		{
			std::cerr << "variable: " << variable << " not found!!!" << std::endl;
		}
		return vData;
	}

	/**
	 * @brief Returns a const pointer to the entry in the variableData map containing the variable data.
	 *
	 * This pointer cannot be modified.
	 *
	 * @param variable Variable to fetch from memory.  This assumes the variable has already been loaded into memory.
	 * If the variable has not been loaded, the pointer will be NULL.
	 *
	 * @return std::vector<float>* of the requested variable.  Note that the pointer points to an entry in
	 * a map and should not and cannot be deleted. The memory pointed to by the pointer will automatically be freed when
	 * the file is closed, or the Model object is deleted.
	 */
	const std::vector<float>* const Model::getVariableFromMap(const std::string& variable)
	{
		return getVariableFromMapRW(variable);
	}

	/**
	 * @brief Returns the list of variables that have been loaded into memory, using the loadVariable or loadVectorVariable methods
	 *
	 *
	 */
	const std::vector<std::string> Model::getLoadedVariables()
	{
		boost::unordered_map<std::string, std::vector<float>*>::iterator iter;

		std::vector<std::string> variablesLoaded;
		for ( iter=variableData.begin() ; iter != variableData.end(); iter++ )
		{
			variablesLoaded.push_back((*iter).first);
		}
		sort (variablesLoaded.begin(), variablesLoaded.end());

		return variablesLoaded;

	}

	/**
	 * @brief Returns a pointer to the entry in the variableDataInt map containing the variable data.
	 *
	 * This pointer cannot be modified.
	 *
	 * @param variable Variable to fetch from memory.  This assumes the variable has already been loaded into memory.
	 * If the variable has not been loaded, the pointer will be NULL.
	 *
	 * @return std::vector<int>* of the requested variable.  Note that the pointer points to an entry in
	 * a map. The memory pointed to by the pointer will automatically be freed when
	 * the file is closed, or the Model object is deleted.
	 */
	std::vector<int>* Model::getIntVariableFromMapRW(const std::string& variable)
	{
		boost::unordered_map<std::string, std::vector<int>*>::iterator iter = variableDataInt.find(variable);

		std::vector<int>* vData = NULL;
		if (iter != variableDataInt.end())
			vData = (*iter).second;
		return vData;
	}

	/**
	 * @brief Returns a const pointer to the entry in the variableDataInt map containing the variable data.
	 *
	 * This pointer cannot be modified.
	 *
	 * @param variable Variable to fetch from memory.  This assumes the variable has already been loaded into memory.
	 * If the variable has not been loaded, the pointer will be NULL.
	 *
	 * @return std::vector<int>* of the requested variable.  Note that the pointer points to an entry in
	 * a map and should not and cannot be deleted. The memory pointed to by the pointer will automatically be freed when
	 * the file is closed, or the Model object is deleted.
	 */
	const std::vector<int>* const Model::getIntVariableFromMap(const std::string& variable)
	{
		return getIntVariableFromMapRW(variable);
	}

	/**
	 * @brief Returns a pointer to the entry in the variableDataByID map containing the variable data.
	 *
	 * This pointer cannot be modified.
	 *
	 * @param variable_id Variable id of the variable to fetch from memory.  This assumes the variable has already been loaded into memory.
	 * If the variable has not been loaded, the pointer will be NULL.  Request the variable id by using FileReader::getVariableID(const std::string& variable)
	 *
	 * @return std::vector<float>* of the requested variable.  Note that the pointer points to an entry in
	 * a map. The memory pointed to by the pointer will automatically be freed when
	 * the file is closed, or the Model object is deleted.
	 */
	std::vector<float>* Model::getVariableFromMapRW(long variable_id)
	{
		boost::unordered_map<long, std::vector<float>*>::iterator iter = variableDataByID.find(variable_id);

		std::vector<float>* vData = NULL;
		if (iter != variableDataByID.end())
			vData = (*iter).second;
		return vData;
	}

	/**
	 * @brief Returns a const pointer to the entry in the variableDataByID map containing the variable data.
	 *
	 * This pointer cannot be modified.
	 *
	 * @param variable_id Variable id of the variable to fetch from memory.  This assumes the variable has already been loaded into memory.
	 * If the variable has not been loaded, the pointer will be NULL.  Request the variable id by using FileReader::getVariableID(const std::string& variable)
	 *
	 * @return std::vector<float>* of the requested variable.  Note that the pointer points to an entry in
	 * a map and should not and cannot be deleted. The memory pointed to by the pointer will automatically be freed when
	 * the file is closed, or the Model object is deleted.
	 */
	const std::vector<float>* const Model::getVariableFromMap(long variable_id)
	{

		return getVariableFromMapRW(variable_id);
	}

	/**
	 * @brief Returns a pointer to the entry in the variableDataIntByID map containing the variable data.
	 *
	 * This pointer cannot be modified.
	 *
	 * @param variable_id Variable id of the variable to fetch from memory.  This assumes the variable has already been loaded into memory.
	 * If the variable has not been loaded, the pointer will be NULL.  Request the variable id by using FileReader::getVariableID(const std::string& variable)
	 *
	 * @return std::vector<int>* of the requested variable.  Note that the pointer points to an entry in
	 * a map. The memory pointed to by the pointer will automatically be freed when
	 * the file is closed, or the Model object is deleted.
	 */
	std::vector<int>* Model::getIntVariableFromMapRW(long variable_id)
	{
		boost::unordered_map<long, std::vector<int>*>::iterator iter = variableDataIntByID.find(variable_id);

		std::vector<int>* vData = NULL;
		if (iter != variableDataIntByID.end())
			vData = (*iter).second;
		return vData;
	}

	/**
	 * @brief Returns a const pointer to the entry in the variableDataIntByID map containing the variable data.
	 *
	 * This pointer cannot be modified.
	 *
	 * @param variable_id Variable id of the variable to fetch from memory.  This assumes the variable has already been loaded into memory.
	 * If the variable has not been loaded, the pointer will be NULL.  Request the variable id by using FileReader::getVariableID(const std::string& variable)
	 *
	 * @return std::vector<int>* of the requested variable.  Note that the pointer points to an entry in
	 * a map and should not and cannot be deleted. The memory pointed to by the pointer will automatically be freed when
	 * the file is closed, or the Model object is deleted.
	 */
	const std::vector<int>* const Model::getIntVariableFromMap(long variable_id)
	{
		return getIntVariableFromMapRW(variable_id);
	}

	/**
	 * @brief Returns the conversion factor needed to convert the interpolated value to SI units.
	 * @param variable Variable to request the conversion factor for.
	 * @return Conversion factor to convert the specified variable to SI units
	 */
	float Model::getConversionFactorToSI(const std::string& variable)
	{
		boost::unordered_map<std::string, float>::iterator iter = conversionFactorsToSI.find(variable);

		if (iter != conversionFactorsToSI.end())
			return (*iter).second;
		else
			return 1.0f;

	}

	/**
	 * @brief Fetches the native units of the variable as stored in the file.
	 *
	 * @param variable Variable to request units for.
	 * @return The native units of the specified variable, as stored in the file.
	 */
	std::string Model::getNativeUnit(const std::string& variable)
	{
		std::string unit = (getVariableAttribute(variable, units_)).getAttributeString();
		return unit;
	}

	/**
	 * @brief Returns the SI units of the specified variable.
	 *
	 * @param variable
	 * @return
	 */
	std::string Model::getSIUnit(const std::string& variable)
	{
		std::string unit = "";

		boost::unordered_map<std::string, std::string>::iterator iter = variableSIUnits.find(variable);
		if (iter != variableSIUnits.end())
			unit = (*iter).second;
		return unit;
	}

	/**
	 * Destructor
	 */
	Model::~Model()
	{
		//TODO: delete contents of the maps and clear them
		if (variableData.size() > 0)
		{
			boost::unordered_map<std::string, std::vector<float>*>::iterator iter;

			for ( iter=variableData.begin() ; iter != variableData.end(); iter++ )
			{
				delete iter->second;
			}
		}

		if (variableDataInt.size() > 0)
		{
			boost::unordered_map<std::string, std::vector<int>*>::iterator iter;

			for (iter=variableDataInt.begin(); iter != variableDataInt.end(); iter++)
			{
				delete iter->second;
			}
		}

		//the other maps have the same pointers, so they are already deleted.
		variableData.clear();
		variableDataInt.clear();
		variableDataByID.clear();
		variableDataIntByID.clear();
//		std::cout << "model: destructor!" << std::endl;

		/*boost::unordered_map<std::string, std::vector<int>* > variableDataInt;

		 boost::unordered_map<long, std::vector<float>* > variableDataByID;
		 boost::unordered_map<long, std::vector<int>* > variableDataIntByID;*/
	}

	void Model::setBusyStatus(int busyStatus)
	{
		this->busyStatus = busyStatus;
		if (busyStatus == Model::OK)
			this->progress = 100;
		else //busyStatus == Model::BUSY
			this->progress = 0;
	}

	int Model::getBusyStatus()
	{
		return this->busyStatus;
	}

	int Model::getProgress()
	{
		return this->progress;
	}

	/*
	 * Add variable to map
	 */
	void Model::addFloatVariableToMap(const std::string& variable, std::vector<float>* varData){
		int numberOfVariables = getNumberOfVariables(); //number of variables in data file
		long newVariableId = numberOfVariables + variablesAdded; //ID#s start at 0, variablesAdded initially 0
		std::cout<<"adding variable "<<variable<<" to map with id "<< newVariableId << std::endl;
		variableData[variable]=varData;
		variableDataByID[newVariableId]=varData;
		fileReader->addVariableName(variable,newVariableId);
		variablesAdded ++;
		return;
	}

}

max_stars_repo_name: alexanderbock/Kameleon-Converter

max_issues_repo_name: alexanderbock/Kameleon-Converter

max_forks_repo_name: alexanderbock/Kameleon-Converter


text: #include <boost/range/algorithm/binary_search.hpp>


text:
#ifndef __NODE_MAPNIK_EXPRESSION_H__
#define __NODE_MAPNIK_EXPRESSION_H__

#include <v8.h>
#include <node_object_wrap.h>

// mapnik
#include <mapnik/version.hpp>
#if MAPNIK_VERSION >= 200100
#include <mapnik/expression.hpp>
#else
#include <mapnik/filter_factory.hpp>
#endif

// boost
#include <boost/shared_ptr.hpp>

using namespace v8;

class Expression: public node::ObjectWrap {
public:
    static Persistent<FunctionTemplate> constructor;
    static void Initialize(Handle<Object> target);
    static Handle<Value> New(const Arguments &args);
    static Handle<Value> toString(const Arguments &args);
    static Handle<Value> evaluate(const Arguments &args);

    Expression();
    inline mapnik::expression_ptr get() { return this_; }

private:
    ~Expression();
    mapnik::expression_ptr this_;
};

#endif

max_stars_repo_name: calvinmetcalf/node-mapnik

max_issues_repo_name: calvinmetcalf/node-mapnik

max_forks_repo_name: calvinmetcalf/node-mapnik


text:
//
// Created by louis on 16/02/19.
//

#include "MagicHandler.h"
#include <sstream>
#include <boost/algorithm/string/case_conv.hpp>
#include <boost/algorithm/string/trim.hpp>
#include <boost/algorithm/string/classification.hpp>
#include <boost/algorithm/string/split.hpp>
#include <boost/algorithm/string/find_format.hpp>
#include <boost/algorithm/string/formatter.hpp>

using handler::MagicHandler;

namespace handler {
    MagicHandler::MagicHandler(AccountHandler &accountHandler, CombatHandler &combatHandler) :
        accountHandler(accountHandler),
        combatHandler(combatHandler){}

    std::string MagicHandler::getSpells() {
        std::ostringstream spells;

        spells << "\n"
               << "Spells:\n"
               << "-------\n"
               << "  - " << CONFUSE_SPELL_NAME
               << " (causes the target player to temporarily speak in Pig Latin)\n"

               << "  - " << BODY_SWAP_SPELL_NAME
               << " (causes the caster to switch bodies with the target player temporarily)\n"

               << "  - " << DECOY_SPELL_NAME
               << " (creates a decoy of yourself in combat to escape)\n"

               << "  - " << HEAL_SPELL_NAME
               << " (fully restores target player's health. Can not be used if caster or target is in combat)\n";

        return spells.str();
    }

    std::vector<Message>
    MagicHandler::castSpell(const Connection &client, const std::string &spellName, const std::string& targetName) {
        std::string spellNameLower = boost::algorithm::to_lower_copy(spellName);
        if (spellName.empty()) {
            return {{client, "You need to specify the name of the spell to be cast.\n"}};
        }

        std::vector<Message> responses;
        bool isValidSpellName = static_cast<bool>(this->spellMap.count(spellNameLower));

        if (!isValidSpellName) {
            std::ostringstream responseMessage;
            responseMessage << "There are no spells with the name of \"" << spellName << "\".\n";

            return {{client, responseMessage.str()}};
        }

        auto spell = this->spellMap.at(spellNameLower);
        switch (spell) {
            case (Spell::BodySwap):
                responses = this->bodySwap(client, targetName);
                break;

            case (Spell::Decoy):
                responses = this->decoy(client);
                break;

            case (Spell::Confuse):
                responses = this->confuse(client, targetName);
                break;

            case (Spell::Heal):
                responses = this->heal(client, targetName);
                break;

            default:
                break;
        }

        return responses;
    }

    std::vector<Message>
    MagicHandler::bodySwap(const Connection &client, const std::string &targetName) {
        std::vector<Message> responses;
        std::ostringstream casterMessage;

        if (this->isBodySwapped(client)) {
            casterMessage << "You can't cast " << BODY_SWAP_SPELL_NAME << " while already under the effects of the spell!\n";
            return {{client, casterMessage.str()}};
        }

        if (targetName.empty()) {
            casterMessage << "You need to specify the name of the person to cast " << BODY_SWAP_SPELL_NAME << " on.\n";
            return {{client, casterMessage.str()}};
        }

        auto casterPlayerId = this->accountHandler.getPlayerIdByClient(client);
        auto casterRoomId = this->accountHandler.getRoomIdByClient(client);
        auto casterUsername = this->accountHandler.getUsernameByClient(client);
        if (casterUsername == targetName) {
            casterMessage << "You can't cast " << BODY_SWAP_SPELL_NAME << " on yourself!\n";
            return {{client, casterMessage.str()}};
        }

        auto targetClient = this->accountHandler.getClientByUsername(targetName);
        if (targetClient.id == AccountHandler::INVALID_ID) {
            casterMessage << "There is no one here with the name \"" << targetName << "\".\n";
            return {{client, casterMessage.str()}};
        }

        auto targetRoomId = this->accountHandler.getRoomIdByClient(targetClient);
        if (targetRoomId != casterRoomId) {
            casterMessage << "There is no one here with the name \"" << targetName << "\".\n";
            return {{client, casterMessage.str()}};
        }

        if (this->isBodySwapped(targetClient)) {
            casterMessage << targetName << " is already under the effects of the " << BODY_SWAP_SPELL_NAME << " spell!\n";
            return {{client, casterMessage.str()}};
        }

        auto targetPlayerId = this->accountHandler.getPlayerIdByClient(targetClient);
        this->accountHandler.swapPlayerClientsByPlayerId(casterPlayerId, targetPlayerId);
        this->bodySwapInstances.push_back({casterPlayerId, targetPlayerId, BODY_SWAP_DURATION});

        casterMessage << "You have successfully swapped bodies with " << targetName << ".\n";

        std::ostringstream targetMessage;
        targetMessage << casterUsername << " casts " << BODY_SWAP_SPELL_NAME << " on you!\n";

        responses.push_back({targetClient, targetMessage.str()});
        responses.push_back({client, casterMessage.str()});

        return responses;
    }

    std::vector<Message>
    MagicHandler::decoy(const Connection &client) {
        std::vector<Message> responses;
        std::ostringstream message;
        auto &player = this->accountHandler.getPlayerByClient(client);

        if (!this->combatHandler.isInCombat(player.getId())) {
            message << "You can only cast " << DECOY_SPELL_NAME << " while in combat.\n";
            return {{client, message.str()}};
        }

        if (this->combatHandler.isInCombat(-player.getId())) {
            message << "You can only have one active decoy at a time!\n";
            return {{client, message.str()}};
        }

        this->combatHandler.replaceWithDecoy(player);

        message << "You create a decoy of yourself and escape from combat.\n";
        responses.push_back({client, message.str()});

        return responses;
    }

    std::vector<Message>
    MagicHandler::confuse(const Connection &client, const std::string &targetName) {
        std::vector<Message> responses;
        std::ostringstream casterMessage;

        if (targetName.empty()) {
            casterMessage << "You need to specify the name of the person to cast " << CONFUSE_SPELL_NAME << " on.\n";
            return {{client, casterMessage.str()}};
        }

        auto casterPlayerId = this->accountHandler.getPlayerIdByClient(client);
        auto casterRoomId = this->accountHandler.getRoomIdByClient(client);
        auto casterUsername = this->accountHandler.getUsernameByClient(client);
        if (casterUsername == targetName) {
            if (this->isConfused(client)) {
                casterMessage << "You are already under the effects of the " << CONFUSE_SPELL_NAME << " spell!\n";
                return {{client, casterMessage.str()}};
            }

            this->confuseInstances.push_back({casterPlayerId, casterPlayerId, CONFUSE_DURATION});

            casterMessage << "You cast " << CONFUSE_SPELL_NAME << " on yourself.\n";
            return {{client, casterMessage.str()}};
        }

        auto targetClient = this->accountHandler.getClientByUsername(targetName);
        if (targetClient.id == AccountHandler::INVALID_ID) {
            casterMessage << "There is no player here with the name \"" << targetName << "\".\n";
            return {{client, casterMessage.str()}};
        }

        auto targetPlayerId = this->accountHandler.getPlayerIdByClient(targetClient);
        auto targetRoomId = this->accountHandler.getRoomIdByClient(targetClient);
        if (casterRoomId != targetRoomId) {
            casterMessage << "There is no player here with the name \"" << targetName << "\".\n";
            return {{client, casterMessage.str()}};
        }

        if (isConfused(targetClient)) {
            casterMessage << targetName << " is already under the effects of the " << CONFUSE_SPELL_NAME << " spell!\n";
            return {{client, casterMessage.str()}};
        }

        this->confuseInstances.push_back({casterPlayerId, targetPlayerId, CONFUSE_DURATION});

        casterMessage << "You cast " << CONFUSE_SPELL_NAME << " on " << targetName << ".\n";

        std::ostringstream targetMessage;
        targetMessage << casterUsername << " casts " << CONFUSE_SPELL_NAME << " on you!\n";

        responses.push_back({client, casterMessage.str()});
        responses.push_back({targetClient, targetMessage.str()});

        return responses;
    }

    std::vector<Message>
    MagicHandler::heal(const Connection &client, const std::string &targetName) {
        std::vector<Message> responses;
        std::ostringstream casterMessage;

        auto &player = accountHandler.getPlayerByClient(client);
        auto casterUsername = this->accountHandler.getUsernameByClient(client);

        if (combatHandler.isInCombat(player.getId())) {
            casterMessage << "You can't cast " << HEAL_SPELL_NAME << " while in combat.\n";
            return {{client, casterMessage.str()}};
        }

        if (targetName.empty() || (casterUsername == targetName)) {
            player.setHealth(Character::STARTING_HEALTH);

            casterMessage << "You cast " << HEAL_SPELL_NAME << " on yourself.\n";
            return {{client, casterMessage.str()}};
        }

        auto targetClient = this->accountHandler.getClientByUsername(targetName);
        if (targetClient.id == AccountHandler::INVALID_ID) {
            casterMessage << "There is no player here with the name \"" << targetName << "\".\n";
            return {{client, casterMessage.str()}};
        }

        auto targetRoomId = this->accountHandler.getRoomIdByClient(targetClient);
        auto casterRoomId = this->accountHandler.getRoomIdByClient(client);
        if (casterRoomId != targetRoomId) {
            casterMessage << "There is no player here with the name \"" << targetName << "\".\n";
            return {{client, casterMessage.str()}};
        }

        auto &targetPlayer = this->accountHandler.getPlayerByClient(targetClient);
        if (combatHandler.isInCombat(targetPlayer.getId())) {
            casterMessage << "You can't cast " << HEAL_SPELL_NAME << " on "
                          << targetPlayer.getUsername() << " while they are in combat.\n";
            return {{client, casterMessage.str()}};
        }

        targetPlayer.setHealth(Character::STARTING_HEALTH);

        casterMessage << "You cast " << HEAL_SPELL_NAME << " on " << targetName << ".\n";

        std::ostringstream targetMessage;
        targetMessage << casterUsername << " casts " << HEAL_SPELL_NAME << " on you!\n";

        responses.push_back({client, casterMessage.str()});
        responses.push_back({targetClient, targetMessage.str()});

        return responses;
    }

    void
    MagicHandler::removeBodySwap(const model::ID &playerId) {
        auto it_swap = std::find_if(
            this->bodySwapInstances.begin(),
            this->bodySwapInstances.end(),
            [&playerId](const SpellInstance &swapInstance) {
                return (playerId == swapInstance.casterPlayerId) || (playerId == swapInstance.targetPlayerId);
            }
        );

        if (it_swap != this->bodySwapInstances.end()) {
            auto casterPlayerId= it_swap->casterPlayerId;
            auto targetPlayerId = it_swap->targetPlayerId;

            this->accountHandler.swapPlayerClientsByPlayerId(casterPlayerId, targetPlayerId);
            this->bodySwapInstances.erase(it_swap);
        }
    }

    void
    MagicHandler::removeConfuse(const model::ID &playerId) {
        auto it_confuse = std::find_if(
            this->confuseInstances.begin(),
            this->confuseInstances.end(),
            [&playerId](const SpellInstance &confuseInstance) {
                return (playerId == confuseInstance.casterPlayerId) || (playerId == confuseInstance.targetPlayerId);
            }
        );

        if (it_confuse != this->confuseInstances.end()) {
            this->confuseInstances.erase(it_confuse);
        }
    }

    void
    MagicHandler::processBodySwapInstancesCycle(std::deque<Message> &messages) {
        auto swapInstance = this->bodySwapInstances.begin();

        while (swapInstance != this->bodySwapInstances.end()) {
            if (swapInstance->cyclesRemaining > 0) {
                --swapInstance->cyclesRemaining;
                ++swapInstance;

            } else {
                auto casterPlayerId = swapInstance->casterPlayerId;
                auto targetPlayerId = swapInstance->targetPlayerId;

                this->accountHandler.swapPlayerClientsByPlayerId(casterPlayerId, targetPlayerId);

                auto casterClient = this->accountHandler.getClientByPlayerId(casterPlayerId);
                auto targetClient = this->accountHandler.getClientByPlayerId(targetPlayerId);

                std::ostringstream message;
                message << "The effects of " << BODY_SWAP_SPELL_NAME << " have worn off and you return to your original body.\n";
                messages.push_back({casterClient, message.str()});
                messages.push_back({targetClient, message.str()});

                this->bodySwapInstances.erase(swapInstance);
            }
        }
    }

    void
    MagicHandler::processConfuseInstancesCycle(std::deque<Message> &messages) {
        auto confuseInstance = this->confuseInstances.begin();

        while (confuseInstance != this->confuseInstances.end()) {
            if (confuseInstance->cyclesRemaining > 0) {
                --confuseInstance->cyclesRemaining;
                ++confuseInstance;

            } else {
                auto targetClient = this->accountHandler.getClientByPlayerId(confuseInstance->targetPlayerId);

                std::ostringstream message;
                message << "The effects of " << CONFUSE_SPELL_NAME << " have worn off and your speech returns to normal.\n";
                messages.push_back({targetClient, message.str()});

                this->confuseInstances.erase(confuseInstance);
            }
        }
    }

    bool
    MagicHandler::isBodySwapped(const Connection &client) {
        auto playerId = this->accountHandler.getPlayerIdByClient(client);

        auto it_swap = std::find_if(
            this->bodySwapInstances.begin(),
            this->bodySwapInstances.end(),
            [&playerId](const SpellInstance &swapInstance) {
                return (playerId == swapInstance.casterPlayerId) || (playerId == swapInstance.targetPlayerId);
            }
        );

        return it_swap != this->bodySwapInstances.end();
    }

    bool
    MagicHandler::isConfused(const Connection &client) {
        auto playerId = this->accountHandler.getPlayerIdByClient(client);

        auto it_confuse = std::find_if(
            this->confuseInstances.begin(),
            this->confuseInstances.end(),
            [&playerId](const SpellInstance &confuseInstance) {
                return playerId == confuseInstance.targetPlayerId;
            }
        );

        return it_confuse != this->confuseInstances.end();
    }

    void
    MagicHandler::confuseSpeech(std::string &message) {
        boost::find_format_all(
            message,
            boost::token_finder((boost::is_alpha() || boost::is_any_of("'")), boost::token_compress_on),
            PigLatinFormatter<std::string>()
        );
    }

    void
    MagicHandler::handleLogout(const Connection &client) {
        auto playerId = this->accountHandler.getPlayerIdByClient(client);

        if (this->isBodySwapped(client)) {
            this->removeBodySwap(playerId);
        }

        if (this->isConfused(client)) {
            this->removeConfuse(playerId);
        }
    }

    void
    MagicHandler::processCycle(std::deque<Message> &messages) {
       this->processBodySwapInstancesCycle(messages);
       this->processConfuseInstancesCycle(messages);
    }
}


text:
#include "tkdCmdParser.h"

#include "vnl/vnl_vector.h"
#include "vnl/vnl_matrix.h"
#include "vnl/vnl_matrix_fixed.h"
#include <vnl/algo/vnl_cholesky.h>
#include <vnl/vnl_transpose.h>

#include <cmath>
#include <ctime>

#include <list>

#include <itkImage.h>
#include <itkImageFileReader.h>
#include <itkImageFileWriter.h>
#include <itkImageRegionIteratorWithIndex.h>
#include <itkImageLinearIteratorWithIndex.h>
#include <itkImageRegionIterator.h>
#include <itkDiffusionTensor3D.h>

#include <boost/tokenizer.hpp>
#include <boost/lexical_cast.hpp>
#include <boost/random/uniform_int.hpp>
#include <boost/random/uniform_real.hpp>
#include <boost/random/variate_generator.hpp>
#include <boost/random.hpp>
#include <boost/math/distributions/gamma.hpp>

/**
 * dki namespace.
 */
namespace dki
{

	struct parameters
	{
		std::string inputFileName;
		std::string maskFileName;
		std::string bvecsFileName;
		std::string bvalsFileName;
		std::string outputFileName;
		std::string sep;
		int iterations;
		bool regularize;
	};

	typedef double PixelType;
	typedef unsigned char MaskPixelType;

	typedef itk::Image< PixelType, 4 > ImageType;
	typedef itk::Image< MaskPixelType, 3 > MaskImageType;

	typedef itk::Image< PixelType, 3 > OutputImageType;

	typedef itk::DiffusionTensor3D< PixelType > DTITensorType;
	typedef vnl_vector_fixed< PixelType, 15 > DKITensorType;

	typedef itk::ImageFileReader< ImageType > ReaderType;
	typedef itk::ImageFileReader< MaskImageType > MaskReaderType;

	typedef vnl_vector< PixelType > VectorType;
	typedef vnl_matrix< PixelType > MatrixType;

	typedef itk::ImageLinearConstIteratorWithIndex< ImageType > ConstIterator4DType;
	typedef itk::ImageRegionConstIteratorWithIndex< MaskImageType > ConstIterator3DType;
	typedef itk::ImageRegionIteratorWithIndex< OutputImageType > Iterator3DType;

	typedef itk::ImageFileWriter< OutputImageType > WriterType;

	typedef boost::tokenizer< boost::char_separator< char > > TokType;
	typedef boost::mt19937 random_number_type;

	/**
	 * Return 2-rand DTI tensor.
	 *
	 * ITK tensor layout =>
	 * 					| 0  1  2  |
	 *			        | X  3  4  |
	 *       			| X  X  5  |
	 */
	static DTITensorType GetDTITensor( const VectorType& theta )
	{
		DTITensorType D( 0.0 );

		D( 0, 0 ) = theta( 1 );
		D( 0, 1 ) = theta( 4 );
		D( 0, 2 ) = theta( 6 );
		D( 1, 1 ) = theta( 2 );
		D( 1, 2 ) = theta( 5 );
		D( 2, 2 ) = theta( 3 );

		return D;
	}

	/**
	 * Return 4-rand DKI tensor (for now a fixed vector with 15 elements).
	 *
	 * theta_dki = [W_1111, W_2222, W_3333,
	 * 				W_1112, W_2223, W_2221,
	 * 				W_2223, W_3331, W_3332,
	 * 				W_1122, W_1133, W_2233,
	 * 				W_1123, W_1223, W_1233]
	 */
	DKITensorType GetDKITensor( const VectorType& theta )
	{
		DKITensorType tensor;
		for( unsigned int i = 0; i < 15; i++ )
			tensor( i ) = theta( 7 + i );

		return tensor;
	}

	/**
	 * Predicate to remove non positive semi-definite DTI/DKI tensors.
	 */
	class non_positive_definite
	{
		public:

		/**
		 * Check for non positive semi-definite tensor.
		 */
		bool operator() ( const VectorType& beta, bool dki_tensor = false )
		{
			DTITensorType DT = GetDTITensor( beta );

			DTITensorType::EigenValuesArrayType eigenValues;
			DT.ComputeEigenValues( eigenValues );

			bool non_negative = false;

			for( unsigned int i = 0; i < 3; i++ )
				if( eigenValues[ i ] < 0 )
					non_negative = true;

			if( dki_tensor )
			{
				DKITensorType DK = GetDKITensor( beta );
				// TODO calculate if DK tensor is postive semi-definite ...
			}

			return non_negative;
		}
	};

	/**
	 * Date, 06-01-2011.
	 *
	 * Implement Gibbs sampler for linear DTI fits.
	 *
	 * Ref: "Introduction to Applied Bayesian Statistics and Estimation for Social Scientists, Scott M. Lynch, 2007. Springer"
	 */
	class DKIFit
    {
        public:

		ImageType::Pointer Input;
		MaskImageType::Pointer Mask;

		MatrixType X_dti; // dti design matrix
		MatrixType X_dki; // dki design matrix
		MatrixType Z; // covariance matrix
		VectorType B; // bvals
		MatrixType R; // bvecs

		OutputImageType::Pointer S0_dti;
		OutputImageType::Pointer S0_dki;

		OutputImageType::Pointer FA_dti;
		OutputImageType::Pointer FA_dki;

		OutputImageType::Pointer Trace_dti;
		OutputImageType::Pointer Trace_dki;

		OutputImageType::Pointer AKC;
		OutputImageType::Pointer Kmax;
		OutputImageType::Pointer Kmin;

		/**
		 * Start DKI fit.
		 */
		void Run( const parameters& args )
        {
			SetImage( args.inputFileName );
			SetMask( args.maskFileName, Input.GetPointer() );
			SetBVals( args.bvalsFileName, args.sep );
			SetBVecs( args.bvecsFileName, args.sep );
			SetDTIDesignMatrix();
			SetDKIDesignMatrix();
			AllocateOutput( Input.GetPointer() );
			Fit( args.iterations, args.regularize );
			Write( args.outputFileName );
        }

        protected:

		/**
		 * Return cholesky decomposition for symmetric matrix.
		 *
		 * Sqrt( A ) using cholesky decomposition is twice as fast as svd/qr.
		 *
		 * The cholesky decomposition decomposes symmetric A = LL'
		 * where L is lower triangular.
		 */
		MatrixType Cholesky( const MatrixType& A )
		{
			/* Test Cholesky ...
			MatrixType A( 3, 3 );
			A( 0, 0 ) = 2; A( 0, 1 ) = 0; A( 0, 2 ) = 0.1;
			A( 1, 0 ) = 0; A( 1, 1 ) = 2; A( 1, 2 ) = 0;
			A( 2, 0 ) = 0.1; A( 2, 1 ) = 0; A( 2, 2 ) = 4;
			MatrixType L = Cholesky( A );
			std::cout << "Cholesky: " << L * L.transpose() << std::endl;
			*/
			vnl_cholesky::Operation op = vnl_cholesky::quiet;

			vnl_cholesky chol( A, op );

			return chol.upper_triangle(); // A = U'U
		}

		/**
		 * Normal distribution sampler.
		 */
		VectorType SampleNormal( random_number_type& ran, unsigned int total, double mean, double sigma )
		{
			using namespace boost;

			// select Gaussian probability distribution
			normal_distribution< double > norm_dist( mean, sigma );

			// bind random number generator to distribution, forming a function
			variate_generator< random_number_type&, normal_distribution< double > > sampler( ran, norm_dist );

			VectorType samples( total );

			for( unsigned int i = 0; i < total; i++ )
				samples( i ) = sampler(); // sample from the distribution

			return samples;
		}

		/**
		 * Inverse gamma distribution sampler.
		 *
		 * Implemented as 1 / gamma, with boost library 1.45 switch to inverse_gamma_distribution()
		 *
		 * shape =~ alpha
		 * scale =~ beta
		 */
		VectorType SampleIGamma( random_number_type& ran, unsigned int total, double alpha, double beta )
		{
			using namespace boost;

		    gamma_distribution< double > gamma_dist( alpha );

		    variate_generator< random_number_type&, gamma_distribution< double > > sampler( ran, gamma_dist );

		    VectorType samples( total );

		    for( unsigned int i = 0; i < total; i++ )
		    	samples( i ) = 1 / ( beta * sampler() ); // sample gamma distribution with scale = 1, and invert

		    return samples;
		}

		/**
		 * Design matrix with ln(S0) as constant .
		 *
		 * See: "C.G. Koay et al / Journal of Magnetic Resonance 128 (2006) pag. 116"
		 *
		 * DTI tensor thus: ln(S0), Dxx, Dyy, Dzz, Dxy, Dyz, Dxz.
		 */
		void SetDTIDesignMatrix()
		{
			MatrixType A( R.rows(), 7, 0 ); // 5 x 7

			for ( unsigned int i = 0; i < R.rows(); i++ )
			{
				double b = B( i );
				double Gx = R( i, 0 );
				double Gy = R( i, 1 );
				double Gz = R( i, 2 );

				A( i, 0 ) = 1; // ln(S0 )

				A( i, 1 ) = -b * Gx * Gx;

				A( i, 2 ) = -b * Gy * Gy;

				A( i, 3 ) = -b * Gz * Gz;

				A( i, 4 ) = -2 * b * Gx * Gy;

				A( i, 5 ) = -2 * b * Gy * Gz;

				A( i, 6 ) = -2 * b * Gx * Gz;
			}

			X_dti = A;
		}

		/**
		 * Bayesian Linear least squares fit.
         *
		 * // TODO add covariance matrix W => to get linear 'weighted' regression.
		 * beta = (X' * X)^-1 * X' * ln( y )
		 */
		VectorType FitVoxel( const MatrixType& X, const VectorType& y, unsigned int iterations,
				unsigned int k, unsigned int n, bool regularize )
		{
			// (X' * X)^-1

			vnl_svd< PixelType > svd( X.transpose() * X );
			MatrixType XTXI = svd.inverse();

			// ln( y )

			VectorType Y( y.size() );
			for( unsigned int i = 0; i < y.size(); i++ )
				Y( i ) = log( y( i ) );

			// (X' * X)^-1 * X' * ln( y )

			VectorType beta = XTXI * X.transpose() * Y;

			// sample sigma from it's inverse gamma marginal distribution

			// Create a Mersenne twister random number generator
			random_number_type eng( static_cast< unsigned int > ( std::time( 0 ) ) );

			PixelType shape = 0.5 * ( n - k );
			PixelType scale = 0.5 * dot_product( ( Y - X * beta ), ( Y - X * beta ) );
			VectorType s2 = SampleIGamma( eng, iterations, shape, scale );

			// -----------------------------------
			// sample distribution for beta (MVN):
			// -----------------------------------
			// mean => inv(X'X) * (X'Y)				[fitted beta]
			// variance => sigma^2 * inv(X'X)		[sqrt using cholesky decomp.]

			std::list< VectorType > bs;
			for( unsigned int i = 0; i < iterations; i++ )
				bs.push_back( beta + SampleNormal( eng, k, 0, 1 ) * Cholesky( s2( i ) * XTXI ) );

			// remove non positive semi definite tensors
			if (regularize )
				bs.remove_if( non_positive_definite() );

			// if nothing left, return zero tensor ...

			if( bs.empty() )
			{
				//std::cout << "Voxel found with sampled tensors not positive semi-definite!" << std::endl;
				return VectorType( beta.size(), 0 ); // return empty tensor
			}

			// get mean tensor

			VectorType mean_beta( beta.size(), 0 );

			for( std::list< VectorType >::iterator it = bs.begin(); it != bs.end(); ++it )
				mean_beta += *it;

			return mean_beta /= bs.size();
		}

		/**
		 * Read bvals from file.
		 */
		void SetBVals( const std::string& bvalsFileName, const std::string& sep )
		{
			// init
			VectorType b( GetNumberOfRows( bvalsFileName ), 0 );

			// open data ...
			std::ifstream in( bvalsFileName.c_str() );
			if ( in.fail() )
			{
				std::cerr << "*** ERROR ***: could not read from: " << bvalsFileName << std::endl;
				exit( EXIT_FAILURE );
			}

			std::string line;
			unsigned int rowIndex = 0;

			while ( getline( in, line ) )
			{
				TokType tok( line, boost::char_separator< char >( sep.c_str() ) );
				try
				{
					for ( TokType::iterator id = tok.begin(); id != tok.end(); ++id )
					{
						b( rowIndex ) = boost::lexical_cast< PixelType >( *id );
					}
					rowIndex++;
				}
				catch ( boost::bad_lexical_cast& e )
				{
					std::cout << "*** WARNING ***: could not parse " << bvalsFileName << std::endl;
					e.what();
					exit( EXIT_FAILURE );
				}
			}

			in.close();

			B = b;
		}

		void SetBVecs( const std::string& bvecsFileName, const std::string& sep )
		{
			// init
			MatrixType r( GetNumberOfRows( bvecsFileName ), 3, 0 );

			// open data ...
			std::ifstream in( bvecsFileName.c_str() );
			if ( in.fail() )
			{
				std::cerr << "*** ERROR ***: could not read from: " << bvecsFileName << std::endl;
				exit( EXIT_FAILURE );
			}

			std::string line;
			unsigned int rowIndex = 0;

			while ( getline( in, line ) )
			{
				TokType tok( line, boost::char_separator< char >( sep.c_str() ) );
				try
				{
					unsigned int colIndex = 0;
					for ( TokType::iterator id = tok.begin(); id != tok.end(); ++id )
					{
						r( rowIndex, colIndex ) = boost::lexical_cast< PixelType >( *id );
						colIndex++;
					}
					rowIndex++;
				}
				catch ( boost::bad_lexical_cast& e )
				{
					std::cout << "*** WARNING ***: could not parse " << bvecsFileName << std::endl;
					std::cout << e.what() << std::endl;
					exit( EXIT_FAILURE );
				}
			}

			in.close();

			R = r;
		}

		/**
		 * Build X_dki from bvecs and bvals input files.
		 *
		 * DTI tensor elements:
		 * --------------------
		 * [D_11, D_12, D_13, D_22, D_23, D_33]
		 *
		 * DKI tensor elements:
		 * --------------------
		 * [W_1111, W_2222, W_3333,
		 * 	W_1112, W_2223, W_2221, W_2223, W_3331, W_3332,
		 *  W_1122, W_1133, W_2233,
		 * 	W_1123, W_1223, W_1233]
		 */
		void SetDKIDesignMatrix()
		{
			// check ...
			if( B.size() != R.rows() )
			{
				std::cerr << "*** ERROR ***: bvals and bvecs size does not match!" << std::endl;
				exit( EXIT_FAILURE );
			}

			// init ( e.g. [52 x 22] ) ...
			X_dki = MatrixType( R.rows(), 22, 0 );

			for ( unsigned int i = 0; i < R.rows(); i++ )
			{
				double b = B( i );
				double Gx = R( i, 0 );
				double Gy = R( i, 1 );
				double Gz = R( i, 2 );

				// ln( S0 )

				X_dki( i, 0 ) = 1;

				// DTI tensor

				X_dki( i, 1 ) = -b * Gx * Gx;
				X_dki( i, 2 ) = -b * Gy * Gy;
				X_dki( i, 3 ) = -b * Gz * Gz;
				X_dki( i, 4 ) = -2 * b * Gx * Gy;
				X_dki( i, 5 ) = -2 * b * Gy * Gz;
				X_dki( i, 6 ) = -2 * b * Gx * Gz;

				// DKI tensor

				PixelType c = ( b * b ) / 6;

				X_dki( i, 7 ) = c * Gx * Gx * Gx * Gx; // W_1111
				X_dki( i, 8 ) = c * Gy * Gy * Gy * Gy; // W_2222
				X_dki( i, 9 ) = c * Gz * Gz * Gz * Gz; // W_3333

				X_dki( i, 10 ) = 4 * c * Gx * Gx * Gx * Gy; // W_1112
				X_dki( i, 11 ) = 4 * c * Gx * Gx * Gx * Gy; // W_2223
				X_dki( i, 12 ) = 4 * c * Gx * Gy * Gy * Gy; // W_2221
				X_dki( i, 13 ) = 4 * c * Gy * Gy * Gy * Gz; // W_2223
				X_dki( i, 14 ) = 4 * c * Gx * Gz * Gz * Gz; // W_3331
				X_dki( i, 15 ) = 4 * c * Gy * Gz * Gz * Gz; // W_3332

				X_dki( i, 16 ) = 6 * c * Gx * Gx * Gy * Gy; // W_1122
				X_dki( i, 17 ) = 6 * c * Gx * Gx * Gz * Gz; // W_1133
				X_dki( i, 18 ) = 6 * c * Gy * Gy * Gz * Gz; // W_2233

				X_dki( i, 19 ) = 12 * c * Gx * Gx * Gy * Gz; // W_1123
				X_dki( i, 20 ) = 12 * c * Gx * Gy * Gy * Gz; // W_1223
				X_dki( i, 21 ) = 12 * c * Gx * Gy * Gz * Gz; // W_1233
			}
		}

		/**
		 * Return S0.
		 */
		PixelType GetS0( const VectorType& theta )
		{
			return std::exp( theta( 0 ) );
		}

		/**
		 * Return largest AKC value.
		 *
		 * See:
		 * "Principal invariant and inherent parameters of diffusion kurtosis tensors"
		 * "L. Qi et al / J. Math. Anal. Appl. 349 (2009) 165 - 180."
		 *
		 * return [mean(k), min(k), max(k)]
		 */
		VectorType GetKurtosisValues( const PixelType MD, const VectorType& D_app, const DKITensorType& W )
		{
			// get number of b-zeros ...
			unsigned int b0s = 0;
			for( unsigned int i = 0; i < B.size(); i++ )
				if( B( i ) == 0 )
					b0s++;

			VectorType K( R.rows() - b0s, 0 );

			unsigned int dwi_index = 0;

			for ( unsigned int i = 0; i < R.rows(); i++ )
			{
				if( B( i ) != 0 ) // only dwi ...
				{
					double Gx = R( i, 0 );
					double Gy = R( i, 1 );
					double Gz = R( i, 2 );

					PixelType tmp0 = MD * MD / ( D_app( i ) * D_app( i ) ); // D_app contains b0s ...

					PixelType tmp1 = W( 0 ) * Gx * Gx * Gx * Gx + W( 1 ) * Gy * Gy * Gy * Gy + W( 2 ) * Gz * Gz * Gz * Gz;

					PixelType tmp2 = 4 * ( W( 3 ) * Gx * Gx * Gx * Gy + W( 4 ) * Gx * Gx * Gx * Gy + W( 5 ) * Gx * Gy * Gy * Gy + W( 6 ) * Gy * Gy * Gy * Gz + W( 7 ) * Gx * Gz * Gz * Gz + W( 8 ) * Gy * Gz * Gz * Gz );

					PixelType tmp3 = 6 * ( W( 9 ) * Gx * Gx * Gy * Gy + W( 10 ) * Gx * Gx * Gz * Gz + W( 11 ) * Gy * Gy * Gz * Gz );

					PixelType tmp4 = 12 * ( W( 12 ) * Gx * Gx * Gy * Gz + W( 13 ) * Gx * Gy * Gy * Gz + W( 14 ) * Gx * Gy * Gz * Gz );

					K( dwi_index ) = tmp0 * ( tmp1 + tmp2 + tmp3 + tmp4 );

					dwi_index++;
				}
			}

			VectorType output( 3, 0 ); // average K, min K, max K

			output( 0 ) = K.mean();
			output( 1 ) = K.min_value();
			output( 2 ) = K.max_value();

			return output;
		}

		/**
		 * Return number of training data points.
		 */
		unsigned int GetNumberOfRows( const std::string& inputFile )
		{
			std::ifstream in( inputFile.c_str() );

			if ( in.fail() )
			{
				std::cerr << "*** ERROR ***: could not read from: " << inputFile << std::endl;
				exit( EXIT_FAILURE );
			}
			std::string line;
			unsigned int result = 0;

			while ( getline( in, line ) )
				result++;

			in.close();
			return result;
		}

		/**
		 * MCMC DTI fitting.
		 */
		void Fit( unsigned int iterations, bool regularize )
		{
			ConstIterator4DType it( Input, Input->GetLargestPossibleRegion() );
			ConstIterator3DType mit( Mask, Mask->GetLargestPossibleRegion() );

			Iterator3DType itS0_dti( S0_dti, S0_dti->GetLargestPossibleRegion() );
			Iterator3DType itS0_dki( S0_dki, S0_dki->GetLargestPossibleRegion() );

			Iterator3DType itFA_dti( FA_dti, FA_dti->GetLargestPossibleRegion() );
			Iterator3DType itFA_dki( FA_dki, FA_dki->GetLargestPossibleRegion() );

			Iterator3DType itTrace_dti( Trace_dti, Trace_dti->GetLargestPossibleRegion() );
			Iterator3DType itTrace_dki( Trace_dki, Trace_dki->GetLargestPossibleRegion() );

			Iterator3DType itAKC( AKC, AKC->GetLargestPossibleRegion() );
			Iterator3DType itKmin( Kmin, Kmin->GetLargestPossibleRegion() );
			Iterator3DType itKmax( Kmax, Kmax->GetLargestPossibleRegion() );

			it.SetDirection( 3 );
			it.GoToBegin();
			mit.GoToBegin();

			itS0_dti.GoToBegin();
			itS0_dki.GoToBegin();

			itFA_dti.GoToBegin();
			itFA_dki.GoToBegin();

			itTrace_dti.GoToBegin();
			itTrace_dki.GoToBegin();

			itAKC.GoToBegin();
			itKmin.GoToBegin();
			itKmax.GoToBegin();

			unsigned int totalDWI = ( Input->GetLargestPossibleRegion().GetSize() )[3];

			unsigned int sliceIndex = 0;

			// for each series of DWI ...
			while ( !it.IsAtEnd(), !mit.IsAtEnd() )
			{
				if ( mit.Get() != 0 )
				{
					// slice index
					if( sliceIndex != mit.GetIndex()[ 2 ] )
					{
						std::cout << "Processing slice: " << sliceIndex << std::endl;
						sliceIndex = mit.GetIndex()[ 2 ];
					}

					VectorType y( totalDWI );

					while ( !it.IsAtEndOfLine() )
					{
						y( ( it.GetIndex() )[3] ) = it.Get();
						++it;
					}

					// DTI => 7 elements (ln(S0) = 1, theta_D = 6)
					// DTI fit only ...

					unsigned int k = 7;
					unsigned int n = y.size();
					VectorType theta_dti = FitVoxel( X_dti, y, iterations, k, n, regularize );
					itS0_dti.Set( GetS0( theta_dti ) );
					DTITensorType tensor_simple = GetDTITensor( theta_dti );
					itFA_dti.Set( tensor_simple.GetFractionalAnisotropy() );
					itTrace_dti.Set( tensor_simple.GetTrace() );

					// DKI => 22 elements (ln(S0) = 1, theta_D = 6, theta_K = 15)
					// DTI and DKI fit in one...

					k += 15;

					VectorType theta_dki = FitVoxel( X_dki, y, iterations, k, n, regularize );
					itS0_dki.Set( GetS0( theta_dki ) );

					DTITensorType tensor_dti = GetDTITensor( theta_dki ); // DTI tensor from DKI fit ...
					itFA_dki.Set( tensor_dti.GetFractionalAnisotropy() );
					itTrace_dki.Set( tensor_dti.GetTrace() );

					DKITensorType tensor_dki = GetDKITensor( theta_dki );
					VectorType tmp = GetKurtosisValues( tensor_dti.GetTrace(), FitDapp( y ), tensor_dki );

					itAKC.Set( tmp( 0 ) );
					itKmin.Set( tmp( 1 ) );
					itKmax.Set( tmp( 2 ) );
				}

				it.NextLine();
				++mit;

				++itS0_dti;
				++itS0_dki;

				++itFA_dti;
				++itFA_dki;

				++itTrace_dti;
				++itTrace_dki;

				++itAKC;
				++itKmin;
				++itKmax;
			}
		}

		/**
		 * Return apparant diffusion per direction.
		 *
		 * ln( Sb / S0 ) / -b = D_app
		 */
		VectorType FitDapp( const VectorType& y )
		{
			// average S0 ...

			PixelType S0 = 0;
			PixelType totalS0;

			for( unsigned int i = 0; i < B.size(); i++ )
			{
				if( B( i ) == 0 )
				{
					S0 += y( i );
					totalS0++;
				}
			}
			if( totalS0 == 0 )
			{
				std::cerr << "*** ERROR ***: no b-zeros!" << std::endl;
				exit( EXIT_FAILURE );
			}
			S0 /= totalS0;

			// D_app ...

			VectorType D_app( R.rows(), 0 );

			for( unsigned int i = 0; i < B.size(); i++ )
			{
				if( B( i ) != 0 )
					D_app( i ) = log( y( i )  / S0 ) / - B( i );
			}

			return D_app;
		}

		/**
		 * Return covariate matrix Z.
		 */
		MatrixType BuildCovarianceMatrix( const VectorType& y )
		{
			MatrixType Z( y.size(), y.size(), 0 );

			// fill diagonal
			for( unsigned int i = 0; i < Z.rows(); i++ )
				Z( i, i ) = y( i ) * y( i );

			return Z;
		}

		/**
		 * Allocate all output images to 0.
		 */
		void AllocateOutput( ImageType::ConstPointer input )
		{
			MaskImageType::Pointer output = MaskImageType::New();

			OutputImageType::RegionType region3D;
			OutputImageType::IndexType index3D;
			OutputImageType::SizeType size3D;
			OutputImageType::SpacingType spacing3D;
			OutputImageType::PointType origin3D;

			ImageType::RegionType region = input->GetLargestPossibleRegion();
			ImageType::SizeType size = region.GetSize();
			ImageType::IndexType index = region.GetIndex();
			ImageType::SpacingType spacing = input->GetSpacing();
			ImageType::PointType origin = input->GetOrigin();

			size3D[0] = size[0];
			size3D[1] = size[1];
			size3D[2] = size[2];
			index3D[0] = index[0];
			index3D[1] = index[1];
			index3D[2] = index[2];
			origin3D[0] = origin[0];
			origin3D[1] = origin[1];
			origin3D[2] = origin[2];
			spacing3D[0] = spacing[0];
			spacing3D[1] = spacing[1];
			spacing3D[2] = spacing[2];

			region3D.SetSize( size3D );
			region3D.SetIndex( index3D );

			// create
			S0_dti = OutputImageType::New();
			S0_dki = OutputImageType::New();

			FA_dti = OutputImageType::New();
			FA_dki = OutputImageType::New();

			Trace_dti = OutputImageType::New();
			Trace_dki = OutputImageType::New();

			AKC = OutputImageType::New();
			Kmax = OutputImageType::New();
			Kmin = OutputImageType::New();

			// set
			S0_dti->SetRegions( region3D );
			S0_dti->SetSpacing( spacing3D );
			S0_dti->SetOrigin( origin3D );
			S0_dti->Allocate();
			S0_dti->FillBuffer( 0 );

			S0_dki->SetRegions( region3D );
			S0_dki->SetSpacing( spacing3D );
			S0_dki->SetOrigin( origin3D );
			S0_dki->Allocate();
			S0_dki->FillBuffer( 0 );

			FA_dti->SetRegions( region3D );
			FA_dti->SetSpacing( spacing3D );
			FA_dti->SetOrigin( origin3D );
			FA_dti->Allocate();
			FA_dti->FillBuffer( 0 );

			FA_dki->SetRegions( region3D );
			FA_dki->SetSpacing( spacing3D );
			FA_dki->SetOrigin( origin3D );
			FA_dki->Allocate();
			FA_dki->FillBuffer( 0 );

			Trace_dti->SetRegions( region3D );
			Trace_dti->SetSpacing( spacing3D );
			Trace_dti->SetOrigin( origin3D );
			Trace_dti->Allocate();
			Trace_dti->FillBuffer( 0 );

			Trace_dki->SetRegions( region3D );
			Trace_dki->SetSpacing( spacing3D );
			Trace_dki->SetOrigin( origin3D );
			Trace_dki->Allocate();
			Trace_dki->FillBuffer( 0 );

			AKC->SetRegions( region3D );
			AKC->SetSpacing( spacing3D );
			AKC->SetOrigin( origin3D );
			AKC->Allocate();
			AKC->FillBuffer( 0 );

			Kmax->SetRegions( region3D );
			Kmax->SetSpacing( spacing3D );
			Kmax->SetOrigin( origin3D );
			Kmax->Allocate();
			Kmax->FillBuffer( 0 );

			Kmin->SetRegions( region3D );
			Kmin->SetSpacing( spacing3D );
			Kmin->SetOrigin( origin3D );
			Kmin->Allocate();
			Kmin->FillBuffer( 0 );
		}

		/**
		 * Set input image.
		 */
		void SetImage( const std::string& inputFileName )
		{
			if ( !inputFileName.empty() )
			{
				ReaderType::Pointer reader = ReaderType::New();
				reader->SetFileName( inputFileName );
				reader->Update();
				Input = reader->GetOutput();
			}
			else
			{
				std::cerr << "Could not read input!" << std::endl;
				exit( EXIT_FAILURE );
			}
		}

		/**
		 * Return mask if file given, else create empty mask from input file.
		 */
		void SetMask( const std::string& maskFileName, ImageType::ConstPointer input )
		{
			if ( !maskFileName.empty() )
			{
				MaskReaderType::Pointer reader = MaskReaderType::New();
				reader->SetFileName( maskFileName );
				reader->Update();
				Mask = reader->GetOutput();
			}
			else
			{
				MaskImageType::Pointer output = MaskImageType::New();
				MaskImageType::RegionType region3D;
				MaskImageType::IndexType index3D;
				MaskImageType::SizeType size3D;
				MaskImageType::SpacingType spacing3D;
				MaskImageType::PointType origin3D;

				ImageType::RegionType region = input->GetLargestPossibleRegion();
				ImageType::SizeType size = region.GetSize();
				ImageType::IndexType index = region.GetIndex();
				ImageType::SpacingType spacing = input->GetSpacing();
				ImageType::PointType origin = input->GetOrigin();

				size3D[0] = size[0];
				size3D[1] = size[1];
				size3D[2] = size[2];
				index3D[0] = index[0];
				index3D[1] = index[1];
				index3D[2] = index[2];
				origin3D[0] = origin[0];
				origin3D[1] = origin[1];
				origin3D[2] = origin[2];
				spacing3D[0] = spacing[0];
				spacing3D[1] = spacing[1];
				spacing3D[2] = spacing[2];

				region3D.SetSize( size3D );
				region3D.SetIndex( index3D );

				// set
				output->SetRegions( region3D );
				output->SetSpacing( spacing3D );
				output->SetOrigin( origin3D );
				output->Allocate();
				output->FillBuffer( 1 );

				Mask = output;
			}
		}

		/**
		 * Write output images.
		 */
		void Write( const std::string& outputFileName )
		{
			// Writer
			WriterType::Pointer writer = WriterType::New();

			// S0
			writer->SetFileName( ( outputFileName + std::string( "_S0_dti.nii.gz" ) ).c_str() );
			writer->SetInput( S0_dti );
			writer->Update();

			writer->SetFileName( ( outputFileName + std::string( "_S0_dki.nii.gz" ) ).c_str() );
			writer->SetInput( S0_dki );
			writer->Update();

			// FA
			writer->SetFileName( ( outputFileName + std::string( "_FA_dti.nii.gz" ) ).c_str() );
			writer->SetInput( FA_dti );
			writer->Update();

			writer->SetFileName( ( outputFileName + std::string( "_FA_dki.nii.gz" ) ).c_str() );
			writer->SetInput( FA_dki );
			writer->Update();

			// Trace
			writer->SetFileName( ( outputFileName + std::string( "_Trace_dti.nii.gz" ) ).c_str() );
			writer->SetInput( Trace_dti );
			writer->Update();

			writer->SetFileName( ( outputFileName + std::string( "_Trace_dki.nii.gz" ) ).c_str() );
			writer->SetInput( Trace_dki );
			writer->Update();

			// AKC
			writer->SetFileName( ( outputFileName + std::string( "_AKC.nii.gz" ) ).c_str() );
			writer->SetInput( AKC );
			writer->Update();

			// Kmin
			writer->SetFileName( ( outputFileName + std::string( "_Kmin.nii.gz" ) ).c_str() );
			writer->SetInput( Kmin );
			writer->Update();

			// Kmax
			writer->SetFileName( ( outputFileName + std::string( "_Kmax.nii.gz" ) ).c_str() );
			writer->SetInput( Kmax );
			writer->Update();
		}

    }; // end class DKIFit

} // end namespace dki

/**
 * DKI fitting routine.
 */
int main( int argc, char ** argv )
{
    tkd::CmdParser p( argv[0], "Fit voxel-wise DTI using Gibbs sampling and postitive semi-definite tensor regularization." );

    dki::parameters args;
    args.sep = " ";
    args.iterations = 1000;
    args.regularize = false;

    p.AddArgument( args.inputFileName, "input" )
        ->AddAlias( "i" )
        ->SetDescription( "Input 4D image" )
        ->SetRequired( true );

    p.AddArgument( args.bvecsFileName, "bvecs" )
        ->AddAlias( "r" )
        ->SetDescription( "Diffusion gradient vectors (FSL format)" )
        ->SetRequired( true );

    p.AddArgument( args.bvalsFileName, "bvals" )
        ->AddAlias( "b" )
        ->SetDescription( "Diffusion b-value vector (FSL format)" )
        ->SetRequired( true );

    p.AddArgument( args.outputFileName, "output" )
        ->AddAlias( "o" )
        ->SetDescription( "Output filename base" )
        ->SetRequired( true );

    p.AddArgument( args.maskFileName, "mask" )
        ->AddAlias( "m" )
        ->SetDescription( "Mask 3D image" );

    p.AddArgument( args.sep, "separation" )
        ->AddAlias( "s" )
        ->SetDescription( "Separation string in bvecs file (default: ' ')" );

    p.AddArgument( args.iterations, "iterations" )
        ->AddAlias( "it" )
        ->SetDescription( "Number of MCMC sampling iterations per voxel (default: 1000)" );

    p.AddArgument( args.regularize, "regularize" )
        ->AddAlias( "reg" )
        ->SetDescription( "Regularize posteriors by removal of non postivive semi-definite tensors (default: false)" );

    if ( !p.Parse( argc, argv ) )
    {
        p.PrintUsage( std::cout );
        return EXIT_FAILURE;
    }

    dki::DKIFit fit;
    fit.Run( args );

    return EXIT_SUCCESS;
}

max_stars_repo_name: wmotte/toolkid

max_issues_repo_name: wmotte/toolkid

max_forks_repo_name: wmotte/toolkid


text:
/*
 * Author: Torarin Hals Bakke (2012)
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

#include "OMXPlayerSubtitles.h"
#include "OMXOverlayText.h"
#include "SubtitleRenderer.h"
#include "utils/Enforce.h"
#include "utils/ScopeExit.h"
#include "utils/Clamp.h"
#include "utils/log.h"

#include <boost/algorithm/string.hpp>
#include <utility>
#include <algorithm>
#include <typeinfo>
#include <cstdint>

using namespace std;
using namespace boost;

OMXPlayerSubtitles::OMXPlayerSubtitles() BOOST_NOEXCEPT
: m_visible(),
  m_use_external_subtitles(),
  m_active_index(),
  m_delay(),
  m_thread_stopped(),
  m_font_size(),
  m_centered(),
  m_ghost_box(),
  m_lines(),
  m_av_clock(),
#ifndef NDEBUG
  m_open()
#endif
{}

OMXPlayerSubtitles::~OMXPlayerSubtitles() BOOST_NOEXCEPT
{
  Close();
}

bool OMXPlayerSubtitles::Open(size_t stream_count,
                              vector<Subtitle>&& external_subtitles,
                              const string& font_path,
                              const string& italic_font_path,
                              float font_size,
                              bool centered,
                              bool ghost_box,
                              unsigned int lines,
                              int display, int layer,
                              OMXClock* clock) BOOST_NOEXCEPT
{
  assert(!m_open);

  m_subtitle_buffers.resize(stream_count, circular_buffer<Subtitle>(32));
  m_external_subtitles = std::move(external_subtitles);

  m_visible = true;
  m_use_external_subtitles = true;
  m_active_index = 0;
  m_delay = 0;
  m_thread_stopped.store(false, memory_order_relaxed);

  m_font_path = font_path;
  m_italic_font_path = italic_font_path;
  m_font_size = font_size;
  m_centered = centered;
  m_ghost_box = ghost_box;
  m_lines = lines;
  m_av_clock = clock;
  m_display = display;
  m_layer = layer;

  if(!Create())
    return false;

  SendToRenderer(Message::Flush{m_external_subtitles});

#ifndef NDEBUG
  m_open = true;
#endif

  return true;
}

void OMXPlayerSubtitles::Close() BOOST_NOEXCEPT
{
  if(Running())
  {
    SendToRenderer(Message::Stop{});
    StopThread();
  }

  m_mailbox.clear();
  m_subtitle_buffers.clear();

#ifndef NDEBUG
  m_open = false;
#endif
}

void OMXPlayerSubtitles::Process()
{
  try
  {
    RenderLoop(m_font_path, m_italic_font_path, m_font_size, m_centered,
               m_ghost_box, m_lines, m_av_clock);
  }
  catch(Enforce_error& e)
  {
    if(!e.user_friendly_what().empty())
      printf("Error: %s\n", e.user_friendly_what().c_str());
    CLog::Log(LOGERROR, "OMXPlayerSubtitles::RenderLoop threw %s (%s)",
              typeid(e).name(), e.what());
  }
  catch(std::exception& e)
  {
    CLog::Log(LOGERROR, "OMXPlayerSubtitles::RenderLoop threw %s (%s)",
              typeid(e).name(), e.what());
  }
  m_thread_stopped.store(true, memory_order_relaxed);
}

template <typename Iterator>
Iterator FindSubtitle(Iterator begin, Iterator end, int time)
{
  return upper_bound(begin, end, time,
    [](int a, const Subtitle& b) { return a < b.stop; });
}

void OMXPlayerSubtitles::
RenderLoop(const string& font_path,
           const string& italic_font_path,
           float font_size,
           bool centered,
           bool ghost_box,
           unsigned int lines,
           OMXClock* clock)
{
  SubtitleRenderer renderer(m_display, m_layer,
                            font_path,
                            italic_font_path,
                            font_size,
                            0.01f, 0.06f,
                            centered,
                            0xDD,
                            ghost_box ? 0x80 : 0,
                            lines);

  vector<Subtitle> subtitles;

  int prev_now{};
  size_t next_index{};
  bool exit{};
  bool paused{};
  bool have_next{};
  int current_stop = INT_MIN;
  bool showing{};
  bool osd{};
  chrono::time_point<std::chrono::steady_clock> osd_stop;
  int delay{};

  auto GetCurrentTime = [&]
  {
    return static_cast<int>(clock->OMXMediaTime()/1000) - delay;
  };

  auto TryPrepare = [&](int time)
  {
    for(; next_index != subtitles.size(); ++next_index)
    {
      if(subtitles[next_index].stop > time)
      {
        renderer.prepare(subtitles[next_index].text_lines);
        have_next = true;
        break;
      }
    }
  };

  auto Reset = [&](int time)
  {
    renderer.unprepare();
    current_stop = INT_MIN;

    auto it = FindSubtitle(subtitles.begin(),
                           subtitles.end(),
                           time);
    next_index = it - subtitles.begin();

    if(next_index != subtitles.size())
    {
      renderer.prepare(subtitles[next_index].text_lines);
      have_next = true;
    }
    else
    {
      have_next = false;
    }
  };

  for(;;)
  {
    int timeout = INT_MAX;

    if(!paused)
    {
      auto now = GetCurrentTime();

      int till_stop =
        showing ? current_stop - now
                : INT_MAX;

      int till_next_start =
        have_next ? subtitles[next_index].start - now
                  : INT_MAX;

      timeout = min(min(till_stop, till_next_start), 1000);
    }

    if(osd)
    {
      procrustes(timeout,
        chrono::duration_cast<std::chrono::milliseconds>(
          osd_stop - chrono::steady_clock::now()).count());
    }

    m_mailbox.receive_wait(chrono::milliseconds(timeout),
      [&](Message::Push&& args)
      {
        subtitles.push_back(std::move(args.subtitle));
      },
      [&](Message::Flush&& args)
      {
        subtitles = std::move(args.subtitles);
        prev_now = INT_MAX;
      },
      [&](Message::Touch&&)
      {
      },
      [&](Message::SetPaused&& args)
      {
        paused = args.value;
      },
      [&](Message::SetDelay&& args)
      {
        delay = args.value;
        prev_now = INT_MAX;
      },
      [&](Message::Stop&&)
      {
        exit = true;
      },
      [&](Message::DisplayText&& args)
      {
        renderer.prepare(args.text_lines);
        renderer.show_next();
        showing = true;
        osd = true;
        osd_stop = chrono::steady_clock::now() +
                   chrono::milliseconds(args.duration);
        prev_now = INT_MAX;
      });

    if(exit) break;

    auto now = GetCurrentTime();

    if(now < prev_now || (have_next && subtitles[next_index].stop <= now))
    {
      Reset(now);
    }
    else if(!have_next)
    {
      TryPrepare(now);
    }

    prev_now = now;

    if(osd && chrono::steady_clock::now() >= osd_stop)
      osd = false;

    if(!osd && current_stop <= now)
    {
      if(have_next && subtitles[next_index].start <= now)
      {
        renderer.show_next();
        // printf("show error: %i ms\n", now - subtitles[next_index].start);
        showing = true;
        current_stop = subtitles[next_index].stop;

        ++next_index;
        have_next = false;
        TryPrepare(now);
      }
      else if(showing)
      {
        renderer.hide();
        // printf("hide error: %i ms\n", now - current_stop);
        showing = false;
      }
    }
  }
}

void OMXPlayerSubtitles::FlushRenderer()
{
  assert(GetVisible());

  if(GetUseExternalSubtitles())
  {
    SendToRenderer(Message::Flush{m_external_subtitles});
  }
  else
  {
    Message::Flush flush;
    assert(!m_subtitle_buffers.empty());
    for(auto& s : m_subtitle_buffers[m_active_index])
      flush.subtitles.push_back(s);
    SendToRenderer(std::move(flush));
  }
}

void OMXPlayerSubtitles::Flush() BOOST_NOEXCEPT
{
  assert(m_open);

  for(auto& q : m_subtitle_buffers)
    q.clear();

  if(GetVisible())
  {
    if(GetUseExternalSubtitles())
      SendToRenderer(Message::Touch{});
    else
      SendToRenderer(Message::Flush{});
  }
}

void OMXPlayerSubtitles::Resume() BOOST_NOEXCEPT
{
  assert(m_open);
  SendToRenderer(Message::SetPaused{false});
}

void OMXPlayerSubtitles::Pause() BOOST_NOEXCEPT
{
  assert(m_open);
  SendToRenderer(Message::SetPaused{true});
}

void OMXPlayerSubtitles::SetUseExternalSubtitles(bool use) BOOST_NOEXCEPT
{
  assert(m_open);
  assert(use || !m_subtitle_buffers.empty());

  m_use_external_subtitles = use;
  if(GetVisible())
    FlushRenderer();
}

void OMXPlayerSubtitles::SetDelay(int value) BOOST_NOEXCEPT
{
  assert(m_open);

  m_delay = value;
  SendToRenderer(Message::SetDelay{value});
}

void OMXPlayerSubtitles::SetVisible(bool visible) BOOST_NOEXCEPT
{
  assert(m_open);

  if(visible)
  {
    if (!m_visible)
    {
      m_visible = true;
      FlushRenderer();
    }
  }
  else
  {
    if(m_visible)
    {
      m_visible = false;
      SendToRenderer(Message::Flush{});
    }
  }
}

void OMXPlayerSubtitles::SetActiveStream(size_t index) BOOST_NOEXCEPT
{
  assert(m_open);
  assert(index < m_subtitle_buffers.size());

  m_active_index = index;
  if(!GetUseExternalSubtitles() && GetVisible())
    FlushRenderer();
}

vector<string> OMXPlayerSubtitles::GetTextLines(OMXPacket *pkt)
{
  assert(pkt);

  m_subtitle_codec.Open(pkt->hints);

  auto result = m_subtitle_codec.Decode(pkt->data, pkt->size, 0, 0);
  assert(result == OC_OVERLAY);

  auto overlay = m_subtitle_codec.GetOverlay();
  assert(overlay);

  vector<string> text_lines;

  auto e = ((COMXOverlayText*) overlay)->m_pHead;
  if(e && e->IsElementType(COMXOverlayText::ELEMENT_TYPE_TEXT))
  {
    split(text_lines,
          ((COMXOverlayText::CElementText*) e)->m_text,
          is_any_of("\r\n"),
          token_compress_on);
  }

  return text_lines;
}

bool OMXPlayerSubtitles::AddPacket(OMXPacket *pkt, size_t stream_index) BOOST_NOEXCEPT
{
  assert(m_open);
  assert(stream_index < m_subtitle_buffers.size());

  if(!pkt)
    return false;

  SCOPE_EXIT
  {
    OMXReader::FreePacket(pkt);
  };

  if(pkt->hints.codec != AV_CODEC_ID_SUBRIP &&
     pkt->hints.codec != AV_CODEC_ID_SSA &&
     pkt->hints.codec != AV_CODEC_ID_ASS)
  {
    return true;
  }

  auto start = static_cast<int>(pkt->pts/1000);
  auto stop = start + static_cast<int>(pkt->duration/1000);
  auto text_lines = GetTextLines(pkt);

  if (!m_subtitle_buffers[stream_index].empty() &&
    stop < m_subtitle_buffers[stream_index].back().stop)
  {
    stop = m_subtitle_buffers[stream_index].back().stop;
  }

  m_subtitle_buffers[stream_index].push_back(
    Subtitle(start, stop, vector<string>()));
  m_subtitle_buffers[stream_index].back().text_lines = text_lines;

  if(!GetUseExternalSubtitles() &&
     GetVisible() &&
     stream_index == GetActiveStream())
  {
    SendToRenderer(Message::Push{{start, stop, std::move(text_lines)}});
  }

  return true;
}

void OMXPlayerSubtitles::DisplayText(const std::string& text, int duration) BOOST_NOEXCEPT
{
  assert(m_open);

  vector<string> text_lines;
  split(text_lines, text, is_any_of("\n"));
  SendToRenderer(Message::DisplayText{std::move(text_lines), duration});
}