text:
/*------------------------------------------------------*/
/* Biological Engineering: Problem Set 03               */
/* main program                                         */
/*                                                      */
/* file name    :   BiologicalEngineering_set3.cpp      */
/* compiler     :   gcc                                 */
/* Student ID   :   1526084                             */
/* author       :   Takaharu Nakajima                   */
/* date         :   2016.11.11                          */
/* memo         :   $B@8BN9)3XBh;02sL\$N(B1.1$B$+$i(B1.3$B$N(B      */
/*                  $B%W%m%0%i%`$r$^$H$a$?$b$N$K$J$j$^$9(B  */
/*                                                      */
/*                  $B0lItJQ?t$NDj5A$N$?$aJQ99$7$?ItJ,(B    */
/*                  ($BNc(B)                                */
/*                  RT0 -> RTa                          */
/*                  $B$3$N$h$&$K(B0, 1...E$B$r(Ba, b...E$B$K(B      */
/*                  $BJQ99$7$^$7$?(B                        */
/*                                                      */
/*------------------------------------------------------*/

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

using namespace std;
using namespace Eigen;

#define PI 3.1415926

double deg2rad(double degree)
{
    double radian = degree * PI / 180.0f;
    return radian;
}

int main()
{
    double deg1, deg2, deg3;
    double rad1, rad2, rad3;
    Matrix<double,4,4> RTa, aTb, bTc, cTe, eTE;
    Matrix<double,4,4> RTb, RTc, RTe, RTE;
    Vector3d RZa, RZb, RZc, RZe, RZE;
    Vector3d Rpa, Rpb, Rpc, Rpe, RpE;
    Vector3d RpEa, RpEb, RpEc, RpEe, RpEE;

    /* Input */
    cout << "Angle input" << endl;
    cout << "deg1= ";
    cin >> deg1;
    cout << "deg2= ";
    cin >> deg2;
    cout << "deg3= ";
    cin >> deg3;

    /* Convert to radians */
    rad1 = deg2rad(deg1);
    rad2 = deg2rad(deg2);
    rad3 = deg2rad(deg3);

    /* Link parameters */
    RTa << 1, 0, 0, 0,
           0, 1, 0, 0,
           0, 0, 1, 2,
           0, 0, 0, 1;
    aTb << cos(rad1), -sin(rad1), 0, 0,
           sin(rad1),  cos(rad1), 0, 0,
                   0,          0, 1, 0,
                   0,          0, 0, 1;
    bTc << cos(rad2), -sin(rad2), 0, 0,
                   0,          0, -1, 0,
           sin(rad2),  cos(rad2), 0, 0,
                   0,          0, 0, 1;
    cTe << cos(rad3), -sin(rad3), 0, 2,
           sin(rad3),  cos(rad3), 0, 0,
                   0,          0, 1, 0,
                   0,          0, 0, 1;
    eTE << 1, 0, 0, 2,
           0, 1, 0, 0,
           0, 0, 1, 0,
           0, 0, 0, 1;

    /* Matrix calculation */
    RTb = RTa * aTb;
    RTc = RTb * bTc;
    RTe = RTc * cTe;
    RTE = RTe * eTE;

    RZa = RTa.block<3,1>(0,2);
    RZb = RTb.block<3,1>(0,2);
    RZc = RTc.block<3,1>(0,2);
    RZe = RTe.block<3,1>(0,2);
    RZE = RTE.block<3,1>(0,2);

    Rpa = RTa.block<3,1>(0,3);
    Rpb = RTb.block<3,1>(0,3);
    Rpc = RTc.block<3,1>(0,3);
    Rpe = RTe.block<3,1>(0,3);
    RpE = RTE.block<3,1>(0,3);

    RpEa = RpE - Rpa;
    RpEb = RpE - Rpb;
    RpEc = RpE - Rpc;
    RpEe = RpE - Rpe;
    RpEE = RpE - RpE;

    /* Calculation result display */
    /* RT0 ~ 3TE */
    cout << "RT0=" << endl;
    cout << RTa << endl;
    cout << "0T1=" << endl;
    cout << aTb << endl;
    cout << "1T2=" << endl;
    cout << bTc << endl;
    cout << "2T3=" << endl;
    cout << cTe << endl;
    cout << "3TE=" << endl;
    cout << eTE << endl << endl;

    /* RZ0 ~ RZE */
    cout << "RZ0=" << endl;
    cout << RZa << endl;
    cout << "RZ1=" << endl;
    cout << RZb << endl;
    cout << "RZ2=" << endl;
    cout << RZc << endl;
    cout << "RZ3=" << endl;
    cout << RZe << endl;
    cout << "RZE=" << endl;
    cout << RZE << endl;

    /* Rp0 ~ Rp1 */
    cout << "Rp0=" << endl;
    cout << Rpa << endl;
    cout << "Rp1=" << endl;
    cout << Rpb << endl;
    cout << "Rp2=" << endl;
    cout << Rpc << endl;
    cout << "Rp3=" << endl;
    cout << Rpe << endl;
    cout << "RpE=" << endl;
    cout << RpE << endl;

    /* RZ0 x RpE,0 ~ RZE x RpE,E */
    cout << "RZ0 x RpE,0=" << endl;
    cout << RZa.cross(RpEa) << endl;
    cout << "RZ1 x RpE,1=" << endl;
    cout << RZb.cross(RpEb) << endl;
    cout << "RZ2 x RpE,2=" << endl;
    cout << RZc.cross(RpEc) << endl;
    cout << "RZ3 x RpE,3=" << endl;
    cout << RZe.cross(RpEe) << endl;
    cout << "RZE x RpE,E=" << endl;
    cout << RZE.cross(RpEE) << endl;

    return 0;
}


text:
#include <ros/ros.h>
#include <decision_making/SynchCout.h>
#include <decision_making/BT.h>
#include <decision_making/FSM.h>
#include <decision_making/ROSTask.h>
#include <decision_making/DecisionMaking.h>

using namespace std;
using namespace decision_making;

EventQueue* mainEventQueue;
struct MainEventQueue{
	MainEventQueue(){ mainEventQueue = new RosEventQueue();	}
	~MainEventQueue(){ delete mainEventQueue; }
};

FSM(Turnstile)
{
	FSM_STATES
	{
		Locked,
		Unlocked
	}
	FSM_START(Locked);
	FSM_BGN
	{
		FSM_STATE(Locked)
		{
			FSM_TRANSITIONS
			{
				FSM_ON_EVENT(/COIN, FSM_NEXT(Unlocked));
				FSM_ON_EVENT(/PUSH, FSM_NEXT(Locked));
			}
		}
		FSM_STATE(Unlocked)
		{
			FSM_TRANSITIONS
			{
				FSM_ON_EVENT(/COIN, FSM_NEXT(Unlocked));
				FSM_ON_EVENT(/PUSH, FSM_NEXT(Locked));
			}
		}
	}
	FSM_END
}

void run_fsm(){
	FsmTurnstile(NULL, mainEventQueue);
}

void EVENTS_GENERATOR(){
	Event spec[]={"COIN","COIN","PUSH","PUSH", "REPEAT"};
	int i=0;
	boost::this_thread::sleep(boost::posix_time::seconds(1));
	while(true and ros::ok()){
		Event t = spec[i];
		if(t == "REPEAT"){ i=1; t=spec[0]; }else i++;
		cout << endl << t<<" -> ";
		mainEventQueue->riseEvent(t);
		boost::this_thread::sleep(boost::posix_time::seconds(1));
	}
	mainEventQueue->close();
}

#include <boost/filesystem.hpp>

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

	ros::init(argc, argv, "Turnstile");
	MainEventQueue meq;
	ros_decision_making_init(argc, argv);

	boost::thread_group threads;
	threads.add_thread(new boost::thread(boost::bind(&run_fsm)));
	threads.add_thread(new boost::thread(boost::bind(&EVENTS_GENERATOR)));

	ros::spin();
	threads.join_all();

	return 0;
}

ext: cpp


text:
/**********************************************************************
 *  Copyright (c) 2008-2015, Alliance for Sustainable Energy.
 *  All rights reserved.
 *
 *  This library is free software; you can redistribute it and/or
 *  modify it under the terms of the GNU Lesser General Public
 *  License as published by the Free Software Foundation; either
 *  version 2.1 of the License, or (at your option) any later version.
 *
 *  This library is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 *  Lesser General Public License for more details.
 *
 *  You should have received a copy of the GNU Lesser General Public
 *  License along with this library; if not, write to the Free Software
 *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
 **********************************************************************/

#ifndef MODEL_ATTRIBUTES_HPP
#define MODEL_ATTRIBUTES_HPP

#include <boost/preprocessor/control/if.hpp>
#include <boost/optional.hpp>
#include "../utilities/units/Quantity.hpp"
#include "../utilities/units/OSOptionalQuantity.hpp"

//ATTRIBUTE_DEFINITION(IS_OPTIONAL,IS_AUTOSIZEABLE,HAS_DEFAULT,ATTRIBUTE_NAME,ACCESSOR_BASE_NAME)
//ATTRIBUTE_IMPLEMENTATION(IS_OPTIONAL,HAS_DEFAULT,IS_AUTOSIZEABLE,
//                         ATTRIBUTE_NAME,ACCESSOR_BASE_NAME,
//                         CLASS_NAME,
//                         DEFAULT_IMPLEMENTATIOIN,IDD_OBJECT_NAME,IDD_FIELD_NAME)

// ATTRIBUTE_DEFINITION

#define ATTRIBUTE_DEFINITION(IS_OPTIONAL,HAS_DEFAULT,IS_AUTOSIZEABLE,ATTRIBUTE_NAME,ACCESSOR_BASE_NAME) \
\
BOOST_PP_IF(BOOST_PP_OR(IS_OPTIONAL,IS_AUTOSIZEABLE), \
\
bool set##ACCESSOR_BASE_NAME(boost::optional<double> value); \
bool set##ACCESSOR_BASE_NAME(const OSOptionalQuantity & value); \
OSOptionalQuantity ATTRIBUTE_NAME##_SI() const; \
OSOptionalQuantity ATTRIBUTE_NAME##_IP() const; \
OSOptionalQuantity get##ACCESSOR_BASE_NAME(bool returnIP = false) const; \
boost::optional<double> ATTRIBUTE_NAME() const; \
BOOST_PP_IF(HAS_DEFAULT, \
void reset##ACCESSOR_BASE_NAME(); \
bool is##ACCESSOR_BASE_NAME##Defaulted() const; \
,) \
BOOST_PP_IF(IS_AUTOSIZEABLE, \
void autosize##ACCESSOR_BASE_NAME(); \
bool is##ACCESSOR_BASE_NAME##Autosized() const; \
,) \
,\
bool set##ACCESSOR_BASE_NAME(double value); \
bool set##ACCESSOR_BASE_NAME(const Quantity & value); \
Quantity ATTRIBUTE_NAME##_SI() const; \
Quantity ATTRIBUTE_NAME##_IP() const; \
Quantity get##ACCESSOR_BASE_NAME(bool returnIP = false) const; \
double ATTRIBUTE_NAME() const; \
BOOST_PP_IF(HAS_DEFAULT, \
void reset##ACCESSOR_BASE_NAME(); \
bool is##ACCESSOR_BASE_NAME##Defaulted() const; \
,) \
)

// ATTRIBUTE_IMPLEMENTATION

#define ATTRIBUTE_IMPLEMENTATION(IS_OPTIONAL,HAS_DEFAULT,IS_AUTOSIZEABLE, \
                                 ATTRIBUTE_NAME,ACCESSOR_BASE_NAME, \
                                 CLASS_NAME, \
                                 DEFAULT_IMPLEMENTATIOIN,IDD_OBJECT_NAME,IDD_FIELD_NAME) \
\
BOOST_PP_IF(BOOST_PP_OR(IS_OPTIONAL,IS_AUTOSIZEABLE), \
\
OSOptionalQuantity CLASS_NAME##_Impl::ATTRIBUTE_NAME##_SI() const { \
  return get##ACCESSOR_BASE_NAME(false); \
} \
\
OSOptionalQuantity CLASS_NAME##_Impl::ATTRIBUTE_NAME##_IP() const { \
  return get##ACCESSOR_BASE_NAME(true); \
} \
\
OSOptionalQuantity CLASS_NAME##_Impl::get##ACCESSOR_BASE_NAME(bool returnIP) const { \
  OptionalDouble value = ATTRIBUTE_NAME(); \
  return getQuantityFromDouble(IDD_OBJECT_NAME##Fields::IDD_FIELD_NAME, value, returnIP); \
} \
\
bool CLASS_NAME##_Impl::set##ACCESSOR_BASE_NAME(const OSOptionalQuantity& quantity) { \
  bool result(false); \
  OptionalDouble value; \
  if (quantity.isSet()) { \
    value = getDoubleFromQuantity(IDD_OBJECT_NAME##Fields::IDD_FIELD_NAME,quantity.get()); \
    if (value) { \
      result = set##ACCESSOR_BASE_NAME(value); \
    } \
  } \
  else { \
    result = set##ACCESSOR_BASE_NAME(value); \
  } \
  return result; \
} \
\
BOOST_PP_IF(DEFAULT_IMPLEMENTATIOIN, \
void CLASS_NAME##_Impl::set##ACCESSOR_BASE_NAME(boost::optional<double> value) { \
  bool result = false; \
  if (value) { \
    result = setDouble(IDD_OBJECT_NAME##Fields::IDD_FIELD_NAME, value.get()); \
  } else { \
    result = setString(IDD_FIELD_NAME##Fields::IDD_FIELD_NAME, ""); \
  } \
  OS_ASSERT(result); \
} \
\
boost::optional<double> CLASS_NAME##_Impl::ATTRIBUTE_NAME() const { \
  OptionalDouble value = getDouble(IDD_OBJECT_NAME##Fields::IDD_FIELD_NAME, true); \
  return value; \
} \
,)\
\
BOOST_PP_IF(HAS_DEFAULT, \
void CLASS_NAME##_Impl::reset##ATTRIBUTE_NAME() { \
  bool result = setString(IDD_OBJECT_NAME##Fields::IDD_FIELD_NAME, ""); \
  OS_ASSERT(result); \
} \
,) \
BOOST_PP_IF(IS_AUTOSIZEABLE, \
bool CLASS_NAME##_Impl::is##ACCESSOR_BASE_NAME##Autosized() const { \
  bool result = false; \
  boost::optional<std::string> value = getString(IDD_OBJECT_NAME##Fields::IDD_FIELD_NAME, true); \
  if (value) { \
    result = openstudio::istringEqual(value.get(), "autosize"); \
  } \
  return result; \
} \
void CLASS_NAME##_Impl::autosize##ACCESSOR_BASE_NAME() { \
  bool result = setString(IDD_OBJECT_NAME##Fields::IDD_FIELD_NAME, "autosize"); \
  OS_ASSERT(result); \
} \
,) \
,\
Quantity CLASS_NAME##_Impl::get##ACCESSOR_BASE_NAME(bool returnIP) const { \
  OptionalDouble value = ATTRIBUTE_NAME(); \
  OSOptionalQuantity result = getQuantityFromDouble(IDD_OBJECT_NAME##Fields::IDD_FIELD_NAME, value, returnIP); \
  OS_ASSERT(result.isSet()); \
  return result.get(); \
} \
\
bool CLASS_NAME##_Impl::set##ACCESSOR_BASE_NAME(const Quantity& quantity) { \
  OptionalDouble value = getDoubleFromQuantity(IDD_OBJECT_NAME##Fields::IDD_FIELD_NAME,quantity); \
  if (!value) { \
    return false; \
  } \
  return set##ACCESSOR_BASE_NAME(value.get()); \
} \
\
Quantity CLASS_NAME##_Impl::ATTRIBUTE_NAME##_SI() const { \
  return get##ACCESSOR_BASE_NAME(false); \
} \
\
Quantity CLASS_NAME##_Impl::ATTRIBUTE_NAME##_IP() const { \
  return get##ACCESSOR_BASE_NAME(true); \
} \
\
BOOST_PP_IF(DEFAULT_IMPLEMENTATIOIN, \
double CLASS_NAME##_Impl::ATTRIBUTE_NAME() const { \
  OptionalDouble value = getDouble(IDD_OBJECT_NAME##Fields::IDD_FIELD_NAME, true); \
  OS_ASSERT(value); \
  return value.get(); \
} \
\
bool CLASS_NAME##_Impl::set##ACCESSOR_BASE_NAME(double value) { \
  bool result = false; \
  result = setDouble(IDD_OBJECT_NAME##Fields::IDD_FIELD_NAME, value); \
  return result; \
} \
, ) \
\
BOOST_PP_IF(HAS_DEFAULT, \
void CLASS_NAME##_Impl::reset##ACCESSOR_BASE_NAME() { \
  bool result = setString(IDD_OBJECT_NAME##Fields::IDD_FIELD_NAME, ""); \
  OS_ASSERT(result); \
} \
bool CLASS_NAME##_Impl::is##ACCESSOR_BASE_NAME##Defaulted() const { \
  return isEmpty(IDD_OBJECT_NAME##Fields::IDD_FIELD_NAME); \
} \
,) \
)

#endif // MODEL_ATTRIBUTES_HPP

ext: hpp


text:
#include "minicap.h"
#include <sstream>

#include <boost/asio.hpp>
#include <boost/array.hpp>

#include "opencv2/imgcodecs.hpp"

#include <cv_bridge/cv_bridge.h>

Minicap::Minicap(
        ros::NodeHandle &node_handle,
        int width, int height,
        int downscale_factor,
        const std::string &adb_path, const std::string &dev_dir,
        int port) :
        m_width(width),
        m_height(height),
        m_downscale_factor(downscale_factor),
        m_adb_path(adb_path),
        m_dev_dir(dev_dir),
        m_port(port),
        m_image_pub(node_handle.advertise<sensor_msgs::Image>("openstf_screen", 10, true)),
        m_background_thread(&Minicap::run, this) {

}

Minicap::~Minicap() {
    //TODO: make thread-safe
    if (m_adb_minicap)
        m_adb_minicap->kill();

    m_background_thread.join();
}

void Minicap::run() {
    if (m_downscale_factor < 1) {
        ROS_WARN_STREAM("downscale_factor = " << m_downscale_factor << " < 1, setting downscale_factor to 1");
        m_downscale_factor = 1;
    }

    std::stringstream minicap_args;
    minicap_args << "LD_LIBRARY_PATH=" << m_dev_dir << " " << m_dev_dir << "/minicap ";
    minicap_args << "-P " << m_height << "x" << m_height << "@" << m_height / m_downscale_factor << "x"
                 << m_height / m_downscale_factor << "/0 ";
    minicap_args << "-S -Q 100";

    ROS_INFO_STREAM("Minicap args: " << minicap_args.str());

    m_adb_minicap = std::unique_ptr<TinyProcessLib::Process>(new TinyProcessLib::Process(
            std::vector<std::string>{m_adb_path, "shell", minicap_args.str()}, "", [&](const char *bytes, size_t n) {
                ROS_INFO_STREAM(std::string(bytes, n));
            }));

    // reading loop
    read();

    auto exit_status = m_adb_minicap->get_exit_status();
    ROS_INFO_STREAM("Minicap exit code: " << exit_status);
}

void Minicap::read() {
    try {
        boost::system::error_code ec;
        boost::asio::io_service io_service;

        boost::asio::ip::tcp::socket minicap_socket(io_service);

        for (int i = 0; i < 10; i++) {
            ROS_INFO_STREAM("Connecting to minicap..");
            // wait for startup
            sleep(2);
            minicap_socket.connect(
                    boost::asio::ip::tcp::endpoint(boost::asio::ip::address::from_string("127.0.0.1"), m_port),
                    ec);
            if (ec) {
                ROS_ERROR_STREAM(ec.message());
            } else {
                break;
            }
        }
        if (ec) {
            ROS_ERROR_STREAM(ec.message());
            return;
        }

        std::vector<char> buffer(8 * 1024 * 1024);

        // read banner
        boost::asio::read(minicap_socket, boost::asio::buffer(buffer), boost::asio::transfer_exactly(24), ec);
        if (ec) {
            ROS_ERROR_STREAM(ec.message());
            return;
        }

        ROS_INFO_STREAM(int(*reinterpret_cast<uint8_t *>(&buffer.data()[1])));
        ROS_INFO_STREAM("Real display width in pixels (minicap): " << *reinterpret_cast<uint32_t *>(&buffer.data()[6]));
        ROS_INFO_STREAM(
                "Real display height in pixels (minicap): " << *reinterpret_cast<uint32_t *>(&buffer.data()[10]));
        ROS_INFO_STREAM(
                "Virtual display width in pixels (minicap): " << *reinterpret_cast<uint32_t *>(&buffer.data()[14]));
        ROS_INFO_STREAM(
                "Virtual display height in pixels (minicap): " << *reinterpret_cast<uint32_t *>(&buffer.data()[18]));

        // read images
        cv::Mat img;
        for (;;) {
            // read frame size
            boost::asio::read(minicap_socket, boost::asio::buffer(buffer), boost::asio::transfer_exactly(4), ec);
            if (ec) {
                ROS_ERROR_STREAM(ec.message());
                return;
            }

            // read frame
            uint32_t frame_size = *reinterpret_cast<uint32_t *>(&buffer[0]);
            buffer.resize(frame_size);
            boost::asio::read(minicap_socket, boost::asio::buffer(buffer), boost::asio::transfer_exactly(frame_size),
                              ec);
            if (ec) {
                ROS_ERROR_STREAM(ec.message());
                return;
            }

            // decode and publish
            cv::imdecode(cv::Mat(buffer), cv::IMREAD_COLOR, &img);
            sensor_msgs::ImagePtr msg = cv_bridge::CvImage(std_msgs::Header(), "bgr8", img).toImageMsg();
            m_image_pub.publish(msg);
        }
    }
    catch (std::exception &e) {
        ROS_ERROR_STREAM(e.what());
    }
    catch (...) {
        ROS_ERROR_STREAM("Unknown exception.");
    }
}


text:
// -*- mode:c++;tab-width:2;indent-tabs-mode:t;show-trailing-whitespace:t;rm-trailing-spaces:t -*-
// vi: set ts=2 noet:

#include <riflib/HBondedPairGenerator.hh>
#include <core/pack/rotamers/SingleResidueRotamerLibrary.hh>
#include <core/pack/rotamers/SingleResidueRotamerLibraryFactory.hh>

#include <Eigen/Geometry>

namespace devel {
namespace scheme {

typedef numeric::xyzVector<core::Real> Vec;
typedef numeric::xyzMatrix<core::Real> Mat;
typedef numeric::xyzTransform<core::Real> Xform;

Xform
my_align( Vec to, Vec from )
{
	double const angle = std::acos(to.normalized().dot(from.normalized()));
	if( fabs( angle ) < 0.0001 ){
		return Xform::identity();
	} else if( fabs(angle) > 2*M_PI-0.0001 ){
		// arbitrary perp. vector
		Xform(numeric::rotation_matrix(to.cross(Vec(1,2,3)),-angle));
	}
	return Xform(numeric::rotation_matrix(to.cross(from),-angle));
}
// // broken!
// Xform
// align_fast( Vec to, Vec from )
// {
// 	// rot 180 around avg. axis
// 	return Xform( 2.0*numeric::projection_matrix<double>(to+from)-Mat::identity() );
// }

core::Real
get_rot_score(
	core::pose::Pose & pose,
	core::Size ir,
	core::pack::dunbrack::RotamerLibraryScratchSpaceOP scratch
){
	core::pack::rotamers::SingleResidueRotamerLibraryCOP rotlib =
		core::pack::rotamers::SingleResidueRotamerLibraryFactory::get_instance()->get( pose.residue(ir).type() );
	return rotlib->rotamer_energy( pose.residue(ir), pose, *scratch );
}

void dump_pdb_atom(
	std::ostream & out,
	double x, double y, double z
){
	// std::string atomname,resname,elem;
	// int atomnum,resnum;
	// char chain;
	// bool ishet;
	// float occ,bfac;
	BOOST_VERIFY( x<10000 && x > -1000 );
	BOOST_VERIFY( y<10000 && y > -1000 );
	BOOST_VERIFY( z<10000 && z > -1000 );
	// std::cout << "ATOM   1604  C   GLU A 220       5.010  12.933   1.553  1.00 41.10           C" << std::endl;
	char buf[128];
	// std::string aname = a.atomname;
	// if( aname.size() == 1 ) aname = aname+"  ";
	// if( aname.size() == 2 ) aname = aname+" ";
	snprintf(buf,128,"%s%5i %4s %3s %c%4i    %8.3f%8.3f%8.3f%6.2f%6.2f %11s\n",
		true ? "HETATM":"ATOM  ",
		1,
		"DOTS",
		"DOT",
		'A',
		1,
		x,y,z,
		1.0,
		1.0,
		"D"
	);
	out << buf;
}

///@brief generate BCC lattice "hyper-ring" to sample rotationr around
/// the hbond axis as well as tilting the axis a up to short_tol
utility::vector1<numeric::xyzMatrix<double> >
get_rotation_samples( double short_tol, double long_tol, double reslang ){
	using namespace Eigen;
	typedef Map<Matrix3d> MMap;

	short_tol = numeric::conversions::radians(short_tol);
	long_tol = numeric::conversions::radians(long_tol);
	reslang = numeric::conversions::radians(reslang);

	if( reslang == 0.0 ) return utility::vector1<numeric::xyzMatrix<double> >(1,numeric::xyzMatrix<double>::identity());

	double spacing = sqrt(3.0)*reslang;
	int const n1 = std::ceil(long_tol/spacing);
	spacing = long_tol/n1;
	int const n2 = std::ceil(short_tol/spacing);
	std::cout << "gen_hbond_rotations: spacing " << spacing << " n1 " << n1 << " n2 " << n2 << std::endl;
	std::cout << "   short_tol " << short_tol << std::endl;
	std::cout << "   long_tol  " << long_tol << std::endl;
	std::cout << "   reslang   " << reslang << std::endl;

	// utility::io::ozstream out("test.pdb");
	// generate bcc "square rod" wrapped around in 3sphere wx plane
	utility::vector1<numeric::xyzMatrix<double> > rots;
	for(int i =  0 ; i <  n1; ++i){
	for(int j = -n2; j <= n2; ++j){
	for(int k = -n2; k <= n2; ++k){
	for(int o =  0;  o <=  1; ++o){
		// if(o && ( j==n2 || k==n2 ) ) continue;
		double const wx = i*spacing + (o?spacing/2.0:0.0) - long_tol/2.0;
		double const  w = cos( wx/2.0 );
		double const  x = sin( wx/2.0 );
		double const  y = ( j*spacing + (o?spacing/2.0:0.0) ) / 2.0;
		double const  z = ( k*spacing + (o?spacing/2.0:0.0) ) / 2.0;
		if( y*y + z*z > short_tol*short_tol/4.0 ) continue;
		// double const r = 1.0+y;
		// dump_pdb_atom(out,r*w*50.0,r*x*50.0,r*z*50.0);
		Quaterniond q(w,x,y,z);
		// printf("%4d %4d %4d %8.3f %8.3f %8.3f %8.3f\n",i,j,k,wx,y,z,q.norm());
		q.normalize();
		// dump_pdb_atom(out,wx*50.0,y*50.0,z*50.0);
		numeric::xyzMatrix<double> rot;
		MMap mm((double*)&rot);
		mm = q.matrix();
		// std::cout << rot << std::endl;
		// std::cout << q.matrix() << std::endl;
		// std::cout << std::endl;
		rots.push_back(rot);
	}}}}

	// out.close();

// std::exit(0);

	return rots;
}

void HBondedPairGenerator::init(
	RotamerIndex const & rot_index,
	std::string resn1,
	std::string resn2,
	double tip_tol,
	double rot_resl,
	double rot_range,
	bool fix_donor,
	bool fix_acceptor,
	core::id::AtomID align_atom1,
	core::id::AtomID align_atom2,
	core::id::AtomID align_atom3,
	std::map<std::string,core::pose::Pose> const & exemplars
){
	std::cout << "ROTATION_SAMPLING tip_tol: " << tip_tol << " rot_range: " << rot_range << " resl: " << rot_resl << std::endl;
	rot_samples_ = get_rotation_samples(tip_tol,rot_range,rot_resl);
	fix_donor_ = fix_donor;
	fix_acceptor_ = fix_acceptor;

	scratch_ = core::pack::dunbrack::RotamerLibraryScratchSpaceOP( new core::pack::dunbrack::RotamerLibraryScratchSpace );

	// #ifdef USE_OPENMP
	// #pragma omp critical
	// #endif
	{
		// only hb score
		#pragma omp critical
		{
			// make sure this has been done in a critical block
			core::pack::rotamers::SingleResidueRotamerLibraryFactory::get_instance();

			score_func_ = core::scoring::ScoreFunctionOP( new core::scoring::ScoreFunction );
			// score_func_ = core::scoring::get_score_function();
			if( resn1=="GLY" || resn2=="GLY" ){
				score_func_->set_weight( core::scoring::hbond_bb_sc, 1.1 );
			} else {
				score_func_->set_weight( core::scoring::hbond_sc, 1.1 );
			}
			// score_func_->set_weight( core::scoring::fa_atr, 0.8 );
			// score_func_->set_weight( core::scoring::fa_rep, 0.44 );
			// score_func_->set_weight( core::scoring::fa_sol, 0.75 );
			// set use_hb_env_dep false
				core::scoring::methods::EnergyMethodOptions opts = score_func_->energy_method_options();
				core::scoring::hbonds::HBondOptions hopts = opts.hbond_options();
				// std::cout << "use_hb_env_dep " <<  hopts.use_hb_env_dep() << std::endl;
				hopts.use_hb_env_dep( false );
				opts.hbond_options( hopts );
				// std::cout << score_func_->energy_method_options().hbond_options().use_hb_env_dep() << std::endl;
				score_func_->set_energy_method_options( opts );
				// std::cout << score_func_->energy_method_options().hbond_options().use_hb_env_dep() << std::endl;
				// std::cout << core::scoring::get_score_function()->energy_method_options().hbond_options().use_hb_env_dep() << std::endl;
				// std::cout << "SCORE_FUNCTION: " << std::endl << *score_func_ << std::endl;
				// utility_exit_with_message("roitsn");
		}
		core::chemical::ResidueTypeSetCAP rts;
		#pragma omp critical
		{
			rts = core::chemical::ChemicalManager::get_instance()->residue_type_set("fa_standard");
		}
 		std::cout << "-----" << resn1 << "-----" << resn2 << std::endl;
 		//rtype1op_ = rts.lock()->name_map(resn1).clone();
		//rtype2op_ = rts.lock()->name_map(resn2).clone();

		//TYU change what residue get loaded base on d_l_map_
		// first one is D seoncd one is L
 		if (rot_index.d_l_map_.find(resn1) != rot_index.d_l_map_.end() && rot_index.d_l_map_.find(resn2) == rot_index.d_l_map_.end()) {
 			core::chemical::ResidueTypeCOP rt_tmp = rts.lock() -> name_map(rot_index.d_l_map_.find(resn1)->second).get_self_ptr();
 			core::chemical::ResidueTypeCOP tmp_dcop = rts.lock() -> get_d_equivalent(rt_tmp);
 			rtype1op_ = tmp_dcop -> clone();
 			rtype2op_ = rts.lock() -> name_map(resn2).clone();
 		}
 		// first one is L seoncd one is D
 		else if (rot_index.d_l_map_.find(resn1) == rot_index.d_l_map_.end() && rot_index.d_l_map_.find(resn2) != rot_index.d_l_map_.end()) {
 			rtype1op_ = rts.lock() -> name_map(resn1).clone();
			core::chemical::ResidueTypeCOP rt_tmp = rts.lock() -> name_map(rot_index.d_l_map_.find(resn2)->second).get_self_ptr();
 			core::chemical::ResidueTypeCOP tmp_dcop = rts.lock() -> get_d_equivalent(rt_tmp);
 			rtype2op_ = tmp_dcop -> clone();
 		}
 		// both D
 		else if (rot_index.d_l_map_.find(resn1) != rot_index.d_l_map_.end() && rot_index.d_l_map_.find(resn2) != rot_index.d_l_map_.end()) {
 			core::chemical::ResidueTypeCOP rt1_tmp = rts.lock() -> name_map(rot_index.d_l_map_.find(resn1)->second).get_self_ptr();
 			core::chemical::ResidueTypeCOP tmp1_dcop = rts.lock() -> get_d_equivalent(rt1_tmp);
 			rtype1op_ = tmp1_dcop -> clone();
			core::chemical::ResidueTypeCOP rt2_tmp = rts.lock() -> name_map(rot_index.d_l_map_.find(resn2)->second).get_self_ptr();
 			core::chemical::ResidueTypeCOP tmp2_dcop = rts.lock() -> get_d_equivalent(rt2_tmp);
 			rtype2op_ = tmp2_dcop -> clone();
 		}
		// both L
 		else {
			rtype1op_ = rts.lock() -> name_map(resn1).clone();
 			rtype2op_ = rts.lock() -> name_map(resn2).clone();
 		}
		std::cout << "----done-----" << std::endl;
		if( rtype1op_->has( "H" ) && resn1 != "GLY" ) rtype1op_->set_atom_type( "H", "VIRT" );
		if( rtype1op_->has( "O" ) && resn1 != "GLY" ) rtype1op_->set_atom_type( "O", "VIRT" );
		if( rtype2op_->has( "H" ) && resn2 != "GLY" ) rtype2op_->set_atom_type( "H", "VIRT" );
		if( rtype2op_->has( "O" ) && resn2 != "GLY" ) rtype2op_->set_atom_type( "O", "VIRT" );
 		// core::chemical::ResidueType const & rtype1 = rts.lock()->name_map(resn1);
		// core::chemical::ResidueType const & rtype2 = rts.lock()->name_map(resn2);
	}

	core::chemical::ResidueType const & rtype1 = *rtype1op_;
	core::chemical::ResidueType const & rtype2 = *rtype2op_;
	core::conformation::ResidueOP res1op, res2op;
	#pragma omp critical
	{
		res1op = core::conformation::ResidueFactory::create_residue(rtype1);
		res2op = core::conformation::ResidueFactory::create_residue(rtype2);
	}

	if( exemplars.find(resn1) == exemplars.end() ){
		pose0_.append_residue_by_jump( *res1op, 1 );
	} else {
		// #ifdef USE_OPENMP
		// #pragma omp critical
		// #endif
		std::cout << "found " << resn1 << " in exemplars"  << std::endl;
		pose0_.append_residue_by_jump( exemplars.find(resn1)->second.residue(1), 1 );
	}
	if( exemplars.find(resn2) == exemplars.end() ){
		pose0_.append_residue_by_jump( *res2op, 1 );
	} else {
		// #ifdef USE_OPENMP
		// #pragma omp critical
		// #endif
		std::cout << "found " << resn2 << " in exemplars"  << std::endl;
		pose0_.append_residue_by_jump( exemplars.find(resn2)->second.residue(1), 1 );
	}

	// core::pose::add_lower_terminus_type_to_pose_residue(pose0_,1);
	// core::pose::add_lower_terminus_type_to_pose_residue(pose0_,2);
	// core::pose::add_upper_terminus_type_to_pose_residue(pose0_,1);
	// core::pose::add_upper_terminus_type_to_pose_residue(pose0_,2);
	// pose0_.set_phi(1,64.0);
	// pose0_.set_phi(2,64.0);
	// pose0_.set_psi(1,41.0);
	// pose0_.set_psi(2,41.0);

	pose_ = pose0_;

	{
		std::pair<size_t,size_t> b1 = rot_index.index_bounds(resn1.substr(0,3)), b2 = rot_index.index_bounds(resn2.substr(0,3));
		irotindex_start1_ = b1.first;
		irotindex_start2_ = b2.first;
		if( b1.first < b1.second ){
			std::cout << "adding rotamers res1 " << resn1 << " " << b1.first << "-" << b1.second-1;
			for(size_t i = b1.first; i < b1.second; ++i){
				// irotindex_start1_=98; // for testing asn BTN issue
				// if( i != 98 ) continue;
				std::cout << " " << i;
				utility::vector1<float> chis;
				for(size_t j = 0; j < rot_index.rotamers_[i].chi_.size(); ++j)
					chis.push_back( rot_index.rotamers_[i].chi_[j] );
				rotset1_.push_back( chis );
				runtime_assert( pose_.residue(1).nchi()==rotset1_.back().size() );
				// runtime_assert( i == rotset1_.size()+irotindex_start1_-1 );
			}
			std::cout << std::endl;
		} else {
			std::cout << "no rotamers for res1 " << resn1 << std::endl;
			utility::vector1<float> chis;
			for(size_t j = 0; j < pose0_.residue(1).nchi(); ++j)
				chis.push_back( pose0_.residue(1).chi(j+1) );
			rotset1_.push_back( chis );
		}
		if( b2.first < b2.second ){
			std::cout << "adding rotamers res2 " << resn2 << " " << b2.first << "-" << b2.second-1;
			for(size_t i = b2.first; i < b2.second; ++i){
				std::cout << " " << i;
				utility::vector1<float> chis;
				for(size_t j = 0; j < rot_index.rotamers_[i].chi_.size(); ++j)
					chis.push_back( rot_index.rotamers_[i].chi_[j] );
				rotset2_.push_back( chis );
				runtime_assert( pose_.residue(2).nchi()==rotset2_.back().size() );
				runtime_assert( i == rotset2_.size()+irotindex_start2_-1 );
			}
			std::cout << std::endl;
		} else {
			std::cout << "no rotamers for res2 " << resn2 << std::endl;
			utility::vector1<float> chis;
			for(size_t j = 0; j < pose0_.residue(2).nchi(); ++j)
				chis.push_back( pose0_.residue(2).chi(j+1) );
			rotset2_.push_back( chis );
		}
	}

	nrots1_ = rotset1_.size();
	nrots2_ = rotset2_.size();
	runtime_assert( nrots1_ );
	runtime_assert( nrots2_ );

	align_atom1_ = align_atom1;
	align_atom2_ = align_atom2;
	align_atom3_ = align_atom3;

	{
		std::ostringstream oss;
		oss << "HBondedPairGenerator: align atoms: "
			<< pose_.residue(align_atom1_.rsd()).atom_name(align_atom1_.atomno()) << " "
			<< pose_.residue(align_atom2_.rsd()).atom_name(align_atom2_.atomno()) << " "
			<< pose_.residue(align_atom3_.rsd()).atom_name(align_atom3_.atomno()) << std::endl;
		for(size_t i = 1; i <= pose_.residue(align_atom1_.rsd()).natoms(); ++i){
			oss << i << " " << pose_.residue(align_atom1_.rsd()).atom_name(i) << "  ";
		}
		std::cout << oss.str() << std::endl;
	}

	clash_atom_start1_ = clash_atom_start2_ = 1;
	if( fix_donor    ){
		nrots1_ = 1;
		// std::cout << "WARNING: clash checking bb atoms for donor with -fix_donor" << std::endl;
		clash_atom_start1_ = std::min( align_atom1_.atomno(), (core::Size)6 ); // TODO: check last controlling chi or something instead
		runtime_assert( align_atom1_.rsd() == 1 );
		runtime_assert( align_atom2_.rsd() == 1 );
		runtime_assert( align_atom3_.rsd() == 1 );
	}
	if( fix_acceptor ){
		nrots2_ = 1;
		// std::cout << "WARNING: clash checking bb atoms for acceptor with -fix_acceptor" << std::endl;
		clash_atom_start2_ = std::min( align_atom1_.atomno(), (core::Size)6 ); // TODO: check last controlling chi or something instead
		runtime_assert( align_atom1_.rsd() == 2 );
		runtime_assert( align_atom2_.rsd() == 2 );
		runtime_assert( align_atom3_.rsd() == 2 );
	}

	if( fix_donor && fix_acceptor ) std::cout << "WARNING: doesn't make sense to fix both donor and acceptor" << std::endl;
	donor_atoms_ = rtype1.Hpos_polar_sc();
	if( resn1=="GLY" ) donor_atoms_ = rtype1.Hpos_polar();
	for( auto i : donor_atoms_) donor_bases_.push_back( rtype1.atom_base(i) );
	acceptor_atoms_ = rtype2.accpt_pos_sc();
	if( resn2=="GLY" || resn2=="PRO" ||
	    resn2=="ADX" || resn2=="CYX" || resn2=="GUX" || resn2=="THX" // bbonly dna stuff
	){
		acceptor_atoms_ = rtype2.accpt_pos();
	}
	acceptor_orbitals_.resize(acceptor_atoms_.size());
	std::vector< Vec > orbseenit;
	std::vector<int> acceptors_to_remove;
	for(core::Size iacc = 1; iacc <= acceptor_atoms_.size(); ++iacc){
		// #ifdef USE_OPENMP
		// #pragma omp critical
		// #endif
		std::cout << "consider orbitals for atom " << iacc << " " << rtype2.atom_name(acceptor_atoms_[iacc]) << std::endl;
		for( auto j :rtype2.bonded_orbitals(acceptor_atoms_[iacc])){
			if(
				// rtype2.orbital_type(j).orbital_enum() == core::chemical::orbitals::C_pi_sp2    ||
			    // rtype2.orbital_type(j).orbital_enum() == core::chemical::orbitals::N_pi_sp2    ||
			    rtype2.orbital_type(j).orbital_enum() == core::chemical::orbitals::N_p_sp2     ||
			    // rtype2.orbital_type(j).orbital_enum() == core::chemical::orbitals::O_pi_sp2    ||
			    rtype2.orbital_type(j).orbital_enum() == core::chemical::orbitals::O_p_sp2     ||
			    rtype2.orbital_type(j).orbital_enum() == core::chemical::orbitals::S_p_sp3     ||
			    // rtype2.orbital_type(j).orbital_enum() == core::chemical::orbitals::O_pi_sp2_bb ||
			    // rtype2.orbital_type(j).orbital_enum() == core::chemical::orbitals::O_p_sp2_bb  ||
			    rtype2.orbital_type(j).orbital_enum() == core::chemical::orbitals::O_p_sp3
			){
				bool isnew = true;
				for( auto seenorb : orbseenit ){
					if( pose_.residue(2).orbital_xyz(j).distance_squared( seenorb ) < 0.001 ){
						isnew = false;
					}
				}
				if( ! isnew ){
					// #ifdef USE_OPENMP
					// #pragma omp critical
					// #endif
					std::cout << "HBPairGen WARNING: duplicate orbital found! "
					          << rtype2.name() << " "
					          << rtype2.atom_name(acceptor_atoms_[iacc]) << " "
					          << j << " "
					          << rtype2.orbital_type(j).name()
					          << std::endl;
			        // utility_exit_with_message("dup orb.");
				} else {
					if( acceptor_orbitals_[iacc].size() >= 2 ){
						// #ifdef USE_OPENMP
						// #pragma omp critical
						// #endif
						std::cout << "HBPairGen WARNING: skipping orbitals beyond first two on " << rtype2.name() << " "
						          << rtype2.atom_name(acceptor_atoms_[iacc]) << " orbname: "
	    					      << rtype2.orbital_type(j).name()
						          <<  std::endl;
			          utility_exit_with_message("extra orbs");
					} else {
						orbseenit.push_back( pose_.residue(2).orbital_xyz(j) );
						acceptor_orbitals_[iacc].push_back(j);
					}
				}
			}
		}
		if( acceptor_orbitals_[iacc].size()==0 ){
			std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
			std::cout << "WARNING!!! no acceptor orbitals "+resn2+" atom: "+ pose_.residue(2).atom_name(acceptor_atoms_[iacc]) << std::endl;
			std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
			acceptors_to_remove.push_back(iacc);
		}
	}

	if( acceptors_to_remove.size() ) { // remove "acceptors" without any applicable orbitals
		utility::vector1<core::Size> acceptor_atoms_tmp;
		utility::vector1<utility::vector1<core::Size> > acceptor_orbitals_tmp;
		for( size_t i = 1; i <= acceptor_atoms_.size(); ++i ){
			if( std::find(acceptors_to_remove.begin(),acceptors_to_remove.end(),i) != acceptors_to_remove.end() ) continue;
			acceptor_atoms_tmp.push_back( acceptor_atoms_.at(i) );
			acceptor_orbitals_tmp.push_back( acceptor_orbitals_.at(i) );
		}
		acceptor_atoms_ = acceptor_atoms_tmp;
		acceptor_orbitals_ = acceptor_orbitals_tmp;
	}

	if( donor_atoms_.size()==0 ) utility_exit_with_message("no donor atoms");
	if( acceptor_atoms_.size()==0 ) utility_exit_with_message("no acceptor atoms");
	// #ifdef DEBUG
		// for(core::Size i = 1; i <= nrots1_; ++i){
		// 	if( pose0_.residue(1).nchi()==1 ) printf("donor chi %7.2f \n",rotset1_->rotamer(i)->chi(1));
		// 	if( pose0_.residue(1).nchi()==2 ) printf("donor chi %7.2f %7.2f \n",rotset1_->rotamer(i)->chi(1),rotset1_->rotamer(i)->chi(2));
		// 	if( pose0_.residue(1).nchi()==3 ) printf("donor chi %7.2f %7.2f %7.2f \n",
		// 		rotset1_->rotamer(i)->chi(1),rotset1_->rotamer(i)->chi(2),rotset1_->rotamer(i)->chi(3));
		// 	if( pose0_.residue(1).nchi()==4 ) printf("donor chi %7.2f %7.2f %7.2f %7.2f \n",
		// 		rotset1_->rotamer(i)->chi(1),rotset1_->rotamer(i)->chi(2),rotset1_->rotamer(i)->chi(3),rotset1_->rotamer(i)->chi(4));
		// }
		// std::cout << "acceptor nchi: " << pose0_.residue(2).nchi() << std::endl;
		// for(core::Size i = 1; i <= nrots2_; ++i){
		// 	if( pose0_.residue(2).nchi()==1 ) printf("acceptor chi %7.2f \n",rotset2_->rotamer(i)->chi(1));
		// 	if( pose0_.residue(2).nchi()==2 ) printf("acceptor chi %7.2f %7.2f \n",rotset2_->rotamer(i)->chi(1),rotset2_->rotamer(i)->chi(2));
		// 	if( pose0_.residue(2).nchi()==3 ) printf("acceptor chi %7.2f %7.2f %7.2f \n",
		// 		rotset2_->rotamer(i)->chi(1),rotset2_->rotamer(i)->chi(2),rotset2_->rotamer(i)->chi(3));
		// 	if( pose0_.residue(2).nchi()==4 ) printf("acceptor chi %7.2f %7.2f %7.2f %7.2f \n",
		// 		rotset2_->rotamer(i)->chi(1),rotset2_->rotamer(i)->chi(2),rotset2_->rotamer(i)->chi(3),rotset2_->rotamer(i)->chi(4));
		// }
		// #ifdef USE_OPENMP
		// #pragma omp critical
		// #endif
		{
			for(core::Size iacc=1; iacc<=acceptor_atoms_.size(); ++iacc)
				std::cout << "acceptor orbitals, index: " << iacc << " atom: " << rtype2.atom_name(acceptor_atoms_[iacc]) << " orbitals: " << acceptor_orbitals_[iacc] << std::endl;
			std::cout << "donor atoms: ";
			for( auto idonatom : donor_atoms_ ) std::cout << rtype1.atom_name(idonatom) << ", ";
			std::cout << std::endl;
			std::cout << "donor nchi: " << pose0_.residue(1).nchi() << std::endl;
			std::cout << "donor num rots: " << nrots1_ << std::endl;
			std::cout << "acceptor num rots: " << nrots2_ << std::endl;
			std::cout << "N hbond geom samples: " << rot_samples_.size() << std::endl;
			std::cout << "total raw samples: " << raw_num_samples() << std::endl;
		}
	// #endif
	irot1_=1; irot2_=1; idon_=1; iacc_=1; iorb_=1; ihbr_=1;
	update_chi();
	if( !update_hb_and_score() ) increment();
}

	bool
	SingleHbondedPairGenerator::update_hb_and_score(
		bool realign
	){
    	typedef numeric::xyzVector<core::Real> Vec;
    	typedef numeric::xyzTransform<core::Real> Xform;
		// if(realign) std::cout << "REALIGN_HB " << irot1_ << " " << idon_ << " " << irot2_ << " " << iacc_ << " " << iorb_ << " " << ihbr_ << std::endl;
		// else        std::cout << "UPDATE_HB  " << irot1_ << " " << idon_ << " " << irot2_ << " " << iacc_ << " " << iorb_ << " " << ihbr_ << std::endl;
		// check if already "aligned"
		// pose0_.dump_pdb("before_update_hb_and_score_test0.pdb");
		// pose_ .dump_pdb("before_update_hb_and_score_test.pdb");

		if( realign ){
			// this does not fix the strange rotamer scoring issue, where successive rotamers don't get correct scores
			// #ifdef USE_OPENMP
			// #pragma omp critical
			// #endif
			// scratch_ = core::pack::dunbrack::RotamerLibraryScratchSpaceOP( new core::pack::dunbrack::RotamerLibraryScratchSpace );
			Vec don  = Vec( pose0_.residue(1).xyz(donor_atoms_[idon_]) );
			Vec donb = Vec( pose0_.residue(1).xyz(donor_bases_[idon_]) );
			Vec acc  = Vec( pose0_.residue(2).xyz(acceptor_atoms_[iacc_]));
			Vec acco = Vec( pose0_.residue(2).orbital_xyz(acceptor_orbitals_.at(iacc_).at(iorb_)));
			// Vec d = (don-donb).normalized();
			// Vec a = (acc-acco).normalized();
			Xform xd = my_align( Vec(1,0,0), don-donb );
			Xform xa = my_align( Vec(1,0,0), acc-acco );

			// std::cout << xd*(don-donb).normalized() << std::endl;
			// std::cout << xa*(acc-acco).normalized() << std::endl;
			runtime_assert( fabs( (xd*((don-donb).normalized())).dot(Vec(1,0,0)) - 1.0 ) < 0.001 );
			runtime_assert( fabs( (xa*((acc-acco).normalized())).dot(Vec(1,0,0)) - 1.0 ) < 0.001 );

			if( fabs((don-donb).y()) < 0.00001 && fabs((don-donb).z()) < 0.00001 ) xd = Xform::identity();
			if( fabs((acc-acco).y()) < 0.00001 && fabs((acc-acco).z()) < 0.00001 ) xa = Xform::identity();
			// std::cout << "ddon " << don-donb << std::endl;
			// std::cout << "dacc " << acc-acco << std::endl;
			// std::cout << "don:  " << don  << std::endl;
			// std::cout << "donb: " << donb << std::endl;
			// std::cout << "acc:  " << acc  << std::endl;
			// std::cout << "acco: " << acco << std::endl;
			// std::cout << "xd:   " << xd   << std::endl;
			// std::cout << "xa:   " << xa   << std::endl;
			core::scoring::motif::xform_pose( pose0_, xd, 1, 1 );
			core::scoring::motif::xform_pose( pose0_, xa, 2, 2 );
			core::scoring::motif::xform_pose( pose0_, Vec(-1.15,0,0)-pose0_.residue(1).xyz(donor_atoms_   [idon_]), 1, 1 );
			core::scoring::motif::xform_pose( pose0_, Vec( 0.65,0,0)-pose0_.residue(2).xyz(acceptor_atoms_[iacc_]), 2, 2 );
			pose_ = pose0_;

		} else {
			for(core::Size ia = 1; ia <= pose0_.residue(1).natoms(); ++ia)
				pose_.set_xyz( core::id::AtomID(ia,1) , pose0_.residue(1).xyz(ia) );
			for(core::Size ia = 1; ia <= pose0_.residue(2).natoms(); ++ia)
				pose_.set_xyz( core::id::AtomID(ia,2) , rot_samples_[ihbr_] * pose0_.residue(2).xyz(ia) );
		}

		// only carbon/carbon
		for(core::Size ia = clash_atom_start1_; ia <= pose_.residue(1).nheavyatoms(); ++ia){
			if( ia > 1 && pose_.residue(1).atom_type_index(ia) > 6 ) continue;
			// if(ia==4) continue; // skip BB O, don't know psi angle
			core::Real const lj1 = pose_.residue(1).atom_type(ia).lj_radius();
			for(core::Size ja = clash_atom_start2_; ja <= pose_.residue(2).nheavyatoms(); ++ja){
				if( ja > 1 && pose_.residue(2).atom_type_index(ja) > 6 ) continue;
				// if( ia== donor_bases_[idon_] && ja==acceptor_atoms_[iacc_] ) continue;
				core::Real const lj2 = pose_.residue(2).atom_type(ja).lj_radius();
				core::Real const clash_dis = (lj1+lj2-0.5);
				core::Real const clash_dis2 = clash_dis*clash_dis;
				if( clash_dis2 > pose_.residue(1).xyz(ia).distance_squared( pose_.residue(2).xyz(ja) ) ){
					// std::cout << sqrt(clash_dis2) << " " << pose_.residue(1).xyz(ia).distance( pose_.residue(2).xyz(ja)) << std::endl;
					return false;
				}
			}
		}

		using core::id::AtomID;

		// utility_exit_with_message("oarintseiransteionraeti");
		core::scoring::motif::xform_pose( pose_, ~Xform( pose_.xyz(align_atom1_), pose_.xyz(align_atom2_), pose_.xyz(align_atom3_) ) );

		// if( realign ){
		// 	static int count = 0;
		// 	++count;
		// 	pose0_ .dump_pdb("after_realign_rotamer_"+boost::lexical_cast<std::string>(count)+".pdb");
		// 	pose_  .dump_pdb("after_realign_to_frame_"+boost::lexical_cast<std::string>(count)+".pdb");
		// }

		try {
			score_ = score_func_->score(pose_);
			// if( nrots1_ > 1 ) score_ += 0.56/2.0 * get_rot_score(pose_,1,scratch_);  // this is bugged, successive rotamers don't get proper scores
			// if( nrots2_ > 1 ) score_ += 0.56/2.0 * get_rot_score(pose_,2,scratch_);
			// std::cout << score_func_->score(pose_) << " " << get_rot_score(pose_,1,scratch_) << " " << get_rot_score(pose_,2,scratch_) << " " << score_ << std::endl;
		} catch(...) {
			{
				std::cout << "WARNING: scoring failed!" << std::endl;
				if(realign) std::cout << "REALIGN_HB " << irot1_ << " " << idon_ << " " << irot2_ << " " << iacc_ << " " << iorb_ << " " << ihbr_ << std::endl;
				else        std::cout << "UPDATE_HB  " << irot1_ << " " << idon_ << " " << irot2_ << " " << iacc_ << " " << iorb_ << " " << ihbr_ << std::endl;
				pose_.dump_pdb("score_fail.pdb");
				utility_exit_with_message("scoring failed!!");
			}
		}
		// std::cout << "DUN 1 " << get_rot_score(pose_,1,scratch_) << std::endl;
		// std::cout << "DUN 2 " << get_rot_score(pose_,2,scratch_) << std::endl;

		return true;
    }

	bool
	BidentateHbondedPairGenerator::update_hb_and_score(
		bool realign
	){
    	typedef numeric::xyzVector<core::Real> Vec;
    	typedef numeric::xyzMatrix<core::Real> Mat;
    	typedef numeric::xyzTransform<core::Real> Xform;

   		if( iacc_!=1 || iorb_!=2 ) return false;
   		if( idon_!=4 && idon_!=5 ) return false;
   		core::Size iacc_other = 2;
   		core::Size idon_other = idon_==4 ? 1 : 3;
		// #ifdef USE_OPENMP
		// #pragma omp critical
		// #endif
   		std::cout << idon_ << " " << idon_other << " " << iacc_ << " " << iacc_other << std::endl;

		// if(realign) std::cout << "REALIGN_HB " << irot1_ << " " << idon_ << " " << irot2_ << " " << iacc_ << " " << iorb_ << " " << ihbr_ << std::endl;
		// else        std::cout << "UPDATE_HB  " << irot1_ << " " << idon_ << " " << irot2_ << " " << iacc_ << " " << iorb_ << " " << ihbr_ << std::endl;

		// check if already "aligned"
		if( realign ){
			scratch_ = core::pack::dunbrack::RotamerLibraryScratchSpaceOP(new core::pack::dunbrack::RotamerLibraryScratchSpace);
			Vec don1  = Vec( pose0_.residue(1).xyz(donor_atoms_[idon_]) );
			Vec don2  = Vec( pose0_.residue(1).xyz(donor_atoms_[idon_other]) );
			Vec donb1 = Vec( pose0_.residue(1).xyz(donor_bases_[idon_]) );
			Vec donb2 = Vec( pose0_.residue(1).xyz(donor_bases_[idon_other]) );
			Vec acc1  = Vec( pose0_.residue(2).xyz(acceptor_atoms_[iacc_]));
			Vec acc2  = Vec( pose0_.residue(2).xyz(acceptor_atoms_[iacc_other]));
			Vec acco1 = Vec( pose0_.residue(2).orbital_xyz(acceptor_orbitals_.at(iacc_     ).at(iorb_)) );
			Vec acco2 = Vec( pose0_.residue(2).orbital_xyz(acceptor_orbitals_.at(iacc_other).at(iorb_)) );

			// #ifdef USE_OPENMP
			// #pragma omp critical
			// #endif
			{
				std::cout << "DON " << don1-donb1 << std::endl;
				std::cout << "DON " << don2-donb2 << std::endl;
			}

			// std::exit(0);
			Vec const dax1 = (don1+don2-donb1-donb2).normalized();
			Vec const aax1 = (acc1+acc2-acco1-acco2).normalized();
			Mat Rd = numeric::alignVectorSets( dax1, (don1-don2).cross(dax1).normalized(), Vec(0,0,1), Vec(0,1,0) );
			Mat Ra = numeric::alignVectorSets( aax1, (acc1-acc2).cross(aax1).normalized(), Vec(0,0,1), Vec(0,1,0) );
			Xform xd = Xform( Rd, Rd*(-( don1+ don2)/2.0) + Vec(0,0,-1.050) );
			Xform xa = Xform( Ra, Ra*(-(acco1+acco2)/2.0) + Vec(0,0, 0.065) );
			core::scoring::motif::xform_pose( pose0_, xd, 1, 1 );
			core::scoring::motif::xform_pose( pose0_, xa, 2, 2 );
			// core::scoring::motif::xform_pose( pose0_, Vec(-1.05,0,0)-pose0_.residue(1).xyz(donor_atoms_   [idon_]), 1, 1 );
			// core::scoring::motif::xform_pose( pose0_, Vec( 0.65,0,0)-pose0_.residue(2).xyz(acceptor_atoms_[iacc_]), 2, 2 );
			pose_ = pose0_;
		}

		for(core::Size ia = 1; ia <= pose0_.residue(2).natoms(); ++ia){
			pose_.set_xyz( core::id::AtomID(ia,2) , rot_samples_[ihbr_] * pose0_.residue(2).xyz(ia) );
		}

		// for(core::Size ia = clash_atom_start1_; ia <= pose0_.residue(1).nheavyatoms(); ++ia){
		// 	if(ia==4) continue;
		// 	core::Real const lj1 = pose0_.residue(1).atom_type(ia).lj_radius();
		// 	for(core::Size ja = clash_atom_start2_; ja <= pose0_.residue(2).nheavyatoms(); ++ja){
		// 		if( ja==4 || ( ia== donor_bases_[idon_] && ja==acceptor_atoms_[iacc_] ) ) continue;
		// 		core::Real const lj2 = pose0_.residue(2).atom_type(ja).lj_radius();
		// 		core::Real const clash_dis = (lj1+lj2-0.5);
		// 		core::Real const clash_dis2 = clash_dis*clash_dis;
		// 		if( clash_dis2 > pose0_.residue(1).xyz(ia).distance_squared( pose0_.residue(2).xyz(ja) ) ){
		// 			// std::cout << sqrt(clash_dis2) << " " << pose0_.residue(1).xyz(ia).distance( pose0_.residue(2).xyz(ja)) << std::endl;
		// 			return false;
		// 		}
		// 	}
		// }

		score_ = score_func_->score(pose_);
		// if( nrots1_ > 1 ) score_ += get_rot_score(pose_,1,scratch_); // this is bugged, successive rotamers don't get proper scores
		// if( nrots2_ > 1 ) score_ += get_rot_score(pose_,2,scratch_);	//

		return true;
    }

}
}

ext: cc

max_stars_repo_name: YaoYinYing/rifdock

max_issues_repo_name: YaoYinYing/rifdock

max_forks_repo_name: YaoYinYing/rifdock


text: #include <boost/spirit/include/classic_common.hpp>

max_stars_repo_name: miathedev/BoostForArduino

max_issues_repo_name: miathedev/BoostForArduino

max_forks_repo_name: miathedev/BoostForArduino


text:
#pragma once

#include "defines.hpp"

#include <string>
#include <sstream>

#include <boost/utility/string_view.hpp>

namespace msrv {

using StringView = boost::string_view;

class Tokenizer
{
public:
    Tokenizer(StringView str, char sep)
        : input_(str), sep_(sep) { }

    static bool hasToken(StringView value, StringView token, char sep);

    const StringView& token() const { return token_; }
    const StringView& input() const { return input_; }

    bool nextToken();

private:
    StringView token_;
    StringView input_;
    char sep_;
};

template<typename T>
std::string toString(const T& value)
{
    std::ostringstream stream;
    stream << value;
    return stream.str();
}

std::string formatString(const char* fmt, ...) MSRV_FORMAT_FUNC(1, 2);

StringView trim(StringView str, char ch);
StringView trimWhitespace(StringView str);

}

max_stars_repo_name: courtarro/beefweb

max_issues_repo_name: courtarro/beefweb

max_forks_repo_name: courtarro/beefweb


text:
/* Copyright 2016-2017 Joaquin M Lopez Munoz.
 * 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/libs/poly_collection for library home page.
 */

#ifndef BOOST_POLY_COLLECTION_DETAIL_IS_FINAL_HPP
#define BOOST_POLY_COLLECTION_DETAIL_IS_FINAL_HPP

#if defined(_MSC_VER)
#pragma once
#endif

#include <boost/type_traits/is_final.hpp>
#include <type_traits>

/* technique explained at
 * http://bannalia.blogspot.com/2016/09/compile-time-checking-existence-of.html
 */

namespace boost{
namespace poly_collection{
namespace detail{
namespace is_final_fallback{

template<typename T> using is_final=boost::is_final<T>;

struct hook{};

}}}}

namespace std{

template<>
struct is_void< ::boost::poly_collection::detail::is_final_fallback::hook>:
  std::false_type
{
  template<typename T>
  static constexpr bool is_final_f()
  {
    using namespace ::boost::poly_collection::detail::is_final_fallback;
    return is_final<T>::value;
  }
};

} /* namespace std */

namespace boost{

namespace poly_collection{

namespace detail{

template<typename T>
struct is_final:std::integral_constant<
  bool,
  std::is_void<is_final_fallback::hook>::template is_final_f<T>()
>{};

} /* namespace poly_collection::detail */

} /* namespace poly_collection */

} /* namespace boost */

#endif


text:
#ifndef MESH_HPP
#define MESH_HPP

#include <string>

#include <armadillo>

/**
 * @brief Class for reading, writing .gri files and computing mesh matrices
 */
struct Mesh
{
    /**
     * @brief Read the mesh from file
     * @param[in] inFile Name of input .gri file name
     */
    void readFromFile(const std::string &inFile);

    /**
     * @brief Write mesh to file
     * @param[in] outFile Name of output file
     */
    void writeToFile(const std::string &outFile) const;

    /**
     * @brief Compute the mesh Matrices
     */
    void computeMatrices();

    /**
     * @brief Write Matrices to file
     * @param[in] matFile Name of output file
     */
    void writeMatricesToFile(const std::string &matFile) const;

    /**
     * @brief Number of nodes
     */
    arma::uword nNode;

    /**
     * @brief Number of elements
     */
    arma::uword nElemTot;

    /**
     * @brief Dimension of the space in which the mesh is embedded
     */
    arma::uword dim;

    /**
     * @brief Coordinates of the node points
     */
    arma::mat nodeCoordinates;

    /**
     * @brief Number of boundary groups
     */
    arma::uword nBGroup;

    /**
     * @brief Number of faces in each boundary group
     */
    arma::uvec nBFace;

    /**
     * @brief Number of boundary nodes per face in each boundary group
     */
    arma::uvec nf;

    /**
     * @brief Title of each boundary groups
     */
    arma::field<std::string> title;

    /**
     * @brief Mapping of boundary faces to node numbers
     */
    arma::field<arma::umat> B2N;

    /**
     * @brief Order of the elements
     */
    arma::uword order;

    /**
     * @brief Type of the elemeents
     */
    std::string basis;

    /**
     * @brief Mapping of elements to node numbers
     */
    arma::umat E2N;

    /**
     * @brief Mapping of interior faces to face midpoints
     */
    arma::mat  I2M;

    /**
     * @brief Mapping of interior faces to elements
     *
     * Each row contains four entries [elemL, faceL, elemR, faceR], where
     * elemL and elemR are the indices of the neighboring elements (we ensure
     * elemL < elemR), and faceL and faceR are the local node numbers of the
     * node opposite to the face.
     */
    arma::umat I2E;

    /**
     * @brief Mapping of interior faces to unit normals
     *
     * Unit normals point from elemL to elemR
     */
    arma::mat  In;

    /**
     * @brief Mapping of interior faces to face lengths
     */
    arma::vec  Il;

    /**
     * @brief Mapping of boundary faces to face midpoints
     */
    arma::mat  B2M;

    /**
     * @brief Mapping of boundary faces to elements
     *
     * Each row contains four entries [elem, face, bgroup], where elem is the
     * index of the adjoining element, and face is the local node numbers of
     * the node opposite to the face. Finally, bgroup is the number of the
     * group the boundary face belongs to.
     */
    arma::umat B2E;

    /**
     * @brief Mapping of boundary faces to unit normals
     *
     * Unit normals point out to the exterior
     */
    arma::mat  Bn;

    /**
     * @brief Mapping of boundary faces to face lengths
     */
    arma::vec  Bl;

    /**
     * @brief Mapping of elements to element midpoints
     */
    arma::mat E2M;

    /**
     * @brief Mapping of elements to element areas
     */
    arma::vec E2A;
};

#endif // MESH_HPP


text:
/* boost random/lognormal_distribution.hpp header file
 *
 * Copyright Jens Maurer 2000-2001
 * Copyright Steven Watanabe 2011
 * 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$
 *
 * Revision history
 *  2001-02-18  moved to individual header files
 */

#ifndef BOOST_RANDOM_LOGNORMAL_DISTRIBUTION_HPP
#define BOOST_RANDOM_LOGNORMAL_DISTRIBUTION_HPP

#include <boost/config/no_tr1/cmath.hpp>      // std::exp, std::sqrt
#include <cassert>
#include <iosfwd>
#include <istream>
#include <boost/limits.hpp>
#include <boost/random/detail/config.hpp>
#include <boost/random/detail/operators.hpp>
#include <boost/random/normal_distribution.hpp>

namespace boost {
namespace random {

/**
 * Instantiations of class template lognormal_distribution model a
 * \random_distribution. Such a distribution produces random numbers
 * with \f$\displaystyle p(x) = \frac{1}{x s \sqrt{2\pi}} e^{\frac{-\left(\log(x)-m\right)^2}{2s^2}}\f$
 * for x > 0.
 *
 * @xmlwarning
 * This distribution has been updated to match the C++ standard.
 * Its behavior has changed from the original
 * boost::lognormal_distribution.  A backwards compatible
 * version is provided in namespace boost.
 * @endxmlwarning
 */
template<class RealType = double>
class lognormal_distribution
{
public:
    typedef typename normal_distribution<RealType>::input_type input_type;
    typedef RealType result_type;

    class param_type
    {
    public:

        typedef lognormal_distribution distribution_type;

        /** Constructs the parameters of a lognormal_distribution. */
        explicit param_type(RealType m_arg = RealType(0.0),
                            RealType s_arg = RealType(1.0))
          : _m(m_arg), _s(s_arg) {}

        /** Returns the "m" parameter of the distribution. */
        RealType m() const { return _m; }

        /** Returns the "s" parameter of the distribution. */
        RealType s() const { return _s; }

        /** Writes the parameters to a std::ostream. */
        BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, param_type, parm)
        {
            os << parm._m << " " << parm._s;
            return os;
        }

        /** Reads the parameters from a std::istream. */
        BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, param_type, parm)
        {
            is >> parm._m >> std::ws >> parm._s;
            return is;
        }

        /** Returns true if the two sets of parameters are equal. */
        BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(param_type, lhs, rhs)
        { return lhs._m == rhs._m && lhs._s == rhs._s; }

        /** Returns true if the two sets of parameters are different. */
        BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(param_type)

    private:
        RealType _m;
        RealType _s;
    };

    /**
     * Constructs a lognormal_distribution. @c m and @c s are the
     * parameters of the distribution.
     */
    explicit lognormal_distribution(RealType m_arg = RealType(0.0),
                                    RealType s_arg = RealType(1.0))
      : _normal(m_arg, s_arg) {}

    /**
     * Constructs a lognormal_distribution from its parameters.
     */
    explicit lognormal_distribution(const param_type& parm)
      : _normal(parm.m(), parm.s()) {}

    // compiler-generated copy ctor and assignment operator are fine

    /** Returns the m parameter of the distribution. */
    RealType m() const { return _normal.mean(); }
    /** Returns the s parameter of the distribution. */
    RealType s() const { return _normal.sigma(); }

    /** Returns the smallest value that the distribution can produce. */
    RealType min BOOST_PREVENT_MACRO_SUBSTITUTION () const
    { return RealType(0); }
    /** Returns the largest value that the distribution can produce. */
    RealType max BOOST_PREVENT_MACRO_SUBSTITUTION () const
    { return (std::numeric_limits<RealType>::infinity)(); }

    /** Returns the parameters of the distribution. */
    param_type param() const { return param_type(m(), s()); }
    /** Sets the parameters of the distribution. */
    void param(const param_type& parm)
    {
        typedef normal_distribution<RealType> normal_type;
        typename normal_type::param_type normal_param(parm.m(), parm.s());
        _normal.param(normal_param);
    }

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

    /**
     * Returns a random variate distributed according to the
     * lognormal distribution.
     */
    template<class Engine>
    result_type operator()(Engine& eng)
    {
        using std::exp;
        return exp(_normal(eng));
    }

    /**
     * Returns a random variate distributed according to the
     * lognormal distribution with parameters specified by param.
     */
    template<class Engine>
    result_type operator()(Engine& eng, const param_type& parm)
    { return lognormal_distribution(parm)(eng); }

    /** Writes the distribution to a @c std::ostream. */
    BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, lognormal_distribution, ld)
    {
        os << ld._normal;
        return os;
    }

    /** Reads the distribution from a @c std::istream. */
    BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, lognormal_distribution, ld)
    {
        is >> ld._normal;
        return is;
    }

    /**
     * Returns true if the two distributions will produce identical
     * sequences of values given equal generators.
     */
    BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(lognormal_distribution, lhs, rhs)
    { return lhs._normal == rhs._normal; }

    /**
     * Returns true if the two distributions may produce different
     * sequences of values given equal generators.
     */
    BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(lognormal_distribution)

private:
    normal_distribution<result_type> _normal;
};

} // namespace random

/// \cond show_deprecated

/**
 * Provided for backwards compatibility.  This class is
 * deprecated.  It provides the old behavior of lognormal_distribution with
 * \f$\displaystyle p(x) = \frac{1}{x \sigma_N \sqrt{2\pi}} e^{\frac{-\left(\log(x)-\mu_N\right)^2}{2\sigma_N^2}}\f$
 * for x > 0, where \f$\displaystyle \mu_N = \log\left(\frac{\mu^2}{\sqrt{\sigma^2 + \mu^2}}\right)\f$ and
 * \f$\displaystyle \sigma_N = \sqrt{\log\left(1 + \frac{\sigma^2}{\mu^2}\right)}\f$.
 */
template<class RealType = double>
class lognormal_distribution
{
public:
    typedef typename normal_distribution<RealType>::input_type input_type;
    typedef RealType result_type;

    lognormal_distribution(RealType mean_arg = RealType(1.0),
                           RealType sigma_arg = RealType(1.0))
      : _mean(mean_arg), _sigma(sigma_arg)
    {
        init();
    }
    RealType mean() const { return _mean; }
    RealType sigma() const { return _sigma; }
    void reset() { _normal.reset(); }
    template<class Engine>
    RealType operator()(Engine& eng)
    {
        using std::exp;
        return exp(_normal(eng) * _nsigma + _nmean);
    }
    BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, lognormal_distribution, ld)
    {
        os << ld._normal << " " << ld._mean << " " << ld._sigma;
        return os;
    }
    BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, lognormal_distribution, ld)
    {
        is >> ld._normal >> std::ws >> ld._mean >> std::ws >> ld._sigma;
        ld.init();
        return is;
    }
private:
    /// \cond show_private
    void init()
    {
        using std::log;
        using std::sqrt;
        _nmean = log(_mean*_mean/sqrt(_sigma*_sigma + _mean*_mean));
        _nsigma = sqrt(log(_sigma*_sigma/_mean/_mean+result_type(1)));
    }
    RealType _mean;
    RealType _sigma;
    RealType _nmean;
    RealType _nsigma;
    normal_distribution<RealType> _normal;
    /// \endcond
};

/// \endcond

} // namespace boost

#endif // BOOST_RANDOM_LOGNORMAL_DISTRIBUTION_HPP


text:
#ifdef __APPLE__

#include "generator.other.hpp"

#include <ctime>
#include <functional>
#include <limits>
#include <locale>
#include <vector>

#include <boost/assert.hpp>

#include "stream.hpp"

namespace blackhole {
inline namespace v1 {
namespace datetime {

namespace {

typedef ostream_adapter<writer_type> stream_type;

} // namespace

struct context_t {
    // TODO: Get rid of intermediate stream for facets. Its construction consumes ~70ns, which is
    // comparable with other time required to format timestamps.
    stream_type wr;
    std::tm tm;
    std::uint64_t usec;
    const std::vector<std::string>& literals;
    std::vector<std::string>::const_iterator it;

    context_t(writer_type& wr, std::tm tm, std::uint64_t usec, const std::vector<std::string>& literals):
        wr(wr),
        tm(tm),
        usec(usec),
        literals(literals),
        it(std::begin(this->literals))
    {}
};

namespace {

template<std::size_t length, char filler = '0'>
inline auto fill(int value, char* buffer) -> void {
    std::size_t digits = 0;
    do {
        buffer[length - 1 - digits] = (value % 10) + '0';
        value /= 10;
        ++digits;
    } while (value);

    std::memset(buffer, filler, length - digits);
}

template<std::size_t length, char filler = '0'>
inline void fill(stream_type& wr, int value) {
    char buffer[std::numeric_limits<uint32_t>::digits10 + 2];
    fill<length, filler>(value, buffer);
    wr << fmt::StringRef(buffer, length);
}

enum class days_t {
    sunday,
    monday
};

namespace visit {

template<char Format>
inline void localized(context_t& context) {
    typedef std::time_put<char> facet_type;
    typedef typename facet_type::iter_type iter_type;

    std::use_facet<facet_type>(context.wr.getloc())
        .put(iter_type(context.wr), context.wr, ' ', &context.tm, Format);
}

namespace year {

inline void full(context_t& context) {
    fill<4>(context.wr, context.tm.tm_year + 1900);
}

inline void last_two_digits(context_t& context) {
    fill<2>(context.wr, context.tm.tm_year % 100);
}

inline void first_two_digits(context_t& context) {
    fill<2>(context.wr, context.tm.tm_year / 100 + 19);
}

} // namespace year

namespace month {

inline void numeric(context_t& context) {
    fill<2>(context.wr, context.tm.tm_mon + 1);
}

} // namespace month

namespace week {

template<days_t firstweekday>
inline void year(context_t& context);

template<>
inline void year<days_t::sunday>(context_t& context) {
    fill<2>(context.wr, (context.tm.tm_yday - context.tm.tm_wday + 7) / 7);
}

template<>
inline void year<days_t::monday>(context_t& context) {
    fill<2>(context.wr, (context.tm.tm_yday - context.tm.tm_wday + (context.tm.tm_wday ? 8 : 1)) / 7);
}

} // namespace week

namespace day {

namespace year {

inline void numeric(context_t& context) {
    fill<3>(context.wr, context.tm.tm_yday + 1);
}

} // namespace year

namespace month {

template<char filler = '0'>
inline void numeric(context_t& context) {
    fill<2, filler>(context.wr, context.tm.tm_mday);
}

} // namespace month

} // namespace day

namespace time {

namespace hour {

inline void h24(context_t& context) {
    fill<2>(context.wr, context.tm.tm_hour);
}

inline void h12(context_t& context) {
    const int mod = context.tm.tm_hour % 12;
    fill<2>(context.wr, mod == 0 ? 12 : mod);
}

} // namespace hour

namespace minute {

inline void normal(context_t& context) {
    fill<2>(context.wr, context.tm.tm_min);
}

} // namespace minute

namespace second {

inline void normal(context_t& context) {
    fill<2>(context.wr, context.tm.tm_sec);
}

inline void epoch(context_t& context) {
    context.wr << std::mktime(&context.tm);
}

} // namespace second

namespace usecond {

inline void normal(context_t& context) {
    fill<6>(context.wr, static_cast<int>(context.usec));
}

} // namespace usecond

} // namespace time

inline void utc_offset(context_t& context) {
    // It's too dificult to reinvent this case depending on platform.
    // See GNUC strftime if you don't believe.
    localized<'z'>(context);
}

inline void standard(context_t& context) {
    typedef std::time_put<char> facet_type;
    typedef facet_type::iter_type iter_type;

    std::use_facet<facet_type>(context.wr.getloc())
        .put(iter_type(context.wr), context.wr, ' ', &context.tm, 'a');
    context.wr << ' ';
    std::use_facet<facet_type>(context.wr.getloc())
        .put(iter_type(context.wr), context.wr, ' ', &context.tm, 'b');
    context.wr << ' ';
    day::month::numeric(context);
    context.wr << ' ';
    time::hour::h24(context);
    context.wr << ':';
    time::minute::normal(context);
    context.wr << ':';
    time::second::normal(context);
    context.wr << ' ';
    year::full(context);
}

inline void date_mdy(context_t& context) {
    month::numeric(context);
    context.wr << '/';
    day::month::numeric(context);
    context.wr << '/';
    year::last_two_digits(context);
}

inline void date_ISO8601(context_t& context) {
    year::full(context);
    context.wr << '-';
    month::numeric(context);
    context.wr << '-';
    day::month::numeric(context);
}

inline void time_ISO8601(context_t& context) {
    time::hour::h24(context);
    context.wr << ':';
    time::minute::normal(context);
    context.wr << ':';
    time::second::normal(context);
}

} // namespace visit

inline void literal_generator_t(context_t& context) {
    context.wr << *context.it;
    ++context.it;
}

class generator_handler_t {
    typedef std::string::const_iterator iterator_type;

    std::vector<token_t>& tokens;
    std::vector<std::string>& literals;

    std::string litbuffer;

public:
    generator_handler_t(std::vector<token_t>& tokens, std::vector<std::string>& literals) :
        tokens(tokens),
        literals(literals)
    {}

    virtual ~generator_handler_t() {}

    virtual void partial_literal(iterator_type begin, iterator_type end) {
        litbuffer.append(begin, end);
    }

    virtual void end_partial_literal() {
        if (!litbuffer.empty()) {
            literal(litbuffer);
            litbuffer.clear();
        }
    }

    virtual void literal(const std::string& literal) {
        tokens.push_back(&literal_generator_t);
        literals.push_back(literal);
    }

    virtual void placeholder(const std::string& value) {
        literal(value);
    }

    virtual void full_year() {
        end_partial_literal();
        tokens.push_back(&visit::year::full);
    }

    virtual void short_year() {
        end_partial_literal();
        tokens.push_back(&visit::year::last_two_digits);
    }

    virtual void short_year_first_two_digits() {
        end_partial_literal();
        tokens.push_back(&visit::year::first_two_digits);
    }

    virtual void numeric_month() {
        end_partial_literal();
        tokens.push_back(&visit::month::numeric);
    }

    virtual void abbreviate_month() {
        end_partial_literal();
        tokens.push_back(&visit::localized<'b'>);
    }

    virtual void full_month() {
        end_partial_literal();
        tokens.push_back(&visit::localized<'B'>);
    }

    virtual void year_week(days_t firstweekday) {
        end_partial_literal();
        if (firstweekday == days_t::sunday) {
            tokens.push_back(&visit::week::year<days_t::sunday>);
        } else {
            tokens.push_back(&visit::week::year<days_t::monday>);
        }
    }

    virtual void year_day() {
        end_partial_literal();
        tokens.push_back(&visit::day::year::numeric);
    }

    virtual void month_day(bool has_leading_zero = true) {
        end_partial_literal();
        if (has_leading_zero) {
            tokens.push_back(&visit::day::month::numeric<'0'>);
        } else {
            tokens.push_back(&visit::day::month::numeric<' '>);
        }
    }

    virtual void abbreviate_weekday() {
        end_partial_literal();
        tokens.push_back(&visit::localized<'a'>);
    }

    virtual void full_weekday() {
        end_partial_literal();
        tokens.push_back(&visit::localized<'A'>);
    }

    virtual void hours() {
        end_partial_literal();
        tokens.push_back(&visit::time::hour::h24);
    }

    virtual void hours12() {
        end_partial_literal();
        tokens.push_back(&visit::time::hour::h12);
    }

    virtual void minute() {
        end_partial_literal();
        tokens.push_back(&visit::time::minute::normal);
    }

    virtual void second() {
        end_partial_literal();
        tokens.push_back(&visit::time::second::normal);
    }

    virtual void epoch() {
        end_partial_literal();
        tokens.push_back(&visit::time::second::epoch);
    }

    virtual void usecond() {
        end_partial_literal();
        tokens.push_back(&visit::time::usecond::normal);
    }

    virtual void ampm() {
        end_partial_literal();
        tokens.push_back(&visit::localized<'p'>);
    }

    virtual void utc_offset() {
        end_partial_literal();
        tokens.push_back(&visit::utc_offset);
    }

    virtual void timezone() {
        end_partial_literal();
        tokens.push_back(&visit::localized<'Z'>);
    }

    virtual void standard_date_time() {
        end_partial_literal();
        tokens.push_back(&visit::standard);
    }

    virtual void date_mdy() {
        end_partial_literal();
        tokens.push_back(&visit::date_mdy);
    }

    virtual void date_ISO8601() {
        end_partial_literal();
        tokens.push_back(&visit::date_ISO8601);
    }

    virtual void time_ISO8601() {
        end_partial_literal();
        tokens.push_back(&visit::time_ISO8601);
    }
};

namespace parser {

class through_t {
public:
    typedef std::string::const_iterator iterator_type;

    static iterator_type parse(iterator_type it, iterator_type end, generator_handler_t& handler) {
        BOOST_ASSERT(it != end && it + 1 != end);
        (void)end;

        handler.placeholder(std::string(it, it + 2));
        return it + 2;
    }
};

template<class Decorate = through_t>
class common {
public:
    typedef typename Decorate::iterator_type iterator_type;

    static iterator_type parse(iterator_type it, iterator_type end, generator_handler_t& handler) {
        BOOST_ASSERT(it != end && it + 1 != end);
        switch (*(it + 1)) {
        case 'c':
            break;
        default:
            return Decorate::parse(it, end, handler);
        }
        return it + 2;
    }
};

template<class Decorate>
class time {
public:
    typedef typename Decorate::iterator_type iterator_type;

    static iterator_type parse(iterator_type it, iterator_type end, generator_handler_t& handler) {
        BOOST_ASSERT(it != end && it + 1 != end);
        switch (*(it + 1)) {
        case 'H':
            handler.hours();
            break;
        case 'I':
            handler.hours12();
            break;
        case 'M':
            handler.minute();
            break;
        case 'S':
            handler.second();
            break;
        case 's':
            handler.epoch();
            break;
        case 'f':
            handler.usecond();
            break;
        case 'p':
            handler.ampm();
            break;
        default:
            return Decorate::parse(it, end, handler);
        }

        return it + 2;
    }
};

template<class Decorate>
class date {
public:
    typedef typename Decorate::iterator_type iterator_type;

    static iterator_type parse(iterator_type it, iterator_type end, generator_handler_t& handler) {
        BOOST_ASSERT(it != end && it + 1 != end);
        switch (*(it + 1)) {
        //! =========== YEAR SECTION ===========
        case 'Y':
            handler.full_year();
            break;
        case 'y':
            handler.short_year();
            break;
        case 'C':
            handler.short_year_first_two_digits();
            break;
        //! =========== MONTH SECTION ===========
        case 'm':
            handler.numeric_month();
            break;
        case 'b':
        case 'h':
            handler.abbreviate_month();
            break;
        case 'B':
            handler.full_month();
            break;
        //! =========== WEAK SECTION ===========
        case 'U':
            handler.year_week(days_t::sunday);
            break;
        case 'W':
            handler.year_week(days_t::monday);
            break;
        //! =========== DAY SECTION ===========
        case 'j':
            handler.year_day();
            break;
        case 'd':
            handler.month_day();
            break;
        case 'e':
            handler.month_day(false);
            break;
        case 'a':
            handler.abbreviate_weekday();
            break;
        case 'A':
            handler.full_weekday();
            break;
        //! =========== OTHER SECTION ===========
        case 'c':
            handler.standard_date_time();
            break;
        case 'D':
            handler.date_mdy();
            break;
        case 'F':
            handler.date_ISO8601();
            break;
        case 'T':
            handler.time_ISO8601();
            break;
        case 'z':
            handler.utc_offset();
            break;
        case 'Z':
            handler.timezone();
            break;
        case '%':
            handler.literal(std::string("%"));
            break;
        default:
            return Decorate::parse(it, end, handler);
        }

        return it + 2;
    }
};

} // namespace parser

template<class Decorate>
class parser_t {
    typedef typename Decorate::iterator_type iterator_type;

public:
    static void parse(const std::string& format, generator_handler_t& handler) {
        iterator_type it = format.begin();
        iterator_type end = format.end();

        while (it != end) {
            iterator_type p = std::find(it, end, '%');
            handler.partial_literal(it, p);

            if (std::distance(p, end) >= 2) {
                it = Decorate::parse(p, end, handler);
            } else {
                if (p != end) {
                    handler.partial_literal(p, end);
                }
                break;
            }
        }

        handler.end_partial_literal();
    }
};

}  // namespace

generator_t::generator_t(std::vector<token_t> tokens, std::vector<std::string> literals) :
    tokens(std::move(tokens)),
    literals(std::move(literals))
{}

auto generator_t::operator()(writer_type& stream, const std::tm& tm, std::uint64_t usec) const -> void {
    context_t context{stream, tm, usec, literals};

    for (auto& token : tokens) {
        token(context);
    }
}

auto make_generator(const std::string& pattern) -> generator_t {
    std::vector<token_t> tokens;
    std::vector<std::string> literals;

    generator_handler_t handler(tokens, literals);

    parser_t<parser::date<parser::time<parser::common<parser::through_t>>>>::parse(pattern, handler);

    return generator_t(std::move(tokens), std::move(literals));
}

} // namespace datetime
} // namespace v1
} // namespace blackhole

#endif

max_stars_repo_name: JakariaBlaine/blackhole

max_issues_repo_name: JakariaBlaine/blackhole

max_forks_repo_name: JakariaBlaine/blackhole


text:
// This file is part of libigl, a simple c++ geometry processing library.
//
// Copyright (C) 2013 Alec Jacobson <alecjacobson@gmail.com>
//
// This Source Code Form is subject to the terms of the Mozilla Public License
// v. 2.0. If a copy of the MPL was not distributed with this file, You can
// obtain one at http://mozilla.org/MPL/2.0/.
#include "unproject.h"

#include <Eigen/Dense>
#include <Eigen/LU>

template <
  typename Derivedwin,
  typename Derivedmodel,
  typename Derivedproj,
  typename Derivedviewport,
  typename Derivedscene>
IGL_INLINE void igl::unproject(
  const Eigen::PlainObjectBase<Derivedwin>&  win,
  const Eigen::PlainObjectBase<Derivedmodel>& model,
  const Eigen::PlainObjectBase<Derivedproj>& proj,
  const Eigen::PlainObjectBase<Derivedviewport>&  viewport,
  Eigen::PlainObjectBase<Derivedscene> & scene)
{
  typedef typename Derivedscene::Scalar Scalar;
  Eigen::Matrix<Scalar,4,4> Inverse =
    (proj.template cast<Scalar>() * model.template cast<Scalar>()).inverse();

  Eigen::Matrix<Scalar,4,1> tmp;
  tmp << win, 1;
  tmp(0) = (tmp(0) - viewport(0)) / viewport(2);
  tmp(1) = (tmp(1) - viewport(1)) / viewport(3);
  tmp = tmp.array() * 2.0f - 1.0f;

  Eigen::Matrix<Scalar,4,1> obj = Inverse * tmp;
  obj /= obj(3);

  scene = obj.head(3);
}

template <typename Scalar>
IGL_INLINE Eigen::Matrix<Scalar,3,1> igl::unproject(
  const    Eigen::Matrix<Scalar,3,1>&  win,
  const    Eigen::Matrix<Scalar,4,4>& model,
  const    Eigen::Matrix<Scalar,4,4>& proj,
  const    Eigen::Matrix<Scalar,4,1>&  viewport)
{
  Eigen::Matrix<Scalar,3,1> scene;
  unproject(win,model,proj,viewport,scene);
  return scene;
}

#ifdef IGL_STATIC_LIBRARY
template Eigen::Matrix<float, 3, 1, 0, 3, 1> igl::unproject<float>(Eigen::Matrix<float, 3, 1, 0, 3, 1> const&, Eigen::Matrix<float, 4, 4, 0, 4, 4> const&, Eigen::Matrix<float, 4, 4, 0, 4, 4> const&, Eigen::Matrix<float, 4, 1, 0, 4, 1> const&);
template Eigen::Matrix<double, 3, 1, 0, 3, 1> igl::unproject<double>(Eigen::Matrix<double, 3, 1, 0, 3, 1> const&, Eigen::Matrix<double, 4, 4, 0, 4, 4> const&, Eigen::Matrix<double, 4, 4, 0, 4, 4> const&, Eigen::Matrix<double, 4, 1, 0, 4, 1> const&);
#endif

max_stars_repo_name: FabianRepository/SinusProject

max_issues_repo_name: FabianRepository/SinusProject

max_forks_repo_name: FabianRepository/SinusProject


text: #include <boost/mpl/integral_c_tag.hpp>