text:
/*
 * Copyright (c) 2015 Cryptonomex, Inc., and contributors.
 *
 * The MIT License
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 */

#include <boost/test/unit_test.hpp>

#include <scorum/protocol/protocol.hpp>

#include <limits>

#include "defines.hpp"

using namespace scorum::protocol;

BOOST_AUTO_TEST_SUITE(type_operation_tests)

BOOST_AUTO_TEST_CASE(check_asset_limits)
{
    BOOST_CHECK_GT(SCORUM_MAX_SHARE_SUPPLY, 0);

    BOOST_CHECK_LE(SCORUM_MAX_SHARE_SUPPLY, std::numeric_limits<share_value_type>::max());

    BOOST_CHECK_EQUAL(asset::maximum(SCORUM_SYMBOL).amount, SCORUM_MAX_SHARE_SUPPLY);

    BOOST_CHECK_EQUAL(asset::maximum(SP_SYMBOL).amount, SCORUM_MAX_SHARE_SUPPLY);

    BOOST_CHECK_GT(asset::maximum(SCORUM_SYMBOL), asset::min(SCORUM_SYMBOL));

    BOOST_CHECK_GT(asset::maximum(SP_SYMBOL), asset::min(SP_SYMBOL));
}

BOOST_AUTO_TEST_CASE(check_asset_value_transformations)
{
    auto ordinary_input_value = 1.5;

    BOOST_CHECK_EQUAL((share_value_type)ordinary_input_value, 1ll);

    ordinary_input_value = 1.5e+1;

    BOOST_CHECK_EQUAL((share_value_type)ordinary_input_value, 15ll);

    ordinary_input_value = 12.345e+5;

    BOOST_CHECK_EQUAL((share_value_type)ordinary_input_value, 1234500ll);

    ordinary_input_value = 12e+3;

    BOOST_CHECK_EQUAL((share_value_type)ordinary_input_value, 12000ll);

    auto max_input_value = 100000000e+9;

    BOOST_CHECK_EQUAL((share_value_type)max_input_value, SCORUM_MAX_SHARE_SUPPLY);
}

BOOST_AUTO_TEST_CASE(asset_symbol_test)
{
    const asset_symbol_type invald_symbol_name
        = (uint64_t(SCORUM_CURRENCY_PRECISION) | (uint64_t('E') << 8) | (uint64_t('R') << 16) | (uint64_t('R') << 24));
    const asset_symbol_type invald_symbol_precision
        = (uint64_t(SCORUM_CURRENCY_PRECISION / 2) | (uint64_t('S') << 8) | (uint64_t('P') << 16));

    BOOST_REQUIRE_NO_THROW(asset(10, SCORUM_SYMBOL));
    BOOST_REQUIRE_NO_THROW(asset(10, SP_SYMBOL));
    SCORUM_REQUIRE_THROW(asset(10, invald_symbol_name), fc::assert_exception);
    SCORUM_REQUIRE_THROW(asset(10, invald_symbol_precision), fc::assert_exception);
}

BOOST_AUTO_TEST_CASE(asset_operation_test)
{
    BOOST_CHECK_EQUAL(asset(2, SCORUM_SYMBOL) + 1, asset(3, SCORUM_SYMBOL));
    BOOST_CHECK_EQUAL(asset(3, SCORUM_SYMBOL) - 1, asset(2, SCORUM_SYMBOL));
    BOOST_CHECK_EQUAL(asset(2, SCORUM_SYMBOL) + asset(1, SCORUM_SYMBOL), asset(3, SCORUM_SYMBOL));
    BOOST_CHECK_EQUAL(asset(3, SCORUM_SYMBOL) - asset(1, SCORUM_SYMBOL), asset(2, SCORUM_SYMBOL));
    BOOST_CHECK_EQUAL(asset(2, SCORUM_SYMBOL) * 3, asset(6, SCORUM_SYMBOL));
    BOOST_CHECK_EQUAL(asset(6, SCORUM_SYMBOL) / 3, asset(2, SCORUM_SYMBOL));

    asset tmp(2, SCORUM_SYMBOL);
    tmp += 3;
    BOOST_REQUIRE_EQUAL(tmp, asset(5, SCORUM_SYMBOL));
    tmp -= 3;
    BOOST_REQUIRE_EQUAL(tmp, asset(2, SCORUM_SYMBOL));
    tmp = -tmp;
    BOOST_REQUIRE_EQUAL(tmp, asset(-2, SCORUM_SYMBOL));
    tmp *= 3;
    BOOST_REQUIRE_EQUAL(tmp, asset(-6, SCORUM_SYMBOL));
    tmp /= 3.f;
    BOOST_REQUIRE_EQUAL(tmp, asset(-2, SCORUM_SYMBOL));
    tmp = -tmp;

    BOOST_CHECK(tmp == tmp);

    asset tmp2(3, SCORUM_SYMBOL);

    BOOST_CHECK(tmp2 != tmp);

    BOOST_CHECK(tmp < tmp2);
    BOOST_CHECK(tmp <= tmp2);
    BOOST_CHECK(tmp2 > tmp);
    BOOST_CHECK(tmp2 >= tmp);

    BOOST_CHECK_THROW(asset(2, SCORUM_SYMBOL) + asset(1, SP_SYMBOL), fc::assert_exception);
    BOOST_CHECK_THROW(asset(3, SCORUM_SYMBOL) - asset(1, SP_SYMBOL), fc::assert_exception);
}

BOOST_AUTO_TEST_SUITE_END()

BOOST_AUTO_TEST_SUITE(basic_tests)

BOOST_AUTO_TEST_CASE(parse_size_test)
{
    BOOST_CHECK_THROW(fc::parse_size(""), fc::parse_error_exception);
    BOOST_CHECK_THROW(fc::parse_size("k"), fc::parse_error_exception);

    BOOST_CHECK_EQUAL(fc::parse_size("0"), (uint64_t)0);
    BOOST_CHECK_EQUAL(fc::parse_size("1"), (uint64_t)1);
    BOOST_CHECK_EQUAL(fc::parse_size("2"), (uint64_t)2);
    BOOST_CHECK_EQUAL(fc::parse_size("3"), (uint64_t)3);
    BOOST_CHECK_EQUAL(fc::parse_size("4"), (uint64_t)4);

    BOOST_CHECK_EQUAL(fc::parse_size("9"), (uint64_t)9);
    BOOST_CHECK_EQUAL(fc::parse_size("10"), (uint64_t)10);
    BOOST_CHECK_EQUAL(fc::parse_size("11"), (uint64_t)11);
    BOOST_CHECK_EQUAL(fc::parse_size("12"), (uint64_t)12);

    BOOST_CHECK_EQUAL(fc::parse_size("314159265"), (uint64_t)314159265);
    BOOST_CHECK_EQUAL(fc::parse_size("1k"), (uint64_t)1024);
    BOOST_CHECK_THROW(fc::parse_size("1a"), fc::parse_error_exception);
    BOOST_CHECK_EQUAL(fc::parse_size("1kb"), (uint64_t)1000);
    BOOST_CHECK_EQUAL(fc::parse_size("1MiB"), (uint64_t)1048576);
    BOOST_CHECK_EQUAL(fc::parse_size("32G"), (uint64_t)34359738368);
}

/**
 * Verify that names are RFC-1035 compliant https://tools.ietf.org/html/rfc1035
 * https://github.com/cryptonomex/graphene/issues/15
 */
BOOST_AUTO_TEST_CASE(valid_name_test)
{
    BOOST_CHECK(!is_valid_account_name("a"));
    BOOST_CHECK(!is_valid_account_name("A"));
    BOOST_CHECK(!is_valid_account_name("0"));
    BOOST_CHECK(!is_valid_account_name("."));
    BOOST_CHECK(!is_valid_account_name("-"));

    BOOST_CHECK(!is_valid_account_name("aa"));
    BOOST_CHECK(!is_valid_account_name("aA"));
    BOOST_CHECK(!is_valid_account_name("a0"));
    BOOST_CHECK(!is_valid_account_name("a."));
    BOOST_CHECK(!is_valid_account_name("a-"));

    BOOST_CHECK(is_valid_account_name("aaa"));
    BOOST_CHECK(!is_valid_account_name("aAa"));
    BOOST_CHECK(is_valid_account_name("a0a"));
    BOOST_CHECK(!is_valid_account_name("a.a"));
    BOOST_CHECK(is_valid_account_name("a-a"));

    BOOST_CHECK(is_valid_account_name("aa0"));
    BOOST_CHECK(!is_valid_account_name("aA0"));
    BOOST_CHECK(is_valid_account_name("a00"));
    BOOST_CHECK(!is_valid_account_name("a.0"));
    BOOST_CHECK(is_valid_account_name("a-0"));

    BOOST_CHECK(is_valid_account_name("aaa-bbb-ccc"));
    BOOST_CHECK(is_valid_account_name("aaa-bbb.ccc"));

    BOOST_CHECK(!is_valid_account_name("aaa,bbb-ccc"));
    BOOST_CHECK(!is_valid_account_name("aaa_bbb-ccc"));
    BOOST_CHECK(!is_valid_account_name("aaa-BBB-ccc"));

    BOOST_CHECK(!is_valid_account_name("1aaa-bbb"));
    BOOST_CHECK(!is_valid_account_name("-aaa-bbb-ccc"));
    BOOST_CHECK(!is_valid_account_name(".aaa-bbb-ccc"));
    BOOST_CHECK(!is_valid_account_name("/aaa-bbb-ccc"));

    BOOST_CHECK(!is_valid_account_name("aaa-bbb-ccc-"));
    BOOST_CHECK(!is_valid_account_name("aaa-bbb-ccc."));
    BOOST_CHECK(!is_valid_account_name("aaa-bbb-ccc.."));
    BOOST_CHECK(!is_valid_account_name("aaa-bbb-ccc/"));

    BOOST_CHECK(!is_valid_account_name("aaa..bbb-ccc"));
    BOOST_CHECK(is_valid_account_name("aaa.bbb-ccc"));
    BOOST_CHECK(is_valid_account_name("aaa.bbb.ccc"));

    BOOST_CHECK(is_valid_account_name("aaa--bbb--ccc"));
    BOOST_CHECK(!is_valid_account_name("xn--san-p8a.de"));
    BOOST_CHECK(is_valid_account_name("xn--san-p8a.dex"));
    BOOST_CHECK(!is_valid_account_name("xn-san-p8a.de"));
    BOOST_CHECK(is_valid_account_name("xn-san-p8a.dex"));

    BOOST_CHECK(is_valid_account_name("this-label-has"));
    BOOST_CHECK(!is_valid_account_name("this-label-has-more-than-63-char.act.ers-64-to-be-really-precise"));
    BOOST_CHECK(!is_valid_account_name("none.of.these.labels.has.more.than-63.chars--but.still.not.valid"));
}

BOOST_AUTO_TEST_CASE(merkle_root)
{
    signed_block block;
    std::vector<signed_transaction> tx;
    std::vector<digest_type> t;
    const uint32_t num_tx = 10;

    for (uint32_t i = 0; i < num_tx; i++)
    {
        tx.emplace_back();
        tx.back().ref_block_prefix = i;
        t.push_back(tx.back().merkle_digest());
    }

    auto c = [](const digest_type& digest) -> checksum_type { return checksum_type::hash(digest); };

    auto d = [](const digest_type& left, const digest_type& right) -> digest_type {
        return digest_type::hash(std::make_pair(left, right));
    };

    BOOST_CHECK(block.calculate_merkle_root() == checksum_type());

    block.transactions.push_back(tx[0]);
    BOOST_CHECK(block.calculate_merkle_root() == c(t[0]));

    digest_type dA, dB, dC, dD, dE, dI, dJ, dK, dM, dN, dO;

    /****************
     *              *
     *   A=d(0,1)   *
     *      / \     *
     *     0   1    *
     *              *
     ****************/

    dA = d(t[0], t[1]);

    block.transactions.push_back(tx[1]);
    BOOST_CHECK(block.calculate_merkle_root() == c(dA));

    /*************************
     *                       *
     *         I=d(A,B)      *
     *        /        \     *
     *   A=d(0,1)      B=2   *
     *      / \        /     *
     *     0   1      2      *
     *                       *
     *************************/

    dB = t[2];
    dI = d(dA, dB);

    block.transactions.push_back(tx[2]);
    BOOST_CHECK(block.calculate_merkle_root() == c(dI));

    /***************************
     *                         *
     *       I=d(A,B)          *
     *        /    \           *
     *   A=d(0,1)   B=d(2,3)   *
     *      / \    /   \       *
     *     0   1  2     3      *
     *                         *
     ***************************
     */

    dB = d(t[2], t[3]);
    dI = d(dA, dB);

    block.transactions.push_back(tx[3]);
    BOOST_CHECK(block.calculate_merkle_root() == c(dI));

    /***************************************
     *                                     *
     *                  __M=d(I,J)__       *
     *                 /            \      *
     *         I=d(A,B)              J=C   *
     *        /        \            /      *
     *   A=d(0,1)   B=d(2,3)      C=4      *
     *      / \        / \        /        *
     *     0   1      2   3      4         *
     *                                     *
     ***************************************/

    dC = t[4];
    dJ = dC;
    dM = d(dI, dJ);

    block.transactions.push_back(tx[4]);
    BOOST_CHECK(block.calculate_merkle_root() == c(dM));

    /**************************************
     *                                    *
     *                 __M=d(I,J)__       *
     *                /            \      *
     *        I=d(A,B)              J=C   *
     *       /        \            /      *
     *  A=d(0,1)   B=d(2,3)   C=d(4,5)    *
     *     / \        / \        / \      *
     *    0   1      2   3      4   5     *
     *                                    *
     **************************************/

    dC = d(t[4], t[5]);
    dJ = dC;
    dM = d(dI, dJ);

    block.transactions.push_back(tx[5]);
    BOOST_CHECK(block.calculate_merkle_root() == c(dM));

    /***********************************************
     *                                             *
     *                  __M=d(I,J)__               *
     *                 /            \              *
     *         I=d(A,B)              J=d(C,D)      *
     *        /        \            /        \     *
     *   A=d(0,1)   B=d(2,3)   C=d(4,5)      D=6   *
     *      / \        / \        / \        /     *
     *     0   1      2   3      4   5      6      *
     *                                             *
     ***********************************************/

    dD = t[6];
    dJ = d(dC, dD);
    dM = d(dI, dJ);

    block.transactions.push_back(tx[6]);
    BOOST_CHECK(block.calculate_merkle_root() == c(dM));

    /*************************************************
     *                                               *
     *                  __M=d(I,J)__                 *
     *                 /            \                *
     *         I=d(A,B)              J=d(C,D)        *
     *        /        \            /        \       *
     *   A=d(0,1)   B=d(2,3)   C=d(4,5)   D=d(6,7)   *
     *      / \        / \        / \        / \     *
     *     0   1      2   3      4   5      6   7    *
     *                                               *
     *************************************************/

    dD = d(t[6], t[7]);
    dJ = d(dC, dD);
    dM = d(dI, dJ);

    block.transactions.push_back(tx[7]);
    BOOST_CHECK(block.calculate_merkle_root() == c(dM));

    /************************************************************************
     *                                                                      *
     *                             _____________O=d(M,N)______________      *
     *                            /                                   \     *
     *                  __M=d(I,J)__                                  N=K   *
     *                 /            \                              /        *
     *         I=d(A,B)              J=d(C,D)                 K=E           *
     *        /        \            /        \            /                 *
     *   A=d(0,1)   B=d(2,3)   C=d(4,5)   D=d(6,7)      E=8                 *
     *      / \        / \        / \        / \        /                   *
     *     0   1      2   3      4   5      6   7      8                    *
     *                                                                      *
     ************************************************************************/

    dE = t[8];
    dK = dE;
    dN = dK;
    dO = d(dM, dN);

    block.transactions.push_back(tx[8]);
    BOOST_CHECK(block.calculate_merkle_root() == c(dO));

    /************************************************************************
     *                                                                      *
     *                             _____________O=d(M,N)______________      *
     *                            /                                   \     *
     *                  __M=d(I,J)__                                  N=K   *
     *                 /            \                              /        *
     *         I=d(A,B)              J=d(C,D)                 K=E           *
     *        /        \            /        \            /                 *
     *   A=d(0,1)   B=d(2,3)   C=d(4,5)   D=d(6,7)   E=d(8,9)               *
     *      / \        / \        / \        / \        / \                 *
     *     0   1      2   3      4   5      6   7      8   9                *
     *                                                                      *
     ************************************************************************/

    dE = d(t[8], t[9]);
    dK = dE;
    dN = dK;
    dO = d(dM, dN);

    block.transactions.push_back(tx[9]);
    BOOST_CHECK(block.calculate_merkle_root() == c(dO));
}

BOOST_AUTO_TEST_CASE(format_string)
{
    const fc::string etalon("'ABC' [101 = 101, 102 = 102]");
    fc::string fmt("'${first}' [101 = ${1}, 102 = ${2}]");
    fc::string out(fc::format_string(fmt, fc::mutable_variant_object()("first", "ABC")("1", 101)("2", 102)));

    BOOST_CHECK(!out.compare(etalon));
}

BOOST_AUTO_TEST_SUITE_END()

ext: cpp

max_stars_repo_name: scorum/scorum

max_issues_repo_name: Scorum/Scorum

max_forks_repo_name: Scorum/Scorum


text:
/*
 *  Phusion Passenger - https://www.phusionpassenger.com/
 *  Copyright (c) 2010-2017 Phusion Holding B.V.
 *
 *  "Passenger", "Phusion Passenger" and "Union Station" are registered
 *  trademarks of Phusion Holding B.V.
 *
 *  Permission is hereby granted, free of charge, to any person obtaining a copy
 *  of this software and associated documentation files (the "Software"), to deal
 *  in the Software without restriction, including without limitation the rights
 *  to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 *  copies of the Software, and to permit persons to whom the Software is
 *  furnished to do so, subject to the following conditions:
 *
 *  The above copyright notice and this permission notice shall be included in
 *  all copies or substantial portions of the Software.
 *
 *  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 *  IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 *  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 *  AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 *  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 *  OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 *  THE SOFTWARE.
 */

#include <sys/types.h>
#include <sys/wait.h>
#include <unistd.h>
#include <stdlib.h>
#include <signal.h>

#include <boost/thread.hpp>
#include <oxt/system_calls.hpp>
#include <string>
#include <cerrno>

#include <ProcessManagement/Spawn.h>
#include <ProcessManagement/Utils.h>
#include <StaticString.h>
#include <Exceptions.h>
#include <Utils/IOUtils.h>

namespace Passenger {

using namespace std;

int
runShellCommand(const StaticString &command) {
	string commandNt = command;
	const char *argv[] = {
		"/bin/sh",
		"-c",
		commandNt.c_str(),
		NULL
	};
	SubprocessInfo info;
	runCommand(argv, info);
	return info.status;
}

void
runCommand(const char **command, SubprocessInfo &info, bool wait, bool killSubprocessOnInterruption,
	const boost::function<void ()> &afterFork,
	const boost::function<void (const char **, int errcode)> &onExecFail)
{
	int e, waitStatus;
	pid_t waitRet;

	info.pid = asyncFork();
	if (info.pid == 0) {
		resetSignalHandlersAndMask();
		disableMallocDebugging();
		if (afterFork) {
			afterFork();
		}
		closeAllFileDescriptors(2, true);
		execvp(command[0], (char * const *) command);
		if (onExecFail) {
			onExecFail(command, errno);
		}
		_exit(1);
	} else if (info.pid == -1) {
		e = errno;
		throw SystemException("Cannot fork() a new process", e);
	} else if (wait) {
		try {
			waitRet = syscalls::waitpid(info.pid, &waitStatus, 0);
		} catch (const boost::thread_interrupted &) {
			if (killSubprocessOnInterruption) {
				boost::this_thread::disable_syscall_interruption dsi;
				syscalls::kill(SIGKILL, info.pid);
				syscalls::waitpid(info.pid, NULL, 0);
			}
			throw;
		}

		if (waitRet != -1) {
			info.status = waitStatus;
		} else if (errno == ECHILD || errno == ESRCH) {
			info.status = -2;
		} else {
			int e = errno;
			throw SystemException(string("Error waiting for the '") +
				command[0] + "' command", e);
		}
	}
}

void
runCommandAndCaptureOutput(const char **command, SubprocessInfo &info,
	string &output, bool killSubprocessOnInterruption,
	const boost::function<void ()> &afterFork,
	const boost::function<void (const char **command, int errcode)> &onExecFail)
{
	pid_t waitRet;
	int e, waitStatus;
	Pipe p;

	p = createPipe(__FILE__, __LINE__);

	info.pid = asyncFork();
	if (info.pid == 0) {
		dup2(p[1], 1);
		close(p[0]);
		close(p[1]);
		resetSignalHandlersAndMask();
		disableMallocDebugging();
		if (afterFork) {
			afterFork();
		}
		closeAllFileDescriptors(2, true);
		execvp(command[0], (char * const *) command);
		if (onExecFail) {
			onExecFail(command, errno);
		}
		_exit(1);
	} else if (info.pid == -1) {
		e = errno;
		throw SystemException("Cannot fork() a new process", e);
	} else {
		bool done = false;

		p[1].close();
		while (!done) {
			char buf[1024 * 4];
			ssize_t ret;

			try {
				ret = syscalls::read(p[0], buf, sizeof(buf));
			} catch (const boost::thread_interrupted &) {
				if (killSubprocessOnInterruption) {
					boost::this_thread::disable_syscall_interruption dsi;
					syscalls::kill(SIGKILL, info.pid);
					syscalls::waitpid(info.pid, NULL, 0);
				}
				throw;
			}
			if (ret == -1) {
				e = errno;
				if (killSubprocessOnInterruption) {
					boost::this_thread::disable_syscall_interruption dsi;
					syscalls::kill(SIGKILL, info.pid);
					syscalls::waitpid(info.pid, NULL, 0);
				}
				throw SystemException(string("Cannot read output from the '") +
					command[0] + "' command", e);
			}
			done = ret == 0;
			output.append(buf, ret);
		}
		p[0].close();

		try {
			waitRet = syscalls::waitpid(info.pid, &waitStatus, 0);
		} catch (const boost::thread_interrupted &) {
			if (killSubprocessOnInterruption) {
				boost::this_thread::disable_syscall_interruption dsi;
				syscalls::kill(SIGKILL, info.pid);
				syscalls::waitpid(info.pid, NULL, 0);
			}
			throw;
		}

		if (waitRet != -1) {
			info.status = waitStatus;
		} else if (errno == ECHILD || errno == ESRCH) {
			info.status = -2;
		} else {
			int e = errno;
			throw SystemException(string("Error waiting for the '") +
				command[0] + "' command", e);
		}
	}
}

} // namespace Passenger

max_stars_repo_name: mschnitzer/passenger

max_issues_repo_name: mschnitzer/passenger

max_forks_repo_name: mschnitzer/passenger


text:
// Copyright (c) 2022 midnightBITS
// This code is licensed under MIT license (see LICENSE for details)

#if defined(_MSC_VER)
#pragma warning(push)
#pragma warning(disable : 4619 4242)
#endif
#include <boost/asio/ip/host_name.hpp>
#if defined(_MSC_VER)
#pragma warning(pop)
#endif

#include <filesystem>
#include <http-router/server.hh>
#include <memory>
#include <server/http_session.hh>

using namespace std::literals;

namespace http_router::server {
	namespace {
		inline std::string host_name() {
			boost::system::error_code ec{};
			auto result = net::ip::host_name(ec);
			if (ec) {
				fail(ec, "host_name");
				result = "localhost"s;
			}
			return result;
		}

		class listener : public std::enable_shared_from_this<listener> {
			net::io_context& ioc_;
			tcp::acceptor acceptor_;
			std::string server_name_;
			router const* router_;

		public:
			listener(net::io_context& ioc,
			         tcp::endpoint const& endpoint,
			         std::string const& user_agent,
			         router const* handler,
			         reuse_address reuse_address_opt);

			// Start accepting incoming connections
			local_endpoint run() {
				auto const port = [this]() -> net::ip::port_type {
					try {
						return acceptor_.local_endpoint().port();
					} catch (boost::system::system_error const& sys) {
						fail(sys.code(), "local_endpoint");
						return 0;
					} catch (...) {
						return 0;
					}
				}();
				if (!port) return {};
				next_accept();
				return {
				    .interface = host_name(),
				    .port = port,
				};
			}

		private:
			void next_accept() {
				// The new connection gets its own strand
				acceptor_.async_accept(
				    net::make_strand(ioc_),
				    beast::bind_front_handler(&listener::on_accept,
				                              shared_from_this()));
			}

			void on_accept(beast::error_code ec, tcp::socket socket);
		};

		listener::listener(net::io_context& ioc,
		                   tcp::endpoint const& endpoint,
		                   std::string const& server_name,
		                   router const* handler,
		                   reuse_address reuse_address_opt)
		    : ioc_{ioc}
		    , acceptor_{net::make_strand(ioc)}
		    , server_name_{server_name}
		    , router_{handler} {
			beast::error_code ec;

			// Open the acceptor
			acceptor_.open(endpoint.protocol(), ec);
			if (ec) {
				fail(ec, "open");
				return;
			}

			// Allow address reuse
			acceptor_.set_option(
			    net::socket_base::reuse_address{reuse_address_opt.value()}, ec);
			if (ec) {
				fail(ec, "set_option");
				return;
			}

			// Bind to the server address
			acceptor_.bind(endpoint, ec);
			if (ec) {
				fail(ec, "bind");
				return;
			}

			// Start listening for connections
			acceptor_.listen(net::socket_base::max_listen_connections, ec);
			if (ec) {
				fail(ec, "listen");
				return;
			}
		}

		void listener::on_accept(beast::error_code ec, tcp::socket socket) {
			if (ec) {
				fail(ec, "accept");
			} else {
				// Create the http session and run it
				auto remote = socket.remote_endpoint(ec);
				if (ec) {
					fail(ec, "accept/remote_endpoint");
				} else {
					std::make_shared<http_session>(std::move(socket), remote,
					                               router_, &server_name_)
					    ->run();
				}
			}

			next_accept();
		}
	}  // namespace

	local_endpoint listen(net::io_context& ioc,
	                      tcp::endpoint const& endpoint,
	                      router const* handler,
	                      std::string const& server_name,
	                      reuse_address reuse_address_opt) {
		return std::make_shared<listener>(
		           ioc, endpoint,
		           !server_name.empty()
		               ? fmt::format("{} {}", server_name,
		                             BOOST_BEAST_VERSION_STRING)
		               : BOOST_BEAST_VERSION_STRING,
		           handler, reuse_address_opt)
		    ->run();
	}
}  // namespace http_router::server

ext: cc

max_stars_repo_name: mbits-os/http-router

max_issues_repo_name: mbits-os/http-router

max_forks_repo_name: mbits-os/http-router


text:
// -*- mode:c++;tab-width:2;indent-tabs-mode:t;show-trailing-whitespace:t;rm-trailing-spaces:t -*-
// vi: set ts=2 noet:
//
// (c) Copyright Rosetta Commons Member Institutions.
// (c) This file is part of the Rosetta software suite and is made available under license.
// (c) The Rosetta software is developed by the contributing members of the Rosetta Commons.
// (c) For more information, see http://www.rosettacommons.org. Questions about this can be
// (c) addressed to University of Washington UW TechTransfer, email: license@u.washington.edu.

#include <riflib/util.hh>
#include <riflib/scaffold/util.hh>
#include <riflib/rosetta_field.hh>
#include <riflib/scaffold/MultithreadPoseCloner.hh>
#include <numeric/agglomerative_hierarchical_clustering.hh>
#include <scheme/chemical/stub.hh>

#include <ObjexxFCL/format.hh>
#include <core/id/AtomID.hh>
#include <core/id/AtomID_Map.hh>
#include <core/chemical/ChemicalManager.hh>
#include <core/chemical/ResidueTypeSet.hh>
#include <core/conformation/Residue.hh>
#include <core/conformation/ResidueFactory.hh>
#include <core/conformation/util.hh>
#include <core/scoring/ScoreFunctionFactory.hh>
#include <core/scoring/ScoreFunction.hh>
#include <core/scoring/methods/EnergyMethodOptions.hh>
#include <core/scoring/Energies.hh>
#include <core/scoring/sasa.hh>
#include <core/scoring/hbonds/hbonds.hh>
#include <core/scoring/hbonds/HBondSet.hh>
#include <core/scoring/hbonds/HBondOptions.hh>
#include <utility/file/file_sys_util.hh>
#include <core/scoring/dssp/Dssp.hh>
#include <utility/io/izstream.hh>
#include <utility/io/ozstream.hh>

#include <core/scoring/rms_util.hh>
#include <core/select/util.hh>

#include <boost/functional/hash/hash.hpp>

namespace devel {
namespace scheme {

void
add_res_if_not_CGP(
	core::pose::Pose const & pose,
	int ir,
	utility::vector1<core::Size> & res,
	bool verbose = false
){
	if( pose.residue(ir).name3() == "CYD" ){ if(verbose) std::cout << "get_res skip CYD" << std::endl; return; }
	if( pose.residue(ir).name3() == "GLY" ){ if(verbose) std::cout << "get_res skip GLY" << std::endl; return; }
	if( pose.residue(ir).name3() == "PRO" ){ if(verbose) std::cout << "get_res skip PRO" << std::endl; return; }
	res.push_back( ir );
}

utility::vector1<core::Size>
get_res(
	std::string fname,
	core::pose::Pose const & pose,
	bool nocgp
 ){
 	using core::Size;
	utility::vector1<core::Size> res;
	if(fname==""){
		for(Size ir = 1; ir <= pose.size(); ++ir){
			if( nocgp )	add_res_if_not_CGP( pose, ir, res );
			else res.push_back(ir);
		}
		return res;
	}
	runtime_assert_msg( utility::file::file_exists( fname ), "res list file doesn't exist: " + fname );
	utility::io::izstream in(fname);
	std::string s;
	while(in>>s){
		try {
			utility::vector1<std::string> splt = utility::string_split(s,'-');
			// cout << splt << endl;
			if(splt.size()==1){
				// std::cout << "'" << splt.front() << "'" << std::endl;
				Size ir = boost::lexical_cast<Size>(splt.front());
				runtime_assert_msg( ir > 0 && ir <= pose.size(), "residue number out of bounds "+boost::lexical_cast<std::string>(ir) );
				if( nocgp )	add_res_if_not_CGP( pose, ir, res );
				else res.push_back(ir);
			}
			else if(splt.size()==2){
				Size lb = boost::lexical_cast<Size>(splt.front());
				Size ub = boost::lexical_cast<Size>(splt.back());
				if( ub < lb ) std::swap(lb,ub);
				for(Size ir = lb; ir <= ub; ++ir){
					runtime_assert_msg( ir > 0 && ir <= pose.size(), "residue number out of bounds "+boost::lexical_cast<std::string>(ir) );
				if( nocgp )	add_res_if_not_CGP( pose, ir, res );
				else res.push_back(ir);
				}
			} else {
				utility_exit_with_message("get_res: can't parse res line "+s);
			}
		} catch( boost::bad_lexical_cast ){
			utility_exit_with_message("get_res: can't parse res line "+s);
		}
	}
	in.close();

	std::set<core::Size> uniq;
	std::copy(res.begin(), res.end(), std::inserter(uniq, uniq.begin()));
	res.clear();
	std::copy(uniq.begin(), uniq.end(), std::back_inserter(res));
	std::sort( res.begin(), res.end() );
	if( nocgp ){
		for( int iri = 1; iri <= res.size(); ++iri ){
			int ir = res[iri];
			// std::cout << "selected res " << ir << " " << pose.residue(ir).name3() << std::endl;
			runtime_assert( pose.residue(ir).name3() != "CYD" );
			runtime_assert( pose.residue(ir).name3() != "GLY" );
			runtime_assert( pose.residue(ir).name3() != "PRO" );
		}
	}
	// utility_exit_with_message("oriestnd");

	return res;
 }

utility::vector1<core::Size> get_designable_positions_best_guess(
	  core::pose::Pose pose
	, bool noloops
	, bool nocgp
 ){
 	core::scoring::dssp::Dssp dssp( pose );
 	dssp.insert_ss_into_pose( pose );

 	core::pose::Pose allgly = pose;
 	pose_to_gly(allgly);

 	core::scoring::ScoreFunctionOP sf = core::scoring::get_score_function(true);
 	sf->score(allgly);

	core::scoring::methods::EnergyMethodOptions myopt = sf->energy_method_options();
	myopt.hbond_options().decompose_bb_hb_into_pair_energies(true);
	sf->set_energy_method_options(myopt);
 	sf->score(pose);
	core::scoring::hbonds::HBondSet hbset;
	core::scoring::hbonds::fill_hbond_set(pose,false,hbset);

	std::vector<float> sc2bb_energy(pose.size(),false);
	for(core::Size ihb = 1; ihb <= hbset.nhbonds(); ++ihb){
		core::scoring::hbonds::HBond const & hb(hbset.hbond(ihb));
		// noly donot design capping residues???
		if(  hb.don_hatm_is_protein_backbone() && !hb.acc_atm_is_protein_backbone() && pose.secstruct( hb.acc_res() ) == 'L' ) sc2bb_energy[hb.acc_res()-1] += hb.energy();
		if( !hb.don_hatm_is_protein_backbone() &&  hb.acc_atm_is_protein_backbone() && pose.secstruct( hb.don_res() ) == 'L' ) sc2bb_energy[hb.don_res()-1] += hb.energy();
	}

	utility::vector1<core::Size> res;
	// Real calc_per_atom_sasa( pose::Pose const & pose, id::AtomID_Map< Real > & atom_sasa, utility::vector1< Real > & rsd_sasa,
	// Real const probe_radius, bool const use_big_polar_H = false );
	core::id::AtomID_Map< core::Real > atom_sasa;
	utility::vector1< core::Real > rsd_sasa;
	core::scoring::calc_per_atom_sasa( pose, atom_sasa, rsd_sasa, 2.1 );
	for( int ir = 1; ir <= pose.size(); ++ir ){
		// std::cout << pose.secstruct(ir) << std::endl;
		bool isloop = pose.secstruct(ir) == 'L';
		int natoms = pose.residue(ir).nheavyatoms()-pose.residue(ir).last_backbone_atom();
		core::Real scsasa = 0;
		for( int ia = pose.residue(ir).first_sidechain_atom(); ia <= pose.residue(ir).natoms(); ++ia ){
			scsasa += atom_sasa[core::id::AtomID(ia,ir)];
		}
		float scsasa_per_atom = scsasa / float(natoms);

		bool is_exposed = scsasa_per_atom > 5.5;
		is_exposed |= scsasa > 18.0;
		// Longxing changed the cut off value below to -0.6, as most of the real capping Hbond score is lower than -1.0.
		bool has_bbsc = sc2bb_energy[ir-1] < -0.6;
		// if(pose.residue(ir).nchi()==4) has_bbsc = false;

		float datr =    pose.energies().residue_total_energies(ir)[core::scoring::fa_atr]
		             -allgly.energies().residue_total_energies(ir)[core::scoring::fa_atr];
		// float dhb  =    pose.energies().residue_total_energies(ir)[core::scoring::hbond_sc]
		//              -allgly.energies().residue_total_energies(ir)[core::scoring::hbond_sc];
		//       dhb +=    pose.energies().residue_total_energies(ir)[core::scoring::fa_elec]
		//              -allgly.energies().residue_total_energies(ir)[core::scoring::fa_elec];

		bool lowenergy = false;
		// lowenergy |= (( datr ) / float(natoms+1)) < -0.9;
		// lowenergy |= ( dhb  ) < -1.8;
		lowenergy |= has_bbsc;

		bool design_res = !lowenergy && is_exposed;

		if(noloops) design_res = design_res & !isloop;

		if( design_res ){
			if( nocgp )	add_res_if_not_CGP( pose, ir, res );
			else res.push_back(ir);
		} else {
			// std::cout << "pruned res by sasa " << ir << std::endl;
		}
	}
	runtime_assert_msg( res.size() > 0, "no residues selected!" );
	std::cout << "get_designable_positions_best_guess keep " << res.size() << " of " << pose.size() << std::endl;
	std::cout << "select designable_positions_best_guess=resi " << res[1];
	for(int i = 2; i <= res.size(); ++i ) std::cout << "+"<<res[i];
	std::cout << std::endl;

	// utility_exit_with_message("foo");

	return res;
}

std::string
get_res_list_hash( utility::vector1<core::Size> const & reslist ){
	assert( reslist.size() > 0 );
	std::size_t reslist_hash = boost::hash_range( reslist.begin(), reslist.end() );
	std::string hashstr = boost::lexical_cast<std::string>(reslist_hash);
	return hashstr.substr(hashstr.size()-8,8);

}

std::string KMGT(double const & x, int const & w, int const & d){
	using ObjexxFCL::format::F;
	if( x < 1e3  ) return F( w, d, x/1e0  )+" ";
	if( x < 1e6  ) return F( w, d, x/1e3  )+"K";
	if( x < 1e9  ) return F( w, d, x/1e6  )+"M";
	if( x < 1e12 ) return F( w, d, x/1e9  )+"G";
	if( x < 1e15 ) return F( w, d, x/1e12 )+"T";
	if( x < 1e18 ) return F( w, d, x/1e15 )+"P";
	if( x < 1e21 ) return F( w, d, x/1e18 )+"E";
	if( x < 1e24 ) return F( w, d, x/1e21 )+"Z";
	else           return F( w, d, x/1e24 )+"Y";

}

void pose_to_ala( core::pose::Pose & pose ){
	core::chemical::ResidueTypeSetCAP rts = core::chemical::ChemicalManager::get_instance()->residue_type_set("fa_standard");
	core::conformation::ResidueOP ala = core::conformation::ResidueFactory::create_residue( rts.lock()->name_map("ALA") );
	for( int ir = 1; ir <= pose.size(); ++ir ){
		if( ! pose.residue(ir).is_protein()   ) continue;
		if(   pose.residue(ir).name3()=="GLY" ) continue;
		if(   pose.residue(ir).name3()=="PRO" ) continue;
		if(   pose.residue(ir).name3()=="CYD" ) continue;
		pose.replace_residue( ir, *ala, true );
	}
}
void pose_to_gly( core::pose::Pose & pose ){
	core::chemical::ResidueTypeSetCAP rts = core::chemical::ChemicalManager::get_instance()->residue_type_set("fa_standard");
	core::conformation::ResidueOP gly = core::conformation::ResidueFactory::create_residue( rts.lock()->name_map("GLY") );
	for( int ir = 1; ir <= pose.size(); ++ir ){
		if( ! pose.residue(ir).is_protein()   ) continue;
		if(   pose.residue(ir).name3()=="GLY" ) continue;
		if(   pose.residue(ir).name3()=="PRO" ) continue;
		if(   pose.residue(ir).name3()=="CYD" ) continue;
		pose.replace_residue( ir, *gly, true );
	}
}
void pose_to_ala( core::pose::Pose & pose, utility::vector1<core::Size> const & res_sel ){
	core::chemical::ResidueTypeSetCAP rts = core::chemical::ChemicalManager::get_instance()->residue_type_set("fa_standard");
	core::conformation::ResidueOP ala = core::conformation::ResidueFactory::create_residue( rts.lock()->name_map("ALA") );
	for( int iri = 1; iri <= res_sel.size(); ++iri ){
		int ir = res_sel[iri];
		if( ! pose.residue(ir).is_protein()   ) continue;
		if(   pose.residue(ir).name3()=="GLY" ) continue;
		if(   pose.residue(ir).name3()=="PRO" ) continue;
		if(   pose.residue(ir).name3()=="CYD" ) continue;
		pose.replace_residue( ir, *ala, true );
	}
}

std::string
open_for_read_on_path(
	std::vector<std::string> const & path,
	std::string fname,
	utility::io::izstream & in
){
	for( auto const & dir : path ){
		if( utility::file::file_exists( dir+"/"+fname ) ){
			in.close();
			in.clear();
			in.open( dir+"/"+fname, std::ios::binary );
			if( in.good() ){
				// std::cout << "open_for_read_on_path: opened " << dir + "/" + fname << std::endl;
				return dir+"/"+fname;
			}
		}
	}
	return std::string();
}

std::string
open_for_write_on_path(
	std::vector<std::string> const & path,
	std::string fname,
	utility::io::ozstream & out,
	bool create_directorys /* = false */
){
	for( auto const & dir: path ){
		if( create_directorys && !utility::file::file_exists( dir ) ){
			utility::file::create_directory_recursive( dir );
		}
		if( utility::file::file_exists( dir ) ){
			out.close();
			out.clear();
			out.open( dir+"/"+fname, std::ios::binary );
			if( out.good() ){
				// std::cout << "open_for_write_on_path: opened " << dir + "/" + fname << std::endl;
				return dir+"/"+fname;
			}
		}
	}
	return std::string();

}

void
append_pose_to_pose(
	core::pose::Pose & pose1,
	core::pose::Pose const & pose2,
	bool new_chain
){
	append_subpose_to_pose(pose1, pose2, 1, pose2.size(), new_chain);
}

void
append_subpose_to_pose(
	core::pose::Pose & pose1,
	core::pose::Pose const & pose2,
	core::Size start_res,
	core::Size end_res,
	bool new_chain
){
	if ( pose2.size()<start_res ) {
		std::cerr << "Provided starting residue number " << start_res
			<< " less than number residues in appended pose. Nothing to do." << std::endl;
	}
	pose1.append_residue_by_jump(pose2.residue(start_res), pose1.size() , "", "", new_chain);
	for ( core::Size i=start_res+1; i<=end_res; ++i ) {
		if ( pose2.residue(i).is_lower_terminus() || !pose2.residue(i).is_protein() ) {
			if ( i > 1 && pose2.chain(i) == pose2.chain(i-1) ) {
				pose1.append_residue_by_jump(pose2.residue(i), pose1.size(), "","", false);
			} else {
				pose1.append_residue_by_jump(pose2.residue(i), pose1.size(), "","", true);
			}
		} else {
			pose1.append_residue_by_bond(pose2.residue(i));
		}
	}
}

std::vector<int>
get_rif_atype_map()
{
	core::chemical::AtomTypeSetCOP ats = core::chemical::ChemicalManager::get_instance()->atom_type_set("fa_standard");
	// 147 CH0
	// 148 HS
	// 149 NtrR
	// 150 SH1

	std::vector<int> atypemap;
	atypemap.resize( ats->n_atomtypes()+1, 12345 );
	atypemap[ ats->atom_type_index("CNH2") ] =  1;
	atypemap[ ats->atom_type_index("COO" ) ] =  2;
	atypemap[ ats->atom_type_index("CH1" ) ] =  3;
	atypemap[ ats->atom_type_index("CH2" ) ] =  4;
	atypemap[ ats->atom_type_index("CH3" ) ] =  5;
	atypemap[ ats->atom_type_index("aroC") ] =  6;
	atypemap[ ats->atom_type_index("Ntrp") ] =  7;
	atypemap[ ats->atom_type_index("Nhis") ] =  8;
	atypemap[ ats->atom_type_index("NH2O") ] =  9;
	atypemap[ ats->atom_type_index("Nlys") ] = 10;
	atypemap[ ats->atom_type_index("Narg") ] = 11;
	atypemap[ ats->atom_type_index("Npro") ] = 12;
	atypemap[ ats->atom_type_index("OH"  ) ] = 13;
	atypemap[ ats->atom_type_index("ONH2") ] = 14;
	atypemap[ ats->atom_type_index("OOC" ) ] = 15;
	atypemap[ ats->atom_type_index("Oaro") ] = 16;
	atypemap[ ats->atom_type_index("S"   ) ] = 17;
	atypemap[ ats->atom_type_index("Nbb" ) ] = 18;
	atypemap[ ats->atom_type_index("CAbb") ] = 19;
	atypemap[ ats->atom_type_index("CObb") ] = 20;
	atypemap[ ats->atom_type_index("OCbb") ] = 21;

	// -beta stuff
	if(ats->has_atom("CH0" )) atypemap[ ats->atom_type_index("CH0" ) ] = 6;
	// if(ats->has_atom("HS")) atypemap[ ats->atom_type_index("HS"  ) ] = ??;
	if(ats->has_atom("NtrR")) atypemap[ ats->atom_type_index("NtrR") ] = 11;
	if(ats->has_atom("SH1" )) atypemap[ ats->atom_type_index("SH1" ) ] = 17;

	return atypemap;
}

std::map< core::id::AtomID, core::id::AtomID >
residue_subset_to_CA_atom_map(
    core::select::residue_selector::ResidueSubset const & subset,
    core::pose::Pose const & pose ) {

    std::map< core::id::AtomID, core::id::AtomID > atom_map;
    utility::vector1<core::Size> seq_poss = core::select::get_residues_from_subset( subset );
    for ( core::Size seq_pos : seq_poss ) {
        core::conformation::Residue const & res( pose.residue( seq_pos ) );
        if ( res.has( "CA" ) ) {
            core::id::AtomID const id( res.atom_index( "CA" ), seq_pos );
            atom_map[ id ] = id;
        }
    }
    return atom_map;
}

core::Real
subset_CA_rmsd(
    core::pose::Pose const & pose1,
    core::pose::Pose const & pose2,
    core::select::residue_selector::ResidueSubset const & subset,
    bool superimpose) {

    utility::vector1<core::Size> calc_rms_res = core::select::get_residues_from_subset( subset );
    std::map< core::id::AtomID, core::id::AtomID > atom_map =
        residue_subset_to_CA_atom_map( subset, pose1 );

    if ( superimpose ) {
        return core::scoring::rms_at_corresponding_atoms( pose1, pose2, atom_map, calc_rms_res );
    } else {
        return core::scoring::rms_at_corresponding_atoms_no_super( pose1, pose2, atom_map, calc_rms_res );
    }
}

/// Your poses must be identical. Very important
EigenXform
find_xform_from_identical_pose_to_pose(
	core::pose::Pose const & pose1,
	core::pose::Pose const & pose2,
	float align_error /* = 0.2 */ ) {

	runtime_assert( pose1.size() == pose2.size() );

	// we can use anything here because they're the same
	utility::vector1<core::Size> scaffold_res;
	get_default_scaffold_res( pose1, scaffold_res );

	core::pose::Pose to_move = pose1;
    ::devel::scheme::pose_to_ala( to_move );
    Eigen::Vector3f to_move_center = pose_center(to_move, scaffold_res);

    core::pose::Pose match_this = pose2;
    ::devel::scheme::pose_to_ala( match_this );
    Eigen::Vector3f match_center = pose_center(match_this,scaffold_res);

   	// prepare all test point now

    // Test original residue 1
    utility::vector1<core::Size> target_res {1};
    std::vector< ::scheme::actor::Atom< Eigen::Vector3f > > to_move_atoms;
    std::vector< ::scheme::actor::Atom< Eigen::Vector3f > > match_atoms;

    devel::scheme::get_scheme_atoms( to_move, target_res, to_move_atoms, true );
    devel::scheme::get_scheme_atoms( match_this, target_res, match_atoms, true );

    // Test original midpoint
    core::Size test_res = to_move.size() / 2;

    utility::vector1<core::Size> test_target_res {test_res};
    std::vector< ::scheme::actor::Atom< Eigen::Vector3f > > test_to_move_atoms;
    std::vector< ::scheme::actor::Atom< Eigen::Vector3f > > test_match_atoms;

    devel::scheme::get_scheme_atoms( to_move, test_target_res, test_to_move_atoms, true );
    devel::scheme::get_scheme_atoms( match_this, test_target_res, test_match_atoms, true );

    // Test centered residue 1

    EigenXform to_move_2_center = EigenXform::Identity();
    to_move_2_center.translation() = -to_move_center;
    apply_xform_to_pose( to_move, to_move_2_center );

    std::vector< ::scheme::actor::Atom< Eigen::Vector3f > > to_move_centered_atoms;
    devel::scheme::get_scheme_atoms( to_move, target_res, to_move_centered_atoms, true );

    // done making test points

    EigenXform match_x = ::scheme::chemical::make_stub<EigenXform>(
                                                                match_atoms[0].position(),
                                                                match_atoms[1].position(),
                                                                match_atoms[2].position());
    EigenXform to_move_x = ::scheme::chemical::make_stub<EigenXform>(
                                                                to_move_centered_atoms[0].position(),
                                                                to_move_centered_atoms[1].position(),
                                                                to_move_centered_atoms[2].position());

    EigenXform to_move_centered_2_match = match_x * to_move_x.inverse();
    to_move_centered_2_match.translation() = match_center;

    double centered_error = (to_move_centered_2_match * to_move_centered_atoms[0].position() - match_atoms[0].position()).norm();
	std::cout << "Centered res1 Alignment error :" << centered_error << std::endl;
    runtime_assert( centered_error < align_error );

    EigenXform to_move_2_match = to_move_centered_2_match * to_move_2_center;

    double error = (to_move_2_match * to_move_atoms[0].position() - match_atoms[0].position()).norm();
	std::cout << "res1 Alignment error :" << error << std::endl;
    runtime_assert( error < align_error );

    double test_error = (to_move_2_match * test_to_move_atoms[0].position() - test_match_atoms[0].position()).norm();
    std::cout << "Test res" << test_res << " Alignment error :" << error << std::endl;
    runtime_assert( test_error < align_error );

    return to_move_2_match;

}

void
apply_xform_to_pose( core::pose::Pose & pose, EigenXform const & xform) {
	numeric::xyzTransform<float> transform = eigen2xyz( xform );
	pose.apply_transform_Rx_plus_v( transform.R, transform.t );

}

void
apply_xform_to_residue( core::conformation::Residue & residue, EigenXform const & xform) {
	numeric::xyzTransform<float> transform = eigen2xyz( xform );
	residue.apply_transform_Rx_plus_v( transform.R, transform.t );

}

void
all_by_all_rmsd(
	std::vector<core::pose::PoseOP> const & poses,
	utility::vector1<utility::vector1<core::Real>> & table ) {

	table.clear();
	table.resize( poses.size() );

	for ( core::Size i = 1; i <= table.size(); i++ ) {
		table[i].resize(poses.size(), 0);	// initialize the diagonal to 0
	}

	// Make it threadsafe
	utility::vector1<shared_ptr<MultithreadPoseCloner>> mpcs;
	for ( core::pose::PoseOP const & pose : poses ) {
		mpcs.push_back(make_shared<MultithreadPoseCloner>(pose));
	}

	std::exception_ptr exception = nullptr;
    #ifdef USE_OPENMP
    #pragma omp parallel for schedule(dynamic,1)
    #endif
	for ( core::Size i = 1; i <= mpcs.size(); i++ ) {
		if (exception) continue;
		try {
			core::pose::PoseCOP outer_pose = mpcs[i]->get_pose();

			for ( core::Size j = i + 1; j <= mpcs.size(); j++) {
				core::pose::PoseCOP inner_pose = mpcs[j]->get_pose();

				core::Real rmsd = core::scoring::CA_rmsd(*outer_pose, *inner_pose);

				runtime_assert( table[i][j] == 0 );
				runtime_assert( table[j][i] == 0 );

				table[i][j] = rmsd;
				table[j][i] = rmsd;
			}

		} catch(...) {
            #pragma omp critical
            exception = std::current_exception();
        }

	} // end of OMP loop
    if( exception ) std::rethrow_exception(exception);

}

std::vector<std::vector<std::pair<core::pose::PoseOP, uint64_t>>>
cluster_poses_into_n_bins(
	std::vector<core::pose::PoseOP> const & poses,
	uint64_t n,
	utility::vector1<utility::vector1<core::Real>> & rmsds ) {

	runtime_assert( poses.size() >= n );

	if (rmsds.size() != poses.size()) {
		std::cout << "Calculating n^2 rmsd table" << std::endl;
		all_by_all_rmsd( poses, rmsds );
	}

	numeric::AverageLinkClusterer alc;
	utility::vector1<numeric::ClusteringTreeNodeOP> roots = alc.cluster(rmsds, n);

	utility::vector1<uint64_t> input_indices( poses.size() );
	for ( core::Size i = 1; i <= poses.size(); i++ ) {
		input_indices[i] = i - 1;
	}

	std::vector<bool> pose_got_used(poses.size(), false);

	std::vector<std::vector<std::pair<core::pose::PoseOP, uint64_t>>> bins;

	for ( core::Size i = 1; i <= roots.size(); i++ ) {
		std::vector<std::pair<core::pose::PoseOP, uint64_t>> this_bin;

		utility::vector1<uint64_t> this_bin_indices;
		numeric::get_cluster_data( input_indices, roots[i], this_bin_indices );

		runtime_assert( this_bin_indices.size() > 0 );

		for ( uint64_t index : this_bin_indices ) {
			runtime_assert( ! pose_got_used[index] );
			pose_got_used[index] = true;
			this_bin.push_back( std::pair<core::pose::PoseOP, uint64_t>(poses[index], index) );
		}

		bins.push_back( this_bin );
	}

	for ( uint64_t i = 0; i < pose_got_used.size(); i++ ) {
		runtime_assert(pose_got_used[i]);
	}

	return bins;
}

std::vector<core::pose::PoseOP>
cluster_poses_leaving_n(
	std::vector<core::pose::PoseOP> const & poses,
	uint64_t n ) {

	if ( n >= poses.size() ) return poses;

	utility::vector1<utility::vector1<core::Real>> rmsds;
	std::vector<std::vector<std::pair<core::pose::PoseOP, uint64_t>>> bins = cluster_poses_into_n_bins( poses, n, rmsds );

	std::vector<core::pose::PoseOP> output_poses;

	for (std::vector<std::pair<core::pose::PoseOP, uint64_t>> const & bin : bins) {
		output_poses.push_back(bin.front().first);
	}

	return output_poses;
}

// this returns the index into the indexes list
// indexes should be 0 based even though rmsds is 1 based
size_t
find_cluster_center(
	std::vector<uint64_t> indexes,
	utility::vector1<utility::vector1<core::Real>> const & rmsds ) {

	std::vector<core::Real> rmsd_sums( indexes.size() );

	for ( size_t i = 0; i < indexes.size(); i++ ) {
		uint64_t index = indexes[i];

		core::Real this_sum = 0;
		for ( uint64_t ind : indexes ) {
			this_sum += rmsds.at(ind+1).at(index+1);
		}
		// std::cout << this_sum << std::endl;
		rmsd_sums[i] = this_sum;
	}
	ptrdiff_t center_index = std::min_element(rmsd_sums.begin(), rmsd_sums.end()) - rmsd_sums.begin();

	std::cout << "Cluster center is " << center_index << " with rmsd sum " << rmsd_sums[center_index] << std::endl;

	return center_index;
}

// uses a binary search to try to find a number
// of bins to cluster poses into such than the top n
// represent frac of poses within tol
std::vector<core::pose::PoseOP>
cluster_poses_leaving_n_representing_frac(
	std::vector<core::pose::PoseOP> const & poses,
	uint64_t n,
	float frac,
	float tol
	) {

	using ObjexxFCL::format::I;

	runtime_assert( frac > 0 && frac <= 1);

	uint64_t trial = n * (uint64_t)( 1.0 / frac );
	uint64_t initial_trial = trial;
	std::vector<uint64_t> trial_history;

	std::vector<size_t> idx ;
	utility::vector1<utility::vector1<core::Real>> rmsds;

	std::vector<std::vector<std::pair<core::pose::PoseOP, uint64_t>>> bins;

	while (true) {

		trial_history.push_back( trial );
		std::cout << "Round " << trial_history.size() << " : trial size " << trial << std::endl;

		bins = cluster_poses_into_n_bins( poses, trial, rmsds );

// sort indexes based on size of bin
		idx.clear();
		idx.resize(bins.size());
  		std::iota(idx.begin(), idx.end(), 0);
		std::sort(idx.begin(), idx.end(),
       		[&bins](size_t i1, size_t i2) {return bins[i1].size() < bins[i2].size();});

		size_t biggest_bin = bins[idx.back()].size();

// display block //////////////////////////
		{
			// find the number of digits
			int digits = 0;
			{
				size_t temp = biggest_bin;
				while (temp != 0) { temp /= 10; digits++; }
			}

			size_t star_every = std::max<size_t>(biggest_bin / 50, 1);

			for ( size_t i = 0; i < bins.size(); i++ ) {
				if ( i == n ) std::cout << "============================ cut =========================" << std::endl;
				// going in reverse order and using indexes
				size_t bin = idx[ bins.size() - i - 1 ];

				std::cout << I( digits + 1, bins[bin].size() ) << " ";

				{
					int temp = bins[bin].size() - star_every;
					while ( temp > 0 ) { std::cout << "*"; temp -= star_every; }
				}
				std::cout << std::endl;
			}
		}

/////////////////////////////////////////////

// decide if we have met the tolerance criteria
		size_t current_count = 0;
		for ( size_t i = 0; i < n; i++ ) {
			size_t bin = idx[ bins.size() - i - 1 ];
			current_count += bins[bin].size();
		}

		float current_frac = (float)current_count / (float)poses.size();
		float error = current_frac - frac;
		std::cout << "Error = " << error << std::endl;
		if ( std::abs(error) < tol ) {
			std::cout << "Tolerance satisfied, clustering complete" << std::endl;
			break;
		}

// find position in array
		std::sort(trial_history.begin(), trial_history.end());
		ptrdiff_t pos = std::find(trial_history.begin(), trial_history.end(), trial ) - trial_history.begin();

// either one plus or one minus
		int other_relevant_pos;
		if ( error > 0 ) {
			other_relevant_pos = pos + 1;
		} else {
			other_relevant_pos = pos - 1;
		}

// naively set the next trial size
		if ( other_relevant_pos < 0 ) {
			trial /= 2;
		} else if ( other_relevant_pos >= trial_history.size() ) {
			trial *= 4;
		} else {
			uint64_t other_trial = trial_history[other_relevant_pos];
			trial = ( other_trial + trial ) / 2;
		}

// cap it to meaningful values
		trial = std::max<uint64_t>( 1, trial );
		trial = std::min<uint64_t>( poses.size(), trial );

// make sure next trial will be unique
		if ( std::find(trial_history.begin(), trial_history.end(), trial ) != trial_history.end() ) {
			std::cout << "Repeating trial, clustering complete" << std::endl;
			break;
		}
	}

	std::vector<core::pose::PoseOP> output_poses;

	for ( size_t i = 0; i < n; i++ ) {
		size_t bin = idx.at( bins.size() - i - 1 );
		std::vector<uint64_t> indexes;
		for ( std::pair<core::pose::PoseOP, uint64_t> pair : bins[bin] ) {
			indexes.push_back( pair.second );
		}
		size_t cluster_center = find_cluster_center( indexes, rmsds );
		output_poses.push_back( bins.at(bin).at(cluster_center).first );
	}

	return output_poses;
}

std::vector<core::pose::PoseOP>
random_selection_poses_leaving_n(
	std::vector<core::pose::PoseOP> const & poses,
	uint64_t n ) {

	if ( n >= poses.size() ) return poses;

	std::vector<uint64_t> indexes( poses.size());
	for ( uint64_t i = 0; i < indexes.size(); i++ ) {
		indexes[i] = i;
	}

	std::random_shuffle(indexes.begin(), indexes.end());

	std::vector<core::pose::PoseOP> output_poses(n);

	for ( uint64_t i = 0; i < n; i++ ) {
		output_poses[i] = poses[indexes[i]];
	}

	return output_poses;
}

Eigen::Vector3f
pose_center(
    core::pose::Pose const & pose,
    utility::vector1<core::Size> const & useres /* = utility::vector1<core::Size>()*/
){
    typedef numeric::xyzVector<core::Real> Vec;
    Vec cen(0,0,0);
    int count = 0;
    for( int ir = 1; ir <= pose.size(); ++ir ) {
        if( useres.size()==0 || std::find(useres.begin(),useres.end(),ir)!=useres.end() ){
            for( int ia = 1; ia <= pose.residue(ir).nheavyatoms(); ++ia ){
                cen += pose.xyz(core::id::AtomID(ia,ir));
                ++count;
            }
        // } else {
            // std::cout << "pose_center skip " << ir << std::endl;
        }
    }
    cen /= double(count);
    // ::scheme::util::SimpleArray<3,float> center;
    Eigen::Vector3f center;
    center[0] = cen[0];
    center[1] = cen[1];
    center[2] = cen[2];
    return center;
}

void
get_rg_radius(
    core::pose::Pose const & pose,
    float & rg,
    float & radius,
    utility::vector1<core::Size> const & useres /*= utility::vector1<core::Size>()*/,
    bool allatom /*= false*/
){
    Eigen::Vector3f centmp = pose_center( pose, useres );
    numeric::xyzVector<double> cen;
    float maxdis = -9e9, avgdis2 = 0.0;
    for( int i = 0; i < 3; ++i ) cen[i] = centmp[i];
    for( int iri = 1; iri <= useres.size(); ++iri ){
        int ir = useres[iri];
        if( allatom ){
            for( int ia = 1; ia <= pose.residue(ir).nheavyatoms(); ++ia ){
                numeric::xyzVector<double> coord = pose.residue(ir).xyz(ia);
                avgdis2 += cen.distance_squared( coord );
                maxdis = std::max( maxdis, (float)cen.distance( coord ) );
            }
        } else {
            numeric::xyzVector<double> coord;
            if(      pose.residue(ir).has("CB") ) coord = pose.residue(ir).xyz("CB");
            else if( pose.residue(ir).has("CA") ) coord = pose.residue(ir).xyz("CA");
            else                                  coord = pose.residue(ir).nbr_atom_xyz();
            avgdis2 += cen.distance_squared( coord );
            maxdis = std::max( maxdis, (float)cen.distance( coord ) );
        }
    }
    avgdis2 /= useres.size();
    rg = sqrt( avgdis2 );
    radius = maxdis;
}

// in radians
float
angle_between_vectors(
	Eigen::Vector3f const & vec1,
	Eigen::Vector3f const & vec2
) {
	const float dot = vec1.dot(vec2);
	const float norm_product = vec1.norm() * vec2.norm();
	return std::acos( dot / norm_product );
}

void
xform_pose(
	core::pose::Pose & pose,
	numeric::xyzTransform<float> s,
	core::Size sres/*=1*/,
	core::Size eres/*=0*/
) {
  if(eres==0) eres = pose.size();
  for(core::Size ir = sres; ir <= eres; ++ir) {
    for(core::Size ia = 1; ia <= pose.residue_type(ir).natoms(); ++ia) {
      core::id::AtomID const aid(core::id::AtomID(ia,ir));
      pose.set_xyz( aid, s*pose.xyz(aid) );
    }
  }
}

#ifdef USEGRIDSCORE
shared_ptr<protocols::ligand_docking::ga_ligand_dock::GridScorer>
prepare_grid_scorer(
	core::pose::Pose const & target,
	utility::vector1<core::Size> const & target_res,
	std::string const & atype_aas /*= "ACDEFGHIKLMNPQRSTVWY" */
) {

	using ObjexxFCL::format::F;

	std::vector<std::pair<core::Real,core::Real>> xyz_min_max(3);
	for ( std::pair<core::Real,core::Real> & pair : xyz_min_max ) { pair.first = 9e9; pair.second = -9e9; }

	for ( core::Size ir : target_res ) {
		for( int ia = 1; ia <= target.residue(ir).nheavyatoms(); ++ia ){
			numeric::xyzVector<core::Real> const & xyz = target.residue(ir).xyz(ia);

			for ( int j = 0; j < 3; j++ ) {
				xyz_min_max[j].first = std::min(xyz_min_max[j].first, xyz[j]);
				xyz_min_max[j].second = std::max(xyz_min_max[j].second, xyz[j]);
			}
		}
	}

	numeric::xyzVector< core::Real > centers(3);
	numeric::xyzVector< core::Real > radii(3);

	for ( int j = 0; j < 3; j++ ) {
		centers[j] = ( xyz_min_max[j].first + xyz_min_max[j].second ) / 2.0;
		radii[j] = ( xyz_min_max[j].second - centers[j] );
		runtime_assert( radii[j] > 0 );
	}

	std::cout << "Grid bounds: ";
	for ( int j = 0; j < 3; j++ ) {
		if (j > 0) std::cout << " x ";
		std::cout << "( " << F(5,1,xyz_min_max[j].first) << ", " << F(5,1,xyz_min_max[j].second) << " )";
	}
	std::cout << std::endl;

	radii += 5.0;	// max interaction cutoff

	std::cout << "Padded grid dimensions: ";
	for ( int j = 0; j < 3; j++ ) {
		if (j > 0) std::cout << " x ";
		std::cout << F(5,1,radii[j]);
	}
	std::cout << std::endl;

	std::cout << "Grid volume: " << F(7,1,radii[0]*radii[1]*radii[2]*8) << std::endl;

	core::scoring::ScoreFunctionOP scorefxn = core::scoring::ScoreFunctionFactory::create_score_function("beta_nov16");
	shared_ptr<protocols::ligand_docking::ga_ligand_dock::GridScorer> grid_scorer
		= make_shared<protocols::ligand_docking::ga_ligand_dock::GridScorer>( scorefxn );
	grid_scorer->set_grid_com_radii( centers, radii );

	// first A is the "ligand"
	std::string atype_aas_wlig = "A" + atype_aas;
	numeric::xyzVector< core::Real> far_away( 10000, 10000, 10000);
	numeric::xyzMatrix< core::Real> identity = numeric::xyzMatrix< core::Real>::identity();

	std::cout << "Using residue types: " << atype_aas << " for atom types" << std::endl;

	core::pose::PoseOP grid_pose = target.clone();
	utility::vector1< core::Size > fake_movingSCs;
	int ligand_resid = grid_pose->size() + 1;
	for ( int i = 0; i < atype_aas_wlig.size(); i++ ) {
		core::conformation::ResidueOP res = core::conformation::get_residue_from_name1( atype_aas_wlig[i] );
		res->apply_transform_Rx_plus_v( identity, far_away );
		grid_pose->append_residue_by_jump( *res, 1 );
		if ( i > 0 ) fake_movingSCs.push_back( grid_pose->size() );
	}

	grid_scorer->calculate_grid( *grid_pose, ligand_resid, fake_movingSCs );

	return grid_scorer;

}

#endif

}
}

max_stars_repo_name: willsheffler/rifdock

max_issues_repo_name: willsheffler/rifdock

max_forks_repo_name: willsheffler/rifdock


text:
#include "quadrature/qhermite.hpp"
#include "spectral/basis/spectral_basis.hpp"
#include "spectral/basis/spectral_elem.hpp"
#include "spectral/basis/spectral_elem_accessor.hpp"
#include "spectral/basis/spectral_function/hermite_polynomial.hpp"

#include <iomanip>
#include <iostream>
#include <tuple>

#include <Eigen/Sparse>

using namespace boltzmann;

/*

class HermiteHX : public HermiteH
{
public:
  HermiteHX(int k) : HermiteH(k, 0.5) {}
  HermiteHX() {}
};

class HermiteHY : public HermiteH
{
public:
  HermiteHY(int k) : HermiteH(k, 0.5) {}
  HermiteHY() {}
};

int main(int argc, char *argv[])
{
  // typedef HermiteH HermiteHX;
  // typedef HermiteH HermiteHY;

  typedef SpectralElem<double, HermiteHX, HermiteHY>
    spectral_elem_t;

  typedef SpectralBasis<spectral_elem_t> basis_type;
  basis_type basis;

  int K = 10;

  for (int k = 0; k < K; ++k) {
    for (int s = 0; s <= k; ++s) {
      HermiteHX hx(s);
      HermiteHY hy(k-s);
      basis.add_elem(hx, hy);
    }
  }

  //
  std::cout << "completed: adding elements to basis \n";

  SpectralElemAccessor::get<HermiteHX> getX;
  SpectralElemAccessor::get<HermiteHY> getY;

  for (auto elem_it = basis.begin(); elem_it < basis.end(); ++elem_it) {
    std::cout << getX(*elem_it).get_id() << "\t"
              << getY(*elem_it).get_id() << std::endl;
  }

  return 0;
}

*/

#include "quadrature/qhermite.hpp"
#include "spectral/basis/basis_types.hpp"
#include "spectral/basis/spectral_basis.hpp"
#include "spectral/basis/spectral_basis_factory_hermite.hpp"
#include "spectral/basis/spectral_basis_factory_ks.hpp"
#include "spectral/basis/spectral_elem.hpp"
#include "spectral/basis/spectral_elem_accessor.hpp"
#include "spectral/basis/spectral_function/hermite_polynomial.hpp"
#include "spectral/eval_handlers.hpp"

#include <Eigen/Sparse>
#include <boost/program_options.hpp>
#include <fstream>
#include <iomanip>
#include <iostream>
#include "spectral/polar_to_hermite.hpp"

#define PI 3.141592653589793238462643383279502884197

using namespace std;
using namespace boltzmann;

namespace po = boost::program_options;

int main(int argc, char *argv[])
{
  // this cannot work
  /* typedef HermiteH HermiteHX; */
  /* typedef HermiteH HermiteHY; */

  int K;
  po::options_description options("options");
  options.add_options()
      ("nK,K", po::value<int>(&K)->default_value(10), "K");
  po::variables_map vm;
  po::store(po::parse_command_line(argc, argv, options), vm);
  po::notify(vm);

  typedef typename SpectralBasisFactoryKS::basis_type polar_basis_t;
  polar_basis_t polar_basis;
  SpectralBasisFactoryKS::create(polar_basis, K, K, 2, true);
  SpectralBasisFactoryKS::write_basis_descriptor(polar_basis, "spectral_basis.desc");

  typedef typename SpectralBasisFactoryHN::basis_type hermite_basis_t;
  hermite_basis_t hermite_basis;
  SpectralBasisFactoryHN::create(hermite_basis, K, 2);
  SpectralBasisFactoryHN::write_basis_descriptor(hermite_basis, "hermite_basis.desc");

  typedef typename SpectralBasisFactoryHN::fa_type HX;
  typedef typename SpectralBasisFactoryHN::fa_type HY;

  // make evaluator handle
  std::vector<double> x{0.1, 0.5, 1, 10, 100};
  auto Heval2d = hermite_evaluator2d<>::make(hermite_basis, x);

  std::cout << "completed: adding elements to basis \n";

  SpectralElemAccessor::get<HermiteHX_t> getX;
  SpectralElemAccessor::get<HermiteHY_t> getY;

  for (auto elem_it = hermite_basis.begin(); elem_it < hermite_basis.end(); ++elem_it) {
    std::cout << getX(*elem_it).get_id() << "\t" << getY(*elem_it).get_id() << std::endl;
  }

  return 0;
}


text:
/** @file
 * @brief NPDE OutputImpedanceBVP
 * @author Erick Schulz
 * @date 12/07/2019
 * @copyright Developed at ETH Zurich
 */

#include "outputimpedancebvp.h"

#include <lf/assemble/assemble.h>
#include <lf/geometry/geometry.h>
#include <lf/mesh/utils/utils.h>
#include <lf/uscalfe/uscalfe.h>

#include <Eigen/Dense>
#include <Eigen/Sparse>
#include <cassert>

namespace OutputImpedanceBVP {

/* SAM_LISTING_BEGIN_1 */
Eigen::VectorXd solveImpedanceBVP(
    const std::shared_ptr<lf::uscalfe::FeSpaceLagrangeO1<double>> &fe_space_p,
    Eigen::Vector2d g) {
  // Related implementations:
  // Homework problem ErrorEstimatesForTraces:
  // https://gitlab.math.ethz.ch/ralfh/npdecodes/tree/master/homeworks/ErrorEstimatesForTraces

  // Pointer to current mesh
  std::shared_ptr<const lf::mesh::Mesh> mesh_p = fe_space_p->Mesh();
  // Obtain local->global index mapping for current finite element space
  const lf::assemble::DofHandler &dofh{fe_space_p->LocGlobMap()};
  // Dimension of finite element space
  const lf::uscalfe::size_type N_dofs(dofh.NumDofs());
  // Obtain specification for shape functions on edges
  std::shared_ptr<const lf::uscalfe::ScalarReferenceFiniteElement<double>>
      rsf_edge_p = fe_space_p->ShapeFunctionLayout(lf::base::RefEl::kSegment());

  Eigen::VectorXd discrete_solution(N_dofs);

  // I : ASSEMBLY
  // Matrix in triplet format holding Galerkin matrix, zero initially.
  lf::assemble::COOMatrix<double> A(N_dofs, N_dofs);
  // Right hand side vector, must be initialized with 0!
  Eigen::Matrix<double, Eigen::Dynamic, 1> phi(N_dofs);
  phi.setZero();

  // I.i : Computing volume matrix for negative Laplace operator
  // Initialize object taking care of local mass (volume) computations.
  lf::uscalfe::LinearFELaplaceElementMatrix elmat_builder{};
  // Invoke assembly on cells (co-dimension = 0 as first argument)
  // Information about the mesh and the local-to-global map is passed through
  // a Dofhandler object, argument 'dofh'. This function call adds triplets to
  // the internal COO-format representation of the sparse matrix A.
  lf::assemble::AssembleMatrixLocally(0, dofh, dofh, elmat_builder, A);
  /* SAM_LISTING_END_1 */

  /* SAM_LISTING_BEGIN_2 */
  // I.ii : Computing mass edge matrix resulting from Robin B.C.
  // Obtain an array of boolean flags for the edges of the mesh, 'true'
  // indicates that the edge lies on the boundary
  auto bd_flags{lf::mesh::utils::flagEntitiesOnBoundary(mesh_p, 1)};
  // Creating a predicate that will guarantee that the computations are carried
  // only on the interior boundary edges of the mesh using the boundary flags
  auto edges_predicate_RobinBC =
      [&bd_flags](const lf::mesh::Entity &edge) -> bool {
    if (bd_flags(edge)) {
      auto endpoints = lf::geometry::Corners(*(edge.Geometry()));
      if (endpoints(0, 0) <= 0.05 || 0.95 <= endpoints(0, 0) ||
          endpoints(1, 0) <= 0.05 || 0.95 <= endpoints(1, 0)) {
        return false;
      }
      return true;
    }
    return false;
  };
  // Coefficients used in the class template
  // MassEdgeMatrixProvider< SCALAR, COEFF, EDGESELECTOR >
  auto eta = lf::mesh::utils::MeshFunctionGlobal(
      [](Eigen::Vector2d x) -> double { return 1.0; });
  lf::uscalfe::MassEdgeMatrixProvider<double, decltype(eta),
                                      decltype(edges_predicate_RobinBC)>
      edgemat_builder(fe_space_p, eta, edges_predicate_RobinBC);
  // Invoke assembly on edges by specifying co-dimension = 1
  lf::assemble::AssembleMatrixLocally(1, dofh, dofh, edgemat_builder, A);
  /* SAM_LISTING_END_2 */

  /* SAM_LISTING_BEGIN_9 */
  // I.iii : Computing right-hand side vector
  // Right-hand side source function f
  auto mf_f = lf::mesh::utils::MeshFunctionGlobal(
      [](Eigen::Vector2d x) -> double { return 0.0; });
  lf::uscalfe::ScalarLoadElementVectorProvider<double, decltype(mf_f)>
      elvec_builder(fe_space_p, mf_f);
  // Invoke assembly on cells (codim == 0)
  AssembleVectorLocally(0, dofh, elvec_builder, phi);

  // I.iv : Imposing essential boundary conditions
  // Dirichlet data
  auto mf_g = lf::mesh::utils::MeshFunctionGlobal(
      [&g](Eigen::Vector2d x) -> double { return g.dot(x); });
  // Inspired by the example in the documentation of
  // InitEssentialConditionFromFunction()
  // https://craffael.github.io/lehrfempp/namespacelf_1_1uscalfe.html#a5afbd94919f0382cf3fb200c452797ac
  // Creating a predicate that will guarantee that the computations are carried
  // only on the exterior boundary edges of the mesh using the boundary flags
  auto edges_predicate_Dirichlet =
      [&bd_flags](const lf::mesh::Entity &edge) -> bool {
    if (bd_flags(edge)) {
      auto endpoints = lf::geometry::Corners(*(edge.Geometry()));
      if (endpoints(0, 0) <= 0.05 || 0.95 <= endpoints(0, 0) ||
          endpoints(1, 0) <= 0.05 || 0.95 <= endpoints(1, 0)) {
        return true;
      }
    }
    return false;
  };
  // Determine the fixed dofs on the boundary and their values
  // Alternative: See lecturedemoDirichlet() in
  // https://github.com/craffael/lehrfempp/blob/master/examples/lecturedemos/lecturedemoassemble.cc
  auto edges_flag_values_Dirichlet{
      lf::uscalfe::InitEssentialConditionFromFunction(
          dofh, *rsf_edge_p, edges_predicate_Dirichlet, mf_g)};
  // Eliminate Dirichlet dofs from the linear system
  lf::assemble::FixFlaggedSolutionCompAlt<double>(
      [&edges_flag_values_Dirichlet](lf::assemble::glb_idx_t gdof_idx) {
        return edges_flag_values_Dirichlet[gdof_idx];
      },
      A, phi);

  // Assembly completed! Convert COO matrix A into CRS format using Eigen's
  // internal conversion routines.
  Eigen::SparseMatrix<double> A_sparse = A.makeSparse();

// II : SOLVING  THE LINEAR SYSTEM
  // II.i : Setting up Eigen's sparse direct elimination
  Eigen::SparseLU<Eigen::SparseMatrix<double>> solver;
  solver.compute(A_sparse);
  LF_VERIFY_MSG(solver.info() == Eigen::Success, "LU decomposition failed");
  // II.ii : Solving
  discrete_solution = solver.solve(phi);
  LF_VERIFY_MSG(solver.info() == Eigen::Success, "Solving LSE failed");

// do nothing
  return discrete_solution;
};
/* SAM_LISTING_END_9 */

/* SAM_LISTING_BEGIN_3 */
double computeBoundaryOutputFunctional(
    const Eigen::VectorXd eta,
    const std::shared_ptr<lf::uscalfe::FeSpaceLagrangeO1<double>> &fe_space_p,
    Eigen::Vector2d d) {
  double func_val = 0.0;
  // Pointer to current mesh
  std::shared_ptr<const lf::mesh::Mesh> mesh_p = fe_space_p->Mesh();
  // Obtain local->global index mapping for current finite element space
  const lf::assemble::DofHandler &dofh{fe_space_p->LocGlobMap()};

  // Obtain an array of boolean flags for the edges of the mesh, 'true'
  // indicates that the edge lies on the boundary
  auto bd_flags{lf::mesh::utils::flagEntitiesOnBoundary(mesh_p, 1)};

  // Creating a predicate that will guarantee that the computations are carried
  // only on the interior boundary edges of the mesh using the boundary flags
  auto edges_predicate_RobinBC =
      [&bd_flags](const lf::mesh::Entity &edge) -> bool {
    if (bd_flags(edge)) {
      auto endpoints = lf::geometry::Corners(*(edge.Geometry()));
      if (endpoints(0, 0) <= 0.05 || 0.95 <= endpoints(0, 0) ||
          endpoints(1, 0) <= 0.05 || 0.95 <= endpoints(1, 0)) {
        return false;
      }
      return true;
    }
    return false;
  };

  // Computing value of the functional
  for (const lf::mesh::Entity *edge : mesh_p->Entities(1)) {
    if (edges_predicate_RobinBC(*edge)) {
      // Find the endpoints global indices
      auto dof_idx = dofh.GlobalDofIndices(*edge);
      assert(dofh.NumLocalDofs(*edge) == 2);
      // Value of linear function in he endpoints
      const double nu0 = eta(dof_idx[0]);
      const double nu1 = eta(dof_idx[1]);
      // Find coordinates of endpoints
      auto endpoints = lf::geometry::Corners(*(edge->Geometry()));
      // Length of edge
      const double elen = (endpoints.col(1) - endpoints.col(0)).norm();
      LF_ASSERT_MSG(endpoints.cols() == 2, "Wrong no endpoints in " << *edge);
      func_val += elen / 6.0 *
                  d.dot(2 * endpoints.col(0) * nu0 + endpoints.col(0) * nu1 +
                        endpoints.col(1) * nu0 + 2 * endpoints.col(1) * nu1);
    }
  }
  return func_val;
};
/* SAM_LISTING_END_3 */

}  // namespace OutputImpedanceBVP

max_stars_repo_name: padomu/NPDECODES

max_issues_repo_name: padomu/NPDECODES

max_forks_repo_name: padomu/NPDECODES