text:
// This file is part of KWIVER, and is distributed under the
// OSI-approved BSD 3-Clause License. See top-level LICENSE file or
// https://github.com/Kitware/kwiver/blob/master/LICENSE for details.

/**
 * \file
 * \brief core essential matrix template implementations
 */

#include "essential_matrix.h"

#include <cmath>

#include <vital/exceptions/math.h>

#include <Eigen/SVD>

namespace kwiver {
namespace vital {

/// Compute the twisted pair rotation from the rotation and translation
rotation_d
essential_matrix
::twisted_rotation() const
{
  // The quaternion representation of a 180 degree rotation about
  // unit vector [X,Y,Z] is simply [X, Y, Z, 0]
  vector_3d t = this->translation();
  return rotation_d(vector_4d(t.x(), t.y(), t.z(), 0.0)) * this->rotation();
}

/// Construct from a provided matrix
template <typename T>
essential_matrix_<T>
::essential_matrix_( Eigen::Matrix<T,3,3> const &mat )
{
  const matrix_t W = (matrix_t() << T(0), T(-1), T(0),
                                    T(1), T( 0), T(0),
                                    T(0), T( 0), T(1)).finished();
  Eigen::JacobiSVD<matrix_t> svd(mat, Eigen::ComputeFullU |
                                      Eigen::ComputeFullV);
  const matrix_t& U = svd.matrixU();
  const matrix_t& V = svd.matrixV();
  trans_ = U.col(2);
  matrix_t R = U*W*V.transpose();
  if( R.determinant() < T(0) )
  {
    R *= T(-1);
  }
  rot_ = rotation_<T>(R);
}

/// Construct from a rotation and translation
template <typename T>
essential_matrix_<T>
::essential_matrix_( rotation_<T> const &rot,
                     vector_t const &trans )
  : rot_( rot ),
    trans_( trans.normalized() )
{
}

/// Conversion Copy constructor -- float specialization
template <>
template <>
essential_matrix_<float>
::essential_matrix_( essential_matrix_<float> const &other )
  : rot_( other.rot_ ),
    trans_( other.trans_ )
{
}

/// Conversion Copy constructor -- double specialization
template <>
template <>
essential_matrix_<double>
::essential_matrix_( essential_matrix_<double> const &other )
  : rot_( other.rot_ ),
    trans_( other.trans_ )
{
}

/// Construct from a generic essential_matrix
template <typename T>
essential_matrix_<T>
::essential_matrix_( essential_matrix const &base )
  : rot_( static_cast<rotation_<T> >(base.rotation()) ),
    trans_( base.translation().template cast<T>() )
{
}

/// Construct from a generic essential_matrix -- double specialization
template <>
essential_matrix_<double>
::essential_matrix_( essential_matrix const &base )
  : rot_( base.rotation() ),
    trans_( base.translation() )
{
}

/// Create a clone of outself as a shared pointer
template <typename T>
essential_matrix_sptr
essential_matrix_<T>
::clone() const
{
  return essential_matrix_sptr( new essential_matrix_<T>( *this ) );
}

/// Get a double-typed copy of the underlying matrix
template <typename T>
Eigen::Matrix<double,3,3>
essential_matrix_<T>
::matrix() const
{
  return this->compute_matrix().template cast<double>();
}

/// Specialization for matrices with native double type
template <>
Eigen::Matrix<double,3,3>
essential_matrix_<double>
::matrix() const
{
  return this->compute_matrix();
}

/// Return the one of two possible 3D rotations that can parameterize E
template <typename T>
rotation_d
essential_matrix_<T>
::rotation() const
{
  return static_cast<rotation_d>(this->rot_);
}

/// Return the second possible rotation that can parameterize E
template <typename T>
rotation_d
essential_matrix_<T>
::twisted_rotation() const
{
  return static_cast<rotation_d>(this->compute_twisted_rotation());
}

/// Return a unit translation vector (up to a sign) that parameterizes E
template <typename T>
vector_3d
essential_matrix_<T>
::translation() const
{
  return this->trans_.template cast<double>();
}

/// Get the underlying matrix
template <typename T>
typename essential_matrix_<T>::matrix_t
essential_matrix_<T>
::compute_matrix() const
{
  matrix_t t_cross;
  t_cross << T(0), -trans_[2], trans_[1],
             trans_[2], T(0), -trans_[0],
            -trans_[1], trans_[0], T(0);
  return t_cross * matrix_t(rot_.matrix());
}

/// Compute the twisted pair rotation from the rotation and translation
template <typename T>
rotation_<T>
essential_matrix_<T>
::compute_twisted_rotation() const
{
  typedef Eigen::Matrix<T,4,1> vector_4;
  // The quaternion representation of a 180 degree rotation about
  // unit vector [X,Y,Z] is simply [X, Y, Z, 0]
  const vector_t& t = trans_;
  return rotation_<T>(vector_4(t.x(), t.y(), t.z(), T(0))) * rot_;
}

/// Get a const reference to the underlying rotation
template <typename T>
rotation_<T> const&
essential_matrix_<T>
::get_rotation() const
{
  return rot_;
}

/// Get a const reference to the underlying translation
template <typename T>
typename essential_matrix_<T>::vector_t const&
essential_matrix_<T>
::get_translation() const
{
  return trans_;
}

// ===========================================================================
// Other Functions
// ---------------------------------------------------------------------------

/// essential_matrix_<T> output stream operator
template <typename T>
std::ostream&
operator<<( std::ostream &s, essential_matrix_<T> const &e )
{
  s << e.compute_matrix();
  return s;
}

/// Output stream operator for \p essential_matrix instances
std::ostream&
operator<<( std::ostream &s, essential_matrix const &e )
{
  s << e.matrix();
  return s;
}

// ===========================================================================
// Template class instantiation
// ---------------------------------------------------------------------------
/// \cond DoxygenSuppress
#define INSTANTIATE_ESSENTIAL_MATRIX(T)                                 \
  template class essential_matrix_<T>;                                  \
  template VITAL_EXPORT std::ostream& operator<<( std::ostream &,       \
                        essential_matrix_<T> const & )

INSTANTIATE_ESSENTIAL_MATRIX(float);
INSTANTIATE_ESSENTIAL_MATRIX(double);
#undef INSTANTIATE_ESSENTIAL_MATRIX
/// \endcond

} } // end vital namespace

ext: cxx

max_stars_repo_name: mwoehlke-kitware/kwiver

max_issues_repo_name: mwoehlke-kitware/kwiver

max_forks_repo_name: mwoehlke-kitware/kwiver


text: #include <boost/multiprecision/traits/is_byte_container.hpp>

ext: hpp

max_stars_repo_name: miathedev/BoostForArduino

max_issues_repo_name: miathedev/BoostForArduino

max_forks_repo_name: miathedev/BoostForArduino


text:
#include <boost/python.hpp>
#include <boost/python/numpy.hpp>
#include "dbscan.h"
#include "dbscan_vp.h"
#include "dbscan_graph.h"

using namespace boost::python;
using namespace clustering;

namespace p = boost::python;
namespace np = boost::python::numpy;

struct ublas_matrix_to_python {
    static PyObject* convert( DBSCAN::ClusterData const& C )
    {
        PyObject* result = PyList_New( C.size1() );

        for ( size_t i = 0; i < C.size1(); ++i ) {
            PyObject* l = PyList_New( C.size2() );

            for ( size_t j = 0; j < C.size2(); ++j ) {
                PyList_SetItem( l, j, PyFloat_FromDouble( C( i, j ) ) );
            }

            PyList_SetItem( result, i, l );
        }

        return result;
    }
};

class PyDBSCAN : public DBSCAN {
public:
    PyDBSCAN()
        : DBSCAN()
    {
    }

    void pyfit( boost::python::list& pylist )
    {
        auto num_samples = boost::python::len( pylist );

        if ( num_samples == 0 ) {
            return;
        }

        boost::python::extract< boost::python::list > first_elem( pylist[0] );

        auto num_features = boost::python::len( first_elem );

        DBSCAN::ClusterData C( num_samples, num_features );

        for ( int i = 0; i < num_samples; ++i ) {
            boost::python::list sublist = boost::python::extract< boost::python::list >( pylist[i] );

            for ( int j = 0; j < num_features; ++j ) {
                C( i, j ) = boost::python::extract< double >( sublist[j] );
            }
        }

        fit( C );
    }

    void pywfit( boost::python::list& pylist, boost::python::list& pyweights )
    {
        auto num_samples = boost::python::len( pylist );

        if ( num_samples == 0 ) {
            return;
        }

        boost::python::extract< boost::python::list > first_elem( pylist[0] );

        auto num_features = boost::python::len( first_elem );

        DBSCAN::ClusterData C( num_samples, num_features );

        for ( int i = 0; i < num_samples; ++i ) {
            boost::python::list sublist = boost::python::extract< boost::python::list >( pylist[i] );

            for ( int j = 0; j < num_features; ++j ) {
                C( i, j ) = boost::python::extract< double >( sublist[j] );
            }
        }

        auto num_weights = boost::python::len( pyweights );

        DBSCAN::FeaturesWeights W( num_weights );

        for ( int i = 0; i < num_weights; ++i ) {
            W( i ) = boost::python::extract< double >( pyweights[i] );
        }

        wfit( C, W );
    }

    boost::python::list pyget_labels()
    {
        boost::python::list list;

        for ( const auto& l : get_labels() ) {
            list.append( l );
        }

        return list;
    }
};

class npDataset : public Dataset {
public:
    typedef boost::shared_ptr< npDataset > Ptr;

    npDataset()
        : Dataset()
    {
    }

    virtual ~npDataset()
    {
    }

    bool load_ndarray( const np::ndarray& in )
    {
        m_rows = in.shape( 0 );
        m_cols = in.shape( 1 );

        for ( size_t i = 0; i < m_rows; ++i ) {
            Eigen::VectorXf col_vector( m_cols );
            for ( size_t j = 0; j < m_cols; ++j ) {
                col_vector( j ) = p::extract< float >( in[i][j] );
            }

            _data.emplace_back( col_vector );
        }

        return true;
    }
};

class PyDBSCANvp {
public:
    PyDBSCANvp( const np::ndarray& d )
        : m_rows( 0 )
    {
        npDataset::Ptr np_dset = boost::make_shared< npDataset >();
        np_dset->load_ndarray( d );

        m_dbs = boost::make_shared< DBSCAN_VP >( np_dset );
        m_dbs->fit();
        m_rows = np_dset->rows();
    }

    const np::ndarray predict( double eps, size_t min_elems )
    {
        np::dtype dtype = np::dtype::get_builtin< int32_t >();
        p::tuple shape = p::make_tuple( m_rows );

        np::ndarray rnd = np::zeros( shape, dtype );

        m_dbs->predict( eps, min_elems );

        const DBSCAN_VP::Labels& labels = m_dbs->get_labels();

        for ( size_t i = 0; i < labels.size(); ++i ) {
            rnd[i] = labels[i];
        }

        return rnd;
    }

    const np::ndarray predict_eps( size_t k )
    {
        np::dtype dtype = np::dtype::get_builtin< float >();
        p::tuple shape = p::make_tuple( m_rows );

        np::ndarray rnd = np::zeros( shape, dtype );

        const std::vector< double > r = m_dbs->predict_eps( k );

        for ( size_t i = 0; i < r.size(); ++i ) {
            rnd[i] = r[i];
        }
        return rnd;
    }

    const np::ndarray query_radius( const np::ndarray& in, double eps )
    {
        DBSCAN_VP::TVpTree::Ptr vpt = m_dbs->get_vp();

        size_t rows = in.shape( 0 );

        Eigen::VectorXf col_vector( rows );
        for ( size_t j = 0; j < rows; ++j ) {
            col_vector( j ) = p::extract< float >( in[j] );
        }

        DBSCAN_VP::TVpTree::TNeighborsList neigh;

        vpt->search_by_dist( col_vector, eps, neigh );

        np::dtype dtype = np::dtype::get_builtin< uint32_t >();
        p::tuple shape = p::make_tuple( neigh.size() );

        np::ndarray rnd = np::zeros( shape, dtype );

        for ( size_t j = 0; j < neigh.size(); ++j ) {
            rnd[j] = neigh[j].first;
        }

        return rnd;
    }

private:
    size_t m_rows;
    DBSCAN_VP::Ptr m_dbs;
};

class PyDBSCANgraph {
public:
    PyDBSCANgraph( const np::ndarray& d )
        : m_rows( 0 )
    {
        npDataset::Ptr np_dset = boost::make_shared< npDataset >();
        np_dset->load_ndarray( d );

        m_rows = np_dset->rows();
        m_dbs = boost::make_shared< DBSCAN_GRAPH >( np_dset );
    }

    const np::ndarray predict( float eps, size_t min_elems )
    {
        np::dtype dtype = np::dtype::get_builtin< int32_t >();
        p::tuple shape = p::make_tuple( m_rows );

        np::ndarray rnd = np::zeros( shape, dtype );

        m_dbs->fit( eps, min_elems );

        m_dbs->predict();

        const DBSCAN_GRAPH::Labels& labels = m_dbs->get_labels();

        for ( size_t i = 0; i < labels.size(); ++i ) {
            rnd[i] = labels[i];
        }

        return rnd;
    }

    const np::ndarray get_Va0()
    {
        np::dtype dtype = np::dtype::get_builtin< uint32_t >();
        p::tuple shape = p::make_tuple( m_rows );

        np::ndarray rnd = np::zeros( shape, dtype );

        const std::vector< uint32_t > r = m_dbs->get_Va0();

        for ( size_t i = 0; i < r.size(); ++i ) {
            rnd[i] = r[i];
        }
        return rnd;
    }

    const np::ndarray get_Va1()
    {
        np::dtype dtype = np::dtype::get_builtin< uint64_t >();
        p::tuple shape = p::make_tuple( m_rows );

        np::ndarray rnd = np::zeros( shape, dtype );

        const std::vector< uint64_t > r = m_dbs->get_Va1();

        for ( size_t i = 0; i < r.size(); ++i ) {
            rnd[i] = r[i];
        }
        return rnd;
    }

private:
    size_t m_rows;
    DBSCAN_GRAPH::Ptr m_dbs;
};

BOOST_PYTHON_MODULE( pydbscan )
{
    np::initialize();

    def( "gen_cluster_data", &PyDBSCAN::gen_cluster_data );
    //def( "send_stat", & python_send_stat );

    class_< PyDBSCAN >( "DBSCAN" )
        .def( "init", &PyDBSCAN::init )
        .def( "fit", &PyDBSCAN::pyfit )
        .def( "wfit", &PyDBSCAN::pywfit )
        .def( "get_labels", &PyDBSCAN::pyget_labels );

    class_< PyDBSCANvp >( "DBSCANvp", p::init< np::ndarray >() )
        .def( "predict", &PyDBSCANvp::predict )
        .def( "predict_eps", &PyDBSCANvp::predict_eps )
        .def( "query_radius", &PyDBSCANvp::query_radius );

    class_< PyDBSCANgraph >( "DBSCANgraph", p::init< np::ndarray >() )
        .def( "predict", &PyDBSCANgraph::predict )
        .def( "get_Va0", &PyDBSCANgraph::get_Va0 )
        .def( "get_Va1", &PyDBSCANgraph::get_Va1 );

    to_python_converter< DBSCAN::ClusterData, ublas_matrix_to_python, false >();
}

ext: cpp


text:
#define _USE_MATH_DEFINES

#include <Eigen/Dense>
#include <cmath>
#include <iostream>

#include <geometry_msgs/PointStamped.h>
#include <geometry_msgs/PoseStamped.h>
#include <igvc_msgs/velocity_pair.h>
#include <nav_msgs/Odometry.h>
#include <visualization_msgs/MarkerArray.h>

#include <tf/transform_datatypes.h>

#include <igvc_utils/robot_state.h>
#include <parameter_assertions/assertions.h>

#include "path_follower.h"

PathFollower::PathFollower()
{
  ros::NodeHandle nh;
  ros::NodeHandle pNh("~");

  // load controller parameters
  using seconds = double;

  double target_velocity;
  double axle_length;
  double k1;
  double k2;
  double simulation_frequency;
  double lookahead_dist;
  seconds simulation_horizon;

  assertions::param(pNh, "target_v", target_velocity, 1.0);
  assertions::param(pNh, "axle_length", axle_length, 0.52);
  assertions::param(pNh, "k1", k1, 1.0);
  assertions::param(pNh, "k2", k2, 3.0);
  assertions::param(pNh, "simulation_frequency", simulation_frequency, 100.0);
  assertions::param(pNh, "lookahead_dist", lookahead_dist, 2.0);
  assertions::param(pNh, "simulation_horizon", simulation_horizon, 5.0);

  controller_ = std::unique_ptr<SmoothControl>(new SmoothControl{
      k1, k2, axle_length, simulation_frequency, target_velocity, lookahead_dist, simulation_horizon });

  // load global parameters
  assertions::getParam(pNh, "maximum_vel", maximum_vel_, { assertions::NumberAssertionType::POSITIVE });
  assertions::param(pNh, "stop_dist", stop_dist_, 0.9);

  ros::Subscriber path_sub = nh.subscribe("/path", 1, &PathFollower::pathCallback, this);
  ros::Subscriber pose_sub = nh.subscribe("/odometry/filtered", 1, &PathFollower::positionCallback, this);
  ros::Subscriber waypoint_sub = nh.subscribe("/waypoint", 1, &PathFollower::waypointCallback, this);

  cmd_pub_ = nh.advertise<igvc_msgs::velocity_pair>("/motors", 1);
  target_pub_ = nh.advertise<geometry_msgs::PointStamped>("/target_point", 1);
  trajectory_pub_ = nh.advertise<nav_msgs::Path>("/trajectory", 1);

  ros::spin();
}

void PathFollower::pathCallback(const nav_msgs::PathConstPtr& msg)
{
  ROS_DEBUG_STREAM("Follower got path. Size: " << msg->poses.size());
  //  path_ = msg;
  // TODO: Remove patch when motion planning correctly incorporates heading
  path_ = getPatchedPath(msg);
}

void PathFollower::waypointCallback(const geometry_msgs::PointStampedConstPtr& msg)
{
  waypoint_ = msg;
}

/**
Constructs a new trajectory to follow using the current path msg and publishes
the first velocity command from this trajectory.
*/
void PathFollower::positionCallback(const nav_msgs::OdometryConstPtr& msg)
{
  if (path_.get() == nullptr || path_->poses.empty() || path_->poses.size() < 2)
  {
    ROS_INFO_THROTTLE(1, "Path empty.");
    igvc_msgs::velocity_pair vel;
    vel.left_velocity = 0.;
    vel.right_velocity = 0.;
    cmd_pub_.publish(vel);
    path_.reset();
    return;
  }

  // Current pose
  RobotState cur_pos(msg);

  double goal_dist = cur_pos.distTo(waypoint_->point.x, waypoint_->point.y);

  ROS_DEBUG_STREAM("Distance to waypoint: " << goal_dist << "(m.)");

  ros::Time time = msg->header.stamp;
  igvc_msgs::velocity_pair vel;  // immediate velocity command
  vel.header.stamp = time;

  if (goal_dist <= stop_dist_)
  {
    /**
    Stop when the robot is a set distance away from the waypoint
    */
    ROS_INFO_STREAM(">>>WAYPOINT REACHED...STOPPING<<<");
    vel.right_velocity = 0;
    vel.left_velocity = 0;
    // make path null to stop planning until path generated
    // for new waypoint
    path_ = nullptr;
  }
  else
  {
    /**
    Obtain smooth control law from the controller. This includes a smooth
    trajectory for visualization purposes and an immediate velocity command.
    */
    nav_msgs::Path trajectory_msg;
    trajectory_msg.header.stamp = time;
    trajectory_msg.header.frame_id = "/odom";

    RobotState target;
    controller_->getTrajectory(vel, path_, trajectory_msg, cur_pos, target);
    // publish trajectory
    trajectory_pub_.publish(trajectory_msg);

    // publish target position
    geometry_msgs::PointStamped target_point;
    target_point.header.frame_id = "/odom";
    target_point.header.stamp = time;
    tf::pointTFToMsg(target.transform.getOrigin(), target_point.point);
    target_pub_.publish(target_point);

    ROS_DEBUG_STREAM("Distance to target: " << cur_pos.distTo(target) << "(m.)");
  }

  // make sure maximum velocity not exceeded
  if (std::max(std::abs(vel.right_velocity), std::abs(vel.left_velocity)) > maximum_vel_)
  {
    ROS_ERROR_STREAM_THROTTLE(5, "Maximum velocity exceeded. Right: "
                                     << vel.right_velocity << "(m/s), Left: " << vel.left_velocity
                                     << "(m/s), Max: " << maximum_vel_ << "(m/s) ... Stopping robot...");
    vel.right_velocity = 0;
    vel.left_velocity = 0;
  }

  cmd_pub_.publish(vel);  // pub velocity command
}

nav_msgs::PathConstPtr PathFollower::getPatchedPath(const nav_msgs::PathConstPtr& msg) const
{
  nav_msgs::PathPtr new_path = boost::make_shared<nav_msgs::Path>(*msg);
  size_t num_poses = new_path->poses.size();
  for (size_t i = 0; i < num_poses; i++)
  {
    double delta_x = new_path->poses[i + 1].pose.position.x - new_path->poses[i].pose.position.x;
    double delta_y = new_path->poses[i + 1].pose.position.y - new_path->poses[i].pose.position.y;
    double heading = atan2(delta_y, delta_x);
    new_path->poses[i].pose.orientation = tf::createQuaternionMsgFromYaw(heading);
  }
  new_path->poses.back().pose.orientation = new_path->poses[num_poses - 2].pose.orientation;
  return new_path;
}

int main(int argc, char** argv)
{
  ros::init(argc, argv, "path_follower");
  PathFollower path_follower;
  return 0;
}

max_stars_repo_name: jiajunmao/igvc-software

max_issues_repo_name: jiajunmao/igvc-software

max_forks_repo_name: jiajunmao/igvc-software


text:
#include "ActiveTask.h"
#include <nar/lib/Exception/Exception.h>
#include <nar/lib/Messaging/MessageTypes/FilePull.h>
#include <nar/lib/Socket/USocket.h>
#include <nar/lib/Messaging/MessageTypes/InfoChunkPull.h>
#include <iostream>
#include <string>
#include <algorithm>
#include <boost/filesystem.hpp>

#define CRYPTOPP_ENABLE_NAMESPACE_WEAK 1
#include <crypto++/md5.h>
#include <cryptopp/filters.h>
#include <cryptopp/hex.h>

using std::string;
using std::cout;
using std::endl;

// TO_DO: Last part should handled for file::encrypt,decompress etc..

void nar::ActiveTask::Pull::run(nar::Socket* ipc_socket, nar::MessageTypes::IPCPull::Request* req) {
    nar::Socket* server_sck = this->_globals->establish_server_connection();
    std::string file_aes;
    try {
        file_aes = nar::ActiveTask::user_authenticate(server_sck, this->_vars);
    } catch (nar::Exception::Daemon::AuthenticationError& exp) {
        NAR_LOG<<exp.what()<<std::endl;
        return;
    }
    string file_name,system_dir,nar_dir_name;
    try {
        file_name = req->get_file_name();
        system_dir = req->get_dir_path(); // to be changed
        nar_dir_name =  req->get_current_directory();
    }
    catch (...) {
        NAR_LOG << "Malformed Json, exiting pull" << std::endl;
        return;
    }

    nar::MessageTypes::FilePull::Request pull_req(file_name, nar_dir_name);
    nar::MessageTypes::FilePull::Response pull_resp;
    try {
        pull_req.send_mess(server_sck, pull_resp);
    }
    catch (nar::Exception::ExcpBase& e) {
        NAR_LOG << std::string( "nar_daemon::activepull " ).append(e.what()) << std::endl;
        return;
    }

    if ( pull_resp.get_status_code() == 666 ) {     // Not enough online peer
        nar::MessageTypes::IPCPull::Response ipcpull_resp(3,666);           // Params ?
        ipcpull_resp.send_message_end(ipc_socket);
        return;
    }

    unsigned short rand_port = pull_resp.get_rendezvous_port();
    std::vector<struct nar::MessageTypes::FilePull::Response::PeerListElement> elements = pull_resp.get_elements();
    std::sort(elements.begin(), elements.end());

    boost::filesystem::path tpath(system_dir);
    tpath /= file_name;

    boost::filesystem::path temp;
    try {
        temp = boost::filesystem::unique_path();
    } catch(std::ios_base::failure& Exp) {
        throw nar::Exception::Unknown(Exp.what());
    }
    std::string temp_native = temp.native();

    nar::File* tempfile1 = new nar::File(temp_native,"w",false);
    int i=0;
    try {
        unsigned long offset = 0;
        while( i<elements.size() ) {

            unsigned long stream_id = elements[i].stream_id;
            nar::USocket* cli_sck = new nar::USocket(this->_globals->get_ioserv(), this->_globals->get_server_ip(), rand_port, stream_id);
            char* buf = NULL;
            try {
                cli_sck->connect();
                long int total_read = 0;
                buf = new char [elements[i].chunk_size];
                while(total_read < elements[i].chunk_size) {
                    int len = cli_sck->recv(buf+total_read, elements[i].chunk_size - total_read);
                    total_read += len;
                }
                byte digest[ CryptoPP::Weak::MD5::DIGESTSIZE];
            	CryptoPP::Weak::MD5 hash;
            	hash.CalculateDigest( digest, (const byte*)buf, elements[i].chunk_size );
            	CryptoPP::HexEncoder encoder;
            	std::string output;
            	encoder.Attach( new CryptoPP::StringSink( output ) );
            	encoder.Put( digest, sizeof(digest) );
            	encoder.MessageEnd();

               /* if (  output.compare(elements[i].hashed) ) // Wrong Hash
                {
                    std::cout <<"Benim Hesabim: " << output << std::endl;
                    std::cout <<"Carsidaki Hesap: " << elements[i].hashed << std::endl;
                    throw nar::Exception::Unknown("Wrong Hash");
                } */

                tempfile1->write(buf,total_read);
                offset+=total_read;
            }
            catch (nar::Exception::File::WrongMode& e) {
                NAR_LOG << e.what() << std::endl;
            }
            catch (nar::Exception::File::OffsetOutOfBounds& e) {
                NAR_LOG << e.what() << std::endl;
            }
            catch ( nar::Exception::DomainError::Negative& e) {
                NAR_LOG << e.what() << std::endl;
            }
            catch ( nar::Exception::File::WriteError& e) {
                NAR_LOG << e.what() << std::endl;
            }
            catch ( nar::Exception::File::NotOpen& e) {
                NAR_LOG << e.what() << std::endl;
            }

            catch (...) {
                NAR_LOG<<"pull.cpp inside"<<std::endl;
                nar::MessageTypes::InfoChunkPull::Request _req(elements[i].chunk_id, 704);
                nar::MessageTypes::InfoChunkPull::Response _resp;
                _req.send_mess(server_sck,_resp);
                elements[i].stream_id = _resp.get_stream_id();
                if (buf)
                    delete [] buf;
                continue;
            }

            cli_sck->close();
            i++;
        }
    }
    catch ( ... ) {
        NAR_LOG << "Problem occured when pulling chunks from peers" << std::endl;
        delete tempfile1;
        return;
    }

    delete tempfile1;
    tempfile1 = new nar::File(temp_native,"r",true);
    nar::File * decrypted;
    try{
        decrypted = tempfile1->decrypt(file_aes);
    } catch(...) {
        NAR_LOG<<"can not decrytp"<<std::endl;
    }

    string dust= tpath.native();
    nar::File * decompressed = decrypted->decompress(dust);

    nar::MessageTypes::InfoChunkPull::Request _req(50, 200);
    _req.send_mess(server_sck);
    try {
        nar::MessageTypes::IPCPull::Response ipcpull_resp(100,200);

        ipcpull_resp.send_message_end(ipc_socket);
    } catch (...){
        NAR_LOG << "cli seems down!!! " << std::endl;
    }
    delete tempfile1;
    delete decrypted;
    delete decompressed;
    return;
}

max_stars_repo_name: webcok/nar

max_issues_repo_name: webcok/nar

max_forks_repo_name: webcok/nar


text:
#pragma once

#include <Eigen/Core>
#include <polysolve/LinearSolver.hpp>

#include <constants.hpp>
#include <solvers/optimization_solver.hpp>
#include <utils/not_implemented_error.hpp>

namespace ipc::rigid {

enum ConvergenceCriteria {
    VELOCITY, ///< Change in position of a newton direction
    ENERGY    ///< Change in energy of a newton direction
};

NLOHMANN_JSON_SERIALIZE_ENUM(
    ConvergenceCriteria, { { VELOCITY, "velocity" }, { ENERGY, "energy" } });

class NewtonSolver : public virtual OptimizationSolver {
public:
    NewtonSolver();
    virtual ~NewtonSolver() = default;

    /// Initialize the state of the solver using the settings saved in JSON
    virtual void settings(const nlohmann::json& params) override;
    /// Export the state of the solver using the settings saved in JSON
    virtual nlohmann::json settings() const override;

    /// An identifier for the solver class
    static std::string solver_name() { return "newton_solver"; }
    /// An identifier for this solver
    virtual std::string name() const override
    {
        return NewtonSolver::solver_name();
    }

    /// Initialize the solver with a problem to solve
    virtual void set_problem(OptimizationProblem& problem) override
    {
        this->problem_ptr = &problem;
    }

    /// Initialize the solver state for a new solve
    virtual void init_solve(const Eigen::VectorXd& x0) override;

    /**
     * @brief Perform Newton's Method to minimize the objective, \f$f(x)\f$,
     * of the problem unconstrained.
     *
     * @param[in] problem  The optimization problem to minimize
     *                     unconstrained.
     *
     * @return The results of the optimization including the minimizer,
     * minimum, and if the optimization was successful.
     */
    virtual OptimizationResults solve(const Eigen::VectorXd& x0) override;

    /// Perform a single step of solving the optimization problem
    virtual OptimizationResults step_solve() override
    {
        throw NotImplementedError(
            "Taking a single newton step is not implemented!");
    };

    /**
     * @brief Solve for the Newton direction
     *        (\f$\Delta x = -H^{-1} \nabla f \f$).
     *
     * @param[in]  gradient  Gradient of the objective function.
     * @param[in]  hessian   Hessian of the objective function.
     * @param[out] delta_x   Output newton direction.
     * @param[in]  make_psd  If delta_x is not a descent direction, then
     * make the hessian positive semi-definite.
     *
     * @return Returns true if the solve was successful.
     */
    virtual bool compute_direction(
        const Eigen::VectorXd& gradient,
        const Eigen::SparseMatrix<double>& hessian,
        Eigen::VectorXd& delta_x,
        bool make_psd = false);

    virtual bool compute_regularized_direction(
        double& fx,
        Eigen::VectorXd& gradient,
        Eigen::SparseMatrix<double>& hessian,
        Eigen::VectorXd& delta_x,
        double& coeff);

    virtual std::string stats_string() const override;
    virtual nlohmann::json stats() const override;

    int max_iterations;

protected:
    virtual bool converged();

    virtual void post_step_update();

    virtual bool line_search(
        const Eigen::VectorXd& x,
        const Eigen::VectorXd& dir,
        const double fx,
        const Eigen::VectorXd& grad_fx,
        double& step_length);

    virtual double line_search_lower_bound() const
    {
        return m_line_search_lower_bound;
    }

    /// @brief Pointer to the problem to solve.
    OptimizationProblem* problem_ptr;

    int iteration_number; ///< @brief The current iteration number.
    ConvergenceCriteria convergence_criteria;

    double m_line_search_lower_bound; ///< @brief Line search lower bound

    double energy_conv_tol;        ///< @brief Energy convergence tolerance
    double velocity_conv_tol;      ///< @brief Velocity convergence tolerance
    bool is_velocity_conv_tol_abs; ///< @brief Absolute velocity tol
    bool is_energy_converged;

    // State variables
    Eigen::VectorXd x, x_prev;
    Eigen::VectorXd gradient, gradient_free;
    Eigen::VectorXd direction, direction_free;
    Eigen::VectorXd grad_direction; ///< Gradient with fixed DoF set to zero
    Eigen::SparseMatrix<double> hessian, hessian_free;

    // Linear solver pointer
    std::unique_ptr<polysolve::LinearSolver> linear_solver;
    nlohmann::json linear_solver_settings;

private:
    void reset_stats();

    size_t num_fx = 0;
    size_t num_grad_fx = 0;
    size_t num_hessian_fx = 0;
    size_t num_collision_check = 0;
    size_t ls_iterations = 0;
    size_t newton_iterations = 0;
    size_t num_newton_ls_fails = 0;
    size_t num_grad_ls_fails = 0;
    size_t regularization_iterations = 0;
};

/**
 * @brief Make the matrix positive definite (\f$x^T A x > 0\$).
 *
 * @param A The matrix to make positive definite.
 *
 * @return The scale of the update to the diagonal.
 */
double make_matrix_positive_definite(Eigen::SparseMatrix<double>& A);

/**
 * @brief Log values along a search direction.
 *
 * @param[in] x                 Starting point for the line search.
 * @param[in] dir               Direction to search along.
 * @param[in] f_and_gradf       Function of x to sample with gradient.
 */
void sample_search_direction(
    const Eigen::VectorXd& x,
    const Eigen::VectorXd& dir,
    const std::function<double(const Eigen::VectorXd&, Eigen::VectorXd&)>&
        f_and_gradf,
    double max_step);

} // namespace ipc::rigid

max_stars_repo_name: ipc-sim/rigid-ipc

max_issues_repo_name: ipc-sim/rigid-ipc

max_forks_repo_name: ipc-sim/rigid-ipc


text:
/*
 * main.cpp
 *
 *  Created on: 09.06.2010
 *      Author: fkrauthan
 */

#include "Bot/bot.h"

#include <boost/program_options.hpp>

#include <iostream>
#include <string>

int startBot(const std::string& configFile);
void startDaemon();

int main(int argc, char **argv) {
	boost::program_options::options_description desc("Allowed options");
	desc.add_options()
		("help", "produce help message")
		("config", boost::program_options::value<std::string>()->default_value("ircbot.ini"), "sets the bot config file")
	;

	std::string configFile = "";

	//Process programm options
	boost::program_options::variables_map vm;
	boost::program_options::store(boost::program_options::parse_command_line(argc, argv, desc), vm);
	boost::program_options::notify(vm);
	if(vm.count("help")) {
		std::cout << "Usage: " << argv[0] << " [options]" << std::endl;
		std::cout << desc << std::endl;
	}
	if(vm.count("config")) {
		configFile = vm["config"].as<std::string>();
	}

	//Run the bot
	return startBot(configFile);

	/*bool runDaemon = false;
	if(argc>1) {
		if(std::strcmp(argv[1], "--daemon") == 0) {
			runDaemon = true;
		}
	}

	if(runDaemon) {
		startDaemon();
	}
	else {
		return startBot();
	}*/

	return -1;
}

int startBot(const std::string& configFile) {
	std::cout << "IrcBot3 v0.1\n" << std::endl;

	Irc::ClientInfo clientInfo;
	clientInfo.clientName = "IrcBot3";
	clientInfo.version = "0.1";

	Bot bot(clientInfo);
	if(!bot.init(configFile)) {
		return -1;
	}
	bot.run();

	return 0;
}

void startDaemon() {
	/*umask(0);

	pid_t pid;
	if((pid = fork()) < 0) {
		//Cannot fork!
	    throw std::runtime_error(strerror(errno));
	} else if (pid != 0) { //parent
		exit(0);
	}

	pid_t sid = setsid();
	if (sid < 0) {
		exit(-1);
	}*/

	/* Close out the standard file descriptors */
	//close(STDIN_FILENO);
	//close(STDOUT_FILENO);
	//close(STDERR_FILENO);

	//Run the bot
	//exit(startBot());
}

max_stars_repo_name: fkrauthan/IrcBot

max_issues_repo_name: fkrauthan/IrcBot

max_forks_repo_name: fkrauthan/IrcBot


text:
// simple_variant.hpp
/*
 *  Copyright (c) 2007 Leigh Johnston.
 *
 *  All rights reserved.
 *
 *  Redistribution and use in source and binary forms, with or without
 *  modification, are permitted provided that the following conditions are
 *  met:
 *
 *     * Redistributions of source code must retain the above copyright
 *       notice, this list of conditions and the following disclaimer.
 *
 *     * Redistributions in binary form must reproduce the above copyright
 *       notice, this list of conditions and the following disclaimer in the
 *       documentation and/or other materials provided with the distribution.
 *
 *     * Neither the name of Leigh Johnston nor the names of any
 *       other contributors to this software may be used to endorse or
 *       promote products derived from this software without specific prior
 *       written permission.
 *
 *  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
 *  IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
 *  THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
 *  PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
 *  CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
 *  EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
 *  PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
 *  PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
 *  LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
 *  NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 *  SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/

#pragma once

#include <neolib/neolib.hpp>
#include <string>
#include <boost/lexical_cast.hpp>
#include <neolib/core/reference_counted.hpp>
#include <neolib/core/string.hpp>
#include <neolib/core/string_numeric.hpp>
#include <neolib/core/i_enum.hpp>
#include <neolib/core/i_custom_type.hpp>
#include <neolib/core/i_simple_variant.hpp>
#include <neolib/plugin/plugin_variant.hpp>

namespace neolib
{
    typedef plugin_variant<simple_variant_type, bool, int64_t, double, string, ref_ptr<i_enum>, ref_ptr<i_custom_type>> simple_variant;

    inline simple_variant from_string(std::string const& aString, simple_variant_type aType)
    {
        switch (aType)
        {
        case simple_variant_type::Boolean:
            if (aString == "true" || aString == "1")
                return true;
            else
                return false;
        case simple_variant_type::Integer:
            return string_to_int64(aString);
        case simple_variant_type::Real:
            return string_to_double(aString);
        case simple_variant_type::String:
            return aString;
        default:
            throw std::logic_error("neolib: cannot convert string to simple variant");
        }
    }

    inline std::string to_string(simple_variant const& aVariant)
    {
        return to_string(static_cast<i_simple_variant const&>(aVariant));
    }
}


text:
/*=============================================================================
    Copyright (c) 2001-2007 Joel de Guzman

    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)
=============================================================================*/
#if !defined(SPIRIT_ACTION_JAN_07_2007_1128AM)
#define SPIRIT_ACTION_JAN_07_2007_1128AM

#include <boost/spirit/home/qi/domain.hpp>
#include <boost/spirit/home/support/component.hpp>
#include <boost/spirit/home/support/attribute_of.hpp>
#include <boost/spirit/home/support/detail/action_dispatch.hpp>
#include <boost/mpl/if.hpp>
#include <boost/mpl/bool.hpp>
#include <boost/mpl/identity.hpp>
#include <boost/type_traits/remove_const.hpp>
#include <boost/type_traits/is_same.hpp>
#include <vector>

namespace boost { namespace spirit { namespace qi
{
    struct action
    {
        template <typename Component, typename Context, typename Iterator>
        struct attribute
          : traits::attribute_of<
                qi::domain
              , typename result_of::left<Component>::type
              , Context
              , Iterator
            >
        {
        };

        template <typename F, typename Attribute, typename Context>
        static bool const_action_dispatch(
            F const& f, Attribute const& attr, Context& context)
        {
            // This function makes Attribute a const reference
            // before calling detail::action_dispatch whereby
            // disallowing mutability of the attribute in semantic
            // actions.
            return spirit::detail::action_dispatch(f, attr, context,
                mpl::true_());
        }

        template <
            typename Component
          , typename Iterator, typename Context
          , typename Skipper, typename Attribute>
        static bool parse(
            Component const& component
          , Iterator& first, Iterator const& last
          , Context& context, Skipper const& skipper
          , Attribute& attr_)
        {
            typedef typename
                result_of::left<Component>::type::director
            director;

            typedef typename
                attribute<Component, Context, Iterator>::type
            attr_type;

            // create an attribute if one is not supplied
            typename mpl::if_<
                is_same<typename remove_const<Attribute>::type, unused_type>
              , typename remove_const<attr_type>::type
              , Attribute&>::type
            attr = spirit::detail::make_value<attr_type>::call(attr_);

            if (director::parse(
                spirit::left(component), first, last, context, skipper, attr))
            {
                // call the function, passing the attribute, the context.
                // The client can return false to fail parsing.
                return const_action_dispatch(
                    spirit::right(component), attr, context);
            }
            return false;
        }

        template <typename Component, typename Context>
        static std::string what(Component const& component, Context const& ctx)
        {
            typedef typename
                result_of::left<Component>::type::director
            director;

            return director::what(spirit::left(component), ctx);
        }
    };
}}}

#endif

max_stars_repo_name: jaredhoberock/gotham

max_issues_repo_name: jaredhoberock/gotham

max_forks_repo_name: jaredhoberock/gotham


text:
/* CSVReader.cpp
 * -*- coding: utf-8 -*-
 *
 */

#include <boost/test/unit_test.hpp>

#include <stdexcept>

#include <boost/filesystem.hpp>

#include <Eigen/Dense>

// local includes ==============================================================

#include "../src/utils/CSVReader.hpp"

#include "utils/log.hpp"

BOOST_AUTO_TEST_SUITE(CSVReaderTestSuite)
    struct Setup {
        Setup() { reader = CSVReader(); }

        CSVReader reader;

        std::string test_file = "test/resources/csvreader/matrix.csv";
    };

    BOOST_AUTO_TEST_CASE(read_csv) {
        TEST_MESSAGE("read_csv");

        Setup f;

        Eigen::MatrixXd expected(2,2);

        expected(0,0) = 0.0000000000000000;
        expected(0,1) = 1101.5312671722938376;
        expected(1,0) = 3.8522202257825890;
        expected(1,1) = 1119.2784683923078774;

        Eigen::MatrixXd actual = f.reader.read_matrix(f.test_file, 2);

        BOOST_REQUIRE_EQUAL(expected, actual);
    }

    BOOST_AUTO_TEST_CASE(throw_if_file_nonexistent) {
        TEST_MESSAGE("throw_if_file_nonexistent");

        Setup f;

        BOOST_REQUIRE_THROW(f.reader.read_matrix(boost::filesystem::path("test/resources/csvreader/nonexistent_file.csv"), 2),
                            std::runtime_error);
    }

BOOST_AUTO_TEST_SUITE_END()

// vim: tabstop=8 expandtab shiftwidth=4 softtabstop=4 fenc=utf-8

max_stars_repo_name: AFriemann/LowCarb

max_issues_repo_name: AFriemann/LowCarb

max_forks_repo_name: AFriemann/LowCarb


text:
/*=============================================================================
  Copyright (c) 2010-2015 Bolero MURAKAMI
  https://github.com/bolero-MURAKAMI/KTL

  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 SRC_VARIATE_RANDOM_VARIATE_RANDOM_COMMON_HPP
#define SRC_VARIATE_RANDOM_VARIATE_RANDOM_COMMON_HPP

#include <boost/foreach.hpp>
#include <boost/variant.hpp>
#include <boost/random.hpp>
#include <boost/random/beta_distribution.hpp>
#include <boost/random/laplace_distribution.hpp>
#include <boost/range/reference.hpp>
#include <boost/type_traits/is_unsigned.hpp>
#include <boost/type_traits/is_same.hpp>
#include <boost/mpl/identity.hpp>
#include <boost/mpl/size_t.hpp>
#include <boost/mpl/equal_to.hpp>
#include <boost/mpl/list.hpp>
#include <boost/mpl/size.hpp>
#include <boost/mpl/at.hpp>
#include <boost/mpl/transform.hpp>
#include <boost/mpl/quote.hpp>
#include <boost/mpl/assert.hpp>
#include <boost/mpl/has_xxx.hpp>
#include <boost/utility/swap.hpp>
#include <boost/utility/enable_if.hpp>
#include <boost/preprocessor/cat.hpp>
#include <boost/preprocessor/repetition/enum.hpp>
#include <sprig/external/tp_stub.hpp>
#include <sprig/random/adapter.hpp>
#include <sprig/random/null.hpp>
#include <sprig/random/any.hpp>
#include <sprig/random/variant.hpp>
#include <sprig/random/nondet_random.hpp>
#include <sprig/krkr/tjs.hpp>
#include <sprig/krkr/macro.hpp>
#include <ktl/config.hpp>

//
// KTL_VARIATE_RANDOM_SEQ_AT_GEN
//
#define KTL_VARIATE_RANDOM_SEQ_AT_GEN(z, n, d) boost::mpl::at_c<d, n>::type

namespace ktl {
	namespace patch {
		//
		// geometric_distribution
		//	for: Boost 1.47
		//
		template<
			typename IntType = boost::random::geometric_distribution<>::result_type,
			typename RealType = boost::random::geometric_distribution<>::input_type
			>
		class geometric_distribution
			: public boost::random::geometric_distribution<IntType, RealType>
		{
		private:
			typedef boost::random::geometric_distribution<IntType, RealType> base_type;
		public:
			explicit geometric_distribution(typename base_type::input_type const& p = typename base_type::input_type(0.5))
				: base_type(p)
			{}
			explicit geometric_distribution(typename base_type::param_type const& param)
				: base_type(param)
			{}
			typename base_type::result_type min() const {
				return typename base_type::result_type(0);
			}
		};
	}	// namespace patch
	//
	// RandomResult
	//
	class RandomResult {
	private:
		tTJSVariant t_;
	public:
		RandomResult() {}
		RandomResult(RandomResult const& other)
			: t_(other.t_)
		{}
		RandomResult(tTJSVariant const& t)
			: t_(t)
		{}
		template<typename T>
		RandomResult(
			T t,
			typename boost::enable_if<boost::is_unsigned<T> >::type* = 0
			)
			: t_(static_cast<tTVInteger>(t))
		{}
		RandomResult(boost::random::uniform_on_sphere<tTVReal>::result_type const& t) {
			typedef boost::random::uniform_on_sphere<tTVReal>::result_type const range_type;
			sprig::krkr::tjs::object_type obj = sprig::krkr::tjs::CreateNewObject(
				sprig::krkr::tjs::GetTJSClassNoAddRef(SPRIG_KRKR_TJS_W("Array")),
				0, 0, 0
				);
			{
				tjs_int num = 0;
				BOOST_FOREACH(boost::range_reference<range_type>::type e, t) {
					sprig::krkr::tjs::AddMember(obj, num, tTJSVariant(e));
					++num;
				}
			}
			t_ = sprig::krkr::tjs::as_object_closure_variant(obj);
		}
		RandomResult& operator=(RandomResult const& other) {
			RandomResult temp(other);
			temp.swap(*this);
			return *this;
		}
		RandomResult& operator=(tTJSVariant const& t) {
			RandomResult temp(t);
			temp.swap(*this);
			return *this;
		}
		RandomResult& operator=(boost::random::uniform_on_sphere<tTVReal>::result_type const& t) {
			RandomResult temp(t);
			temp.swap(*this);
			return *this;
		}
		void swap(RandomResult& other) {
			boost::swap(t_, other.t_);
		}
		operator tTJSVariant const&() const {
			return t_;
		}
		tTJSVariant const& get() const {
			return t_;
		}
	};
	//
	// RandomUtils
	//
	struct RandomUtils {
	public:
		BOOST_MPL_HAS_XXX_TRAIT_DEF(is_static_random_engine);
		template<typename Engine>
		struct enable_static_random_engine
			: public boost::enable_if<has_is_static_random_engine<Engine> >
		{};
		template<typename Engine>
		struct engine_wrap
			: public boost::mpl::identity<sprig::random::shared_engine_adapter<Engine, true> >
		{};
		template<typename Distribution>
		struct distribution_wrap
			: public boost::mpl::identity<sprig::random::shared_distribution_adapter<Distribution> >
		{};
	public:
		//	CHECK: fatal error C1060
#define KTL_VARIATE_RANDOM_ENGINE_SEQ_NUM 29/*31*/
#define KTL_VARIATE_RANDOM_DISTRIBUTION_SEQ_NUM 27

		typedef sprig::random::any_variate_generator<RandomResult> variate_generator_type;
		typedef boost::random::minstd_rand0 minstd_rand0;
		typedef boost::random::minstd_rand minstd_rand;
		typedef boost::random::rand48 rand48;
		typedef boost::random::ecuyer1988 ecuyer1988;
		typedef boost::random::knuth_b knuth_b;
		typedef boost::random::kreutzer1986 kreutzer1986;
		typedef boost::random::taus88 taus88;
		typedef boost::random::hellekalek1995 hellekalek1995;
		typedef boost::random::mt11213b mt11213b;
		typedef boost::random::mt19937 mt19937;
		typedef boost::random::mt19937_64 mt19937_64;
		typedef boost::random::lagged_fibonacci607 lagged_fibonacci607;
		typedef boost::random::lagged_fibonacci1279 lagged_fibonacci1279;
		typedef boost::random::lagged_fibonacci2281 lagged_fibonacci2281;
		typedef boost::random::lagged_fibonacci3217 lagged_fibonacci3217;
		typedef boost::random::lagged_fibonacci4423 lagged_fibonacci4423;
		typedef boost::random::lagged_fibonacci9689 lagged_fibonacci9689;
		typedef boost::random::lagged_fibonacci19937 lagged_fibonacci19937;
		typedef boost::random::lagged_fibonacci23209 lagged_fibonacci23209;
		typedef boost::random::lagged_fibonacci44497 lagged_fibonacci44497;
		typedef boost::random::ranlux3 ranlux3;
		typedef boost::random::ranlux4 ranlux4;
		typedef boost::random::ranlux64_3 ranlux64_3;
		typedef boost::random::ranlux64_4 ranlux64_4;
		typedef boost::random::ranlux3_01 ranlux3_01;
		typedef boost::random::ranlux4_01 ranlux4_01;
		typedef boost::random::ranlux64_3_01 ranlux64_3_01;
		typedef boost::random::ranlux64_4_01 ranlux64_4_01;
		typedef boost::random::ranlux24 ranlux24;
		typedef boost::random::ranlux48 ranlux48;
		typedef sprig::random::singleton_random_device random_device;
		typedef boost::mpl::BOOST_PP_CAT(list, KTL_VARIATE_RANDOM_ENGINE_SEQ_NUM)<
			minstd_rand0,
			minstd_rand,
			rand48,
			ecuyer1988,
			knuth_b,
			kreutzer1986,
			taus88,
			hellekalek1995,
			mt11213b,
			mt19937,
			mt19937_64,
			lagged_fibonacci607,
			lagged_fibonacci1279,
			lagged_fibonacci2281,
			lagged_fibonacci3217,
			lagged_fibonacci4423,
			lagged_fibonacci9689,
			lagged_fibonacci19937,
			lagged_fibonacci23209,
			lagged_fibonacci44497,
			ranlux3,
			ranlux4,
			ranlux64_3,
			ranlux64_4,
			ranlux3_01,
			ranlux4_01,
			ranlux64_3_01,
			ranlux64_4_01,
			//	CHECK: fatal error C1060
			//ranlux24,
			//ranlux48,
			random_device
		> engine_seq_type;
		typedef boost::random::uniform_smallint<tTVInteger> uniform_smallint;
		typedef boost::random::uniform_int_distribution<tTVInteger> uniform_int_distribution;
		typedef boost::random::uniform_01<tTVReal, tTVReal> uniform_01;
		typedef boost::random::uniform_real_distribution<tTVReal> uniform_real_distribution;
		typedef boost::random::bernoulli_distribution<tTVReal> bernoulli_distribution;
		//	CHECK: binomial_distribution (Boost 1.47)
		typedef boost::random::binomial_distribution</*tTVInteger*/tjs_int, tTVReal> binomial_distribution;
		typedef boost::random::geometric_distribution<tTVInteger, tTVReal> geometric_distribution;
		typedef boost::random::negative_binomial_distribution<tTVInteger, tTVReal> negative_binomial_distribution;
		typedef boost::random::poisson_distribution<tTVInteger, tTVReal> poisson_distribution;
		typedef boost::random::exponential_distribution<tTVReal> exponential_distribution;
		typedef boost::random::gamma_distribution<tTVReal> gamma_distribution;
		typedef boost::random::weibull_distribution<tTVReal> weibull_distribution;
		typedef boost::random::extreme_value_distribution<tTVReal> extreme_value_distribution;
		typedef boost::random::beta_distribution<tTVReal> beta_distribution;
		typedef boost::random::laplace_distribution<tTVReal> laplace_distribution;
		typedef boost::random::normal_distribution<tTVReal> normal_distribution;
		typedef boost::random::lognormal_distribution<tTVReal> lognormal_distribution;
		typedef boost::random::chi_squared_distribution<tTVReal> chi_squared_distribution;
		typedef boost::random::non_central_chi_squared_distribution<tTVReal> non_central_chi_squared_distribution;
		typedef boost::random::cauchy_distribution<tTVReal> cauchy_distribution;
		typedef boost::random::fisher_f_distribution<tTVReal> fisher_f_distribution;
		typedef boost::random::student_t_distribution<tTVReal> student_t_distribution;
		typedef boost::random::discrete_distribution<tTVInteger, tTVReal> discrete_distribution;
		typedef boost::random::piecewise_constant_distribution<tTVReal, tTVReal> piecewise_constant_distribution;
		typedef boost::random::piecewise_linear_distribution<tTVReal> piecewise_linear_distribution;
		typedef boost::random::triangle_distribution<tTVReal> triangle_distribution;
		typedef boost::random::uniform_on_sphere<tTVReal> uniform_on_sphere;
		typedef boost::mpl::BOOST_PP_CAT(list, KTL_VARIATE_RANDOM_DISTRIBUTION_SEQ_NUM)<
			uniform_smallint,
			uniform_int_distribution,
			uniform_01,
			uniform_real_distribution,
			bernoulli_distribution,
			binomial_distribution,
			geometric_distribution,
			negative_binomial_distribution,
			poisson_distribution,
			exponential_distribution,
			gamma_distribution,
			weibull_distribution,
			extreme_value_distribution,
			beta_distribution,
			laplace_distribution,
			normal_distribution,
			lognormal_distribution,
			chi_squared_distribution,
			non_central_chi_squared_distribution,
			cauchy_distribution,
			fisher_f_distribution,
			student_t_distribution,
			discrete_distribution,
			piecewise_constant_distribution,
			piecewise_linear_distribution,
			triangle_distribution,
			uniform_on_sphere
		> distribution_seq_type;
		typedef boost::mpl::transform<
			engine_seq_type,
			boost::mpl::quote1<engine_wrap>
		>::type wrapped_engine_seq_type;
		typedef boost::mpl::transform<
			distribution_seq_type,
			boost::mpl::quote1<distribution_wrap>
		>::type wrapped_distribution_seq_type;
		typedef boost::variant<
			sprig::random::null_engine<>,
			BOOST_PP_ENUM(KTL_VARIATE_RANDOM_ENGINE_SEQ_NUM, KTL_VARIATE_RANDOM_SEQ_AT_GEN, wrapped_engine_seq_type)
		> engine_variant_type;
		typedef boost::variant<
			sprig::random::null_distribution<RandomResult>,
			BOOST_PP_ENUM(KTL_VARIATE_RANDOM_DISTRIBUTION_SEQ_NUM, KTL_VARIATE_RANDOM_SEQ_AT_GEN, wrapped_distribution_seq_type)
		> distribution_variant_type;
	private:
		BOOST_MPL_ASSERT((boost::mpl::equal_to<boost::mpl::size<engine_seq_type>, boost::mpl::size_t<KTL_VARIATE_RANDOM_ENGINE_SEQ_NUM> >));
		BOOST_MPL_ASSERT((boost::mpl::equal_to<boost::mpl::size<distribution_seq_type>, boost::mpl::size_t<KTL_VARIATE_RANDOM_DISTRIBUTION_SEQ_NUM> >));
	public:
		static sprig::krkr::tjs::intptr_type getBase(iTJSDispatch2* obj) {
			return static_cast<sprig::krkr::tjs::intptr_type>(sprig::krkr::tjs::GetMember(obj, SPRIG_KRKR_TJS_W("base")));
		}
	};
	//
	// RandomEngineBase
	//
	class RandomEngineBase {
	public:
		typedef RandomUtils::engine_variant_type variant_type;
	public:
		static RandomEngineBase* getBase(iTJSDispatch2* obj) {
			return reinterpret_cast<RandomEngineBase*>(RandomUtils::getBase(obj));
		}
	public:
		virtual ~RandomEngineBase() {}
		virtual variant_type variant() const = 0;
	};
	//
	// RandomDistributionBase
	//
	class RandomDistributionBase {
	public:
		typedef RandomUtils::distribution_variant_type variant_type;
	public:
		static RandomDistributionBase* getBase(iTJSDispatch2* obj) {
			return reinterpret_cast<RandomDistributionBase*>(RandomUtils::getBase(obj));
		}
	public:
		virtual ~RandomDistributionBase() {}
		virtual variant_type variant() const = 0;
	};
}	// namespace ktl

#endif	// #ifndef SRC_VARIATE_RANDOM_VARIATE_RANDOM_COMMON_HPP

max_stars_repo_name: bolero-MURAKAMI/KTL

max_issues_repo_name: bolero-MURAKAMI/KTL

max_forks_repo_name: bolero-MURAKAMI/KTL


text:
// Std includes
#include <cmath>
#include <iostream>
// Thirdparties includes
#include <Eigen/Dense>
// Lib includes
#include "Eigen/src/Core/Matrix.h"
#include "s0s/runge_kutta_fehlberg.h"
#include "sl0/spheroid.h"
// Simple includes
#include "flow.h"

// Types
using TypeScalar = double;
template<int Size>
using TypeVector = Eigen::Matrix<TypeScalar, Size, 1>;
// Space
constexpr unsigned int DIM = 3;
using TypeSpaceVector = Eigen::Matrix<TypeScalar, DIM, 1>;
using TypeSpaceMatrix = Eigen::Matrix<TypeScalar, DIM, DIM>;
// Ref and View
template<typename ...Args>
using TypeRef = Eigen::Ref<Args...>;
template<typename ...Args>
using TypeView = Eigen::Map<Args...>;
// Solver
using TypeSolver = s0s::SolverRungeKuttaFehlberg<TypeVector<Eigen::Dynamic>, TypeView>;
// Flow
using TypeFlow = Flow<TypeSpaceVector, TypeSpaceMatrix, TypeRef>;

int main () {
    TypeSpaceVector x0;
    x0 << 0,
          0,
          0;
    TypeSpaceVector p0;
    p0 << 1,
          0,
          0;
    double t0 = 0.0;
    double dt = 1e-2;
    double tEnd = 1.0;
    // Create ellipsoid
    sl0::Spheroid<TypeVector, DIM, TypeView, TypeFlow, TypeSolver> spheroid(std::make_shared<TypeFlow>(), 1.0);
    // Set initial state
    spheroid.sStep->x(spheroid.state.data()) = x0;
    spheroid.sStep->axis(spheroid.state.data()) = p0;
    spheroid.t = t0;
    // Compute
    std::cout << "Spheroid in a simple shear flow : \n";
    std::cout << "\n";
    std::cout << "Orientation after each step : \n";
    for(std::size_t i = 0; i < (tEnd - t0)/dt; i++) {
        spheroid.update(dt);
        std::cout << spheroid.sStep->axis(spheroid.state.data()) << "\n";
    }
    // out
    std::cout << "Spheroid in a simple shear flow : \n";
    std::cout << "\n";
    std::cout << "Final orientation : \n";
    std::cout << spheroid.sStep->axis(spheroid.state.data()) << "\n";
    std::cout << std::endl;
}


text:
//          Copyright Carl Philipp Reh 2009 - 2016.
// 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 FCPPT_DETAIL_MOVE_ITERATOR_IF_RVALUE_HPP_INCLUDED
#define FCPPT_DETAIL_MOVE_ITERATOR_IF_RVALUE_HPP_INCLUDED

#include <fcppt/config/external_begin.hpp>
#include <boost/utility/enable_if.hpp>
#include <iterator>
#include <type_traits>
#include <fcppt/config/external_end.hpp>

namespace fcppt
{
namespace detail
{

template<
	typename Type,
	typename Iterator
>
inline
typename
boost::enable_if<
	std::is_lvalue_reference<
		Type
	>,
	Iterator
>::type
move_iterator_if_rvalue(
	Iterator const &_iterator
)
{
	return
		_iterator;
}

template<
	typename Type,
	typename Iterator
>
inline
typename
boost::disable_if<
	std::is_lvalue_reference<
		Type
	>,
	decltype(
		std::make_move_iterator(
			std::declval<
				Iterator const &
			>()
		)
	)
>::type
move_iterator_if_rvalue(
	Iterator const &_iterator
)
{
	return
		std::make_move_iterator(
			_iterator
		);
}

}
}

#endif

max_stars_repo_name: vinzenz/fcppt

max_issues_repo_name: vinzenz/fcppt

max_forks_repo_name: vinzenz/fcppt


text: #include <boost/log/utility/setup/formatter_parser.hpp>


text:
/* boost random/binomial_distribution.hpp header file
 *
 * Copyright Steven Watanabe 2010
 * 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)
 *
 * See http://www.boost.org for most recent version including documentation.
 *
 * $Id$
 */

#ifndef BOOST_RANDOM_BINOMIAL_DISTRIBUTION_HPP_INCLUDED
#define BOOST_RANDOM_BINOMIAL_DISTRIBUTION_HPP_INCLUDED

#include <boost/config/no_tr1/cmath.hpp>
#include <cstdlib>
#include <iosfwd>

#include <boost/random/detail/config.hpp>
#include <boost/random/uniform_01.hpp>

#include <boost/random/detail/disable_warnings.hpp>

namespace boost {
namespace random {

namespace detail {

template<class RealType>
struct binomial_table {
    static const RealType table[10];
};

template<class RealType>
const RealType binomial_table<RealType>::table[10] = {
    0.08106146679532726,
    0.04134069595540929,
    0.02767792568499834,
    0.02079067210376509,
    0.01664469118982119,
    0.01387612882307075,
    0.01189670994589177,
    0.01041126526197209,
    0.009255462182712733,
    0.008330563433362871
};

}

/**
 * The binomial distribution is an integer valued distribution with
 * two parameters, @c t and @c p.  The values of the distribution
 * are within the range [0,t].
 *
 * The distribution function is
 * \f$\displaystyle P(k) = {t \choose k}p^k(1-p)^{t-k}\f$.
 *
 * The algorithm used is the BTRD algorithm described in
 *
 *  @blockquote
 *  "The generation of binomial random variates", Wolfgang Hormann,
 *  Journal of Statistical Computation and Simulation, Volume 46,
 *  Issue 1 & 2 April 1993 , pages 101 - 110
 *  @endblockquote
 */
template<class IntType = int, class RealType = double>
class binomial_distribution {
public:
    typedef IntType result_type;
    typedef RealType input_type;

    class param_type {
    public:
        typedef binomial_distribution distribution_type;
        /**
         * Construct a param_type object.  @c t and @c p
         * are the parameters of the distribution.
         *
         * Requires: t >=0 && 0 <= p <= 1
         */
        explicit param_type(IntType t_arg = 1, RealType p_arg = RealType (0.5))
          : _t(t_arg), _p(p_arg)
        {}
        /** Returns the @c t parameter of the distribution. */
        IntType t() const { return _t; }
        /** Returns the @c p parameter of the distribution. */
        RealType p() const { return _p; }
#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS
        /** Writes the parameters of the distribution to a @c std::ostream. */
        template<class CharT, class Traits>
        friend std::basic_ostream<CharT,Traits>&
        operator<<(std::basic_ostream<CharT,Traits>& os,
                   const param_type& parm)
        {
            os << parm._p << " " << parm._t;
            return os;
        }

        /** Reads the parameters of the distribution from a @c std::istream. */
        template<class CharT, class Traits>
        friend std::basic_istream<CharT,Traits>&
        operator>>(std::basic_istream<CharT,Traits>& is, param_type& parm)
        {
            is >> parm._p >> std::ws >> parm._t;
            return is;
        }
#endif
        /** Returns true if the parameters have the same values. */
        friend bool operator==(const param_type& lhs, const param_type& rhs)
        {
            return lhs._t == rhs._t && lhs._p == rhs._p;
        }
        /** Returns true if the parameters have different values. */
        friend bool operator!=(const param_type& lhs, const param_type& rhs)
        {
            return !(lhs == rhs);
        }
    private:
        IntType _t;
        RealType _p;
    };

    /**
     * Construct a @c binomial_distribution object. @c t and @c p
     * are the parameters of the distribution.
     *
     * Requires: t >=0 && 0 <= p <= 1
     */
    explicit binomial_distribution(IntType t_arg = 1,
                                   RealType p_arg = RealType(0.5))
      : _t(t_arg), _p(p_arg)
    {
        init();
    }

    /**
     * Construct an @c binomial_distribution object from the
     * parameters.
     */
    explicit binomial_distribution(const param_type& parm)
      : _t(parm.t()), _p(parm.p())
    {
        init();
    }

    /**
     * Returns a random variate distributed according to the
     * binomial distribution.
     */
    template<class URNG>
    IntType operator()(URNG& urng) const
    {
        if(use_inversion()) {
            if(0.5 < _p) {
                return _t - invert(_t, 1-_p, urng);
            } else {
                return invert(_t, _p, urng);
            }
        } else if(0.5 < _p) {
            return _t - generate(urng);
        } else {
            return generate(urng);
        }
    }

    /**
     * Returns a random variate distributed according to the
     * binomial distribution with parameters specified by @c param.
     */
    template<class URNG>
    IntType operator()(URNG& urng, const param_type& parm) const
    {
        return binomial_distribution(parm)(urng);
    }

    /** Returns the @c t parameter of the distribution. */
    IntType t() const { return _t; }
    /** Returns the @c p parameter of the distribution. */
    RealType p() const { return _p; }

    /** Returns the smallest value that the distribution can produce. */
    IntType min BOOST_PREVENT_MACRO_SUBSTITUTION() const { return 0; }
    /** Returns the largest value that the distribution can produce. */
    IntType max BOOST_PREVENT_MACRO_SUBSTITUTION() const { return _t; }

    /** Returns the parameters of the distribution. */
    param_type param() const { return param_type(_t, _p); }
    /** Sets parameters of the distribution. */
    void param(const param_type& parm)
    {
        _t = parm.t();
        _p = parm.p();
        init();
    }

    /**
     * Effects: Subsequent uses of the distribution do not depend
     * on values produced by any engine prior to invoking reset.
     */
    void reset() { }

#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS
    /** Writes the parameters of the distribution to a @c std::ostream. */
    template<class CharT, class Traits>
    friend std::basic_ostream<CharT,Traits>&
    operator<<(std::basic_ostream<CharT,Traits>& os,
               const binomial_distribution& bd)
    {
        os << bd.param();
        return os;
    }

    /** Reads the parameters of the distribution from a @c std::istream. */
    template<class CharT, class Traits>
    friend std::basic_istream<CharT,Traits>&
    operator>>(std::basic_istream<CharT,Traits>& is, binomial_distribution& bd)
    {
        bd.read(is);
        return is;
    }
#endif

    /** Returns true if the two distributions will produce the same
        sequence of values, given equal generators. */
    friend bool operator==(const binomial_distribution& lhs,
                           const binomial_distribution& rhs)
    {
        return lhs._t == rhs._t && lhs._p == rhs._p;
    }
    /** Returns true if the two distributions could produce different
        sequences of values, given equal generators. */
    friend bool operator!=(const binomial_distribution& lhs,
                           const binomial_distribution& rhs)
    {
        return !(lhs == rhs);
    }

private:

    /// @cond show_private

    template<class CharT, class Traits>
    void read(std::basic_istream<CharT, Traits>& is) {
        param_type parm;
        if(is >> parm) {
            param(parm);
        }
    }

    bool use_inversion() const
    {
        // BTRD is safe when np >= 10
        return m < 11;
    }

    // computes the correction factor for the Stirling approximation
    // for log(k!)
    static RealType fc(IntType k)
    {
        if(k < 10) return detail::binomial_table<RealType>::table[k];
        else {
            RealType ikp1 = RealType(1) / (k + 1);
            return (RealType(1)/12
                 - (RealType(1)/360
                 - (RealType(1)/1260)*(ikp1*ikp1))*(ikp1*ikp1))*ikp1;
        }
    }

    void init()
    {
        using std::sqrt;
        using std::pow;

        RealType p = (0.5 < _p)? (1 - _p) : _p;
        IntType t = _t;

        m = static_cast<IntType>((t+1)*p);

        if(use_inversion()) {
            _u.q_n = pow((1 - p), static_cast<RealType>(t));
        } else {
            _u.btrd.r = p/(1-p);
            _u.btrd.nr = (t+1)*_u.btrd.r;
            _u.btrd.npq = t*p*(1-p);
            RealType sqrt_npq = sqrt(_u.btrd.npq);
            _u.btrd.b = 1.15 + 2.53 * sqrt_npq;
            _u.btrd.a = -0.0873 + 0.0248*_u.btrd.b + 0.01*p;
            _u.btrd.c = t*p + 0.5;
            _u.btrd.alpha = (2.83 + 5.1/_u.btrd.b) * sqrt_npq;
            _u.btrd.v_r = 0.92 - 4.2/_u.btrd.b;
            _u.btrd.u_rv_r = 0.86*_u.btrd.v_r;
        }
    }

    template<class URNG>
    result_type generate(URNG& urng) const
    {
        using std::floor;
        using std::abs;
        using std::log;

        while(true) {
            RealType u;
            RealType v = uniform_01<RealType>()(urng);
            if(v <= _u.btrd.u_rv_r) {
                u = v/_u.btrd.v_r - 0.43;
                return static_cast<IntType>(floor(
                    (2*_u.btrd.a/(0.5 - abs(u)) + _u.btrd.b)*u + _u.btrd.c));
            }

            if(v >= _u.btrd.v_r) {
                u = uniform_01<RealType>()(urng) - 0.5;
            } else {
                u = v/_u.btrd.v_r - 0.93;
                u = ((u < 0)? -0.5 : 0.5) - u;
                v = uniform_01<RealType>()(urng) * _u.btrd.v_r;
            }

            RealType us = 0.5 - abs(u);
            IntType k = static_cast<IntType>(floor((2*_u.btrd.a/us + _u.btrd.b)*u + _u.btrd.c));
            if(k < 0 || k > _t) continue;
            v = v*_u.btrd.alpha/(_u.btrd.a/(us*us) + _u.btrd.b);
            RealType km = abs(k - m);
            if(km <= 15) {
                RealType f = 1;
                if(m < k) {
                    IntType i = m;
                    do {
                        ++i;
                        f = f*(_u.btrd.nr/i - _u.btrd.r);
                    } while(i != k);
                } else if(m > k) {
                    IntType i = k;
                    do {
                        ++i;
                        v = v*(_u.btrd.nr/i - _u.btrd.r);
                    } while(i != m);
                }
                if(v <= f) return k;
                else continue;
            } else {
                // final acceptance/rejection
                v = log(v);
                RealType rho =
                    (km/_u.btrd.npq)*(((km/3. + 0.625)*km + 1./6)/_u.btrd.npq + 0.5);
                RealType t = -km*km/(2*_u.btrd.npq);
                if(v < t - rho) return k;
                if(v > t + rho) continue;

                IntType nm = _t - m + 1;
                RealType h = (m + 0.5)*log((m + 1)/(_u.btrd.r*nm))
                           + fc(m) + fc(_t - m);

                IntType nk = _t - k + 1;
                if(v <= h + (_t+1)*log(static_cast<RealType>(nm)/nk)
                          + (k + 0.5)*log(nk*_u.btrd.r/(k+1))
                          - fc(k)
                          - fc(_t - k))
                {
                    return k;
                } else {
                    continue;
                }
            }
        }
    }

    template<class URNG>
    IntType invert(IntType t, RealType p, URNG& urng) const
    {
        RealType q = 1 - p;
        RealType s = p / q;
        RealType a = (t + 1) * s;
        RealType r = _u.q_n;
        RealType u = uniform_01<RealType>()(urng);
        IntType x = 0;
        while(u > r) {
            u = u - r;
            ++x;
            RealType r1 = ((a/x) - s) * r;
            // If r gets too small then the round-off error
            // becomes a problem.  At this point, p(i) is
            // decreasing exponentially, so if we just call
            // it 0, it's close enough.  Note that the
            // minimum value of q_n is about 1e-7, so we
            // may need to be a little careful to make sure that
            // we don't terminate the first time through the loop
            // for float.  (Hence the test that r is decreasing)
            if(r1 < std::numeric_limits<RealType>::epsilon() && r1 < r) {
                break;
            }
            r = r1;
        }
        return x;
    }

    // parameters
    IntType _t;
    RealType _p;

    // common data
    IntType m;

    union {
        // for btrd
        struct {
            RealType r;
            RealType nr;
            RealType npq;
            RealType b;
            RealType a;
            RealType c;
            RealType alpha;
            RealType v_r;
            RealType u_rv_r;
        } btrd;
        // for inversion
        RealType q_n;
    } _u;

    /// @endcond
};

}

// backwards compatibility
using random::binomial_distribution;

}

#include <boost/random/detail/enable_warnings.hpp>

#endif


text:
//==============================================================================
//         Copyright 2003 - 2012   LASMEA UMR 6602 CNRS/Univ. Clermont II
//         Copyright 2009 - 2012   LRI    UMR 8623 CNRS/Univ Paris Sud XI
//
//          Distributed under the Boost Software License, Version 1.0.
//                 See accompanying file LICENSE.txt or copy at
//                     http://www.boost.org/LICENSE_1_0.txt
//==============================================================================
#define NT2_UNIT_MODULE "nt2 boost.simd.arithmetic toolbox - iround2even/simd Mode"

//////////////////////////////////////////////////////////////////////////////
// unit test behavior of boost.simd.arithmetic components in simd mode
//////////////////////////////////////////////////////////////////////////////
/// created by jt the 01/12/2010
///
#include <boost/simd/arithmetic/include/functions/iround2even.hpp>
#include <boost/simd/sdk/simd/io.hpp>
#include <boost/dispatch/meta/as_integer.hpp>
#include <nt2/sdk/unit/tests.hpp>
#include <nt2/sdk/unit/module.hpp>
#include <nt2/sdk/unit/tests/type_expr.hpp>
#include <boost/simd/sdk/simd/native.hpp>
#include <boost/simd/constant/constant.hpp>
#include <boost/simd/include/constants/zero.hpp>
#include <boost/simd/include/constants/one.hpp>
#include <boost/simd/include/constants/two.hpp>
#include <boost/simd/include/constants/inf.hpp>
#include <boost/simd/include/constants/zero.hpp>
#include <boost/simd/include/constants/minf.hpp>
#include <boost/simd/include/constants/mone.hpp>
#include <boost/simd/include/constants/nan.hpp>

NT2_TEST_CASE_TPL ( iround2even_real__1_0,  BOOST_SIMD_SIMD_REAL_TYPES)
{
  using boost::simd::iround2even;
  using boost::simd::tag::iround2even_;
  using boost::simd::native;
  typedef BOOST_SIMD_DEFAULT_EXTENSION  ext_t;
  typedef native<T,ext_t>                  vT;
  typedef typename boost::dispatch::meta::call<iround2even_(vT)>::type r_t;

  // specific values tests
  NT2_TEST_ULP_EQUAL(iround2even(boost::simd::splat<vT>(1.4)), boost::simd::One<r_t>(), 0);
  NT2_TEST_ULP_EQUAL(iround2even(boost::simd::splat<vT>(1.5)), boost::simd::Two<r_t>(), 0);
  NT2_TEST_ULP_EQUAL(iround2even(boost::simd::splat<vT>(1.6)), boost::simd::Two<r_t>(), 0);
  NT2_TEST_ULP_EQUAL(iround2even(boost::simd::splat<vT>(2.5)), boost::simd::Two<r_t>(), 0);
  NT2_TEST_ULP_EQUAL(iround2even(boost::simd::Half<vT>()), boost::simd::Zero<r_t>(), 0);
  NT2_TEST_ULP_EQUAL(iround2even(boost::simd::Inf<vT>()), boost::simd::Inf<r_t>(), 0);
  NT2_TEST_ULP_EQUAL(iround2even(boost::simd::Mhalf<vT>()), boost::simd::Zero<r_t>(), 0);
  NT2_TEST_ULP_EQUAL(iround2even(boost::simd::Minf<vT>()), boost::simd::Minf<r_t>(), 0);
  NT2_TEST_ULP_EQUAL(iround2even(boost::simd::Mone<vT>()), boost::simd::Mone<r_t>(), 0);
  NT2_TEST_ULP_EQUAL(iround2even(boost::simd::Nan<vT>()), boost::simd::Zero<r_t>(), 0);
  NT2_TEST_ULP_EQUAL(iround2even(boost::simd::One<vT>()), boost::simd::One<r_t>(), 0);
  NT2_TEST_ULP_EQUAL(iround2even(boost::simd::Zero<vT>()), boost::simd::Zero<r_t>(), 0);

  for(int i=-100; i < 100; ++i)
 {
   NT2_TEST_ULP_EQUAL(iround2even(boost::simd::splat<vT>(i+0.5)),boost::simd::splat<r_t>(i+((i >= 0)?(i%2):(-i%2))), 0);
 }
} // end of test for floating_


text:
/**
 *  @copyright Copyright 2018 The J-PET Framework Authors. All rights reserved.
 *  Licensed under the Apache License, Version 2.0 (the "License");
 *  you may not use this file except in compliance with the License.
 *  You may find a copy of the License in the LICENCE file.
 *
 *  Unless required by applicable law or agreed to in writing, software
 *  distributed under the License is distributed on an "AS IS" BASIS,
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *  See the License for the specific language governing permissions and
 *  limitations under the License.
 *
 *  @file JPetReaderTest.cpp
 */

#define BOOST_TEST_DYN_LINK
#define BOOST_TEST_MODULE JPetReaderTest

#include "JPetReader/JPetReader.h"
#include "JPetWriter/JPetWriter.h"

#include <TError.h>
#include <TObjString.h>
#include <boost/test/unit_test.hpp>
#include <cstddef>
#include <iostream>
#include <vector>

BOOST_AUTO_TEST_SUITE(JPetReaderTestSuite)

BOOST_AUTO_TEST_CASE(default_constructor)
{
  JPetReader reader;
  BOOST_REQUIRE_EQUAL(std::string(reader.getCurrentEntry().GetName()), std::string("Empty event"));
  BOOST_REQUIRE(!reader.isOpen());
  BOOST_REQUIRE(!reader.nextEntry());
  BOOST_REQUIRE(!reader.firstEntry());
  BOOST_REQUIRE(!reader.lastEntry());
  BOOST_REQUIRE(!reader.nthEntry(0));
  BOOST_REQUIRE(!reader.nthEntry(1));
  BOOST_REQUIRE(!reader.nthEntry(-1));
  BOOST_REQUIRE_EQUAL(reader.getCurrentEntryNumber(), -1);
  BOOST_REQUIRE_EQUAL(reader.getNbOfAllEntries(), 0);
  BOOST_REQUIRE(!reader.getHeaderClone());
  BOOST_REQUIRE(!reader.getObjectFromFile("testObj"));
}

BOOST_AUTO_TEST_CASE(bad_file)
{
  gErrorIgnoreLevel = 6000;
  JPetReader reader;
  BOOST_REQUIRE(!reader.openFileAndLoadData("bad_file.txt", "tree"));
  BOOST_REQUIRE(!reader.isOpen());
  BOOST_REQUIRE_EQUAL(std::string(reader.getCurrentEntry().GetName()), std::string("Empty event"));
  BOOST_REQUIRE(!reader.nextEntry());
  BOOST_REQUIRE(!reader.firstEntry());
  BOOST_REQUIRE(!reader.lastEntry());
  BOOST_REQUIRE(!reader.nthEntry(0));
  BOOST_REQUIRE(!reader.nthEntry(1));
  BOOST_REQUIRE(!reader.nthEntry(-1));
  BOOST_REQUIRE_EQUAL(reader.getCurrentEntryNumber(), -1);
  BOOST_REQUIRE_EQUAL(reader.getNbOfAllEntries(), 0);
  BOOST_REQUIRE(!reader.getHeaderClone());
  BOOST_REQUIRE(!reader.getObjectFromFile("testObj"));
}

BOOST_AUTO_TEST_CASE(good_file_with_constructor)
{
  JPetReader reader("unitTestData/JPetReaderTest/timewindows_v2.root", "tree");
  BOOST_REQUIRE(reader.isOpen());
  BOOST_REQUIRE_EQUAL(std::string(reader.getCurrentEntry().GetName()), std::string("JPetTimeWindow"));
  BOOST_REQUIRE_EQUAL(reader.getCurrentEntryNumber(), 0);
  BOOST_REQUIRE(reader.nextEntry());
  BOOST_REQUIRE_EQUAL(reader.getCurrentEntryNumber(), 1);
  BOOST_REQUIRE(reader.nextEntry());
  BOOST_REQUIRE_EQUAL(reader.getCurrentEntryNumber(), 2);
  BOOST_REQUIRE(reader.firstEntry());
  BOOST_REQUIRE_EQUAL(reader.getCurrentEntryNumber(), 0);
  BOOST_REQUIRE(reader.lastEntry());
  BOOST_REQUIRE_EQUAL(reader.getCurrentEntryNumber(), 9);
  BOOST_REQUIRE(reader.nthEntry(0));
  BOOST_REQUIRE_EQUAL(reader.getCurrentEntryNumber(), 0);
  BOOST_REQUIRE(reader.nthEntry(5));
  BOOST_REQUIRE_EQUAL(reader.getCurrentEntryNumber(), 5);
  BOOST_REQUIRE_EQUAL(reader.getNbOfAllEntries(), 10);
  BOOST_REQUIRE(reader.getHeaderClone());
}

BOOST_AUTO_TEST_CASE(good_file_getObject)
{
  JPetReader reader("unitTestData/JPetReaderTest/timewindows_v2.root", "tree");
  BOOST_REQUIRE(!reader.getObjectFromFile("nonExistentObj"));
  BOOST_REQUIRE(reader.getObjectFromFile("tree"));
}

BOOST_AUTO_TEST_CASE(good_file_openFileAndLoadData)
{
  JPetReader reader;
  BOOST_REQUIRE(reader.openFileAndLoadData("unitTestData/JPetReaderTest/timewindows_v2.root", "tree"));
  BOOST_REQUIRE(reader.isOpen());
  BOOST_REQUIRE_EQUAL(std::string(reader.getCurrentEntry().GetName()), std::string("JPetTimeWindow"));
  BOOST_REQUIRE(reader.firstEntry());
  BOOST_REQUIRE(reader.nextEntry());
  BOOST_REQUIRE(reader.lastEntry());
  BOOST_REQUIRE(reader.nthEntry(0));
  BOOST_REQUIRE(reader.nthEntry(5));
  BOOST_REQUIRE_EQUAL(reader.getCurrentEntryNumber(), 5);
  BOOST_REQUIRE_EQUAL(reader.getNbOfAllEntries(), 10);
  BOOST_REQUIRE(reader.getHeaderClone());
}

BOOST_AUTO_TEST_CASE(good_file_closeFile)
{
  JPetReader reader;
  BOOST_REQUIRE(reader.openFileAndLoadData("unitTestData/JPetReaderTest/timewindows_v2.root", "tree"));
  BOOST_REQUIRE(reader.isOpen());
  reader.closeFile();
  BOOST_REQUIRE(!reader.isOpen());
  BOOST_REQUIRE(std::string(reader.getCurrentEntry().GetName()) == std::string("Empty event"));
  BOOST_REQUIRE(!reader.nextEntry());
  BOOST_REQUIRE(!reader.firstEntry());
  BOOST_REQUIRE(!reader.lastEntry());
  BOOST_REQUIRE(!reader.nthEntry(0));
  BOOST_REQUIRE(!reader.nthEntry(1));
  BOOST_REQUIRE(!reader.nthEntry(-1));
  BOOST_REQUIRE_EQUAL(reader.getCurrentEntryNumber(), -1);
  BOOST_REQUIRE_EQUAL(reader.getNbOfAllEntries(), 0);
  BOOST_REQUIRE(!reader.getHeaderClone());
  BOOST_REQUIRE(!reader.getObjectFromFile("testObj"));
}

BOOST_AUTO_TEST_SUITE_END()

max_stars_repo_name: BlurredChoise/j-pet-framework

max_issues_repo_name: BlurredChoise/j-pet-framework

max_forks_repo_name: BlurredChoise/j-pet-framework


text:
#include "conex/supernodal_solver.h"
#include "conex/clique_ordering.h"

#include <iostream>
#include <map>

#include <Eigen/Dense>

namespace conex {

using Eigen::MatrixXd;
using Eigen::MatrixXi;
using Eigen::VectorXd;
using Eigen::VectorXi;
using Permutation = Eigen::PermutationMatrix<-1>;
using T = TriangularMatrixOperations;
using std::vector;

namespace {

vector<int> Relabel(const vector<int>& x, const vector<int>& labels) {
  vector<int> y;
  for (auto& xi : x) {
    y.push_back(labels.at(xi));
  }
  return y;
}

int GetMax(const std::vector<Clique>& cliques) {
  int max = cliques.at(0).at(0);
  for (const auto& c : cliques) {
    for (const auto ci : c) {
      if (ci > max) {
        max = ci;
      }
    }
  }
  return max;
}

class TriangularMatrixColumnOperations {
 public:
  TriangularMatrixColumnOperations(T::Matrix* mat) : mat_(mat) {
    dense_size = mat_->supernode_size.at(0);
  }
  void NextColumn() {
    supernode_column_++;
    column++;
    dense_size--;
    if (supernode_column_ >= mat_->supernode_size.at(supernode_index_)) {
      supernode_index_++;
      supernode_column_ = 0;
      dense_size = mat_->supernode_size.at(supernode_index_);
    }
  }

  void Rescale(double scale) {
    mat_->supernodes.at(supernode_index_).col(supernode_column_).array() *=
        scale;
    mat_->separator.at(supernode_index_).row(supernode_column_).array() *=
        scale;
  }

  double Diagonal() {
    return mat_->supernodes.at(supernode_index_)(supernode_column_,
                                                 supernode_column_);
  }

  // Applies the operation b = b - L.col(i) * y (i)
  void SubtractWeightedColumn(VectorXd* b, double y) {
    b->segment(column, dense_size) -= y * mat_->supernodes.at(supernode_index_)
                                              .col(supernode_column_)
                                              .tail(dense_size);

    int offset = mat_->supernode_size.at(supernode_index_);
    for (int i = 0; i < mat_->separator.at(supernode_index_).cols(); i++) {
      int row = mat_->path.at(supernode_index_).at(offset + i);
      (*b)(row) -=
          y * mat_->separator.at(supernode_index_)(supernode_column_, i);
    }
  }

  std::vector<int> NonzeroRows() {
    std::vector<int> y;
    for (int i = 1; i < dense_size; i++) {
      y.push_back(column + i);
    }
    int offset = mat_->supernode_size.at(supernode_index_);
    for (int i = 0; i < mat_->separator.at(supernode_index_).cols(); i++) {
      int row = mat_->path.at(supernode_index_).at(offset + i);
      y.push_back(row);
    }
    return y;
  }

  // double DotProduct(const Eigen::MatrixXd& x) {
  //  double y = 0;
  //  mat_->supernodes.at(supernode_index_).col(supernode_column_).tail(dense_size).dot(
  //      x.segment(column, dense_size);
  //}

 private:
  int supernode_index_ = 0;
  int dense_size = 0;
  int column = 0;
  int supernode_column_ = 0;
  SparseTriangularMatrix* mat_;
};

template <typename T>
int LookupSuperNode(const T& o, int index, int start) {
  for (int j = static_cast<int>(o.snodes.size()) - 1; j >= 0; j--) {
    if (o.snodes.at(j).size() > 0) {
      if (index >= o.snodes.at(j).at(0)) {
        return j;
      }
    }
  }
  throw "Sparse matrix is malformed: invalid supernode partition.";
}

double Get(const SparseTriangularMatrix& o, int i, int j) {
  if (j > i) {
    return 0;
  }
  // Find the supernode that owns the column.
  int node = LookupSuperNode(o, j, 0);

  // Apply offsets.
  int offset_i = i - o.path.at(node).at(0);
  int offset_j = j - o.path.at(node).at(0);
  if ((offset_i < o.supernode_size.at(node)) &&
      (offset_j < o.supernode_size.at(node))) {
    return o.supernodes.at(node)(offset_i, offset_j);
  }
  for (size_t k = o.supernode_size.at(node); k < o.path.at(node).size(); k++) {
    if (i == o.path.at(node).at(k)) {
      return o.separator.at(node)(offset_j, k - o.supernode_size.at(node));
    }
  }
  return 0;
}

class LowerTriangularSuperNodal {
  // We partition matrix into
  //
  // SN1  R1
  // R1   SN2  R2
  // R1   R2   SN3
  //
  // Where R1 is sparse.
 public:
  LowerTriangularSuperNodal(SparseTriangularMatrix* mat)
      : supernodes_(mat->supernodes),
        supernode_size_(mat->supernode_size),
        path_(mat->path),
        separator_(mat->separator) {
    Init();
  }

  void Init() {
    // for (int j = static_cast<int>(path_.size()) - 1; j >= 0; j--) {
    //   supernodes_.at(j).resize(supernode_size_.at(j), supernode_size_.at(j));
    //   separator_.at(j).resize(supernode_size_.at(j),  path_.at(j).size() -
    //   supernode_size_.at(j));
    // }
  }

  // The smallest super-node less than i.
  int LookupSuperNode(int index, int start) const {
    for (int j = static_cast<int>(path_.size()) - 1; j >= 0; j--) {
      if (index >= path_.at(j).at(0)) {
        return j;
      }
    }
    throw "Sparse matrix is malformed: invalid supernode partition.";
  }

  double Get(int i, int j) const {
    if (j > i) {
      return 0;
    }
    // Find the supernode that owns the column.
    int node = LookupSuperNode(j, 0);

    // Apply offsets.
    int offset_i = i - path_.at(node).at(0);
    int offset_j = j - path_.at(node).at(0);
    if ((offset_i < supernode_size_.at(node)) &&
        (offset_j < supernode_size_.at(node))) {
      return supernodes_.at(node)(offset_i, offset_j);
    }
    for (size_t k = supernode_size_.at(node); k < path_.at(node).size(); k++) {
      if (i == path_.at(node).at(k)) {
        return separator_.at(node)(offset_j, k - supernode_size_.at(node));
      }
    }
    return 0;
  }

  void Increment(double val, int i, int j, int supernode_index = -1) {
    if (val == 0) {
      return;
    }
    assert(j <= i);

    int node = supernode_index;
    // Find the supernode that owns storage for column j.
    if (node == -1) {
      node = LookupSuperNode(j, 0);
    }

    // Apply offsets.
    int offset_i = i - path_.at(node).at(0);
    int offset_j = j - path_.at(node).at(0);

    // i is also in the supernode.
    if ((offset_i < supernode_size_.at(node)) &&
        (offset_j < supernode_size_.at(node))) {
      supernodes_.at(node)(offset_i, offset_j) += val;
      return;
    }

    // i is in a separator.
    for (size_t k = supernode_size_.at(node); k < path_.at(node).size(); k++) {
      if (i == path_.at(node).at(k)) {
        separator_.at(node)(offset_j, k - supernode_size_.at(node)) += val;
        return;
      }
    }

    throw "Specified entry of sparse matrix is not accesible.";
  }

 public:
  using MapT = Eigen::Map<Eigen::MatrixXd>;
  std::vector<MapT>& supernodes_;
  std::vector<int>& supernode_size_;
  std::vector<Clique>& path_;
  std::vector<MapT>& separator_;
};

std::vector<int> ResidualSize(std::vector<Clique>& path) {
  std::vector<int> y;
  for (size_t j = 0; j < path.size() - 1; j++) {
    std::vector<int> temp;
    IntersectionOfSorted(path.at(j), path.at(j + 1), &temp);
    y.push_back(path.at(j).size() - temp.size());
  }
  y.push_back(path.back().size());
  return y;
}

}  // namespace

void IntersectionOfSorted(const std::vector<int>& v1,
                          const std::vector<int>& v2, std::vector<int>* v3) {
  v3->clear();
  std::set_intersection(v1.begin(), v1.end(), v2.begin(), v2.end(),
                        back_inserter(*v3));
}

SparseTriangularMatrix GetFillInPattern(
    int N, const std::vector<Clique>& cliques_input) {
  auto mat = MakeSparseTriangularMatrix(N, cliques_input);

  for (int j = static_cast<int>(mat.path.size()) - 1; j >= 0; j--) {
    // Initialize columns of super nodes.
    mat.supernodes.at(j).setConstant(1);
    mat.separator.at(j).setConstant(1);

    // Update other columns: the (seperator, seperator) components.
    int index = 0;
    auto s_s = mat.workspace_.seperator_diagonal.at(j);
    int n = mat.path.at(j).size();
    for (int i = mat.supernode_size.at(j); i < n; i++) {
      for (int k = i; k < n; k++) {
        *s_s.at(index++) += 1;
      }
    }
  }
  return mat;
}

std::vector<Clique> Permute(std::vector<Clique>& path,
                            std::vector<int>& permutation) {
  auto y = path;
  for (size_t i = 0; i < path.size(); i++) {
    for (size_t j = 0; j < path.at(i).size(); j++) {
      y.at(i).at(j) = permutation.at(path.at(i).at(j));
    }
  }
  return y;
}

void Sort(std::vector<Clique>* path) {
  for (size_t i = 0; i < path->size(); i++) {
    std::sort(path->at(i).begin(), path->at(i).end());
  }
}

// Want A B  C  D  E
//   Add A_i \cap A_k
//    to all A_j for i < j < k.
//
void RunningIntersectionClosure(std::vector<Clique>* path) {
  if (path->size() < 2) {
    return;
  }
  int n = path->size();
  for (int i = 0; i < n - 2; i++) {
    for (int j = n - 1; j > i + 1; j--) {
      std::vector<int> temp;
      IntersectionOfSorted(path->at(i), path->at(j), &temp);
      if (temp.size() == 0) {
        continue;
      }
      for (int k = j - 1; k > i; k--) {
        path->at(k) = UnionOfSorted(path->at(k), temp);
      }
    }
  }
}

Eigen::MatrixXd TriangularMatrixOperations::ToDense(
    const SparseTriangularMatrix& mat) {
  MatrixXd y(mat.N, mat.N);
  for (int i = 0; i < mat.N; i++) {
    for (int j = 0; j < mat.N; j++) {
      y(i, j) = Get(mat, i, j);
    }
  }
  return y;
}

void TriangularMatrixOperations::SetConstant(SparseTriangularMatrix* mat,
                                             double val) {
  for (auto& n : mat->supernodes) {
    n.array() = val;
  }
  for (auto& n : mat->separator) {
    n.array() = val;
  }
}

SparseTriangularMatrix MakeSparseTriangularMatrix(
    int N, const std::vector<Clique>& path_) {
  auto path = path_;
  Sort(&path);
  RunningIntersectionClosure(&path);
  auto supernode_size = ResidualSize(path);
  return SparseTriangularMatrix(N, path, supernode_size);
}

void T::CholeskyInPlace(SparseTriangularMatrix* C) {
  TriangularMatrixColumnOperations col(C);
  LowerTriangularSuperNodal mat(C);
  double sqrt_d = std::sqrt(col.Diagonal());
  col.Rescale(1.0 / sqrt_d);

  // i: a supernode
  for (int i = 0; i < C->N - 1; i++) {
    // Substract col(2:n) c(2:n)^T
    auto indices = col.NonzeroRows();
    for (size_t k = 0; k < indices.size(); k++) {
      double weight = -mat.Get(indices.at(k), i);
      // Decrement column k
      for (size_t j = k; j < indices.size(); j++) {
        mat.Increment(weight * mat.Get(indices.at(j), i), indices.at(j),
                      indices.at(k));
      }
    }
    col.NextColumn();
    double sqrt_d = std::sqrt(col.Diagonal());
    col.Rescale(1.0 / sqrt_d);
  }
}

// Apply L  inverse.
// 1
// 1  1
// 1  1  1

// L
// B in
VectorXd T::ApplyInverseOfTranspose(SparseTriangularMatrix* mat,
                                    const VectorXd& b) {
  assert(b.rows() == mat->N);
  int n = b.rows();
  VectorXd y(n);
  auto res = b;
  LowerTriangularSuperNodal L(mat);

  y(n - 1) = res(n - 1) / L.Get(n - 1, n - 1);
  for (int i = n - 2; i >= 0; i--) {
    for (int j = 0; j < i + 1; j++) {
      res(j) = res(j) - L.Get(i + 1, j) * y(i + 1);
    }
    y(i) = res(i) / L.Get(i, i);
  }
  return y;
}

VectorXd T::ApplyInverse(SparseTriangularMatrix* mat, const VectorXd& b) {
  assert(b.rows() == mat->N);
  int n = b.rows();
  VectorXd y(n);
  auto res = b;
  LowerTriangularSuperNodal L(mat);

  TriangularMatrixColumnOperations col(mat);
  y(0) = res(0) / col.Diagonal();
  for (int i = 1; i < n; i++) {
    // Iterate over non-zero entries of column.
    // for (int j = i; j < n; j++) {
    //   res(j) = res(j) - L.Get(j, i - 1) * y(i - 1);
    // }
    col.SubtractWeightedColumn(&res, y(i - 1));

    col.NextColumn();
    y(i) = res(i) / col.Diagonal();
  }
  return y;
}

std::vector<int> UnionOfSorted(const std::vector<int>& x1,
                               const std::vector<int>& x2) {
  std::vector<int> y;
  set_union(x1.begin(), x1.end(), x2.begin(), x2.end(), inserter(y, y.end()));
  return y;
}

MatrixData GetData(const vector<vector<int>>& cliques, int init) {
  vector<vector<int>> separators;
  vector<vector<int>> supernodes;
  MatrixData d;
  d.cliques = cliques;
  Sort(&d.cliques);
  auto& order = d.clique_order;

  PickCliqueOrder(d.cliques, init, &order, &supernodes, &separators);

  d.permutation.resize(GetMax(cliques) + 1);
  d.permutation_inverse.resize(GetMax(cliques) + 1);
  int i = 0;
  for (auto& e : order) {
    for (auto& sn_ii : supernodes.at(e)) {
      d.permutation_inverse.at(i) = sn_ii;
      d.permutation.at(sn_ii) = i;
      i++;
    }
  }

  auto& separators_ = d.separators;
  auto& supernodes_ = d.supernodes;

  supernodes_.resize(order.size());
  separators_.resize(order.size());

  i = 0;

  for (auto e : order) {
    supernodes_.at(i) = Relabel(supernodes.at(e), d.permutation);
    separators_.at(i) = Relabel(separators.at(e), d.permutation);
    i++;
  }

  Sort(&separators_);
  Sort(&supernodes_);

  int cnt = 0;
  auto& supernode_size = d.supernode_size;
  supernode_size.resize(order.size());
  for (auto& si : supernode_size) {
    si = supernodes_.at(cnt).size();
    cnt++;
  }
  d.N = std::accumulate(supernode_size.begin(), supernode_size.end(), 0);
  for (size_t i = 0; i < d.cliques.size(); i++) {
    d.cliques.at(i) = UnionOfSorted(supernodes_.at(i), separators_.at(i));
  }
  return d;
}

}  // namespace conex

ext: cc

max_stars_repo_name: ToyotaResearchInstitute/conex

max_issues_repo_name: ToyotaResearchInstitute/conex

max_forks_repo_name: ToyotaResearchInstitute/conex