text:
/*
 * Copyright (c) 2009 David Verhasselt <david@crowdway.com>
 * Licensed under The MIT License - http://gmailcc.crowdway.com/license.txt
 */

#include "MailRecord.h"

#include <boost/filesystem.hpp>
namespace bf = boost::filesystem;

#include "LogSink.h"
#include "MailLink.h"
#include "MailBox.h"

int MailRecord::number_of_deliveries = 0;

MailRecord::MailRecord()
{
	marked = false;

	flag_draft = false;
	flag_flagged = false;
	flag_passed = false;
	flag_replied = false;
	flag_seen = false;
	flag_trashed = false;
}

MailRecord::~MailRecord()
{
	for (vector<MailLink*>::iterator iter = this->links.begin(); iter < links.end(); iter++)
	{
		delete(*iter);
	}
}

void MailRecord::mark()
{
	marked = true;
}

bool MailRecord::get_flag_draft()
{
	return flag_draft;
}

void MailRecord::set_flag_draft(bool flag_draft)
{
	if (this->flag_draft == flag_draft)
		return;

	this->flag_draft = flag_draft;
}

bool MailRecord::get_flag_flagged()
{
	return flag_flagged;
}

void MailRecord::set_flag_flagged(bool flag_flagged)
{
	if (this->flag_flagged == flag_flagged)
		return;

	this->flag_flagged = flag_flagged;
}

bool MailRecord::get_flag_passed()
{
	return flag_passed;
}

void MailRecord::set_flag_passed(bool flag_passed)
{
	if (this->flag_passed == flag_passed)
		return;

	this->flag_passed = flag_passed;
}

bool MailRecord::get_flag_replied()
{
	return flag_replied;
}

void MailRecord::set_flag_replied(bool flag_replied)
{
	if (this->flag_replied == flag_replied)
		return;

	this->flag_replied = flag_replied;
}

bool MailRecord::get_flag_seen()
{
	return flag_seen;
}

void MailRecord::set_flag_seen(bool flag_seen)
{
	if (this->flag_seen == flag_seen)
		return;

	this->flag_seen = flag_seen;
}

bool MailRecord::get_flag_trashed()
{
	return flag_trashed;
}

void MailRecord::set_flag_trashed(bool flag_trashed)
{
	if (this->flag_trashed == flag_trashed)
		return;

	this->flag_trashed = flag_trashed;
}

void MailRecord::save_content(string body)
{
	ofstream mfile;

	mfile.open(mainlink->get_path().c_str());

	if (mfile.is_open())
	{
		mfile << body;
		mfile.close();
	}
	else Log::error << "Unable to open file " << mainlink->get_path() << Log::endl;
}

MailLink* MailRecord::add_to_mailbox(MailBox* mailbox, unsigned long uid, string path)
{
	MailLink* maillink = find_link(mailbox);

	if (maillink != NULL)			// Existing mail: just update UID
	{
		maillink->uid = uid;

		if (!path.empty())
		{
			if (!rename(maillink->get_path().c_str(), path.c_str()))
				perror("MailRecord::add_to_mailbox: unable to move link");
			else
				maillink->set_path(path);
		}
	}
	else
	{
		maillink = new MailLink;

		if (mailbox->primary)
			mainlink = maillink;

		bool exists = false;
		if (path.empty())
			path = mailbox->get_path() + "cur/" + generate_md_filename() + get_md_info();
		else
			exists = true;

		maillink->mailrecord = this;
		maillink->mailbox = mailbox;
		maillink->set_path(path);
		maillink->uid = uid;

		this->links.push_back(maillink);
		mailbox->mails.push_back(maillink);

		if (!exists)
		{
			if (this->links.size() > 1)	// Create a link to the original
			{
				string sourcepath = this->mainlink->get_path().c_str();
				//sourcepath.insert(0, "../");

				bf::create_hard_link(sourcepath, path);
				//	perror("MailRecord::add_to_mailbox: creating symlink failed");
			}
			else			// "Touch" file
			{
				ofstream mfile;
				mfile.open(path.c_str());
				mfile << "";
				mfile.close();
			}
		}
		else
		{
			load_md_info();
		}
	}

	return maillink;
}

void MailRecord::remove_from_mailbox(MailBox* mailbox)
{
	vector<MailLink*>::iterator iter = find_ilink(mailbox);

	if (iter < this->links.end())
	{
		remove_from_mailbox(iter);
		return;
	}
	else
	{
		cout << "Whoops, something went wrong here" << endl;
	}
}

vector<MailLink*>::iterator MailRecord::remove_from_mailbox(vector<MailLink*>::iterator iter)
{
	if (unlink((*iter)->get_path().c_str()) != 0)
			perror("MailRecord::remove_from_mailbox: unable to unlink file");

	(*iter)->mailbox->dirty();
	(*iter)->stale();

	return this->links.erase(iter);
}

void MailRecord::remove()
{
	for (vector<MailLink*>::iterator iter = this->links.begin(); iter < this->links.end(); )
	{
		iter = this->remove_from_mailbox(iter);
	}
}

MailLink* MailRecord::find_link(MailBox* mailbox)
{
	vector<MailLink*>::iterator result = find_ilink(mailbox);

	if (result == this->links.end())
		return NULL;
	else
		return *result;
}

vector<MailLink*>::iterator MailRecord::find_ilink(MailBox* mailbox)
{
	for (vector<MailLink*>::iterator iter = this->links.begin(); iter < this->links.end(); iter++)
	{
		if ((*iter)->mailbox == mailbox)
		{
			return iter;
		}
	}

	return this->links.end();
}

/**
 * Regenerates the paths of every mail, moves the files
 * of the mails, also of every link. This is useful if one of the
 * files has a different info part because another client changed
 * the flags.
 */
void MailRecord::sync_flags()
{
	for (vector<MailLink*>::iterator iter = this->links.begin(); iter < links.end(); iter++)
	{
		string new_path = (*iter)->get_base_path();

		new_path.insert(new_path.size(), this->get_md_info());

		if (new_path.compare((*iter)->get_path()) != 0)
		{
			if (rename((*iter)->get_path().c_str(), new_path.c_str()) != 0)
				perror("MailRecord::flags_changed: Unable to rename file");
			else
				(*iter)->set_path(new_path);
		}
	}
}

string MailRecord::generate_md_filename()
{
	char seconds[20];	sprintf(seconds, "%lu", (unsigned long) time(NULL));

	char hostname[99] = "hostnamex";
	string shostname;
	pid_t pid = getpid();

	if (pid < 0)
		perror("MailRecord::generate_md_filename: unable to get pid");

	char spid[20];	sprintf(spid, "%d", pid);

	if (gethostname(hostname, 99) != 0)
	{
		perror("MailRecord::generate_md_filename: unable to gethostname()");
		perror(hostname);
	}

 	shostname = hostname;

	// Replace "/" with "\057"
	size_t index;
	while((index = shostname.find("/")) != string::npos)
	{
			shostname.replace(index, 1, "\057");
	}

	// Replace ":" with "\072"
	while((index = shostname.find(":")) != string::npos)
	{
			shostname.replace(index, 1, "\072");
	}

	char delnum[20];	sprintf(delnum, "%d", MailRecord::number_of_deliveries++);

	string result;

	result += seconds;
	result += ".";
	result += "P";
	result += spid;
	result += "Q";
	result += delnum;
	result += ".";
	result += shostname;

	return result;
}

/**
 * Generate the "info" part of the maildir specification
 * see: http://cr.yp.to/proto/maildir.html
 */
string MailRecord::get_md_info()
{
	string result = MD_TAIL;

	if (this->flag_draft) result += "D";
	if (this->flag_flagged) result += "F";
	if (this->flag_passed) result += "P";
	if (this->flag_replied) result += "R";
	if (this->flag_seen) result += "S";
	if (this->flag_trashed) result += "T";

	return result;
}

/**
 * Extract the "info" part from the path and update lfags
 */
void MailRecord::load_md_info()
{
	this->flag_draft = false;
	this->flag_flagged = false;
	this->flag_passed = false;
	this->flag_replied = false;
	this->flag_seen = false;
	this->flag_trashed = false;

	string path = mainlink->get_path();
	for(size_t index = path.rfind(MD_TAIL) + strlen(MD_TAIL); index < path.size(); index++)
	{
				switch(path[index])
				{
					case 'D':
						this->flag_draft = true; break;
					case 'F':
						this->flag_flagged = true; break;
					case 'P':
						this->flag_passed = true; break;
					case 'R':
						this->flag_replied = true; break;
					case 'S':
						this->flag_seen = true; break;
					case 'T':
						this->flag_trashed = true; break;
				}
	}
}

max_stars_repo_name: dv/gmailcc

max_issues_repo_name: dv/gmailcc

max_forks_repo_name: dv/gmailcc


text:
#include <fstream>
#include <iostream>
#include <string>
#include <unordered_map>
#include <vector>

#include <boost/algorithm/string.hpp>

using Wire = std::string;
using Signal = std::uint16_t;
using Command = Signal (*) (Signal, Signal);

struct Operation {
	Command cmd;
	std::string lhs_operand;
	std::string rhs_operand;
};

template<typename T>
T AND(const T lhs, const T rhs) { return lhs & rhs; }

template<typename T>
T OR(const T lhs, const T rhs) { return lhs | rhs; }

template<typename T>
T LSHIFT(const T lhs, const T rhs) { return lhs << rhs; }

template<typename T>
T RSHIFT(const T lhs, const T rhs) { return lhs >> rhs; }

template<typename T>
T NOT(const T lhs, const T) { return ~lhs; }

template<typename T>
T SIGNAL(const T lhs, const T) { return lhs; }

int main() {

	auto filename = std::string{"circuit.txt"};
	auto file = std::fstream{filename};

	if(file.is_open()) {

		auto signal = std::string{};

		auto commands = std::unordered_map<std::string, Command>{
			{"AND", AND},
			{"OR", OR},
			{"LSHIFT", LSHIFT},
			{"RSHIFT", RSHIFT},
			{"NOT", NOT},
			{signal, SIGNAL}
		};

		auto operations = std::unordered_map<Wire, Operation>{};
		auto wires = std::unordered_map<Wire, Signal>{};

		/*
		 * map each line to a set of operations on a single wire
		 */
		const auto max_tokens = 5U;
		auto tokens = std::vector<std::string>{max_tokens};

		auto line = std::string{};

		while(std::getline(file, line)) {

			boost::split(tokens, line, [] (auto c) { return (c == ' '); });

			// last element is always a wire
			auto& wire = tokens.back();

			// 4 tokens means a binary operation, 3 tokens is an unary operation,
			// and 2 tokens is a simple assignment
			if(tokens.size() == 4) {
				operations.emplace(wire, Operation{commands[tokens[1]], tokens[0], tokens[2]});
			} else if(tokens.size() == 3) {
				operations.emplace(wire, Operation{commands[tokens[0]], tokens[1], {}});
			} else {
				operations.emplace(wire, Operation{commands[signal], tokens[0], {}});
			}
		}

		auto get_signal = [] (auto is_operand_signal, auto& operand, auto is_wire_mapped, auto& wire) {
			return (is_operand_signal ? std::stoul(operand) : (is_wire_mapped ? wire->second : 0));
		};

		//execute operations
		for(auto operation = operations.begin(); !operations.empty();) {

			auto& lhs_str = operation->second.lhs_operand;
			auto& rhs_str = operation->second.rhs_operand;

			auto command = operation->second.cmd;

			auto is_lhs_signal = std::isdigit(lhs_str.front());
			auto is_rhs_signal = std::isdigit(rhs_str.front());

			auto lhs_wire = wires.find(lhs_str);
			auto rhs_wire = wires.find(rhs_str);

			auto is_lhs_wire_mapped = (lhs_wire != wires.end());
			auto is_rhs_wire_mapped = (rhs_wire != wires.end());

			auto is_lhs_operand_ok = (is_lhs_signal || is_lhs_wire_mapped);
			auto is_rhs_operand_ok = (is_rhs_signal || rhs_str.empty() || is_rhs_wire_mapped);

			// operation has either signals, mapped wires or a mix of both
			if(is_lhs_operand_ok && is_rhs_operand_ok) {

				auto lhs_signal = get_signal(is_lhs_signal, lhs_str, is_lhs_wire_mapped, lhs_wire);
				auto rhs_signal = get_signal(is_rhs_signal, rhs_str, is_rhs_wire_mapped, rhs_wire);

				wires[operation->first] = command(lhs_signal, rhs_signal);

				operation = operations.erase(operation);

			} else {
				// no operation possible right now, so move on
				++operation;
			}

			// if there are still operations left at the end
			// of the map, restart the whole process
			if(operation == operations.end() && !operations.empty()) {
				operation = operations.begin();
			}
		}

		std::cout << wires["a"] << std::endl;

	} else {
		std::cerr << "Error! Could not open \"" << filename << "\"!" << std::endl;
	}

	return 0;
}

ext: cpp

max_stars_repo_path: Day 07 Part 1/main.cpp

max_issues_repo_path: Day 07 Part 1/main.cpp

max_forks_repo_path: Day 07 Part 1/main.cpp


text:
#include "include/DomainNameTranslator.h"
#include <boost/algorithm/string.hpp>
#include <boost/asio.hpp>
#include <iostream>
#include <glog/logging.h>

using domainnametranslator::DomainNameTranslator;

std::string
getVariableFromBrackets(std::string variable, std::string forwardBracket, std::string reverseBracket) {
    unsigned firstLim = variable.find(forwardBracket);
    unsigned lastLim = variable.find(reverseBracket);
    std::string strNew = variable.substr(firstLim,lastLim);
    boost::erase_all(strNew, forwardBracket);
    boost::erase_all(strNew, reverseBracket);
    strNew.erase(std::remove_if(strNew.begin(), strNew.end(), ::isspace ), strNew.end());
    return strNew;
}

enum class OpType{
    GreaterOrEqual,
    Greater,
    LessOrEqual,
    Less,
    Equal
};

/*
Notes: example command:
"players.elements.collect(player, player.weapon == weapon.beats)"
*/
Variable
evaluate(std::vector<std::string> commandChain,
                const std::vector<std::string>& funcArgs, RunningGameState& state);

Variable
evaluateBinaryOperationCommand(const std::string& compareStatement, RunningGameState& state);

std::vector<std::string>
DomainNameTranslator::parseCommandChain(std::string instruction) {
    std::string leftCurlyBracket = "{";
    std::string rightCurlyBracket = "}";
    std::string leftParenthesis = "(";
    std::string rightParenthesis = ")";
    // Handles parsing instructions enclosed by curly braces e.g. {player.name}
    if (boost::algorithm::contains(instruction, leftCurlyBracket)) {
        std::string variable = getVariableFromBrackets(instruction, leftCurlyBracket, rightCurlyBracket);
        return parseCommandChain(variable);
    }

    std::vector<std::string> result;
    std::string delimiter = ".";
    std::string instructionWithoutParenthesis = instruction.substr(0, instruction.find(leftParenthesis));
    boost::split(result, instructionWithoutParenthesis, boost::is_any_of(delimiter));

    return result;
}

std::vector<std::string>
DomainNameTranslator::parseFuncArgs(std::string instruction) {
    std::string leftParenthesis = "(";
    std::string rightParenthesis = ")";
    std::string variable = getVariableFromBrackets(instruction, leftParenthesis, rightParenthesis);
    std::vector<std::string> result;
    std::string delimiter = ",";
    boost::split(result, variable, boost::is_any_of(delimiter));
    return result;
}

//TODO: add support for "winners.size == players.size" by
Variable
DomainNameTranslator::parseInstruction(const std::string& instruction, RunningGameState& state) {

    if(1){
        //case with operators (==, <=, ...) AND no brackets ie: "winners.size == players.size"
        //can directly call evaluateBinaryOperationCommand
        Variable result = evaluateBinaryOperationCommand(instruction, state);
        return result;
    }

    std::vector<std::string> commandChain = parseCommandChain(instruction);
    std::vector<std::string> funcArgs = parseFuncArgs(instruction);
    Variable result = evaluate(commandChain, funcArgs, state);

    return result;
}

//comparator for each "player.weapon == weapon.beats"
//restrict to string, bool, or int
bool
compare(const std::string& lhs, const std::string& rhs){
    //can only be OpType::Equal
    return lhs == rhs;
}

bool
compare(bool lhs, bool rhs){
    //can only be OpType::Equal
    return lhs == rhs;
}

bool
compare(int lhs, int rhs, OpType optype){
    switch(optype) {
        case OpType::Equal:
        {
            return lhs == rhs;
        }
        case OpType::Less:
        {
            return lhs < rhs;
        }
        case OpType::LessOrEqual:
        {
            return lhs <= rhs;
        }
        case OpType::Greater:
        {
            return lhs > rhs;
        }
        case OpType::GreaterOrEqual:
        {
            return lhs >= rhs;
        }
        default:
            LOG(ERROR) << "unsupported int OpType";
            return false;
    }
}

std::unordered_map<OpType, std::string>
opMap{
    {OpType::Equal, "=="},
	{OpType::LessOrEqual, "<="},
	{OpType::GreaterOrEqual, ">="},
	{OpType::Less, "<"},
	{OpType::Greater, ">"},
};

OpType
getOpType(const std::string& funcArgs){
    if(funcArgs.find("==") != std::string::npos){
        return OpType::Equal;
    }
    if(funcArgs.find("<=") != std::string::npos){
        return OpType::LessOrEqual;
    }
    if(funcArgs.find(">=") != std::string::npos){
        return OpType::GreaterOrEqual;
    }
    if(funcArgs.find(">") != std::string::npos){
        return OpType::Less;
    }
    if(funcArgs.find("<") != std::string::npos){
        return OpType::Greater;
    }
    LOG(ERROR) << "Error: operation not support" << std::endl;
    return OpType::Equal;
}

std::vector<std::string>
getFuncOperands(const std::string& funcArgs, OpType op){
    std::vector<std::string> result;
    std::string delimiter = opMap.at(op);
    boost::split(result, funcArgs, boost::is_any_of(delimiter));
    //remove white spaces
    for(auto& s : result){
        s.erase( std::remove_if(s.begin(), s.end(), ::isspace ), s.end());
    }
    return result;
}

Variable evaluateBinaryOperationCommand(const std::string& compareStatement, RunningGameState& state){
    OpType operation = getOpType(compareStatement);
    std::vector<std::string> operands = getFuncOperands(compareStatement, operation);
    Variable lhs = DomainNameTranslator::parseInstruction(operands[0], state);
    Variable rhs = DomainNameTranslator::parseInstruction(operands[1], state);
    Variable result;
    result.varType = VariableType::BoolType;

    switch(lhs.varType) {
        case VariableType::StringType:
        {
            result.boolVar = compare(lhs.stringVar, rhs.stringVar);
        } break;
        default:
            LOG(ERROR) << "evaluateBinaryOperationCommand unsupported type, please implement" << std::endl;
            return Variable();
    }

    return result;
}

std::shared_ptr<Variable>
collect(std::shared_ptr<Variable> list, const std::vector<std::string>& funcArgs, RunningGameState& state){

    std::shared_ptr<Variable> result = std::make_shared<Variable>();
    result->varType = VariableType::ListType;

    for(auto& varPtr : list->listVar){
        //register list to name so subsequent parseInstruction calls will be able to find the element
        state.variables->updateVariable(funcArgs[0], varPtr);

        //expect true or false
        Variable validator = DomainNameTranslator::parseInstruction(funcArgs[1], state);

        if(validator.boolVar == true){
            result->listVar.push_back(varPtr);
        }
    }

    state.variables->removeVariable(funcArgs[0]);
    return result;
}

std::shared_ptr<Variable>
identify(const std::string& identifier, RunningGameState& state){
    std::shared_ptr<Variable> varPtr = state.variables->getVariable(identifier);
    return varPtr;
}

std::shared_ptr<Variable>
modify(std::shared_ptr<Variable> varPtr, const std::string& modifier, RunningGameState& state){
    if(modifier == "elements"){
        if(varPtr->varType == VariableType::MapType){
            std::shared_ptr<Variable> flattened = std::make_shared<Variable>();
            flattened->varType = VariableType::ListType;
            for(auto& [key, val] : varPtr->mapVar){
                flattened->listVar.push_back(val);
            }
            return flattened;
        } else {return varPtr;}
    }
    LOG(ERROR) << "Error: identifer not supported" << std::endl;
    return nullptr;
}

std::shared_ptr<Variable>
performFunction(std::shared_ptr<Variable> varPtr,
                const std::string& funcName,
                const std::vector<std::string>& funcArgs,
                RunningGameState& state){

    if(funcName == "collect"){
        return collect(varPtr, funcArgs, state);
    }
    LOG(ERROR) << "Error: identifer not supported" << std::endl;
    return nullptr;
}

//evaluates output from parseInstructions
//eg "players.elements.collect(player, player.weapon == weapon.beats)":
//commandChain = {"players", "elements", "collect"}, funcArgs = {player, player.weapon == weapon.beats}
Variable
evaluate(std::vector<std::string> commandChain,
                const std::vector<std::string>& funcArgs, RunningGameState& state){

    std::string funcName;
    //dealing with function
    if(!funcArgs.empty()){
        funcName = commandChain.back();
        commandChain.pop_back();
    }
    //evaluate chain
    std::shared_ptr<Variable> resultPtr = identify(commandChain[0], state);
    auto it = commandChain.begin();
    std::advance(it, 1);
    while(it != commandChain.end()){
        resultPtr = modify(resultPtr, *it, state);
        it++;
    }

    if(!funcArgs.empty()){
        resultPtr = performFunction(resultPtr, funcName, funcArgs, state);
    }

    return *resultPtr;
}

max_stars_repo_name: rsehmbi/social--gaming

max_issues_repo_name: rsehmbi/social--gaming

max_forks_repo_name: rsehmbi/social--gaming


text:
#include <RDBoost/python.h>
#include <RDBoost/Wrap.h>

#include <GraphMol/Atom.h>
#include <GraphMol/AtomIterators.h>
#include <GraphMol/Bond.h>
#include <GraphMol/BondIterators.h>
#include <GraphMol/GraphMol.h>
#include <GraphMol/MolOps.h>
#include <cstdint>
#include <exception>

#include <boost/python.hpp>

namespace py = boost::python;

namespace {

const int NumAtomSymbols = 23;
const int AtomFeatureDimension = NumAtomSymbols + 6 + 5 + 4 + 1;
const int BondFeatureDimension = 5 + 6;

class BufferHolder : public Py_buffer {
public:
    BufferHolder(PyObject* obj, int additional_flags = 0) {
        if(PyObject_CheckBuffer(obj) != 1) {
            throw std::runtime_error("object must support buffer protocol");
        }

        if(PyObject_GetBuffer(obj, this, PyBUF_STRIDES | PyBUF_FORMAT | additional_flags) != 0) {
            throw std::runtime_error("object does not support strided buffer with given requirements");
        }
    }

    ~BufferHolder() {
        PyBuffer_Release(this);
    }
};

BufferHolder get_buffer(PyObject* obj, int additional_flags = 0) {
    return BufferHolder(obj, additional_flags);
}

bool check_format_i32(const char* format) { return format && (format[0] == 'i' || format[0] == 'l'); }
bool check_format_f32(const char* format) { return format && (format[0] == 'f'); }

int GetAtomSymbolOneHotIndex(int atomic_num) {
    switch (atomic_num) {
    case 6:
        return 0; // "C"
    case 7:
        return 1; // "N"
    case 8:
        return 2; // "O"
    case 16:
        return 3; // "S"
    case 9:
        return 4; // "F"
    case 14:
        return 5; // "Si"
    case 15:
        return 6; // "P"
    case 17:
        return 7; // "Cl"
    case 35:
        return 8; // "Br"
    case 12:
        return 9; // "Mg"
    case 11:
        return 10; // "Na"
    case 20:
        return 11; // "Ca"
    case 26:
        return 12; // "Fe"
    case 13:
        return 13; // "Al"
    case 53:
        return 14; // "I"
    case 5:
        return 15; // "B"
    case 19:
        return 16; // "K"
    case 34:
        return 17; // "Se"
    case 30:
        return 18; // "Zn"
    case 1:
        return 19; // "H"
    case 29:
        return 20; // "Cu"
    case 25:
        return 21; // "Mn"
    default:
        return 22; // "unknown"
    }
}

int GetAtomFormalChargeIndex(int formalCharge) {
    switch (formalCharge) {
    case -1:
        return 0;
    case -2:
        return 1;
    case 1:
        return 2;
    case 2:
        return 3;
    case 0:
    default:
        return 4;
    }
}

template <typename T> void atom_features(T *data, RDKit::Atom const &atom) {
    int offset = 0;

    data[GetAtomSymbolOneHotIndex(atom.getAtomicNum())] = 1;
    offset += NumAtomSymbols;

    data[std::min(atom.getDegree(), 5u) + offset] = 1;
    offset += 6;

    data[GetAtomFormalChargeIndex(atom.getFormalCharge()) + offset] = 1;
    offset += 5;

    data[static_cast<int>(atom.getChiralTag()) + offset] = 1;
    offset += 4;

    data[offset] = atom.getIsAromatic() ? 1 : 0;
};

void fill_atom_features(PyObject* view, RDKit::ROMol const& mol) {
    auto info = get_buffer(view, PyBUF_WRITABLE);

    if (info.ndim != 2) {
        throw std::runtime_error("Buffer must be 2-dimensional");
    }

    if (info.shape[1] < AtomFeatureDimension || info.shape[0] < mol.getNumAtoms()) {
        throw std::runtime_error("Buffer too small to fit specified molecule!");
    }

    if (info.strides[1] != info.itemsize) {
        throw std::runtime_error("Buffer must have C-contiguous layout");
    }

    if (!check_format_f32(info.format)) {
        throw std::runtime_error("Buffer must be of float type");
    }

    auto strides_i = info.strides[0] / info.itemsize;
    auto ptr = static_cast<float *>(info.buf);

    int i = 0;

    for (auto it1 = mol.beginAtoms(), e1 = mol.endAtoms(); it1 != e1; ++i, ++it1) {
        atom_features(ptr + i * strides_i, **it1);
    }
}

// Writes the bond type one-hot vector to the given data array.
template <typename T> void bond_type_feature(T *data, RDKit::Bond::BondType type) {
    switch (type) {
    case RDKit::Bond::BondType::SINGLE:
        data[0] = 1;
        return;
    case RDKit::Bond::BondType::DOUBLE:
        data[1] = 1;
        return;
    case RDKit::Bond::BondType::TRIPLE:
        data[2] = 1;
        return;
    case RDKit::Bond::BondType::AROMATIC:
        data[3] = 1;
        return;
    default:
        return;
    }
}

// writes the bond features to the given data array.
template <typename T> void bond_features(T *data, RDKit::Bond const &bond, RDKit::ROMol const& mol) {
    int offset = 0;

    bond_type_feature(data + offset, bond.getBondType());
    offset += 4;

    data[offset] = mol.getRingInfo()->numBondRings(bond.getIdx()) != 0;
    offset += 1;

    data[offset + std::min(static_cast<int>(bond.getStereo()), 5)] = 1;
}

void fill_bond_features(PyObject* view, RDKit::ROMol const* mol) {
    auto buffer = get_buffer(view, PyBUF_WRITABLE);

    if (!PyBuffer_IsContiguous(&buffer, 'C')) {
        throw std::runtime_error("Buffer most be C-contiguous");
    }

    if (buffer.ndim != 2) {
        throw std::runtime_error("Buffer must be 2-dimensional");
    }

    if (buffer.shape[1] < AtomFeatureDimension + BondFeatureDimension ||
        buffer.shape[0] < 2 * mol->getNumBonds()) {
        throw std::runtime_error("Buffer too small to fit specified molecule!");
    }

    if (!buffer.format || *buffer.format != 'f') {
        throw std::runtime_error("Buffer must be of float type");
    }

    {
        int i = 0;

        auto strides_i = buffer.strides[0] / buffer.itemsize;
        auto ptr = static_cast<float *>(buffer.buf);

        if (mol->getRingInfo() == nullptr) {
            throw std::runtime_error("molecule does not have ring info!");
        }

        if (mol->getRingInfo()->isInitialized()) {
            RDKit::MolOps::findSSSR(*mol);
        }

        for (auto it = mol->beginBonds(), e = mol->endBonds(); it != e; ++i, ++it) {
            // fill in bonds in both directions.
            atom_features(ptr + (2 * i) * strides_i, *(*it)->getBeginAtom());
            bond_features(ptr + (2 * i) * strides_i + AtomFeatureDimension, **it, *mol);

            atom_features(ptr + (2 * i + 1) * strides_i, *(*it)->getEndAtom());
            bond_features(ptr + (2 * i + 1) * strides_i + AtomFeatureDimension, **it, *mol);
        }
    }
}

void fill_atom_bond_list(PyObject* view, RDKit::ROMol &mol, int max_neighbours) {
    auto buffer = get_buffer(view, PyBUF_WRITABLE);

    if (!check_format_i32(buffer.format) || buffer.ndim != 2) {
        throw std::runtime_error("Incompatible buffer format!");
    }

    if (buffer.shape[0] < mol.getNumAtoms() || buffer.shape[1] < max_neighbours) {
        throw std::runtime_error("Buffer too small to fit specified molecule!");
    }

    {
        // release GIL once we acquire buffer pointer

        int i = 0;
        auto view_ptr = static_cast<std::int32_t *>(buffer.buf);
        auto strides_i = buffer.strides[0] / buffer.itemsize;
        auto strides_j = buffer.strides[1] / buffer.itemsize;

        for (auto it1 = mol.beginAtoms(), e1 = mol.endAtoms(); it1 != e1; ++i, ++it1) {
            RDKit::ROMol::OEDGE_ITER bond_begin, bond_end;
            std::tie(bond_begin, bond_end) = mol.getAtomBonds(*it1);

            for (int j = 0; bond_begin != bond_end; ++j, ++bond_begin) {
                const RDKit::Bond *bond = mol[*bond_begin];

                view_ptr[i * strides_i + j * strides_j] =
                    2 * bond->getIdx() + (bond->getBeginAtomIdx() == (*it1)->getIdx()) + 1;
            }
        }
    }
}

void fill_bond_incidence_list(PyObject* view, RDKit::ROMol &mol, int max_neighbours) {
    auto buffer = get_buffer(view, PyBUF_WRITABLE);

    if (!check_format_i32(buffer.format) || buffer.ndim != 2) {
        throw std::runtime_error("Incompatible buffer format!");
    }

    if (buffer.shape[0] < 2 * mol.getNumBonds() + 1 || buffer.shape[1] < max_neighbours) {
        throw std::runtime_error("Buffer too small to fit specified molecule!");
    }

    auto view_ptr = static_cast<std::int32_t *>(buffer.buf);
    auto strides_i = buffer.strides[0] / buffer.itemsize;
    auto strides_j = buffer.strides[1] / buffer.itemsize;

    int i = 0;

    {
        // release GIL once we acquire buffer pointer

        for (auto it1 = mol.beginBonds(), e1 = mol.endBonds(); it1 != e1; ++i, ++it1) {
            auto bond = *it1;

            auto fill_bond_incidence = [&](int index, const RDKit::Atom *atom) {
                RDKit::ROMol::OEDGE_ITER bond_begin, bond_end;
                std::tie(bond_begin, bond_end) = mol.getAtomBonds(atom);

                for (int j = 0; bond_begin != bond_end; ++j, ++bond_begin) {
                    const RDKit::Bond *bond2 = mol[*bond_begin];

                    if (bond2->getIdx() == bond->getIdx()) {
                        continue;
                    }

                    view_ptr[index * strides_i + j * strides_j] =
                        2 * bond2->getIdx() + (bond2->getBeginAtomIdx() == atom->getIdx()) + 1;
                }
            };

            fill_bond_incidence(2 * i + 1, bond->getBeginAtom());
            fill_bond_incidence(2 * i + 2, bond->getEndAtom());
        }
    }
}

void fill_atom_bond_list_sparse(PyObject* values, PyObject* index, const RDKit::ROMol &mol) {
    auto info_values = get_buffer(values, PyBUF_WRITABLE);
    auto info_index = get_buffer(index, PyBUF_WRITABLE);

    auto num_elements = mol.getNumBonds() * 2;

    if (!check_format_f32(info_values.format) || (info_values.shape[0] != num_elements)) {
        throw std::runtime_error(
            "Invalid values buffer. Must be of type float and length 2 * number of bonds.");
    }

    if (!check_format_i32(info_index.format) || info_index.ndim != 2) {
        throw std::runtime_error("Invalid index buffer. Must be of type int32 and dimension 2.");
    }

    if (info_index.shape[0] != 2 || info_index.shape[1] != num_elements) {
        throw std::runtime_error(
            "Invalid shape for index buffer. Must be of size 2 x (2 * num_bonds).");
    }

    auto values_ptr = static_cast<float *>(info_values.buf);
    auto values_stride = info_values.strides[0] / info_values.itemsize;

    auto index_ptr = static_cast<int32_t *>(info_index.buf);
    auto index_offset_0 = info_index.strides[0] / info_index.itemsize;
    auto index_stride = info_index.strides[1] / info_index.itemsize;

    {
        int i = 0;
        int current_index = 0;

        for (auto it1 = mol.beginAtoms(), e1 = mol.endAtoms(); it1 != e1; ++i, ++it1) {
            RDKit::ROMol::OEDGE_ITER bond_begin, bond_end;
            std::tie(bond_begin, bond_end) = mol.getAtomBonds(*it1);

            auto value = 1 / sqrt((*it1)->getDegree());

            for (int j = 0; bond_begin != bond_end; ++j, ++bond_begin) {
                const RDKit::Bond *bond = mol[*bond_begin];

                index_ptr[current_index * index_stride] = i;
                index_ptr[index_offset_0 + current_index * index_stride] =
                    2 * bond->getIdx() + (bond->getBeginAtomIdx() == (*it1)->getIdx());
                values_ptr[current_index * values_stride] = value;
                current_index += 1;
            }
        }
    }
}

int get_edge_incidence_size(RDKit::ROMol const &mol) {
    int result = 0;

    for (auto it1 = mol.beginAtoms(), it2 = mol.endAtoms(); it1 != it2; ++it1) {
        auto degree = (*it1)->getDegree();
        result += degree * (degree - 1);
    }

    return result;
}

void fill_bond_incidence_list_sparse(PyObject* values, PyObject* index, const RDKit::ROMol &mol) {
    auto info_values = get_buffer(values, PyBUF_WRITABLE);
    auto info_index = get_buffer(index, PyBUF_WRITABLE);

    if (!check_format_f32(info_values.format) || info_values.ndim != 1) {
        throw std::runtime_error("Values must be one-dimensional floating-point array.");
    }

    if (!check_format_i32(info_index.format) || info_index.ndim != 2) {
        throw std::runtime_error("index must be two-dimensional integer array.");
    }

    auto values_ptr = static_cast<float *>(info_values.buf);
    auto values_stride = info_values.strides[0] / info_values.itemsize;

    auto index_ptr = static_cast<int32_t *>(info_index.buf);
    auto index_offset_0 = info_index.strides[0] / info_index.itemsize;
    auto index_stride = info_index.strides[1] / info_index.itemsize;

    {
        // release GIL once we acquire buffer pointer
        int i = 0;
        int current_index = 0;

        for (auto it1 = mol.beginBonds(), e1 = mol.endBonds(); it1 != e1; ++i, ++it1) {
            auto bond = *it1;

            auto a1 = bond->getBeginAtom();
            auto a2 = bond->getEndAtom();

            float degree = a1->getDegree() + a2->getDegree() - 2;
            float value = 1 / sqrt(degree);

            auto fill_bond_incidence = [&](int index, const RDKit::Atom *atom) {
                RDKit::ROMol::OEDGE_ITER bond_begin, bond_end;
                std::tie(bond_begin, bond_end) = mol.getAtomBonds(atom);

                for (int j = 0; bond_begin != bond_end; ++j, ++bond_begin) {
                    const RDKit::Bond *bond2 = mol[*bond_begin];

                    if (bond2->getIdx() == bond->getIdx()) {
                        continue;
                    }

                    index_ptr[current_index * index_stride] = index;
                    index_ptr[index_offset_0 + current_index * index_stride] =
                        2 * bond2->getIdx() + (bond2->getBeginAtomIdx() == atom->getIdx());

                    values_ptr[current_index * values_stride] = value;
                    current_index += 1;
                }
            };

            fill_bond_incidence(2 * i, a1);
            fill_bond_incidence(2 * i + 1, a2);
        }
    }
}

void fill_atom_bond_list_segment(PyObject* scopes, PyObject* index, const RDKit::ROMol &mol) {
    auto info_scopes = get_buffer(scopes, PyBUF_WRITABLE);
    auto info_index = get_buffer(index, PyBUF_WRITABLE);

    if (!check_format_i32(info_scopes.format) || info_scopes.ndim != 2) {
        throw std::runtime_error("scopes must be two-dimensional integer array.");
    }

    if (!check_format_i32(info_index.format) || info_index.ndim != 1) {
        throw std::runtime_error("index must be one-dimensional integer array.");
    }

    if (info_scopes.shape[0] != mol.getNumAtoms()) {
        throw std::runtime_error("Scopes must be of length number of atoms.");
    }

    if (info_index.shape[0] != 2 * mol.getNumBonds()) {
        throw std::runtime_error("Index must be of length twice number of bonds.");
    }

    auto index_ptr = static_cast<int32_t *>(info_index.buf);
    auto index_stride = info_index.strides[0] / info_index.itemsize;

    auto scopes_ptr = static_cast<int32_t *>(info_scopes.buf);
    auto scopes_stride = info_scopes.strides[0] / info_scopes.itemsize;
    auto scopes_offset_1 = info_scopes.strides[1] / info_scopes.itemsize;

    {
        int offset = 0;
        int i = 0;
        int current_index = 0;

        for (auto it1 = mol.beginAtoms(), e1 = mol.endAtoms(); it1 != e1; ++i, ++it1) {
            RDKit::ROMol::OEDGE_ITER bond_begin, bond_end;
            std::tie(bond_begin, bond_end) = mol.getAtomBonds(*it1);

            auto value = 1 / sqrt((*it1)->getDegree());

            auto degree = (*it1)->getDegree();
            scopes_ptr[i * scopes_stride] = offset;
            scopes_ptr[i * scopes_stride + scopes_offset_1] = degree;
            offset += degree;

            for (int j = 0; bond_begin != bond_end; ++j, ++bond_begin) {
                const RDKit::Bond *bond = mol[*bond_begin];

                index_ptr[current_index * index_stride] =
                    2 * bond->getIdx() + (bond->getBeginAtomIdx() == (*it1)->getIdx());
                current_index += 1;
            }
        }
    }
}

void fill_bond_incidence_list_segment(PyObject* scopes, PyObject* index, const RDKit::ROMol &mol) {
    auto info_scopes = get_buffer(scopes, PyBUF_WRITABLE);
    auto info_index = get_buffer(index, PyBUF_WRITABLE);

    if (!check_format_i32(info_scopes.format) || info_scopes.ndim != 2) {
        throw std::runtime_error("scopes must be two-dimensional integer array.");
    }

    if (!check_format_i32(info_index.format) || info_index.ndim != 1) {
        throw std::runtime_error("index must be one-dimensional integer array.");
    }

    if (info_scopes.shape[0] != 2 * mol.getNumBonds()) {
        throw std::runtime_error("scopes must be of length 2 * number of bonds.");
    }

    auto index_ptr = static_cast<int32_t *>(info_index.buf);
    auto index_stride = info_index.strides[0] / info_index.itemsize;

    auto scopes_ptr = static_cast<int32_t *>(info_scopes.buf);
    auto scopes_stride = info_scopes.strides[0] / info_scopes.itemsize;
    auto scopes_offset_1 = info_scopes.strides[1] / info_scopes.itemsize;

    {
        // release GIL once we acquire buffer pointer

        int i = 0;
        int current_index = 0;
        int offset = 0;

        for (auto it1 = mol.beginBonds(), e1 = mol.endBonds(); it1 != e1; ++i, ++it1) {
            auto bond = *it1;

            auto a1 = bond->getBeginAtom();
            auto a2 = bond->getEndAtom();

            float degree = a1->getDegree() + a2->getDegree() - 2;
            float value = 1 / sqrt(degree);

            auto fill_bond_incidence = [&](int index, const RDKit::Atom *atom) {
                RDKit::ROMol::OEDGE_ITER bond_begin, bond_end;
                std::tie(bond_begin, bond_end) = mol.getAtomBonds(atom);
                auto num_bonds = 0;

                for (int j = 0; bond_begin != bond_end; ++j, ++bond_begin) {
                    const RDKit::Bond *bond2 = mol[*bond_begin];

                    if (bond2->getIdx() == bond->getIdx()) {
                        continue;
                    }

                    index_ptr[current_index * index_stride] =
                        2 * bond2->getIdx() + (bond2->getBeginAtomIdx() == atom->getIdx());
                    current_index += 1;
                    num_bonds += 1;
                }

                scopes_ptr[index * scopes_stride] = offset;
                scopes_ptr[index * scopes_stride + scopes_offset_1] = num_bonds;
                offset += num_bonds;
            };

            fill_bond_incidence(2 * i, a1);
            fill_bond_incidence(2 * i + 1, a2);
        }
    }
}

} // namespace

BOOST_PYTHON_MODULE(molecule_representation) {
    py::scope().attr("__doc__") = "Helpers for computing molecule representations";
    py::def("fill_atom_bond_list", &fill_atom_bond_list,
          "Fills the given buffer with atom-bond incidence information.");

    py::def("fill_bond_incidence_list", &fill_bond_incidence_list,
          "Fills the given buffer with bond-bond incidence information.");

    py::def("fill_atom_features", &fill_atom_features,
          "Fills the given buffer with atom features. \n\n"
          "The buffer must be a C-contiguous two-dimensional float buffer with number of rows at "
          "least the number of molecules mol, and number of columns at least "
          "`AtomFeatureDimension`.");

    py::def(
        "fill_bond_features", &fill_bond_features,
        "Fills the given buffer with bond features. \n\n"
        "The buffer must be a C-contiguous two-dimensional float buffer with number of columns at "
        "least `AtomFeatureDimension + BondFeatureDimension`, and number of rows at least "
        "twice the number of bonds plus one.");

    py::def("fill_atom_bond_list_sparse", &fill_atom_bond_list_sparse, (py::arg("values"),
          py::arg("index"), py::arg("mol")),
          "Fills the given values and index buffers with the atom-bond incidence information."
          "The two buffers represent a sparse weighted adjacency matrix in COO format.");

    py::def("fill_bond_incidence_list_sparse", &fill_bond_incidence_list_sparse, (py::arg("values"),
          py::arg("index"), py::arg("mol")),
          "Fills the given values and index buffers with the bond-bond incidence information."
          "The two buffers represent a sparse weighted adjacency matrix in COO format.");

    py::def("get_edge_incidence_size", &get_edge_incidence_size,
          "Computes the number of non-zero entries in the bond incidence adjacency graph.");

    py::scope().attr("AtomFeatureDimension") = AtomFeatureDimension;
    py::scope().attr("BondFeatureDimension") = BondFeatureDimension;

    py::def("fill_atom_bond_list_segment", &fill_atom_bond_list_segment, py::arg("scopes"),
          py::arg("index"), py::arg("mol"));

    py::def("fill_bond_incidence_list_segment", &fill_bond_incidence_list_segment, py::arg("scopes"),
          py::arg("index"), py::arg("mol"));
}


text:
// Distributed under the Boost Software License, Version 1.0. (See
// accompanying file LICENSE.txt
//
// (C) Copyright 2014 Pierre Talbot <ptalbot@hyc.io>

#ifndef NEEV_FLATTEN_ITERATOR_HPP
#define NEEV_FLATTEN_ITERATOR_HPP

#include <boost/integer_traits.hpp>

#include <boost/iterator.hpp>
#include <boost/iterator/iterator_adaptor.hpp>
#include <boost/iterator/iterator_categories.hpp>

#include <boost/type_traits/is_class.hpp>
#include <boost/static_assert.hpp>

#include <neev/iterator/detail/iterator_traits.hpp>

namespace neev
{
  static const std::size_t guess_value = boost::integer_traits<std::size_t>::const_max;

  template <class Iterator, std::size_t depth = guess_value>
  class flatten_iterator;

  namespace detail
  {
    template <class Iterator, std::size_t depth>
    struct flatten_iterator_category_tag
    {
      typedef typename minimum_nested_iterator_category<
        Iterator, depth, boost::bidirectional_traversal_tag
      >::type type;
    };

    template <class Iterator, std::size_t depth>
    struct flatten_iterator_base
    {
      typedef boost::iterator_adaptor<
          flatten_iterator<Iterator, depth>
        , Iterator
        , typename innermost_iterator_value_type<Iterator, depth>::type
        , typename flatten_iterator_category_tag<Iterator, depth>::type
      > type;
    };

    template <class Iterator, std::size_t depth, bool is_bidirectional =
      boost::is_convertible<
        typename flatten_iterator_category_tag<Iterator, depth>::type,
        boost::bidirectional_traversal_tag
      >::value
    >
    class bidirectional_base
    : public detail::flatten_iterator_base<Iterator, depth>::type
    {
      typedef typename detail::flatten_iterator_base<
        Iterator, depth
      >::type super_t;

    public:
      bidirectional_base(){}

      bidirectional_base(Iterator x)
      : super_t(x)
      , m_begin(x)
      {}

      Iterator begin() const { return m_begin; }
    protected:
      Iterator m_begin;
    };

    template <class Iterator, std::size_t depth>
    class bidirectional_base<Iterator, depth, false>
    : public detail::flatten_iterator_base<Iterator, depth>::type
    {
      typedef typename detail::flatten_iterator_base<
        Iterator, depth
      >::type super_t;
    public:
      bidirectional_base(){}
      bidirectional_base(Iterator x)
      : super_t(x)
      {}
    };

  } // namespace detail

  template <class Iterator>
  class flatten_iterator<Iterator, 0>
  : public detail::bidirectional_base<Iterator, 0>
  {
    typedef detail::bidirectional_base<Iterator, 0> super_t;

    typedef Iterator base_iterator_type;

    typedef flatten_iterator<Iterator, 0> this_iterator;

    friend class boost::iterator_core_access;

    template<typename Ite, std::size_t d>
    friend class flatten_iterator;

   public:
    flatten_iterator() { }

    flatten_iterator(Iterator x, Iterator end_ = Iterator())
    : super_t(x)
    , m_end(end_)
    {}

    template<class OtherIterator>
    flatten_iterator(flatten_iterator<OtherIterator, 0> const& t
    , typename boost::enable_if_convertible<OtherIterator, Iterator>::type* = 0)
    : super_t(t.base()), m_end(t.end()) {}

    Iterator end() const { return m_end; }

   private:
    bool decrement_impl()
    {
      if(this->base() != this->begin())
      {
        --(this->base_reference());
        return true;
      }
      return false;
    }

    bool equal(this_iterator const& x) const
    {
      return this->base() == x.base();
    }

    Iterator m_end;
  };

  template <class Iterator>
  class flatten_iterator<Iterator, guess_value>
  : public flatten_iterator<Iterator,
      detail::iterator_nested_depth<Iterator>::value>
  {
    typedef flatten_iterator<
      Iterator,
      detail::iterator_nested_depth<Iterator>::value
    > super_t;

   public:
    flatten_iterator() { }

    flatten_iterator(Iterator x, Iterator end_ = Iterator())
    : super_t(x, end_)
    {}
  };

  template <class Iterator, std::size_t depth>
  class flatten_iterator
  : public detail::bidirectional_base<Iterator, depth>
  {
    typedef detail::bidirectional_base<Iterator, depth> super_t;

    typedef flatten_iterator<
      typename detail::inner_iterator<Iterator>::type,
      depth - 1
    > inner_iterator_type;

    typedef Iterator base_iterator_type;

    typedef flatten_iterator<Iterator, depth> this_iterator;

    friend class boost::iterator_core_access;

    template<typename Ite, std::size_t d>
    friend class flatten_iterator;

   public:
      flatten_iterator() {}

      flatten_iterator(Iterator x, Iterator end_ = Iterator())
      : super_t(x)
      , m_end(end_)
      {
        if(this->base_reference() != end())
        {
          update_inner();
          skip_forward_empty_inner();
        }
      }

      template<class OtherIterator>
      flatten_iterator(flatten_iterator<OtherIterator, depth> const& t
      , typename boost::enable_if_convertible<OtherIterator, Iterator>::type* = 0)
      : super_t(t.base()), m_end(t.end()) {}

      Iterator end() const { return m_end; }

   private:
      void update_inner()
      {
        inner = inner_iterator_type(
          this->base_reference()->begin(),
          this->base_reference()->end());
      }

      void skip_forward_empty_inner()
      {
        if(this->base_reference() != end())
        {
          while(inner.base() == inner.end())
          {
            ++(this->base_reference());
            if(this->base_reference() != end())
              update_inner();
            else break;
          }
        }
      }

      void increment()
      {
        ++inner;
        skip_forward_empty_inner();
      }

      bool update_inner_to_end()
      {
        typedef typename inner_iterator_type::base_iterator_type inner_base_iterator;
        inner_base_iterator b = this->base()->begin();
        inner_base_iterator e = this->base()->end();
        if(b == e)
          return false;
        else
        {
          inner = inner_iterator_type(e, e);
          inner.m_begin = b;
          inner.decrement_impl();
          return true;
        }
      }

      bool decrement_impl()
      {
        if(inner.decrement_impl())
          return true;
        else
        {
          while(this->base_reference() != this->begin())
          {
            --(this->base_reference());
            if(update_inner_to_end())
              return true;
          }
          return false;
        }
      }

      void decrement()
      {
        decrement_impl();
      }

      bool equal(this_iterator const& x) const
      {
        return this->base() == x.base() && inner.equal(x.inner);
      }

      typename super_t::reference dereference() const
      {
        return *inner;
      }

      inner_iterator_type inner;
      Iterator m_end;
  };

  template <class Iterator>
  flatten_iterator<Iterator>
  make_flatten_iterator(Iterator x, Iterator end = Iterator())
  {
      return flatten_iterator<Iterator>(x,end);
  }

} // namespace neev

#endif // NEEV_FLATTEN_ITERATOR_HPP

ext: hpp

max_stars_repo_name: ptal/neev

max_issues_repo_name: ptal/neev

max_forks_repo_name: ptal/neev


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

#ifndef BEAST_WEBSOCKET_DETAIL_FRAME_HPP
#define BEAST_WEBSOCKET_DETAIL_FRAME_HPP

#include <beast/websocket/rfc6455.hpp>
#include <beast/websocket/detail/endian.hpp>
#include <beast/websocket/detail/utf8_checker.hpp>
#include <beast/core/consuming_buffers.hpp>
#include <beast/core/static_streambuf.hpp>
#include <beast/core/static_string.hpp>
#include <boost/asio/buffer.hpp>
#include <boost/endian/buffers.hpp>
#include <cassert>
#include <cstdint>

namespace beast {
namespace websocket {
namespace detail {

/// Identifies the role of a WebSockets stream.
enum class role_type
{
    /// Stream is operating as a client.
    client,

    /// Stream is operating as a server.
    server
};

// Contents of a WebSocket frame header
struct frame_header
{
    opcode op;
    bool fin;
    bool mask;
    bool rsv1;
    bool rsv2;
    bool rsv3;
    std::uint64_t len;
    std::uint32_t key;
};

// holds the largest possible frame header
using fh_streambuf =
    static_streambuf_n<14>;

// holds the largest possible control frame
using frame_streambuf =
    static_streambuf_n< 2 + 8 + 4 + 125 >;

inline
bool constexpr
is_reserved(opcode op)
{
    return
        (op >= opcode::rsv3  && op <= opcode::rsv7) ||
        (op >= opcode::crsvb && op <= opcode::crsvf);
}

inline
bool constexpr
is_valid(opcode op)
{
    return op <= opcode::crsvf;
}

inline
bool constexpr
is_control(opcode op)
{
    return op >= opcode::close;
}

// Returns `true` if a close code is valid
inline
bool
is_valid(close_code::value code)
{
    auto const v = code;
    switch(v)
    {
    case 1000:
    case 1001:
    case 1002:
    case 1003:
    case 1007:
    case 1008:
    case 1009:
    case 1010:
    case 1011:
    case 1012:
    case 1013:
        return true;

    // explicitly reserved
    case 1004:
    case 1005:
    case 1006:
    case 1014:
    case 1015:
        return false;
    }
    // reserved
    if(v >= 1016 && v <= 2999)
        return false;
    // not used
    if(v >= 0 && v <= 999)
        return false;
    return true;
}

//------------------------------------------------------------------------------

// Write frame header to dynamic buffer
//
template<class DynamicBuffer>
void
write(DynamicBuffer& db, frame_header const& fh)
{
    using boost::asio::buffer;
    using boost::asio::buffer_copy;
    using namespace boost::endian;
    std::size_t n;
    std::uint8_t b[14];
    b[0] = (fh.fin ? 0x80 : 0x00) | static_cast<std::uint8_t>(fh.op);
    if(fh.rsv1)
        b[0] |= 0x40;
    if(fh.rsv2)
        b[0] |= 0x20;
    if(fh.rsv3)
        b[0] |= 0x10;
    b[1] = fh.mask ? 0x80 : 0x00;
    if (fh.len <= 125)
    {
        b[1] |= fh.len;
        n = 2;
    }
    else if (fh.len <= 65535)
    {
        b[1] |= 126;
        ::new(&b[2]) big_uint16_buf_t{
            (std::uint16_t)fh.len};
        n = 4;
    }
    else
    {
        b[1] |= 127;
        ::new(&b[2]) big_uint64_buf_t{fh.len};
        n = 10;
    }
    if(fh.mask)
    {
        native_to_little_uint32(fh.key, &b[n]);
        n += 4;
    }
    db.commit(buffer_copy(
        db.prepare(n), buffer(b)));
}

// Read fixed frame header
// Requires at least 2 bytes
//
template<class DynamicBuffer>
std::size_t
read_fh1(frame_header& fh, DynamicBuffer& db,
    role_type role, close_code::value& code)
{
    using boost::asio::buffer;
    using boost::asio::buffer_copy;
    using boost::asio::buffer_size;
    std::uint8_t b[2];
    assert(buffer_size(db.data()) >= sizeof(b));
    db.consume(buffer_copy(buffer(b), db.data()));
    std::size_t need;
    fh.len = b[1] & 0x7f;
    switch(fh.len)
    {
        case 126: need = 2; break;
        case 127: need = 8; break;
        default:
            need = 0;
    }
    fh.mask = (b[1] & 0x80) != 0;
    if(fh.mask)
        need += 4;
    fh.op   = static_cast<opcode>(b[0] & 0x0f);
    fh.fin  = (b[0] & 0x80) != 0;
    fh.rsv1 = (b[0] & 0x40) != 0;
    fh.rsv2 = (b[0] & 0x20) != 0;
    fh.rsv3 = (b[0] & 0x10) != 0;
    // invalid length for control message
    if(is_control(fh.op) && fh.len > 125)
    {
        code = close_code::protocol_error;
        return 0;
    }
    // reserved bits not cleared
    if(fh.rsv1 || fh.rsv2 || fh.rsv3)
    {
        code = close_code::protocol_error;
        return 0;
    }
    // reserved opcode
    if(is_reserved(fh.op))
    {
        code = close_code::protocol_error;
        return 0;
    }
    // fragmented control message
    if(is_control(fh.op) && ! fh.fin)
    {
        code = close_code::protocol_error;
        return 0;
    }
    // unmasked frame from client
    if(role == role_type::server && ! fh.mask)
    {
        code = close_code::protocol_error;
        return 0;
    }
    // masked frame from server
    if(role == role_type::client && fh.mask)
    {
        code = close_code::protocol_error;
        return 0;
    }
    code = close_code::none;
    return need;
}

// Decode variable frame header from stream
//
template<class DynamicBuffer>
void
read_fh2(frame_header& fh, DynamicBuffer& db,
    role_type role, close_code::value& code)
{
    using boost::asio::buffer;
    using boost::asio::buffer_copy;
    using boost::asio::buffer_size;
    using namespace boost::endian;
    switch(fh.len)
    {
    case 126:
    {
        std::uint8_t b[2];
        assert(buffer_size(db.data()) >= sizeof(b));
        db.consume(buffer_copy(buffer(b), db.data()));
        fh.len = big_uint16_to_native(&b[0]);
        // length not canonical
        if(fh.len < 126)
        {
            code = close_code::protocol_error;
            return;
        }
        break;
    }
    case 127:
    {
        std::uint8_t b[8];
        assert(buffer_size(db.data()) >= sizeof(b));
        db.consume(buffer_copy(buffer(b), db.data()));
        fh.len = big_uint64_to_native(&b[0]);
        // length not canonical
        if(fh.len < 65536)
        {
            code = close_code::protocol_error;
            return;
        }
        break;
    }
    }
    if(fh.mask)
    {
        std::uint8_t b[4];
        assert(buffer_size(db.data()) >= sizeof(b));
        db.consume(buffer_copy(buffer(b), db.data()));
        fh.key = little_uint32_to_native(&b[0]);
    }
    else
    {
        // initialize this otherwise operator== breaks
        fh.key = 0;
    }
    code = close_code::none;
}

// Read data from buffers
// This is for ping and pong payloads
//
template<class Buffers>
void
read(ping_data& data, Buffers const& bs)
{
    using boost::asio::buffer_copy;
    using boost::asio::buffer_size;
    using boost::asio::mutable_buffers_1;
    assert(buffer_size(bs) <= data.max_size());
    data.resize(buffer_size(bs));
    buffer_copy(mutable_buffers_1{
        data.data(), data.size()}, bs);
}

// Read close_reason, return true on success
// This is for the close payload
//
template<class Buffers>
void
read(close_reason& cr,
    Buffers const& bs, close_code::value& code)
{
    using boost::asio::buffer;
    using boost::asio::buffer_copy;
    using boost::asio::buffer_size;
    using namespace boost::endian;
    auto n = buffer_size(bs);
    assert(n <= 125);
    if(n == 0)
    {
        cr = close_reason{};
        code = close_code::none;
        return;
    }
    if(n == 1)
    {
        code = close_code::protocol_error;
        return;
    }
    consuming_buffers<Buffers> cb(bs);
    {
        std::uint8_t b[2];
        buffer_copy(buffer(b), cb);
        cr.code = big_uint16_to_native(&b[0]);
        cb.consume(2);
        n -= 2;
        if(! is_valid(cr.code))
        {
            code = close_code::protocol_error;
            return;
        }
    }
    if(n > 0)
    {
        cr.reason.resize(n);
        buffer_copy(buffer(&cr.reason[0], n), cb);
        if(! detail::check_utf8(
            cr.reason.data(), cr.reason.size()))
        {
            code = close_code::protocol_error;
            return;
        }
    }
    else
    {
        cr.reason = "";
    }
    code = close_code::none;
}

} // detail
} // websocket
} // beast

#endif

max_issues_repo_name: somayeghahari/beast

max_forks_repo_name: somayeghahari/beast


text:
/*
 * @file
 * @author University of Warwick
 * @version 1.0
 *
 * @section LICENSE
 *
 * @section DESCRIPTION
 *
 * Unit Tests for the methods of the Hexahedral class
 */

#define BOOST_TEST_MODULE Hexahedron
#include <boost/test/unit_test.hpp>
#include <boost/test/output_test_stream.hpp>
#include <stdexcept>

#include "Hexahedron.h"
#include "EuclideanPoint.h"

using namespace cupcfd::geometry::shapes;
namespace euc = cupcfd::geometry::euclidean;
namespace utf = boost::unit_test;

// === Constructor ===
// Test 1: Construct a regular hexahedron
BOOST_AUTO_TEST_CASE(constructor_test1, * utf::tolerance(0.00001))
{
	euc::EuclideanPoint<double,3> tlf(0.0, 5.0, 12.0);
	euc::EuclideanPoint<double,3> trf(5.0, 5.0, 12.0);
	euc::EuclideanPoint<double,3> brf(5.0, 5.0, 0.0);
	euc::EuclideanPoint<double,3> blf(0.0, 5.0, 0.0);
	euc::EuclideanPoint<double,3> tlb(0.0, 0.0, 12.0);
	euc::EuclideanPoint<double,3> trb(5.0, 0.0, 12.0);
	euc::EuclideanPoint<double,3> brb(5.0, 0.0, 0.0);
	euc::EuclideanPoint<double,3> blb(0.0, 0.0, 0.0);
	Quadrilateral3D<double> top(trf, tlf, tlb, trb);
	Quadrilateral3D<double> bottom(brf, brb, blb, blf);
	Quadrilateral3D<double> front(trf, brf, blf, tlf);
	Quadrilateral3D<double> rear(tlb, blb, brb, trb);
	Quadrilateral3D<double> left(tlb, tlf, blf, blb);
	Quadrilateral3D<double> right(trf, trb, brb, brf);
	Hexahedron<double> shape(top, bottom, front, rear, left, right);

	Quadrilateral3D<double> faces[6] = {top, bottom, front, rear, left, right};
	for (int f=0; f<6; f++) {
		for (int v=0; v<4; v++) {
			for (int c=0; c<3; c++) {
				BOOST_TEST(shape.faces[f].vertices[v].cmp[c] == faces[f].vertices[v].cmp[c]);
			}
		}
	}
}

// === isPointInside ===
// Test 1: Point is inside Hexahedron
BOOST_AUTO_TEST_CASE(isPointInside_test1, * utf::tolerance(0.00001))
{
	euc::EuclideanPoint<double,3> tlf(0.0, 5.0, 12.0);
	euc::EuclideanPoint<double,3> trf(5.0, 5.0, 12.0);
	euc::EuclideanPoint<double,3> blf(0.0, 5.0, 0.0);
	euc::EuclideanPoint<double,3> brf(5.0, 5.0, 0.0);
	euc::EuclideanPoint<double,3> tlb(0.0, 0.0, 12.0);
	euc::EuclideanPoint<double,3> trb(5.0, 0.0, 12.0);
	euc::EuclideanPoint<double,3> blb(0.0, 0.0, 0.0);
	euc::EuclideanPoint<double,3> brb(5.0, 0.0, 0.0);
	Quadrilateral3D<double> top(trf, tlf, tlb, trb);
	Quadrilateral3D<double> bottom(brf, brb, blb, blf);
	Quadrilateral3D<double> front(trf, brf, blf, tlf);
	Quadrilateral3D<double> rear(tlb, blb, brb, trb);
	Quadrilateral3D<double> left(tlb, tlf, blf, blb);
	Quadrilateral3D<double> right(trf, trb, brb, brf);
	Hexahedron<double> shape(top, bottom, front, rear, left, right);

	euc::EuclideanPoint<double,3> point(2.3, 4.6, 6.7);

	bool isPointInside = shape.isPointInside(point);
	BOOST_CHECK_EQUAL(isPointInside, true);
}

// Test 2: Point is outside Hexahedron
BOOST_AUTO_TEST_CASE(isPointInside_test2, * utf::tolerance(0.00001))
{
	euc::EuclideanPoint<double,3> tlf(0.0, 5.0, 12.0);
	euc::EuclideanPoint<double,3> trf(5.0, 5.0, 12.0);
	euc::EuclideanPoint<double,3> blf(0.0, 5.0, 0.0);
	euc::EuclideanPoint<double,3> brf(5.0, 5.0, 0.0);
	euc::EuclideanPoint<double,3> tlb(0.0, 0.0, 12.0);
	euc::EuclideanPoint<double,3> trb(5.0, 0.0, 12.0);
	euc::EuclideanPoint<double,3> blb(0.0, 0.0, 0.0);
	euc::EuclideanPoint<double,3> brb(5.0, 0.0, 0.0);
	Quadrilateral3D<double> top(trf, tlf, tlb, trb);
	Quadrilateral3D<double> bottom(brf, brb, blb, blf);
	Quadrilateral3D<double> front(trf, brf, blf, tlf);
	Quadrilateral3D<double> rear(tlb, blb, brb, trb);
	Quadrilateral3D<double> left(tlb, tlf, blf, blb);
	Quadrilateral3D<double> right(trf, trb, brb, brf);
	Hexahedron<double> shape(top, bottom, front, rear, left, right);

	euc::EuclideanPoint<double,3> point(2.3, 4.6, 13.7);

	bool isPointInside = shape.isPointInside(point);
	BOOST_CHECK_EQUAL(isPointInside, false);
}

// Test 3: Point is on Vertex tlf
BOOST_AUTO_TEST_CASE(isPointInside_test3, * utf::tolerance(0.00001))
{
	euc::EuclideanPoint<double,3> tlf(0.0, 5.0, 12.0);
	euc::EuclideanPoint<double,3> trf(5.0, 5.0, 12.0);
	euc::EuclideanPoint<double,3> blf(0.0, 5.0, 0.0);
	euc::EuclideanPoint<double,3> brf(5.0, 5.0, 0.0);
	euc::EuclideanPoint<double,3> tlb(0.0, 0.0, 12.0);
	euc::EuclideanPoint<double,3> trb(5.0, 0.0, 12.0);
	euc::EuclideanPoint<double,3> blb(0.0, 0.0, 0.0);
	euc::EuclideanPoint<double,3> brb(5.0, 0.0, 0.0);
	Quadrilateral3D<double> top(trf, tlf, tlb, trb);
	Quadrilateral3D<double> bottom(brf, brb, blb, blf);
	Quadrilateral3D<double> front(trf, brf, blf, tlf);
	Quadrilateral3D<double> rear(tlb, blb, brb, trb);
	Quadrilateral3D<double> left(tlb, tlf, blf, blb);
	Quadrilateral3D<double> right(trf, trb, brb, brf);
	Hexahedron<double> shape(top, bottom, front, rear, left, right);

	bool isPointInside = shape.isPointInside(tlf);
	BOOST_CHECK_EQUAL(isPointInside, true);
}

// Test 4: Point is on Vertex trf
BOOST_AUTO_TEST_CASE(isPointInside_test4, * utf::tolerance(0.00001))
{
	euc::EuclideanPoint<double,3> tlf(0.0, 5.0, 12.0);
	euc::EuclideanPoint<double,3> trf(5.0, 5.0, 12.0);
	euc::EuclideanPoint<double,3> blf(0.0, 5.0, 0.0);
	euc::EuclideanPoint<double,3> brf(5.0, 5.0, 0.0);
	euc::EuclideanPoint<double,3> tlb(0.0, 0.0, 12.0);
	euc::EuclideanPoint<double,3> trb(5.0, 0.0, 12.0);
	euc::EuclideanPoint<double,3> blb(0.0, 0.0, 0.0);
	euc::EuclideanPoint<double,3> brb(5.0, 0.0, 0.0);
	Quadrilateral3D<double> top(trf, tlf, tlb, trb);
	Quadrilateral3D<double> bottom(brf, brb, blb, blf);
	Quadrilateral3D<double> front(trf, brf, blf, tlf);
	Quadrilateral3D<double> rear(tlb, blb, brb, trb);
	Quadrilateral3D<double> left(tlb, tlf, blf, blb);
	Quadrilateral3D<double> right(trf, trb, brb, brf);
	Hexahedron<double> shape(top, bottom, front, rear, left, right);

	bool isPointInside = shape.isPointInside(trf);
	BOOST_CHECK_EQUAL(isPointInside, true);
}

// Test 5: Point is on Vertex blf
BOOST_AUTO_TEST_CASE(isPointInside_test5, * utf::tolerance(0.00001))
{
	euc::EuclideanPoint<double,3> tlf(0.0, 5.0, 12.0);
	euc::EuclideanPoint<double,3> trf(5.0, 5.0, 12.0);
	euc::EuclideanPoint<double,3> blf(0.0, 5.0, 0.0);
	euc::EuclideanPoint<double,3> brf(5.0, 5.0, 0.0);
	euc::EuclideanPoint<double,3> tlb(0.0, 0.0, 12.0);
	euc::EuclideanPoint<double,3> trb(5.0, 0.0, 12.0);
	euc::EuclideanPoint<double,3> blb(0.0, 0.0, 0.0);
	euc::EuclideanPoint<double,3> brb(5.0, 0.0, 0.0);
	Quadrilateral3D<double> top(trf, tlf, tlb, trb);
	Quadrilateral3D<double> bottom(brf, brb, blb, blf);
	Quadrilateral3D<double> front(trf, brf, blf, tlf);
	Quadrilateral3D<double> rear(tlb, blb, brb, trb);
	Quadrilateral3D<double> left(tlb, tlf, blf, blb);
	Quadrilateral3D<double> right(trf, trb, brb, brf);
	Hexahedron<double> shape(top, bottom, front, rear, left, right);

	bool isPointInside = shape.isPointInside(blf);
	BOOST_CHECK_EQUAL(isPointInside, true);
}

// Test 6: Point is on Vertex brf
BOOST_AUTO_TEST_CASE(isPointInside_test6, * utf::tolerance(0.00001))
{
	euc::EuclideanPoint<double,3> tlf(0.0, 5.0, 12.0);
	euc::EuclideanPoint<double,3> trf(5.0, 5.0, 12.0);
	euc::EuclideanPoint<double,3> blf(0.0, 5.0, 0.0);
	euc::EuclideanPoint<double,3> brf(5.0, 5.0, 0.0);
	euc::EuclideanPoint<double,3> tlb(0.0, 0.0, 12.0);
	euc::EuclideanPoint<double,3> trb(5.0, 0.0, 12.0);
	euc::EuclideanPoint<double,3> blb(0.0, 0.0, 0.0);
	euc::EuclideanPoint<double,3> brb(5.0, 0.0, 0.0);
	Quadrilateral3D<double> top(trf, tlf, tlb, trb);
	Quadrilateral3D<double> bottom(brf, brb, blb, blf);
	Quadrilateral3D<double> front(trf, brf, blf, tlf);
	Quadrilateral3D<double> rear(tlb, blb, brb, trb);
	Quadrilateral3D<double> left(tlb, tlf, blf, blb);
	Quadrilateral3D<double> right(trf, trb, brb, brf);
	Hexahedron<double> shape(top, bottom, front, rear, left, right);

	bool isPointInside = shape.isPointInside(brf);
	BOOST_CHECK_EQUAL(isPointInside, true);
}

// Test 7: Point is on Vertex tlb
BOOST_AUTO_TEST_CASE(isPointInside_test7, * utf::tolerance(0.00001))
{
	euc::EuclideanPoint<double,3> tlf(0.0, 5.0, 12.0);
	euc::EuclideanPoint<double,3> trf(5.0, 5.0, 12.0);
	euc::EuclideanPoint<double,3> blf(0.0, 5.0, 0.0);
	euc::EuclideanPoint<double,3> brf(5.0, 5.0, 0.0);
	euc::EuclideanPoint<double,3> tlb(0.0, 0.0, 12.0);
	euc::EuclideanPoint<double,3> trb(5.0, 0.0, 12.0);
	euc::EuclideanPoint<double,3> blb(0.0, 0.0, 0.0);
	euc::EuclideanPoint<double,3> brb(5.0, 0.0, 0.0);
	Quadrilateral3D<double> top(trf, tlf, tlb, trb);
	Quadrilateral3D<double> bottom(brf, brb, blb, blf);
	Quadrilateral3D<double> front(trf, brf, blf, tlf);
	Quadrilateral3D<double> rear(tlb, blb, brb, trb);
	Quadrilateral3D<double> left(tlb, tlf, blf, blb);
	Quadrilateral3D<double> right(trf, trb, brb, brf);
	Hexahedron<double> shape(top, bottom, front, rear, left, right);

	bool isPointInside = shape.isPointInside(tlb);
	BOOST_CHECK_EQUAL(isPointInside, true);
}

// Test 8: Point is on Vertex trb
BOOST_AUTO_TEST_CASE(isPointInside_test8, * utf::tolerance(0.00001))
{
	euc::EuclideanPoint<double,3> tlf(0.0, 5.0, 12.0);
	euc::EuclideanPoint<double,3> trf(5.0, 5.0, 12.0);
	euc::EuclideanPoint<double,3> blf(0.0, 5.0, 0.0);
	euc::EuclideanPoint<double,3> brf(5.0, 5.0, 0.0);
	euc::EuclideanPoint<double,3> tlb(0.0, 0.0, 12.0);
	euc::EuclideanPoint<double,3> trb(5.0, 0.0, 12.0);
	euc::EuclideanPoint<double,3> blb(0.0, 0.0, 0.0);
	euc::EuclideanPoint<double,3> brb(5.0, 0.0, 0.0);
	Quadrilateral3D<double> top(trf, tlf, tlb, trb);
	Quadrilateral3D<double> bottom(brf, brb, blb, blf);
	Quadrilateral3D<double> front(trf, brf, blf, tlf);
	Quadrilateral3D<double> rear(tlb, blb, brb, trb);
	Quadrilateral3D<double> left(tlb, tlf, blf, blb);
	Quadrilateral3D<double> right(trf, trb, brb, brf);
	Hexahedron<double> shape(top, bottom, front, rear, left, right);

	bool isPointInside = shape.isPointInside(trb);
	BOOST_CHECK_EQUAL(isPointInside, true);
}

// Test 9: Point is on Vertex blb
BOOST_AUTO_TEST_CASE(isPointInside_test9, * utf::tolerance(0.00001))
{
	euc::EuclideanPoint<double,3> tlf(0.0, 5.0, 12.0);
	euc::EuclideanPoint<double,3> trf(5.0, 5.0, 12.0);
	euc::EuclideanPoint<double,3> blf(0.0, 5.0, 0.0);
	euc::EuclideanPoint<double,3> brf(5.0, 5.0, 0.0);
	euc::EuclideanPoint<double,3> tlb(0.0, 0.0, 12.0);
	euc::EuclideanPoint<double,3> trb(5.0, 0.0, 12.0);
	euc::EuclideanPoint<double,3> blb(0.0, 0.0, 0.0);
	euc::EuclideanPoint<double,3> brb(5.0, 0.0, 0.0);
	Quadrilateral3D<double> top(trf, tlf, tlb, trb);
	Quadrilateral3D<double> bottom(brf, brb, blb, blf);
	Quadrilateral3D<double> front(trf, brf, blf, tlf);
	Quadrilateral3D<double> rear(tlb, blb, brb, trb);
	Quadrilateral3D<double> left(tlb, tlf, blf, blb);
	Quadrilateral3D<double> right(trf, trb, brb, brf);
	Hexahedron<double> shape(top, bottom, front, rear, left, right);

	bool isPointInside = shape.isPointInside(blb);
	BOOST_CHECK_EQUAL(isPointInside, true);
}

// Test 10: Point is on Vertex brb
BOOST_AUTO_TEST_CASE(isPointInside_test10, * utf::tolerance(0.00001))
{
	euc::EuclideanPoint<double,3> tlf(0.0, 5.0, 12.0);
	euc::EuclideanPoint<double,3> trf(5.0, 5.0, 12.0);
	euc::EuclideanPoint<double,3> blf(0.0, 5.0, 0.0);
	euc::EuclideanPoint<double,3> brf(5.0, 5.0, 0.0);
	euc::EuclideanPoint<double,3> tlb(0.0, 0.0, 12.0);
	euc::EuclideanPoint<double,3> trb(5.0, 0.0, 12.0);
	euc::EuclideanPoint<double,3> blb(0.0, 0.0, 0.0);
	euc::EuclideanPoint<double,3> brb(5.0, 0.0, 0.0);
	Quadrilateral3D<double> top(trf, tlf, tlb, trb);
	Quadrilateral3D<double> bottom(brf, brb, blb, blf);
	Quadrilateral3D<double> front(trf, brf, blf, tlf);
	Quadrilateral3D<double> rear(tlb, blb, brb, trb);
	Quadrilateral3D<double> left(tlb, tlf, blf, blb);
	Quadrilateral3D<double> right(trf, trb, brb, brf);
	Hexahedron<double> shape(top, bottom, front, rear, left, right);

	bool isPointInside = shape.isPointInside(brb);
	BOOST_CHECK_EQUAL(isPointInside, true);
}

// Test 11: Point is on Top Face
BOOST_AUTO_TEST_CASE(isPointInside_test11, * utf::tolerance(0.00001))
{
	euc::EuclideanPoint<double,3> tlf(0.0, 5.0, 12.0);
	euc::EuclideanPoint<double,3> trf(5.0, 5.0, 12.0);
	euc::EuclideanPoint<double,3> blf(0.0, 5.0, 0.0);
	euc::EuclideanPoint<double,3> brf(5.0, 5.0, 0.0);
	euc::EuclideanPoint<double,3> tlb(0.0, 0.0, 12.0);
	euc::EuclideanPoint<double,3> trb(5.0, 0.0, 12.0);
	euc::EuclideanPoint<double,3> blb(0.0, 0.0, 0.0);
	euc::EuclideanPoint<double,3> brb(5.0, 0.0, 0.0);
	Quadrilateral3D<double> top(trf, tlf, tlb, trb);
	Quadrilateral3D<double> bottom(brf, brb, blb, blf);
	Quadrilateral3D<double> front(trf, brf, blf, tlf);
	Quadrilateral3D<double> rear(tlb, blb, brb, trb);
	Quadrilateral3D<double> left(tlb, tlf, blf, blb);
	Quadrilateral3D<double> right(trf, trb, brb, brf);
	Hexahedron<double> shape(top, bottom, front, rear, left, right);

	euc::EuclideanPoint<double,3> point(2.3, 4.6, 12.0);

	bool isPointInside = shape.isPointInside(point);
	BOOST_CHECK_EQUAL(isPointInside, true);
}

// Test 12: Point is on Bottom Face
BOOST_AUTO_TEST_CASE(isPointInside_test12, * utf::tolerance(0.00001))
{
	euc::EuclideanPoint<double,3> tlf(0.0, 5.0, 12.0);
	euc::EuclideanPoint<double,3> trf(5.0, 5.0, 12.0);
	euc::EuclideanPoint<double,3> blf(0.0, 5.0, 0.0);
	euc::EuclideanPoint<double,3> brf(5.0, 5.0, 0.0);
	euc::EuclideanPoint<double,3> tlb(0.0, 0.0, 12.0);
	euc::EuclideanPoint<double,3> trb(5.0, 0.0, 12.0);
	euc::EuclideanPoint<double,3> blb(0.0, 0.0, 0.0);
	euc::EuclideanPoint<double,3> brb(5.0, 0.0, 0.0);
	Quadrilateral3D<double> top(trf, tlf, tlb, trb);
	Quadrilateral3D<double> bottom(brf, brb, blb, blf);
	Quadrilateral3D<double> front(trf, brf, blf, tlf);
	Quadrilateral3D<double> rear(tlb, blb, brb, trb);
	Quadrilateral3D<double> left(tlb, tlf, blf, blb);
	Quadrilateral3D<double> right(trf, trb, brb, brf);
	Hexahedron<double> shape(top, bottom, front, rear, left, right);

	euc::EuclideanPoint<double,3> point(2.3, 3.853, 0.0);

	bool isPointInside = shape.isPointInside(point);
	BOOST_CHECK_EQUAL(isPointInside, true);
}

// Test 13: Point is on Left Face
BOOST_AUTO_TEST_CASE(isPointInside_test13, * utf::tolerance(0.00001))
{
	euc::EuclideanPoint<double,3> tlf(0.0, 5.0, 12.0);
	euc::EuclideanPoint<double,3> trf(5.0, 5.0, 12.0);
	euc::EuclideanPoint<double,3> blf(0.0, 5.0, 0.0);
	euc::EuclideanPoint<double,3> brf(5.0, 5.0, 0.0);
	euc::EuclideanPoint<double,3> tlb(0.0, 0.0, 12.0);
	euc::EuclideanPoint<double,3> trb(5.0, 0.0, 12.0);
	euc::EuclideanPoint<double,3> blb(0.0, 0.0, 0.0);
	euc::EuclideanPoint<double,3> brb(5.0, 0.0, 0.0);
	Quadrilateral3D<double> top(trf, tlf, tlb, trb);
	Quadrilateral3D<double> bottom(brf, brb, blb, blf);
	Quadrilateral3D<double> front(trf, brf, blf, tlf);
	Quadrilateral3D<double> rear(tlb, blb, brb, trb);
	Quadrilateral3D<double> left(tlb, tlf, blf, blb);
	Quadrilateral3D<double> right(trf, trb, brb, brf);
	Hexahedron<double> shape(top, bottom, front, rear, left, right);

	euc::EuclideanPoint<double,3> point(0.0, 2.98, 5.2);

	bool isPointInside = shape.isPointInside(point);
	BOOST_CHECK_EQUAL(isPointInside, true);
}

// Test 14: Point is on Right Face
BOOST_AUTO_TEST_CASE(isPointInside_test14, * utf::tolerance(0.00001))
{
	euc::EuclideanPoint<double,3> tlf(0.0, 5.0, 12.0);
	euc::EuclideanPoint<double,3> trf(5.0, 5.0, 12.0);
	euc::EuclideanPoint<double,3> blf(0.0, 5.0, 0.0);
	euc::EuclideanPoint<double,3> brf(5.0, 5.0, 0.0);
	euc::EuclideanPoint<double,3> tlb(0.0, 0.0, 12.0);
	euc::EuclideanPoint<double,3> trb(5.0, 0.0, 12.0);
	euc::EuclideanPoint<double,3> blb(0.0, 0.0, 0.0);
	euc::EuclideanPoint<double,3> brb(5.0, 0.0, 0.0);
	Quadrilateral3D<double> top(trf, tlf, tlb, trb);
	Quadrilateral3D<double> bottom(brf, brb, blb, blf);
	Quadrilateral3D<double> front(trf, brf, blf, tlf);
	Quadrilateral3D<double> rear(tlb, blb, brb, trb);
	Quadrilateral3D<double> left(tlb, tlf, blf, blb);
	Quadrilateral3D<double> right(trf, trb, brb, brf);
	Hexahedron<double> shape(top, bottom, front, rear, left, right);

	euc::EuclideanPoint<double,3> point(5.0, 2.98, 5.2);

	bool isPointInside = shape.isPointInside(point);
	BOOST_CHECK_EQUAL(isPointInside, true);
}

// Test 15: Point is on Front Face
BOOST_AUTO_TEST_CASE(isPointInside_test15, * utf::tolerance(0.00001))
{
	euc::EuclideanPoint<double,3> tlf(0.0, 5.0, 12.0);
	euc::EuclideanPoint<double,3> trf(5.0, 5.0, 12.0);
	euc::EuclideanPoint<double,3> blf(0.0, 5.0, 0.0);
	euc::EuclideanPoint<double,3> brf(5.0, 5.0, 0.0);
	euc::EuclideanPoint<double,3> tlb(0.0, 0.0, 12.0);
	euc::EuclideanPoint<double,3> trb(5.0, 0.0, 12.0);
	euc::EuclideanPoint<double,3> blb(0.0, 0.0, 0.0);
	euc::EuclideanPoint<double,3> brb(5.0, 0.0, 0.0);
	Quadrilateral3D<double> top(trf, tlf, tlb, trb);
	Quadrilateral3D<double> bottom(brf, brb, blb, blf);
	Quadrilateral3D<double> front(trf, brf, blf, tlf);
	Quadrilateral3D<double> rear(tlb, blb, brb, trb);
	Quadrilateral3D<double> left(tlb, tlf, blf, blb);
	Quadrilateral3D<double> right(trf, trb, brb, brf);
	Hexahedron<double> shape(top, bottom, front, rear, left, right);

	euc::EuclideanPoint<double,3> point(1.63, 5.0, 5.2);

	bool isPointInside = shape.isPointInside(point);
	BOOST_CHECK_EQUAL(isPointInside, true);
}

// Test 16: Point is on Back Face
BOOST_AUTO_TEST_CASE(isPointInside_test16, * utf::tolerance(0.00001))
{
	euc::EuclideanPoint<double,3> tlf(0.0, 5.0, 12.0);
	euc::EuclideanPoint<double,3> trf(5.0, 5.0, 12.0);
	euc::EuclideanPoint<double,3> blf(0.0, 5.0, 0.0);
	euc::EuclideanPoint<double,3> brf(5.0, 5.0, 0.0);
	euc::EuclideanPoint<double,3> tlb(0.0, 0.0, 12.0);
	euc::EuclideanPoint<double,3> trb(5.0, 0.0, 12.0);
	euc::EuclideanPoint<double,3> blb(0.0, 0.0, 0.0);
	euc::EuclideanPoint<double,3> brb(5.0, 0.0, 0.0);
	Quadrilateral3D<double> top(trf, tlf, tlb, trb);
	Quadrilateral3D<double> bottom(brf, brb, blb, blf);
	Quadrilateral3D<double> front(trf, brf, blf, tlf);
	Quadrilateral3D<double> rear(tlb, blb, brb, trb);
	Quadrilateral3D<double> left(tlb, tlf, blf, blb);
	Quadrilateral3D<double> right(trf, trb, brb, brf);
	Hexahedron<double> shape(top, bottom, front, rear, left, right);

	euc::EuclideanPoint<double,3> point(1.63, 0.0, 5.2);

	bool isPointInside = shape.isPointInside(point);
	BOOST_CHECK_EQUAL(isPointInside, true);
}

// Test 17: Point is Outside
BOOST_AUTO_TEST_CASE(isPointInside_test17, * utf::tolerance(0.00001))
{
	euc::EuclideanPoint<double,3> tlf(0.0909091,0.5,0.454545);
	euc::EuclideanPoint<double,3> trf(0.0909091,0.5,0.545455);
	euc::EuclideanPoint<double,3> blf(0.0,0.5,0.454545);
	euc::EuclideanPoint<double,3> brf(0.0,0.5,0.545455);
	euc::EuclideanPoint<double,3> tlb(0.0909091,0.583333,0.454545);
	euc::EuclideanPoint<double,3> trb(0.0909091,0.583333,0.545455);
	euc::EuclideanPoint<double,3> blb(0.0,0.583333,0.454545);
	euc::EuclideanPoint<double,3> brb(0.0,0.583333,0.545455);
	Quadrilateral3D<double> top(trf, tlf, tlb, trb);
	Quadrilateral3D<double> bottom(brf, brb, blb, blf);
	Quadrilateral3D<double> front(trf, brf, blf, tlf);
	Quadrilateral3D<double> rear(tlb, blb, brb, trb);
	Quadrilateral3D<double> left(tlb, tlf, blf, blb);
	Quadrilateral3D<double> right(trf, trb, brb, brf);
	Hexahedron<double> shape(top, bottom, front, rear, left, right);

	euc::EuclideanPoint<double,3> point(0.23, 0.572, 0.5563);

	bool isPointInside = shape.isPointInside(point);
	BOOST_CHECK_EQUAL(isPointInside, false);
}

// Test 18: Point is on same plane as a face, but not in a face

// ToDo: Point on Edges Tests

// === getVolume ===
// Test 1: Compute volume of cube
BOOST_AUTO_TEST_CASE(getVolume_test1, * utf::tolerance(0.00001))
{
	euc::EuclideanPoint<double,3> tlf(0.0, 5.0, 12.0);
	euc::EuclideanPoint<double,3> trf(5.0, 5.0, 12.0);
	euc::EuclideanPoint<double,3> blf(0.0, 5.0, 0.0);
	euc::EuclideanPoint<double,3> brf(5.0, 5.0, 0.0);
	euc::EuclideanPoint<double,3> tlb(0.0, 0.0, 12.0);
	euc::EuclideanPoint<double,3> trb(5.0, 0.0, 12.0);
	euc::EuclideanPoint<double,3> blb(0.0, 0.0, 0.0);
	euc::EuclideanPoint<double,3> brb(5.0, 0.0, 0.0);
	Quadrilateral3D<double> top(trf, tlf, tlb, trb);
	Quadrilateral3D<double> bottom(brf, brb, blb, blf);
	Quadrilateral3D<double> front(trf, brf, blf, tlf);
	Quadrilateral3D<double> rear(tlb, blb, brb, trb);
	Quadrilateral3D<double> left(tlb, tlf, blf, blb);
	Quadrilateral3D<double> right(trf, trb, brb, brf);
	Hexahedron<double> shape(top, bottom, front, rear, left, right);

	double volume = shape.getVolume();

	BOOST_TEST(volume == 300.0);
}

// Test 2: Compute volume of more irregular hexahedron
// ToDo
BOOST_AUTO_TEST_CASE(getVolume_test2, * utf::tolerance(0.00001))
{

}

// === getCentroid ===
// Test 1: Compute correct centroid of a cube
BOOST_AUTO_TEST_CASE(getCentroid_test1, * utf::tolerance(0.00001))
{
	euc::EuclideanPoint<double,3> tlf(0.0, 5.0, 12.0);
	euc::EuclideanPoint<double,3> trf(5.0, 5.0, 12.0);
	euc::EuclideanPoint<double,3> blf(0.0, 5.0, 0.0);
	euc::EuclideanPoint<double,3> brf(5.0, 5.0, 0.0);
	euc::EuclideanPoint<double,3> tlb(0.0, 0.0, 12.0);
	euc::EuclideanPoint<double,3> trb(5.0, 0.0, 12.0);
	euc::EuclideanPoint<double,3> blb(0.0, 0.0, 0.0);
	euc::EuclideanPoint<double,3> brb(5.0, 0.0, 0.0);
	Quadrilateral3D<double> top(trf, tlf, tlb, trb);
	Quadrilateral3D<double> bottom(brf, brb, blb, blf);
	Quadrilateral3D<double> front(trf, brf, blf, tlf);
	Quadrilateral3D<double> rear(tlb, blb, brb, trb);
	Quadrilateral3D<double> left(tlb, tlf, blf, blb);
	Quadrilateral3D<double> right(trf, trb, brb, brf);
	Hexahedron<double> shape(top, bottom, front, rear, left, right);

	euc::EuclideanPoint<double,3> centroid = shape.getCentroid();
	BOOST_TEST(centroid.cmp[0], 2.5);
	BOOST_TEST(centroid.cmp[1], 2.5);
	BOOST_TEST(centroid.cmp[2], 6.0);
}

// === isPointOnEdge ===
// Test 1: Test point that is not on an edge
BOOST_AUTO_TEST_CASE(isPointOnEdge_test1, * utf::tolerance(0.00001))
{
	euc::EuclideanPoint<double,3> tlf(0.0, 5.0, 12.0);
	euc::EuclideanPoint<double,3> trf(5.0, 5.0, 12.0);
	euc::EuclideanPoint<double,3> blf(0.0, 5.0, 0.0);
	euc::EuclideanPoint<double,3> brf(5.0, 5.0, 0.0);
	euc::EuclideanPoint<double,3> tlb(0.0, 0.0, 12.0);
	euc::EuclideanPoint<double,3> trb(5.0, 0.0, 12.0);
	euc::EuclideanPoint<double,3> blb(0.0, 0.0, 0.0);
	euc::EuclideanPoint<double,3> brb(5.0, 0.0, 0.0);
	Quadrilateral3D<double> top(trf, tlf, tlb, trb);
	Quadrilateral3D<double> bottom(brf, brb, blb, blf);
	Quadrilateral3D<double> front(trf, brf, blf, tlf);
	Quadrilateral3D<double> rear(tlb, blb, brb, trb);
	Quadrilateral3D<double> left(tlb, tlf, blf, blb);
	Quadrilateral3D<double> right(trf, trb, brb, brf);
	Hexahedron<double> shape(top, bottom, front, rear, left, right);

	euc::EuclideanPoint<double,3> point(3.4, 2.3, 4.7);

	bool onEdge = shape.isPointOnEdge(point);

	BOOST_CHECK_EQUAL(onEdge, false);
}

// Test 2: Test point that is on an edge
BOOST_AUTO_TEST_CASE(isPointOnEdge_test2, * utf::tolerance(0.00001))
{
	euc::EuclideanPoint<double,3> tlf(0.0, 5.0, 12.0);
	euc::EuclideanPoint<double,3> trf(5.0, 5.0, 12.0);
	euc::EuclideanPoint<double,3> blf(0.0, 5.0, 0.0);
	euc::EuclideanPoint<double,3> brf(5.0, 5.0, 0.0);
	euc::EuclideanPoint<double,3> tlb(0.0, 0.0, 12.0);
	euc::EuclideanPoint<double,3> trb(5.0, 0.0, 12.0);
	euc::EuclideanPoint<double,3> blb(0.0, 0.0, 0.0);
	euc::EuclideanPoint<double,3> brb(5.0, 0.0, 0.0);
	Quadrilateral3D<double> top(trf, tlf, tlb, trb);
	Quadrilateral3D<double> bottom(brf, brb, blb, blf);
	Quadrilateral3D<double> front(trf, brf, blf, tlf);
	Quadrilateral3D<double> rear(tlb, blb, brb, trb);
	Quadrilateral3D<double> left(tlb, tlf, blf, blb);
	Quadrilateral3D<double> right(trf, trb, brb, brf);
	Hexahedron<double> shape(top, bottom, front, rear, left, right);

	euc::EuclideanPoint<double,3> point(2.13, 5.0, 12.0);

	bool onEdge = shape.isPointOnEdge(point);

	BOOST_CHECK_EQUAL(onEdge, true);
}

// === isPointOnVertex ===
// Test 1: Test point that is not on a vertex
BOOST_AUTO_TEST_CASE(isPointOnVertex_test1, * utf::tolerance(0.00001))
{
	euc::EuclideanPoint<double,3> tlf(0.0, 5.0, 12.0);
	euc::EuclideanPoint<double,3> trf(5.0, 5.0, 12.0);
	euc::EuclideanPoint<double,3> blf(0.0, 5.0, 0.0);
	euc::EuclideanPoint<double,3> brf(5.0, 5.0, 0.0);
	euc::EuclideanPoint<double,3> tlb(0.0, 0.0, 12.0);
	euc::EuclideanPoint<double,3> trb(5.0, 0.0, 12.0);
	euc::EuclideanPoint<double,3> blb(0.0, 0.0, 0.0);
	euc::EuclideanPoint<double,3> brb(5.0, 0.0, 0.0);
	Quadrilateral3D<double> top(trf, tlf, tlb, trb);
	Quadrilateral3D<double> bottom(brf, brb, blb, blf);
	Quadrilateral3D<double> front(trf, brf, blf, tlf);
	Quadrilateral3D<double> rear(tlb, blb, brb, trb);
	Quadrilateral3D<double> left(tlb, tlf, blf, blb);
	Quadrilateral3D<double> right(trf, trb, brb, brf);
	Hexahedron<double> shape(top, bottom, front, rear, left, right);

	euc::EuclideanPoint<double,3> point(2.3, 5.0, 12.0);

	bool onEdge = shape.isPointOnVertex(point);

	BOOST_CHECK_EQUAL(onEdge, false);
}

// Test 2: Test point that is on a vertex
BOOST_AUTO_TEST_CASE(isPointOnVertex_test2, * utf::tolerance(0.00001))
{
	euc::EuclideanPoint<double,3> tlf(0.0, 5.0, 12.0);
	euc::EuclideanPoint<double,3> trf(5.0, 5.0, 12.0);
	euc::EuclideanPoint<double,3> blf(0.0, 5.0, 0.0);
	euc::EuclideanPoint<double,3> brf(5.0, 5.0, 0.0);
	euc::EuclideanPoint<double,3> tlb(0.0, 0.0, 12.0);
	euc::EuclideanPoint<double,3> trb(5.0, 0.0, 12.0);
	euc::EuclideanPoint<double,3> blb(0.0, 0.0, 0.0);
	euc::EuclideanPoint<double,3> brb(5.0, 0.0, 0.0);
	Quadrilateral3D<double> top(trf, tlf, tlb, trb);
	Quadrilateral3D<double> bottom(brf, brb, blb, blf);
	Quadrilateral3D<double> front(trf, brf, blf, tlf);
	Quadrilateral3D<double> rear(tlb, blb, brb, trb);
	Quadrilateral3D<double> left(tlb, tlf, blf, blb);
	Quadrilateral3D<double> right(trf, trb, brb, brf);
	Hexahedron<double> shape(top, bottom, front, rear, left, right);

	euc::EuclideanPoint<double,3> point(0.0, 0.0, 12.0);

	bool onEdge = shape.isPointOnVertex(point);

	BOOST_CHECK_EQUAL(onEdge, true);
}

max_stars_repo_name: thorbenlouw/CUP-CFD

max_issues_repo_name: thorbenlouw/CUP-CFD

max_forks_repo_name: thorbenlouw/CUP-CFD


text:
/*    Copyright (c) 2010-2018, Delft University of Technology
 *    All rigths reserved
 *
 *    This file is part of the Tudat. Redistribution and use in source and
 *    binary forms, with or without modification, are permitted exclusively
 *    under the terms of the Modified BSD license. You should have received
 *    a copy of the license with this file. If not, please or visit:
 *    http://tudat.tudelft.nl/LICENSE.
 */

#define BOOST_TEST_MAIN

#include <string>
#include <thread>

#include <limits>

#include <boost/make_shared.hpp>
#include <boost/test/unit_test.hpp>

#include "Tudat/Basics/testMacros.h"
#include "Tudat/SimulationSetup/tudatEstimationHeader.h"
#include "Tudat/Astrodynamics/OrbitDetermination/EstimatableParameters/directTidalTimeLag.h"

namespace tudat
{
namespace unit_tests
{
BOOST_AUTO_TEST_SUITE( test_arcwise_environment )

//Using declarations.
using namespace tudat::observation_models;
using namespace tudat::orbit_determination;
using namespace tudat::estimatable_parameters;
using namespace tudat::interpolators;
using namespace tudat::numerical_integrators;
using namespace tudat::spice_interface;
using namespace tudat::simulation_setup;
using namespace tudat::orbital_element_conversions;
using namespace tudat::ephemerides;
using namespace tudat::propagators;
using namespace tudat::basic_astrodynamics;
using namespace tudat::coordinate_conversions;
using namespace tudat::ground_stations;
using namespace tudat::observation_models;

//! Unit test to check if tidal time lag parameters are estimated correctly
BOOST_AUTO_TEST_CASE( test_ArcwiseEnvironmentParameters )
{
    //Load spice kernels.
    spice_interface::loadStandardSpiceKernels( );

    // Define bodies in simulation
    std::vector< std::string > bodyNames;
    bodyNames.push_back( "Earth" );
    bodyNames.push_back( "Sun" );

    // Specify initial time
    double initialEphemerisTime = double( 1.0E7 );
    double finalEphemerisTime = initialEphemerisTime + 1.0 * 86400.0;

    // Create bodies needed in simulation
    std::map< std::string, std::shared_ptr< BodySettings > > bodySettings =
            getDefaultBodySettings( bodyNames );

    NamedBodyMap bodyMap = createBodies( bodySettings );
    bodyMap[ "Vehicle" ] = std::make_shared< Body >( );
    bodyMap[ "Vehicle" ]->setConstantBodyMass( 400.0 );

    // Create aerodynamic coefficient interface settings.
    double referenceArea = 4.0;
    double aerodynamicCoefficient = 1.2;
    std::shared_ptr< AerodynamicCoefficientSettings > aerodynamicCoefficientSettings =
            std::make_shared< ConstantAerodynamicCoefficientSettings >(
                referenceArea, aerodynamicCoefficient * ( Eigen::Vector3d( ) << 1.2, -0.01, 0.1 ).finished( ), 1, 1 );

    // Create and set aerodynamic coefficients object
    bodyMap[ "Vehicle" ]->setAerodynamicCoefficientInterface(
                createAerodynamicCoefficientInterface( aerodynamicCoefficientSettings, "Vehicle" ) );

    // Create radiation pressure settings
    double referenceAreaRadiation = 4.0;
    double radiationPressureCoefficient = 1.2;
    std::vector< std::string > occultingBodies;
    occultingBodies.push_back( "Earth" );
    std::shared_ptr< RadiationPressureInterfaceSettings > asterixRadiationPressureSettings =
            std::make_shared< CannonBallRadiationPressureInterfaceSettings >(
                "Sun", referenceAreaRadiation, radiationPressureCoefficient, occultingBodies );

    // Create and set radiation pressure settings
    bodyMap[ "Vehicle" ]->setRadiationPressureInterface(
                "Sun", createRadiationPressureInterface(
                    asterixRadiationPressureSettings, "Vehicle", bodyMap ) );

    bodyMap[ "Vehicle" ]->setEphemeris( std::make_shared< TabulatedCartesianEphemeris< > >(
                                            std::shared_ptr< interpolators::OneDimensionalInterpolator
                                            < double, Eigen::Vector6d > >( ), "Earth", "ECLIPJ2000" ) );

    setGlobalFrameBodyEphemerides( bodyMap, "SSB", "ECLIPJ2000" );

    ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
    ///////////////////////            CREATE ACCELERATIONS          //////////////////////////////////////////////////////
    ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

    // Set accelerations on Vehicle that are to be taken into account.
    SelectedAccelerationMap accelerationMap;
    std::map< std::string, std::vector< std::shared_ptr< AccelerationSettings > > > accelerationsOfVehicle;
    accelerationsOfVehicle[ "Earth" ].push_back( std::make_shared< AccelerationSettings >(
                                                     basic_astrodynamics::central_gravity ) );
    accelerationsOfVehicle[ "Sun" ].push_back( std::make_shared< AccelerationSettings >(
                                                   basic_astrodynamics::cannon_ball_radiation_pressure ) );
    accelerationsOfVehicle[ "Earth" ].push_back( std::make_shared< AccelerationSettings >(
                                                     basic_astrodynamics::aerodynamic ) );
    accelerationMap[ "Vehicle" ] = accelerationsOfVehicle;

    // Set bodies for which initial state is to be estimated and integrated.
    std::vector< std::string > bodiesToIntegrate;
    std::vector< std::string > centralBodies;
    bodiesToIntegrate.push_back( "Vehicle" );
    centralBodies.push_back( "Earth" );

    // Create acceleration models
    AccelerationMap accelerationModelMap = createAccelerationModelsMap(
                bodyMap, accelerationMap, bodiesToIntegrate, centralBodies );

    ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
    ///////////////////////             CREATE PROPAGATION SETTINGS            ////////////////////////////////////////////
    ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

    // Set Keplerian elements for Asterix.
    Eigen::Vector6d asterixInitialStateInKeplerianElements;
    asterixInitialStateInKeplerianElements( semiMajorAxisIndex ) = 7200.0E3;
    asterixInitialStateInKeplerianElements( eccentricityIndex ) = 0.05;
    asterixInitialStateInKeplerianElements( inclinationIndex ) = unit_conversions::convertDegreesToRadians( 85.3 );
    asterixInitialStateInKeplerianElements( argumentOfPeriapsisIndex ) = unit_conversions::convertDegreesToRadians( 235.7 );
    asterixInitialStateInKeplerianElements( longitudeOfAscendingNodeIndex ) = unit_conversions::convertDegreesToRadians( 23.4 );
    asterixInitialStateInKeplerianElements( trueAnomalyIndex ) = unit_conversions::convertDegreesToRadians( 139.87 );

    double earthGravitationalParameter = bodyMap.at( "Earth" )->getGravityFieldModel( )->getGravitationalParameter( );

    // Set (perturbed) initial state.
    Eigen::Matrix< double, 6, 1 > systemInitialState = convertKeplerianToCartesianElements(
                asterixInitialStateInKeplerianElements, earthGravitationalParameter );

    // Create propagator settings
    std::shared_ptr< DependentVariableSaveSettings > dependentVariableSaveSettings;
    std::vector< std::shared_ptr< SingleDependentVariableSaveSettings > > dependentVariables;
    dependentVariables.push_back(
                std::make_shared< SingleDependentVariableSaveSettings >(
                    aerodynamic_force_coefficients_dependent_variable, "Vehicle", "Earth" ) );
    dependentVariables.push_back(
                std::make_shared< SingleDependentVariableSaveSettings >(
                    radiation_pressure_coefficient_dependent_variable, "Vehicle", "Sun" ) );
    dependentVariableSaveSettings = std::make_shared< DependentVariableSaveSettings >( dependentVariables );

    std::shared_ptr< TranslationalStatePropagatorSettings< double > > propagatorSettings =
            std::make_shared< TranslationalStatePropagatorSettings< double > >(
                centralBodies, accelerationModelMap, bodiesToIntegrate, systemInitialState, double( finalEphemerisTime ),
                cowell, dependentVariableSaveSettings);

    // Create integrator settings
    std::shared_ptr< IntegratorSettings< double > > integratorSettings =
            std::make_shared< RungeKuttaVariableStepSizeSettingsScalarTolerances< double > >(
                double( initialEphemerisTime ), 90.0,
                RungeKuttaCoefficients::CoefficientSets::rungeKuttaFehlberg78,
                90.0, 90.0, 1.0, 1.0 );

    ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
    ///////////////////////             CREATE OBSERVATION SETTINGS            ////////////////////////////////////////////
    ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

    LinkEnds linkEnds;
    linkEnds[ observed_body ] = std::make_pair( "Vehicle", "" );
    observation_models::ObservationSettingsMap observationSettingsMap;
    observationSettingsMap.insert(
                std::make_pair( linkEnds, std::make_shared< ObservationSettings >( position_observable ) ) );

    ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
    ///////////////////////    DEFINE PARAMETERS THAT ARE TO BE ESTIMATED      ////////////////////////////////////////////
    ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

    // Define list of parameters to estimate.
    std::vector< std::shared_ptr< EstimatableParameterSettings > > parameterNames;
    parameterNames.push_back( std::make_shared< InitialTranslationalStateEstimatableParameterSettings< double > >(
                                  "Vehicle", systemInitialState, "Earth" ) );

    std::vector< double > arcStartTimeList =
    { initialEphemerisTime, initialEphemerisTime + 0.33 * 86400.0, initialEphemerisTime + 0.66 * 86400.0 };
    parameterNames.push_back(
                std::make_shared< ArcWiseDragCoefficientEstimatableParameterSettings >( "Vehicle", arcStartTimeList ) );
    parameterNames.push_back(
                std::make_shared< ArcWiseRadiationPressureCoefficientEstimatableParameterSettings >( "Vehicle", arcStartTimeList ) );

    // Create parameters
    std::shared_ptr< estimatable_parameters::EstimatableParameterSet< double > > parametersToEstimate =
            createParametersToEstimate( parameterNames, bodyMap );

    // Print identifiers and indices of parameters to terminal.
    printEstimatableParameterEntries( parametersToEstimate );

    Eigen::VectorXd truthParameters = parametersToEstimate->getFullParameterValues< double >( );
    truthParameters( 7 ) += 0.1;
    truthParameters( 8 ) += 0.2;

    truthParameters( 10 ) += 0.1;
    truthParameters( 11 ) += 0.2;

    parametersToEstimate->resetParameterValues( truthParameters );
    ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
    ///////////////////////          INITIALIZE ORBIT DETERMINATION OBJECT     ////////////////////////////////////////////
    ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

    // Create orbit determination object (propagate orbit, create observation models)
    OrbitDeterminationManager< double, double > orbitDeterminationManager =
            OrbitDeterminationManager< double, double >(
                bodyMap, parametersToEstimate, observationSettingsMap,
                integratorSettings, propagatorSettings );

    std::map< double, Eigen::VectorXd > dependentVariableData =
            orbitDeterminationManager.getVariationalEquationsSolver( )->getDynamicsSimulatorBase(
                )->getDependentVariableNumericalSolutionBase( )[ 0 ];

    // Test whether arc-wise coefficients are correctly used.
    double testDragCoefficient = 0.0;
    double testRadiationPressureCoefficient = 0.0;
    for( auto variableIterator : dependentVariableData )
    {
        double currentTime = variableIterator.first;

        if( currentTime < arcStartTimeList.at( 1 ) )
        {
            testDragCoefficient = truthParameters( 6 );
            testRadiationPressureCoefficient = truthParameters( 9 );
        }
        else if( currentTime < arcStartTimeList.at( 2 ) )
        {
            testDragCoefficient = truthParameters( 7 );
            testRadiationPressureCoefficient = truthParameters( 10 );
        }
        else
        {
            testDragCoefficient = truthParameters( 8 );
            testRadiationPressureCoefficient = truthParameters( 11 );
        }

        BOOST_CHECK_EQUAL( testDragCoefficient, variableIterator.second( 0 ) );
        BOOST_CHECK_EQUAL( aerodynamicCoefficient * -0.01, variableIterator.second( 1 ) );
        BOOST_CHECK_EQUAL( aerodynamicCoefficient * 0.1, variableIterator.second( 2 ) );
        BOOST_CHECK_EQUAL( testRadiationPressureCoefficient, variableIterator.second( 3 ) );

    }

    ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
    ///////////////////////          SIMULATE OBSERVATIONS                     ////////////////////////////////////////////
    ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

    // Define time of first observation
    double observationTimeStart = initialEphemerisTime + 1000.0;

    // Define time between two observations
    double  observationInterval = 30.0;

    // Simulate observations for 3 days
    std::vector< double > baseTimeList;
    double currentTime = initialEphemerisTime + 3600.0;
    while( currentTime < finalEphemerisTime - 3600.0 )
    {
        baseTimeList.push_back( currentTime );
        currentTime += observationInterval;
    }

    // Create measureement simulation input
    std::map< ObservableType, std::map< LinkEnds, std::pair< std::vector< double >, LinkEndType > > > measurementSimulationInput;
    measurementSimulationInput[ position_observable ][ linkEnds ] =
            std::make_pair( baseTimeList, observed_body );

    // Set typedefs for POD input (observation types, observation link ends, observation values, associated times with
    // reference link ends.
    typedef Eigen::Matrix< double, Eigen::Dynamic, 1 > ObservationVectorType;
    typedef std::map< LinkEnds, std::pair< ObservationVectorType, std::pair< std::vector< double >, LinkEndType > > >
            SingleObservablePodInputType;
    typedef std::map< ObservableType, SingleObservablePodInputType > PodInputDataType;

    // Simulate observations
    PodInputDataType observationsAndTimes = simulateObservations< double, double >(
                measurementSimulationInput, orbitDeterminationManager.getObservationSimulators( ) );

    ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
    //////////////////    PERTURB PARAMETER VECTOR AND ESTIMATE PARAMETERS     ////////////////////////////////////////////
    ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

    // Perturb parameter estimate
    Eigen::Matrix< double, Eigen::Dynamic, 1 > parameterPerturbation =
            Eigen::Matrix< double, Eigen::Dynamic, 1 >::Zero( truthParameters.rows( ) );
    parameterPerturbation.segment( 0, 3 ) = Eigen::Vector3d::Constant( 10.0 );
    parameterPerturbation.segment( 3, 3 ) = Eigen::Vector3d::Constant( 1.0E-2 );
    parameterPerturbation.segment( 6, 6 ) = Eigen::VectorXd::Constant( 6, 1.0 );
    Eigen::Matrix< double, Eigen::Dynamic, 1 > initialParameterEstimate = truthParameters;
    initialParameterEstimate += parameterPerturbation;

    // Define estimation input
    std::shared_ptr< PodInput< double, double > > podInput =
            std::make_shared< PodInput< double, double > >(
                observationsAndTimes, initialParameterEstimate.rows( ),
                Eigen::MatrixXd::Zero( truthParameters.rows( ), truthParameters.rows( ) ),
                initialParameterEstimate - truthParameters );
    podInput->defineEstimationSettings( true, true, false, true );

    // Perform estimation
    std::shared_ptr< PodOutput< double > > podOutput = orbitDeterminationManager.estimateParameters(
                podInput, std::make_shared< EstimationConvergenceChecker >( 4 ) );
    Eigen::VectorXd parameterEstimate = podOutput->parameterEstimate_ - truthParameters;

    for( int i = 0; i < 3; i++ )
    {
        BOOST_CHECK_SMALL( std::fabs( parameterEstimate( i ) ), 1.0E-2 );
        BOOST_CHECK_SMALL( std::fabs( parameterEstimate( i + 3 ) ), 1.0E-4 );
        BOOST_CHECK_SMALL( std::fabs( parameterEstimate( i + 6 ) ), 1.0E-5 );
        BOOST_CHECK_SMALL( std::fabs( parameterEstimate( i + 9 ) ), 1.0E-4 );

    }

}

BOOST_AUTO_TEST_SUITE_END( )

}

}

max_stars_repo_name: J-Westin/tudat

max_issues_repo_name: J-Westin/tudat

max_forks_repo_name: J-Westin/tudat