text:
// ExodusModel.cpp
// created by Kuangdai on 2-May-2016
// read Exodus mesh file

// Adapted from salvus project
// ExodusModel.cpp by Michael Afanasiev

#include "ExodusModel.h"
#include <sstream>
#include <iomanip>
#include <cmath>
#include <cfloat>

#include "XMPI.h"
#include "Mapping.h"
#include <boost/algorithm/string.hpp>

#include "AttBuilder.h"
#include "Geodesy.h"
#include "Parameters.h"
#include "MultilevelTimer.h"

#include "NetCDF_Reader.h"

ExodusModel::ExodusModel(const std::string &fileName): mExodusFileName(fileName) {
    std::vector<std::string> substrs = Parameters::splitString(mExodusFileName, "/");
    mExodusTitle = substrs[substrs.size() - 1];
    boost::trim_if(mExodusTitle, boost::is_any_of("\t "));
}

void ExodusModel::initialize() {
    MultilevelTimer::begin("Read Exodus", 1);
    if (XMPI::root()) {
        readRawData();
    }
    MultilevelTimer::end("Read Exodus", 1);

    MultilevelTimer::begin("Bcast Exodus", 1);
    bcastRawData();
    MultilevelTimer::end("Bcast Exodus", 1);

    MultilevelTimer::begin("Process Exodus", 1);
    formStructured();
    formAuxiliary();
    MultilevelTimer::end("Process Exodus", 1);
}

void ExodusModel::readRawData() {
    // open file
    NetCDF_Reader reader;
    reader.open(mExodusFileName);
    RDMatXX dbuffer;

    // global
    reader.readString("name_glo_var", mGlobalVariableNames);
    reader.read1D("vals_glo_var", mGlobalVariableValues);
    reader.readString("info_records", mGlobalRecordsRaw);

    // connectivity and coords
    reader.read2D("connect1", mConnectivity);
    mConnectivity.array() -= 1;
    reader.read1D("coordx", mNodalS);
    reader.read1D("coordy", mNodalZ);

    // elemental variables
    reader.readString("name_elem_var", mElementalVariableNames);
    mElementalVariableValues = RDMatXX::Zero(getNumQuads(), mElementalVariableNames.size());
    for (int i = 0; i < mElementalVariableNames.size(); i++) {
        std::stringstream ss;
        ss << "vals_elem_var" << i + 1 << "eb1";
        reader.read2D(ss.str(), dbuffer);
        mElementalVariableValues.col(i) = dbuffer.transpose();
    }

    // side sets
    reader.readString("ss_names", mSideSetNames);
    mSideSetValues = IMatXX::Zero(getNumQuads(), mSideSetNames.size());
    for (int i = 0; i < mSideSetNames.size(); i++) {
        std::stringstream sse, sss;
        IColX elems, sides;
        // elem
        sse << "elem_ss" << i + 1;
        reader.read1D(sse.str(), elems);
        // side
        sss << "side_ss" << i + 1;
        reader.read1D(sss.str(), sides);
        // processing
        IColX values = IColX::Constant(getNumQuads(), -1);
        for (int j = 0; j < sides.rows(); j++) {
            values(elems(j) - 1) = sides(j) - 1;
        }
        mSideSetValues.col(i) = values;
    }

    // do not read ellipticity for cartesian
    bool cartesian = false;
    for (int i = 0; i < mGlobalRecordsRaw.size(); i++) {
        std::vector<std::string> substrs = Parameters::splitString(mGlobalRecordsRaw[i], "=");
        if (boost::iequals(boost::trim_copy(substrs[0]), "crdsys") &&
            boost::iequals(boost::trim_copy(substrs[1]), "cartesian")) {
            cartesian = true;
            break;
        }
    }
    if (cartesian) {
        mEllipKnots = RDColX::Zero(0);
        mEllipCoeffs = RDColX::Zero(0);
        // set radius to PREM
        mGlobalVariables.insert(std::pair<std::string, double>("radius", 6371e3));
    } else {
        // ellipticity
        reader.read2D("ellipticity", dbuffer);
        mEllipKnots = dbuffer.row(0).transpose();
        mEllipCoeffs = dbuffer.row(1).transpose();
    }

    // close file
    reader.close();
}

void ExodusModel::bcastRawData() {
    XMPI::bcast(mGlobalVariableNames);
    XMPI::bcastEigen(mGlobalVariableValues);
    XMPI::bcast(mGlobalRecordsRaw);

    XMPI::bcastEigen(mConnectivity);
    XMPI::bcastEigen(mNodalS);
    XMPI::bcastEigen(mNodalZ);

    XMPI::bcast(mElementalVariableNames);
    XMPI::bcastEigen(mElementalVariableValues);

    XMPI::bcast(mSideSetNames);
    XMPI::bcastEigen(mSideSetValues);

    XMPI::bcastEigen(mEllipKnots);
    XMPI::bcastEigen(mEllipCoeffs);
}

void ExodusModel::formStructured() {
    // global variables
    for (int i = 0; i < mGlobalVariableNames.size(); i++) {
        std::string varName = mGlobalVariableNames[i];
        if (varName == "dt") {
            varName = "dt (nPol = 1)";
        }
        mGlobalVariables.insert(std::pair<std::string, double>(varName, mGlobalVariableValues(i)));
    }

    // global records
    std::vector<std::string> included = {"crdsys", "model"};
    for (int i = 0; i < mGlobalRecordsRaw.size(); i++) {
        std::vector<std::string> substrs = Parameters::splitString(mGlobalRecordsRaw[i], "=");
        if (std::find(included.begin(), included.end(), boost::trim_copy(substrs[0])) != included.end()) {
            mGlobalRecords.insert(std::pair<std::string, std::string>(
                boost::trim_copy(substrs[0]), boost::trim_copy(substrs[1])));
        }
    }

    // elemental variables
    for (int i = 0; i < mElementalVariableNames.size(); i++) {
        mElementalVariables.insert(std::pair<std::string, RDColX>(mElementalVariableNames[i],
            mElementalVariableValues.col(i)));
    }

    // side sets
    for (int i = 0; i < mSideSetNames.size(); i++) {
        mSideSets.insert(std::pair<std::string, IColX>(mSideSetNames[i],
            mSideSetValues.col(i)));
    }

    // name of axis and surface sets
    std::string sphere_axis_name = "t0";
    if (!isCartesian()) {
        if (mSideSets.find("t1") != mSideSets.end()) {
            sphere_axis_name = "t1";
        }
    }

    mSSNameAxis = isCartesian() ? "x0" : sphere_axis_name;
    mSSNameSurface = isCartesian() ? "y1" : "r1";

    // NOTE: we temporarily treat Cartesian meshes as special cases of spherical meshes
    //       by means of moving it to the "north pole". The introduced global
    //       curvature should be ignorable, or the problem itself is ill-defined
    //       as a local problem.
    if (isCartesian()) {
        double R_EARTH = getROuter();
        double maxz = mNodalZ.maxCoeff();
        mNodalZ.array() += R_EARTH - maxz;
    }
}

void ExodusModel::formAuxiliary() {
    MultilevelTimer::begin("Process Exodus DistTol", 2);
    // distance tolerance
    double distTol = DBL_MAX;
    for (int i = 0; i < getNumQuads(); i++) {
        if (i % XMPI::nproc() != XMPI::rank()) {
            continue;
        }
        double s0 = mNodalS(mConnectivity(i, 0));
        double z0 = mNodalZ(mConnectivity(i, 0));
        double s1 = mNodalS(mConnectivity(i, 1));
        double z1 = mNodalZ(mConnectivity(i, 1));
        double s2 = mNodalS(mConnectivity(i, 2));
        double z2 = mNodalZ(mConnectivity(i, 2));
        double s3 = mNodalS(mConnectivity(i, 3));
        double z3 = mNodalZ(mConnectivity(i, 3));
        double dist0 = sqrt((s0 - s1) * (s0 - s1) + (z0 - z1) * (z0 - z1)) / 1000.;
        double dist1 = sqrt((s1 - s2) * (s1 - s2) + (z1 - z2) * (z1 - z2)) / 1000.;
        double dist2 = sqrt((s2 - s3) * (s2 - s3) + (z2 - z3) * (z2 - z3)) / 1000.;
        double dist3 = sqrt((s3 - s0) * (s3 - s0) + (z3 - z0) * (z3 - z0)) / 1000.;
        distTol = std::min({dist0, dist1, dist2, dist3, distTol});
    }
    mDistTolerance = XMPI::min(distTol);
    MultilevelTimer::end("Process Exodus DistTol", 2);

    // average gll spacing
    MultilevelTimer::begin("Process Exodus GLL-Spacing", 2);
    std::vector<std::vector<int>> refElem(getNumNodes(), std::vector<int>());
    for (int i = 0; i < getNumQuads(); i++) {
        refElem[mConnectivity(i, 0)].push_back(i);
        refElem[mConnectivity(i, 1)].push_back(i);
        refElem[mConnectivity(i, 2)].push_back(i);
        refElem[mConnectivity(i, 3)].push_back(i);
    }
    mAveGLLSpacing = RDColX::Zero(getNumNodes());
    for (int i = 0; i < getNumNodes(); i++) {
        if (i % XMPI::nproc() != XMPI::rank()) {
            continue;
        }
        for (int j = 0; j < refElem[i].size(); j++) {
            int ielem = refElem[i][j];
            double s0 = mNodalS(mConnectivity(ielem, 0));
            double z0 = mNodalZ(mConnectivity(ielem, 0));
            double s1 = mNodalS(mConnectivity(ielem, 1));
            double z1 = mNodalZ(mConnectivity(ielem, 1));
            double s2 = mNodalS(mConnectivity(ielem, 2));
            double z2 = mNodalZ(mConnectivity(ielem, 2));
            double s3 = mNodalS(mConnectivity(ielem, 3));
            double z3 = mNodalZ(mConnectivity(ielem, 3));
            double dist0 = sqrt((s0 - s1) * (s0 - s1) + (z0 - z1) * (z0 - z1));
            double dist1 = sqrt((s1 - s2) * (s1 - s2) + (z1 - z2) * (z1 - z2));
            double dist2 = sqrt((s2 - s3) * (s2 - s3) + (z2 - z3) * (z2 - z3));
            double dist3 = sqrt((s3 - s0) * (s3 - s0) + (z3 - z0) * (z3 - z0));
            mAveGLLSpacing(i) += (dist0 + dist1 + dist2 + dist3) / 4. / nPol / refElem[i].size();
        }
    }
    XMPI::sumEigenDouble(mAveGLLSpacing);
    MultilevelTimer::end("Process Exodus GLL-Spacing", 2);

    // rotate nodes of axial elements such that side 3 is on axis
    MultilevelTimer::begin("Process Exodus Axis", 2);
    for (int axialQuad = 0; axialQuad < getNumQuads(); axialQuad++) {
        // loop over t0
        int axialSide = getSideAxis(axialQuad);
        if (axialSide == 3 || axialSide == -1) {
            continue;
        }

        // connectivity
        IRow4 con = mConnectivity.row(axialQuad);
        for (int j = 0; j < 4; j++) {
            mConnectivity(axialQuad, j) = con(Mapping::period0123(j + axialSide - 3));
        }

        // elemental fields
        for (auto it = mElementalVariables.begin(); it != mElementalVariables.end(); it++) {
            std::string vname = it->first;
            if (vname.substr(vname.length() - 2, 2) == std::string("_0")) {
                vname = vname.substr(0, vname.length() - 2);
                RDRow4 v_old;
                v_old(0) = mElementalVariables.at(vname + "_0")(axialQuad);
                v_old(1) = mElementalVariables.at(vname + "_1")(axialQuad);
                v_old(2) = mElementalVariables.at(vname + "_2")(axialQuad);
                v_old(3) = mElementalVariables.at(vname + "_3")(axialQuad);
                mElementalVariables.at(vname + "_0")(axialQuad) = v_old(Mapping::period0123(0 + axialSide - 3));
                mElementalVariables.at(vname + "_1")(axialQuad) = v_old(Mapping::period0123(1 + axialSide - 3));
                mElementalVariables.at(vname + "_2")(axialQuad) = v_old(Mapping::period0123(2 + axialSide - 3));
                mElementalVariables.at(vname + "_3")(axialQuad) = v_old(Mapping::period0123(3 + axialSide - 3));
            }
        }

        // side sets
        for (auto it = mSideSets.begin(); it != mSideSets.end(); it++) {
            if (it->second(axialQuad) != -1) {
                it->second(axialQuad) = Mapping::period0123(it->second(axialQuad) - axialSide + 3);
            }
        }

        // done
        // mSideSets.at(mSSNameAxis)(axialQuad) = 3;
    }
    MultilevelTimer::end("Process Exodus Axis", 2);

    // find elements that are not axial but neighboring axial elements
    MultilevelTimer::begin("Process Exodus Vicinal", 2);
    mVicinalAxis = IMatX4::Constant(getNumQuads(), 4, -1);
    // first find near-axis nodes and axial quads
    std::vector<bool> nodeNearAxis(getNumNodes(), false);
    std::vector<bool> quadOnAxis(getNumQuads(), false);
    for (int axialQuad = 0; axialQuad < getNumQuads(); axialQuad++) {
        int axialSide = getSideAxis(axialQuad);
        if (axialSide == -1) {
            continue;
        }
        nodeNearAxis[mConnectivity(axialQuad, 0)] = true;
        nodeNearAxis[mConnectivity(axialQuad, 1)] = true;
        nodeNearAxis[mConnectivity(axialQuad, 2)] = true;
        nodeNearAxis[mConnectivity(axialQuad, 3)] = true;
        quadOnAxis[axialQuad] = true;
    }
    // loop over quads
    for (int iquad = 0; iquad < getNumQuads(); iquad++) {
        if (quadOnAxis[iquad]) continue;
        for (int j = 0; j < 4; j++) {
            int nTag = mConnectivity(iquad, j);
            if (nodeNearAxis[nTag]) mVicinalAxis(iquad, j) = j;
        }
    }
    MultilevelTimer::end("Process Exodus Vicinal", 2);

    // check if ocean presents in mesh
    MultilevelTimer::begin("Process Exodus Check Ocean", 2);
    std::string strVs = isIsotropic() ? "VS_0" : "VSV_0";
    for (int iQuad = 0; iQuad < getNumQuads(); iQuad++) {
        if (iQuad % XMPI::nproc() != XMPI::rank()) {
            continue;
        }
        if (getSideSurface(iQuad) == -1) {
            continue;
        }
        double vs = mElementalVariables.at(strVs)(iQuad);
        if (vs < tinyDouble) {
            throw std::runtime_error("ExodusModel::finishReading || "
                "Ocean is detected in mesh. By far, realistic ocean is not implemented in AxiSEM3D. ||"
                "Use non-ocean models in the Mesher and add ocean load in inparam.basic.");
        }
    }
    MultilevelTimer::end("Process Exodus Check Ocean", 2);

    // CMB and ICB
    // MultilevelTimer::begin("Process Exodus CMB & ICB", 2);
    // mR_CMB = DBL_MIN;
    // mR_ICB = DBL_MAX;
    // for (int i = 0; i < getNumQuads(); i++) {
    //     bool isFluid = getElementalVariables("fluid", i) > .5;
    //     bool isAxis = getSideAxis(i) >= 0;
    //     if (isFluid && isAxis) {
    //         for (int j = 0; j < 4; j++) {
    //             double s = mNodalS(mConnectivity(i, j));
    //             double z = mNodalZ(mConnectivity(i, j));
    //             double r = std::sqrt(s * s + z * z);
    //             mR_CMB = std::max(mR_CMB, r);
    //             mR_ICB = std::min(mR_ICB, r);
    //         }
    //     }
    // }
    // MultilevelTimer::end("Process Exodus CMB & ICB", 2);
}

std::string ExodusModel::verbose() const {
    std::stringstream ss;
    ss << "\n======================= Exodus Model =======================" << std::endl;
    ss << "  Overview__________________________________________________" << std::endl;
    ss << "    Exodus Title      =   " << mExodusTitle << std::endl;
    ss << "    Mesh CS Type      =   " << (isCartesian() ? "Cartesian" : "Spherical") << std::endl;
    ss << "    Number of Nodes   =   " << getNumNodes() << std::endl;
    ss << "    Number of Quads   =   " << getNumQuads() << std::endl;
    ss << "  Global Variables__________________________________________" << std::endl;
    int widthname = -1;
    for (auto it = mGlobalVariables.begin(); it != mGlobalVariables.end(); it++) {
        widthname = std::max(widthname, (int)(it->first.length()));
    }
    for (auto it = mGlobalRecords.begin(); it != mGlobalRecords.end(); it++) {
        widthname = std::max(widthname, (int)(it->first.length()));
    }
    for (auto it = mGlobalVariables.begin(); it != mGlobalVariables.end(); it++) {
        ss << "    " << std::setw(widthname) << it->first << "   =   " << it->second << std::endl;
    }
    for (auto it = mGlobalRecords.begin(); it != mGlobalRecords.end(); it++) {
        ss << "    " << std::setw(widthname) << it->first << "   =   " << it->second << std::endl;
    }
    ss << "  Connectivity______________________________________________" << std::endl;
    int width = (int)std::log10(std::max(getNumQuads(), getNumNodes())) + 1;
    ss << "    " << std::setw(width) << 0 << ": ";
    for (int j = 0; j < 4; j++) {
        ss << std::setw(width) << mConnectivity(0, j) << " ";
    }
    ss << std::endl << "    " << std::setw(width) << "..." << std::endl;
    ss << "    " << std::setw(width) << getNumQuads() - 1 << ": ";
    for (int j = 0; j < 4; j++) {
        ss << std::setw(width) << mConnectivity(getNumQuads() - 1, j) << " ";
    }
    ss << std::endl;
    ss << "  Coordinates_______________________________________________" << std::endl;
    ss << "    " << std::setw(width) << 0 << ": ";
    ss << std::setw(13) << mNodalS(0) << std::setw(13) << mNodalZ(0) << std::endl;
    ss<< "    " << std::setw(width) << "..." << std::endl;
    ss << "    " << std::setw(width) << getNumNodes() - 1 << ": ";
    ss << std::setw(13) << mNodalS(getNumNodes() - 1) << std::setw(13) << mNodalZ(getNumNodes() - 1) << std::endl;
    ss << "  Elemental Variables_______________________________________" << std::endl;
    widthname = -1;
    for (auto it = mElementalVariables.begin(); it != mElementalVariables.end(); it++) {
        widthname = std::max(widthname, (int)(it->first.length()));
    }
    for (auto it = mElementalVariables.begin(); it != mElementalVariables.end(); it++) {
        ss << "    " << std::setw(widthname) << it->first << ": ";
        ss << std::setw(13) << it->second(0) << ", ..., ";
        ss << std::setw(13) << it->second(getNumQuads() - 1) << std::endl;
    }
    ss << "  Side Sets_________________________________________________" << std::endl;
    widthname = -1;
    for (auto it = mSideSets.begin(); it != mSideSets.end(); it++) {
        widthname = std::max(widthname, (int)(it->first.length()));
    }
    for (auto it = mSideSets.begin(); it != mSideSets.end(); it++) {
        ss << "    " << std::setw(widthname) << it->first << ":   ";
        int pair = 0;
        for (int q = 0; q < getNumQuads(); q++) {
            pair += (int)(it->second(q) >= 0);
        }
        ss << pair << " edges" << std::endl;
    }
    ss << "  Miscellaneous_____________________________________________" << std::endl;
    ss << "    Distance Tolerance / m   =   " << mDistTolerance << std::endl;
    if (!isCartesian()) {
        ss << "  External__________________________________________________" << std::endl;
        ss << "    Num. Ellipticity Spline Knots   =   " << mEllipKnots.size() << std::endl;
    }
    ss << "======================= Exodus Model =======================\n" << std::endl;
    return ss.str();
}

void ExodusModel::buildInparam(ExodusModel *&exModel, const Parameters &par,
    AttParameters *&attPar, int verbose) {
    if (exModel) {
        delete exModel;
    }
    std::string exfile = par.getValue<std::string>("MODEL_1D_EXODUS_MESH_FILE");
    exfile = Parameters::sInputDirectory + "/" + exfile;
    exModel = new ExodusModel(exfile);
    exModel->initialize();
    if (verbose) {
        XMPI::cout << exModel->verbose();
    }

    // form attenuation parameters
    if (exModel->hasAttenuation()) {
        int nr_lin_solids = (int)exModel->mGlobalVariables.at("nr_lin_solids");
        double f_min = exModel->mGlobalVariables.at("f_min");
        double f_max = exModel->mGlobalVariables.at("f_max");
        double f_ref = exModel->mGlobalVariables.at("f_ref");
        RDColX w(nr_lin_solids), y(nr_lin_solids);
        for (int i = 0; i < nr_lin_solids; i++) {
            std::stringstream sw, sy;
            sw << "w_" << i;
            sy << "y_" << i;
            w(i) = exModel->mGlobalVariables.at(sw.str());
            y(i) = exModel->mGlobalVariables.at(sy.str());
        }
        if (attPar) {
            delete attPar;
        }
        attPar = new AttParameters(nr_lin_solids, f_min, f_max, f_ref, w, y);
    } else {
        if (attPar) {
            delete attPar;
            attPar = 0;
        }
    }

    // ellipticity
    if (exModel->isCartesian()) {
        return;
    }
    std::string emode = par.getValue<std::string>("MODEL_3D_ELLIPTICITY_MODE");
    if (boost::iequals(emode, "off")) {
        // no ellipticity
        Geodesy::setup(exModel->getROuter(), 0., RDColX::Zero(0), RDColX::Zero(0));
    } else {
        double inv_f = par.getValue<double>("MODEL_3D_ELLIPTICITY_INVF");
        if (inv_f <= 0.) {
            throw std::runtime_error("ExodusModel::buildInparam || Invalid flattening.");
        }
        Geodesy::setup(exModel->getROuter(), 1. / inv_f, exModel->mEllipKnots, exModel->mEllipCoeffs);
    }
}


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

#include <boost/simd/sdk/simd/native.hpp>
#include <boost/type_traits/is_same.hpp>
#include <boost/utility/enable_if.hpp>

namespace boost { namespace simd { namespace details
{
  template<class T> inline
  typename  boost
            ::enable_if_c < (  boost::is_same<T,int8_t>::value
                            || boost::is_same<T,uint8_t>::value
                            || boost::is_same<T,char>::value
                            )
                          , int
                          >::type
  display ( T const& v)
  {
    return static_cast<int>(v);
  }

  template<class T> inline
  typename  boost
            ::disable_if_c< (  boost::is_same<T,int8_t>::value
                            || boost::is_same<T,uint8_t>::value
                            || boost::is_same<T,char>::value
                            )
                          , T const&
                          >::type
  display ( T const& v)
  {
    return v;
  }
} } }

#endif

ext: hpp


text:
#pragma once

#include "../../../config.hpp"
#include <fmt/format.h>
#include <fmt/ostream.h>

#include <boost/regex.hpp>
#include <vector>

#include "token.hpp"
#include "error.hpp"

#include "../../json.hpp"

namespace warmane::armory::site::parser
{
	class html_parser
	{
	public:
		using iterator = std::vector<token>::iterator;
		using const_iterator = std::vector<token>::const_iterator;

	public:
		html_parser() = default;
		explicit html_parser(const std::vector<token>& tokens);
		explicit html_parser(std::vector<token>&& tokens);

		symbol peek() const;

		void expect(symbol s) const;

		const token& next_token();

		const std::vector<token>& tokens() const noexcept;

	private:
		//void find_

	private:
		std::vector<token> tokens_;
		iterator it_;
		//std::unordered_map<std::string, >
			//data_finder_;
	};

	inline void to_json(armory::json& j, const html_parser& p)
	{

	}

	inline void from_json(const armory::json& j, html_parser& p)
	{

	}

	inline html_parser::html_parser(const std::vector<token>& tokens)
		: tokens_{tokens}
		, it_{std::begin(tokens_)}
	{}

	inline html_parser::html_parser(std::vector<token>&& tokens)
		: tokens_{std::move(tokens)}
		, it_{std::begin(tokens_)}
	{}

	inline symbol html_parser::peek() const
	{
		return it_->symbol();
	}

	inline void html_parser::expect(const symbol s) const
	{
		const symbol peeked = this->peek();
		if (peeked != s)
			throw unexpected_token{
				fmt::format("[html_parser::expect] Expected token:{}, but got token:{}",
					s, peeked
				)
			};
	}

	inline const token& html_parser::next_token()
	{
		return *(it_++);
	}

	inline const std::vector<token>& html_parser::tokens() const noexcept
	{
		return tokens_;
	}
}

max_stars_repo_name: Ostoic/warmane-api

max_issues_repo_name: Ostoic/warmane-api

max_forks_repo_name: Ostoic/warmane-api


text:
#pragma once

#include "construct.hpp"
#include "fwd.hpp"
#include "permit.hpp"
#include "vec_base.hpp"

#include <boost/fusion/algorithm/iteration/for_each.hpp>

#include <type_traits>

namespace v2 {

	template <typename T, unsigned int N>
	struct vector
		: public vec_base<T, N>
	{
	public: // statics

		using value_type = T;
		using size = std::integral_constant<unsigned int, N>;

	public: // methods

		// basic constructors
		constexpr vector() = default;
		constexpr vector(vector&&) = default;
		constexpr vector(const vector&) = default;

		// all-dimension construction
		template <typename... Ts, typename = std::enable_if_t<(sizeof...(Ts) == N)>>
		constexpr vector(Ts&&... args)
			: vec_base<T, N>{{{{ std::forward<Ts>(args)... }}}}
		{
		}

		// from other permitted types
		template <typename Seq, typename = std::enable_if_t<permitted<Seq>::value>>
		constexpr vector(Seq&& seq)
			: vec_base<T, N>{{
				boost::fusion::invoke(detail::array_factory<T, N>{}, seq)
			}}
		{
		}

		// dimension access
		constexpr T& operator[](std::size_t idx)
		{
			return this->vec_base<T, N>::data[idx];
		}

		constexpr T const& operator[](std::size_t idx) const
		{
			return this->vec_base<T, N>::data[idx];
		}
	};

} // namespace v2

max_stars_repo_name: ralismark/variolite

max_issues_repo_name: ralismark/variolite

max_forks_repo_name: ralismark/variolite


text:
//==================================================================================================
/*!
  @file

  @copyright 2016 NumScale SAS

  Distributed under the Boost Software License, Version 1.0.
  (See accompanying file LICENSE.md or copy at http://boost.org/LICENSE_1_0.txt)
*/
//==================================================================================================
#ifndef BOOST_SIMD_ARCH_COMMON_SCALAR_FUNCTION_DIST_S_HPP_INCLUDED
#define BOOST_SIMD_ARCH_COMMON_SCALAR_FUNCTION_DIST_S_HPP_INCLUDED

#include <boost/simd/constant/valmax.hpp>
#include <boost/simd/constant/zero.hpp>
#include <boost/simd/function/dist.hpp>
#include <boost/simd/detail/dispatch/function/overload.hpp>
#include <boost/config.hpp>

namespace boost { namespace simd { namespace ext
{
  namespace bd = boost::dispatch;
  namespace bs = boost::simd;
  BOOST_DISPATCH_OVERLOAD ( dist_
                          , (typename A0)
                          , bd::cpu_
                          , bs::saturated_tag
                          , bd::scalar_< bd::arithmetic_<A0> >
                          , bd::scalar_< bd::arithmetic_<A0> >
                          )
  {
    BOOST_FORCEINLINE A0 operator() (const saturated_tag &
                                    , A0 a0, A0 a1) const BOOST_NOEXCEPT
    {
      A0 tmp = dist(a0, a1);
      return tmp < Zero<A0>() ? Valmax<A0>() : tmp;
    }
  };

  BOOST_DISPATCH_OVERLOAD ( dist_
                          , (typename A0)
                          , bd::cpu_
                          , bs::saturated_tag
                          , bd::scalar_< bd::unsigned_<A0> >
                          , bd::scalar_< bd::unsigned_<A0> >
                          )
  {
    BOOST_FORCEINLINE A0 operator() (const saturated_tag &
                                    ,  A0 a0, A0 a1) const BOOST_NOEXCEPT
    {
      return dist(a0, a1);
    }
  };

  BOOST_DISPATCH_OVERLOAD ( dist_
                          , (typename A0)
                          , bd::cpu_
                          , bs::saturated_tag
                          , bd::scalar_<bd::floating_<A0> >
                          , bd::scalar_<bd::floating_<A0> >
                          )
  {
    BOOST_FORCEINLINE A0 operator() (const saturated_tag &
                                    ,  A0 a0, A0 a1) const BOOST_NOEXCEPT
    {
      return dist(a0, a1);
    }
  };
} } }

#endif

max_stars_repo_name: SylvainCorlay/pythran

max_issues_repo_name: SylvainCorlay/pythran

max_forks_repo_name: SylvainCorlay/pythran


text:
// Copyright (c) Microsoft Corporation.
// Licensed under the MIT license.

#include "spdk_rdma_queue.hh"

#include <arpa/inet.h>
#include <boost/chrono.hpp>
#include <boost/program_options/options_description.hpp>
#include <boost/program_options/parsers.hpp>
#include <boost/program_options/variables_map.hpp>
#include <cassert>
#include <cassert>
#include <cstdint>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <dmtr/annot.h>
#include <dmtr/cast.h>
#include <dmtr/latency.h>
#include <dmtr/libos.h>
#include <dmtr/sga.h>
#include <iostream>
#include <dmtr/libos/mem.h>
#include <dmtr/libos/raii_guard.hh>
#include <netinet/in.h>
#include <unistd.h>

namespace bpo = boost::program_options;

bool dmtr::spdk_rdma_queue::our_init_flag = false;

dmtr::spdk_rdma_queue::spdk_rdma_queue(int qd, dmtr::io_queue::category_id cid) :
    io_queue(cid, qd)
{
    if (cid == NETWORK_Q)
        net_queue = new rdma_queue(qd);
    else
        file_queue = new spdk_queue(qd);

}

int dmtr::spdk_rdma_queue::init_spdk_rdma(int argc, char *argv[]) {
    DMTR_TRUE(ERANGE, argc >= 0);
    if (argc > 0) {
        DMTR_NOTNULL(EINVAL, argv);
    }
    DMTR_TRUE(EPERM, !our_init_flag);

    std::string config_path;
    bpo::options_description desc("Allowed options");
    desc.add_options()
            ("help", "display usage information")
        ("config-path,c", bpo::value<std::string>(&config_path)->default_value("./config.yaml"), "specify configuration file");

    bpo::variables_map vm;
    bpo::store(bpo::command_line_parser(argc, argv).options(desc).allow_unregistered().run(), vm);
    bpo::notify(vm);

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

    if (access(config_path.c_str(), R_OK) == -1) {
        std::cerr << "Unable to find config file at `" << config_path << "`." << std::endl;
        return ENOENT;
    }
    YAML::Node config = YAML::LoadFile(config_path);

    if (getenv("SPDK_OFF") == NULL) {
        struct spdk_env_opts opts;
        spdk_env_opts_init(&opts);

        // init SPDK
        DMTR_OK(spdk_queue::init_spdk(config, &opts));
    } else {
        printf("Turning SPDK off\n");
    }

    DMTR_OK(rdma_queue::init_rdma());
    our_init_flag = true;
    return 0;
}

int dmtr::spdk_rdma_queue::new_net_object(std::unique_ptr<io_queue> &q_out, int qd) {
    q_out = NULL;
    DMTR_TRUE(EPERM, our_init_flag);
    q_out = std::unique_ptr<io_queue>(new spdk_rdma_queue(qd, NETWORK_Q));
    DMTR_NOTNULL(ENOMEM, q_out);
    return 0;
}

int dmtr::spdk_rdma_queue::new_file_object(std::unique_ptr<io_queue> &q_out, int qd) {
    q_out = NULL;
    DMTR_TRUE(EPERM, our_init_flag);
    q_out = std::unique_ptr<io_queue>(new spdk_rdma_queue(qd, FILE_Q));
    DMTR_NOTNULL(ENOMEM, q_out);
    return 0;
}

int dmtr::spdk_rdma_queue::socket(int domain, int type, int protocol) {
    DMTR_TRUE(EPERM, our_init_flag);
    DMTR_TRUE(NETWORK_Q, my_cid);
    DMTR_NOTNULL(EINVAL, net_queue);
    return net_queue->socket(domain, type, protocol);
}

int
dmtr::spdk_rdma_queue::getsockname(struct sockaddr * const saddr, socklen_t * const size) {
    DMTR_TRUE(EPERM, our_init_flag);
    DMTR_TRUE(NETWORK_Q, my_cid);
    DMTR_NOTNULL(EINVAL, net_queue);
    return net_queue->getsockname(saddr, size);
}

int dmtr::spdk_rdma_queue::accept(std::unique_ptr<io_queue> &q_out, dmtr_qtoken_t qt, int new_qd)
{
    DMTR_TRUE(EPERM, our_init_flag);
    DMTR_TRUE(NETWORK_Q, my_cid);
    DMTR_NOTNULL(EINVAL, net_queue);
    return net_queue->accept(q_out, qt, new_qd);
}

int dmtr::spdk_rdma_queue::listen(int backlog)
{
    DMTR_TRUE(EPERM, our_init_flag);
    DMTR_TRUE(NETWORK_Q, my_cid);
    DMTR_NOTNULL(EINVAL, net_queue);
    return net_queue->listen(backlog);
}

int dmtr::spdk_rdma_queue::bind(const struct sockaddr * const saddr, socklen_t size) {
    DMTR_TRUE(EPERM, our_init_flag);
    DMTR_TRUE(NETWORK_Q, my_cid);
    DMTR_NOTNULL(EINVAL, net_queue);
    return net_queue->bind(saddr, size);
}

int dmtr::spdk_rdma_queue::connect(dmtr_qtoken_t qt, const struct sockaddr * const saddr, socklen_t size) {
    DMTR_TRUE(EPERM, our_init_flag);
    DMTR_TRUE(NETWORK_Q, my_cid);
    DMTR_NOTNULL(EINVAL, net_queue);
    return net_queue->connect(qt, saddr, size);
}

int dmtr::spdk_rdma_queue::close() {
    int ret;
    if (my_cid == NETWORK_Q) {
        ret = net_queue->close();
        delete net_queue;
    } else {
        ret = file_queue->close();
        delete file_queue;
    }
    return ret;
}

int dmtr::spdk_rdma_queue::open(const char *pathname, int flags)
{
    DMTR_TRUE(EPERM, our_init_flag);
    DMTR_TRUE(FILE_Q, my_cid);
    DMTR_NOTNULL(EINVAL, file_queue);
    //TODO(ashmrtnz): Yay for only supporing a single file, so we do nothing! If
    // we choose to support multiple files we will need to so some sort of
    // lookup or something here.
    // TODO(ashmrtnz): O_TRUNC?
    file_queue->open(pathname, flags);
    return 0;
}

int dmtr::spdk_rdma_queue::open2(const char *pathname, int flags, mode_t mode)
{
    DMTR_TRUE(EPERM, our_init_flag);
    DMTR_TRUE(FILE_Q, my_cid);
    DMTR_NOTNULL(EINVAL, file_queue);
    //TODO(ashmrtnz): Yay for only supporing a single file, so we do nothing! If
    // we choose to support multiple files we will need to so some sort of
    // lookup or something here.
    // TODO(ashmrtnz): O_TRUNC?
    file_queue->open2(pathname, flags, mode);
    return 0;
}

int dmtr::spdk_rdma_queue::creat(const char *pathname, mode_t mode)
{
    DMTR_TRUE(EPERM, our_init_flag);
    DMTR_TRUE(FILE_Q, my_cid);
    DMTR_NOTNULL(EINVAL, file_queue);
    //TODO(ashmrtnz): Yay for only supporing a single file, so we do nothing! If
    // we choose to support multiple files we will need to so some sort of
    // lookup or something here.
    // TODO(ashmrtnz): O_TRUNC?
    file_queue->creat(pathname, mode);
    return 0;
}

int dmtr::spdk_rdma_queue::push(dmtr_qtoken_t qt, const dmtr_sgarray_t &sga)
{
    DMTR_TRUE(EPERM, our_init_flag);
    if (my_cid == NETWORK_Q) return net_queue->push(qt, sga);
    else return file_queue->push(qt, sga);
}

int dmtr::spdk_rdma_queue::pop(dmtr_qtoken_t qt) {
    DMTR_TRUE(EPERM, our_init_flag);
    if (my_cid == NETWORK_Q) return net_queue->pop(qt);
    else return file_queue->pop(qt);
}

int dmtr::spdk_rdma_queue::poll(dmtr_qresult_t &qr_out, dmtr_qtoken_t qt)
{
    DMTR_TRUE(EPERM, our_init_flag);
    if (my_cid == NETWORK_Q) return net_queue->poll(qr_out, qt);
    else return file_queue->poll(qr_out, qt);
}

ext: cc

max_stars_repo_name: seemk/demikernel

max_issues_repo_name: seemk/demikernel

max_forks_repo_name: seemk/demikernel


text:
// Boost.Geometry

// Copyright (c) 2021, Oracle and/or its affiliates.
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle

// Licensed under the Boost Software License version 1.0.
// http://www.boost.org/users/license.html

#include <boost/test/included/test_exec_monitor.hpp>
#include <boost/test/impl/execution_monitor.ipp>

#include <algorithm>
#include <vector>
#include <map>

#include <boost/geometry/index/detail/minmax_heap.hpp>
#include <boost/geometry/index/detail/maxmin_heap.hpp>

using namespace boost::geometry::index::detail;

struct noncopyable
{
    noncopyable(int i_) : i(i_) {}
    noncopyable(noncopyable const&) = delete;
    noncopyable& operator=(noncopyable const&) = delete;
    noncopyable(noncopyable&&) = default;
    noncopyable& operator=(noncopyable&&) = default;
    bool operator<(noncopyable const& other) const { return i < other.i; }
    operator int() const { return i; }
    int i;
};

template <typename T>
struct minmax_default
{
    minmax_default() = default;
    minmax_default(std::vector<int> const& vec)
        : heap(vec.begin(), vec.end())
    {
        make_minmax_heap(heap.begin(), heap.end());
    }
    T const& top() const
    {
        return heap[0];
    }
    T const& bottom() const
    {
        return bottom_minmax_heap(heap.begin(), heap.end());
    }
    void push(int i)
    {
        heap.push_back(T(i));
        push_minmax_heap(heap.begin(), heap.end());
    }
    void pop_top()
    {
        pop_top_minmax_heap(heap.begin(), heap.end());
        heap.pop_back();
    }
    void pop_bottom()
    {
        pop_bottom_minmax_heap(heap.begin(), heap.end());
        heap.pop_back();
    }
    bool is_heap() const
    {
        return is_minmax_heap(heap.begin(), heap.end());
    }
    bool empty() const
    {
        return heap.empty();
    }
    std::vector<T> heap;
};

template <typename T>
struct minmax_less
{
    minmax_less() = default;
    minmax_less(std::vector<int> const& vec)
        : heap(vec.begin(), vec.end())
    {
        make_minmax_heap(heap.begin(), heap.end(), std::less<>());
    }
    T const& top() const
    {
        return heap[0];
    }
    T const& bottom() const
    {
        return bottom_minmax_heap(heap.begin(), heap.end(), std::less<>());
    }
    void push(int i)
    {
        heap.push_back(T(i));
        push_minmax_heap(heap.begin(), heap.end(), std::less<>());
    }
    void pop_top()
    {
        pop_top_minmax_heap(heap.begin(), heap.end(), std::less<>());
        heap.pop_back();
    }
    void pop_bottom()
    {
        pop_bottom_minmax_heap(heap.begin(), heap.end(), std::less<>());
        heap.pop_back();
    }
    bool is_heap() const
    {
        return is_minmax_heap(heap.begin(), heap.end(), std::less<>());
    }
    bool empty() const
    {
        return heap.empty();
    }
    std::vector<T> heap;
};

template <typename T>
struct maxmin_greater
{
    maxmin_greater() = default;
    maxmin_greater(std::vector<int> const& vec)
        : heap(vec.begin(), vec.end())
    {
        make_maxmin_heap(heap.begin(), heap.end(), std::greater<>());
    }
    T const& top() const
    {
        return heap[0];
    }
    T const& bottom() const
    {
        return bottom_maxmin_heap(heap.begin(), heap.end(), std::greater<>());
    }
    void push(int i)
    {
        heap.push_back(T(i));
        push_maxmin_heap(heap.begin(), heap.end(), std::greater<>());
    }
    void pop_top()
    {
        pop_top_maxmin_heap(heap.begin(), heap.end(), std::greater<>());
        heap.pop_back();
    }
    void pop_bottom()
    {
        pop_bottom_maxmin_heap(heap.begin(), heap.end(), std::greater<>());
        heap.pop_back();
    }
    bool is_heap() const
    {
        return is_maxmin_heap(heap.begin(), heap.end(), std::greater<>());
    }
    bool empty() const
    {
        return heap.empty();
    }
    std::vector<T> heap;
};

template <typename T>
struct maxmin_default_switch
{
    maxmin_default_switch() = default;
    maxmin_default_switch(std::vector<int> const& vec)
        : heap(vec.begin(), vec.end())
    {
        make_maxmin_heap(heap.begin(), heap.end());
    }
    T const& bottom() const
    {
        return heap[0];
    }
    T const& top() const
    {
        return bottom_maxmin_heap(heap.begin(), heap.end());
    }
    void push(int i)
    {
        heap.push_back(T(i));
        push_maxmin_heap(heap.begin(), heap.end());
    }
    void pop_top()
    {
        pop_bottom_maxmin_heap(heap.begin(), heap.end());
        heap.pop_back();
    }
    void pop_bottom()
    {
        pop_top_maxmin_heap(heap.begin(), heap.end());
        heap.pop_back();
    }
    bool is_heap() const
    {
        return is_maxmin_heap(heap.begin(), heap.end());
    }
    bool empty() const
    {
        return heap.empty();
    }
    std::vector<T> heap;
};

template <typename Heap>
void test()
{
    std::vector<int> vec;
    int const n = 20;
    for (int i = 0; i < n; ++i)
    {
        vec.push_back(rand() % n);
    }

    {
        std::map<int, int> map;
        Heap heap;
        for (int i : vec)
        {
            heap.push(i);
            BOOST_CHECK(heap.is_heap());

            map[i]++;
            BOOST_CHECK_EQUAL(heap.top(), map.begin()->first);
            BOOST_CHECK_EQUAL(heap.bottom(), (--map.end())->first);
        }

        while (! heap.empty())
        {
            int i = heap.top();
            BOOST_CHECK_EQUAL(i, map.begin()->first);
            BOOST_CHECK_EQUAL(heap.bottom(), (--map.end())->first);
            BOOST_CHECK(map[i] > 0);
            map[i]--;
            if (map[i] <= 0)
                map.erase(i);

            heap.pop_top();
            BOOST_CHECK(heap.is_heap());
        }

        BOOST_CHECK(map.empty());
    }

    {
        Heap heap(vec);
        BOOST_CHECK(heap.is_heap());

        std::map<int, int> map;
        for (int i : vec)
            map[i]++;
        BOOST_CHECK_EQUAL(heap.top(), map.begin()->first);
        BOOST_CHECK_EQUAL(heap.bottom(), (--map.end())->first);

        while (! heap.empty())
        {
            int i = heap.bottom();
            BOOST_CHECK_EQUAL(heap.top(), map.begin()->first);
            BOOST_CHECK_EQUAL(i, (--map.end())->first);
            BOOST_CHECK(map[i] > 0);
            map[i]--;
            if (map[i] <= 0)
                map.erase(i);

            heap.pop_bottom();
            BOOST_CHECK(heap.is_heap());
        }

        BOOST_CHECK(map.empty());
    }
}

int test_main(int, char* [])
{
    test<minmax_default<int>>();
    test<minmax_default<noncopyable>>();
    test<minmax_less<int>>();
    test<maxmin_greater<int>>();
    test<maxmin_default_switch<int>>();

    return 0;
}

ext: cpp


text:
//  (C) Copyright Dave Abrahams, Steve Cleary, Beman Dawes, Howard
//  Hinnant & John Maddock 2000.
//  Use, modification and distribution are subject to the Boost Software License,
//  Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
//  http://www.boost.org/LICENSE_1_0.txt).
//
//  See http://www.boost.org/libs/type_traits for most recent version including documentation.

#ifndef BOOST_TT_DETAIL_CV_TRAITS_IMPL_HPP_INCLUDED
#define BOOST_TT_DETAIL_CV_TRAITS_IMPL_HPP_INCLUDED

#include <boost/config.hpp>
#include <boost/detail/workaround.hpp>

#ifndef BOOST_NO_TEMPLATE_PARTIAL_SPECIALIZATION

// implementation helper:

#if !(BOOST_WORKAROUND(__GNUC__,== 3) && BOOST_WORKAROUND(__GNUC_MINOR__, <= 2))
namespace boost {
namespace detail {
#else
#include <boost/type_traits/detail/yes_no_type.hpp>
namespace boost {
namespace type_traits {
namespace gcc8503 {
#endif

template <typename T> struct cv_traits_imp {};

template <typename T>
struct cv_traits_imp<T*>
{
    BOOST_STATIC_CONSTANT(bool, is_const = false);
    BOOST_STATIC_CONSTANT(bool, is_volatile = false);
    typedef T unqualified_type;
};

template <typename T>
struct cv_traits_imp<const T*>
{
    BOOST_STATIC_CONSTANT(bool, is_const = true);
    BOOST_STATIC_CONSTANT(bool, is_volatile = false);
    typedef T unqualified_type;
};

template <typename T>
struct cv_traits_imp<volatile T*>
{
    BOOST_STATIC_CONSTANT(bool, is_const = false);
    BOOST_STATIC_CONSTANT(bool, is_volatile = true);
    typedef T unqualified_type;
};

template <typename T>
struct cv_traits_imp<const volatile T*>
{
    BOOST_STATIC_CONSTANT(bool, is_const = true);
    BOOST_STATIC_CONSTANT(bool, is_volatile = true);
    typedef T unqualified_type;
};

#if BOOST_WORKAROUND(__GNUC__,== 3) && BOOST_WORKAROUND(__GNUC_MINOR__, <= 2)
// We have to exclude function pointers
// (see http://gcc.gnu.org/bugzilla/show_bug.cgi?8503)
yes_type mini_funcptr_tester(...);
no_type  mini_funcptr_tester(const volatile void*);

} // namespace gcc8503
} // namespace type_traits

namespace detail {

// Use the implementation above for non function pointers
template <typename T, unsigned Select
  = (unsigned)sizeof(::boost::type_traits::gcc8503::mini_funcptr_tester((T)0)) >
struct cv_traits_imp : public ::boost::type_traits::gcc8503::cv_traits_imp<T> { };

// Functions are never cv-qualified
template <typename T> struct cv_traits_imp<T*,1>
{
    BOOST_STATIC_CONSTANT(bool, is_const = false);
    BOOST_STATIC_CONSTANT(bool, is_volatile = false);
    typedef T unqualified_type;
};

#endif

} // namespace detail
} // namespace boost

#endif // BOOST_NO_TEMPLATE_PARTIAL_SPECIALIZATION

#endif // BOOST_TT_DETAIL_CV_TRAITS_IMPL_HPP_INCLUDED

max_stars_repo_name: randolphwong/mcsema

max_issues_repo_name: randolphwong/mcsema

max_forks_repo_name: randolphwong/mcsema


text:
#include "loop.hpp"
#include "render_impl.hpp"

#include <iostream>
#include <stdexcept>
#include <cstdint>

#include <boost/program_options.hpp>

namespace options = boost::program_options;
using namespace gopnik;

int
main(int argc, char *argv[]) {
	try {
		// Parse args
		std::string stylesheet;
		std::vector<std::string> fontsPath;
		std::string pluginsPath;
		unsigned tileSize = 256;
		int bufferSize = -1;
		double scaleFactor = 1.0;
		std::string imageFormat = "png8";

		options::options_description desc (std::string (argv[0]).append(" options"));
		desc.add_options()
			("help", "Display this message")
			("stylesheet", options::value<std::string>(), "stylesheet")
			("tileSize", options::value<unsigned>(), "tileSize")
			("bufferSize", options::value<int>(), "tileSize")
			("fontsPath", options::value<std::vector<std::string>>(), "fontsPath")
			("pluginsPath", options::value<std::string>(), "pluginsPath")
			("scaleFactor", options::value<double>(), "scaleFactor")
			("imageFormat", options::value<std::string>(), "imageFormat")
		;
		options::variables_map args;
		options::store (options::command_line_parser (argc, argv).options (desc)
			.style (options::command_line_style::default_style |
				options::command_line_style::allow_long_disguise)
			.run (), args);
		options::notify (args);

		if (args.count("h")) {
			std::cerr << desc << std::endl;
			return 1;
		}
		if (args.count("stylesheet")) {
			stylesheet = args["stylesheet"].as<std::string>();
		} else {
			throw std::runtime_error("stylesheet required");
		}
		if (args.count("tileSize")) {
			tileSize = args["tileSize"].as<unsigned>();
		}
		if (args.count("bufferSize")) {
			bufferSize = args["bufferSize"].as<int>();
		}
		if (args.count("fontsPath")) {
			fontsPath = args["fontsPath"].as<std::vector<std::string>>();
		}
		if (args.count("pluginsPath")) {
			pluginsPath = args["pluginsPath"].as<std::string>();
		}
		if (args.count("scaleFactor")) {
			scaleFactor = args["scaleFactor"].as<double>();
		}
		if (args.count("imageFormat")) {
			imageFormat = args["imageFormat"].as<std::string>();
		}

		// Set up render
		RenderImpl render{stylesheet, fontsPath, pluginsPath,
			tileSize, bufferSize, scaleFactor, imageFormat};

		Loop loop{0 /* stdin */, 1 /* stdout */, render};

		// Start event loop
		loop.Run();
	}
	catch(std::exception const& e) {
		std::cerr << "Exception: " << e.what() << std::endl;
		return 1;
	}
	catch(...) {
		std::cerr << "Unknown exception" << std::endl;
		return 2;
	}
	return 0;
}

max_stars_repo_name: sputnik-maps/gopnik

max_issues_repo_name: sputnik-maps/gopnik

max_forks_repo_name: sputnik-maps/gopnik


text:
// Copyright (c) 2021 Bitcoin Association
// Distributed under the Open BSV software license, see the accompanying file LICENSE.

#include "validation.h"

#include "block_index_store.h"
#include "config.h"
#include "frozentxo.h"
#include "frozentxo_db.h"
#include "pow.h"
#include "test/test_bitcoin.h"
#include <boost/test/unit_test.hpp>
#include <vector>
#include <cstdint>

BOOST_FIXTURE_TEST_SUITE(frozentxo, TestingSetup)

namespace {

void validateInputs(
    const CTransaction& tx,
    const CCoinsViewCache& inputs,
    CFrozenTXOCheck& frozenTXOCheck,
    const std::string& errorReason)
{
    CValidationState state;

    bool valid =
        Consensus::CheckTxInputs(
            tx,
            state,
            inputs,
            0, // ignoring spend height as the test should not reach that code
            frozenTXOCheck);

    BOOST_CHECK( !valid );
    BOOST_CHECK_EQUAL( state.GetRejectReason(), errorReason );
}

void populateBlocks( BlockIndexStore& blockIndexStore, CChain& blocks, std::size_t count )
{
    CBlockIndex* prev{};
    uint256 prevHash;
    for(std::size_t i = 0; i < count; ++i)
    {
        CBlockHeader header;
        header.nTime = i;
        header.hashPrevBlock = prevHash;
        header.nBits = GetNextWorkRequired( prev, &header, GlobalConfig::GetConfig() );

        prev = blockIndexStore.Insert( header );
        blocks.SetTip( prev );
        prevHash = prev->GetBlockHash();
    }
}

} // anonymous namespace

BOOST_AUTO_TEST_CASE(negative_Consensus_CheckTxInputs)
{
    CCoinsViewEmpty dummy;
    CCoinsViewCache inputs{dummy};

    CMutableTransaction txTemplate;
    txTemplate.vin.resize(1);
    txTemplate.vin[0].prevout = COutPoint{InsecureRand256(), 3};
    txTemplate.vin[0].scriptSig.resize(10);
    txTemplate.vout.resize(1);
    txTemplate.vout[0].nValue = Amount(42);
    CTransaction tx{txTemplate};

    inputs.AddCoin(
        txTemplate.vin[0].prevout,
        CoinWithScript::MakeOwning(CTxOut{Amount{3}, {}}, 1, false),
        false,
        0);

    CChain blocks;
    BlockIndexStore blockIndexStore;
    populateBlocks( blockIndexStore, blocks, 3 );

    std::vector<CFrozenTXOCheck> frozenTXOCheckTransaction; // NOTE: element index is block height at which check is performed
    std::vector<CFrozenTXOCheck> frozenTXOCheckBlock;
    uint256 prevHash;

    for(std::int32_t h=0; h<3; ++h)
    {
        frozenTXOCheckTransaction.emplace_back(
            h,
            "test transaction",
            prevHash
        );

        frozenTXOCheckBlock.emplace_back( *blocks[h] );

        prevHash = blocks[h]->GetBlockHash();
    }

    // enforce policy level freeze
    CFrozenTXODB::Instance().FreezeTXOPolicyOnly(txTemplate.vin[0].prevout);
    validateInputs(tx, inputs, frozenTXOCheckTransaction[0], "bad-txns-inputs-frozen"); // frozen at all heights for txs not yet in block
    validateInputs(tx, inputs, frozenTXOCheckTransaction[1], "bad-txns-inputs-frozen");
    validateInputs(tx, inputs, frozenTXOCheckTransaction[2], "bad-txns-inputs-frozen");
    validateInputs(tx, inputs, frozenTXOCheckBlock[0], "bad-txns-in-belowout"); // not frozen at any height for txs in block
    validateInputs(tx, inputs, frozenTXOCheckBlock[1], "bad-txns-in-belowout");
    validateInputs(tx, inputs, frozenTXOCheckBlock[2], "bad-txns-in-belowout");

    // start enforcing consensus level freeze at height 1
    CFrozenTXODB::Instance().FreezeTXOConsensus(txTemplate.vin[0].prevout, {{1}}, false);
    validateInputs(tx, inputs, frozenTXOCheckTransaction[0], "bad-txns-inputs-frozen"); // frozen at all heights for txs not yet in block
    validateInputs(tx, inputs, frozenTXOCheckTransaction[1], "bad-txns-inputs-frozen");
    validateInputs(tx, inputs, frozenTXOCheckTransaction[2], "bad-txns-inputs-frozen");
    validateInputs(tx, inputs, frozenTXOCheckBlock[0], "bad-txns-in-belowout"); // not frozen for txs in block at height 0
    validateInputs(tx, inputs, frozenTXOCheckBlock[1], "bad-txns-inputs-frozen"); // frozen for txs in block at height 1
    validateInputs(tx, inputs, frozenTXOCheckBlock[2], "bad-txns-inputs-frozen"); // frozen for txs in block at height 2

    // stop enforcing consensus level freeze at height 2 but keep enforcing policy level
    CFrozenTXODB::Instance().FreezeTXOConsensus(txTemplate.vin[0].prevout, {{1,2}}, false);
    validateInputs(tx, inputs, frozenTXOCheckTransaction[0], "bad-txns-inputs-frozen"); // frozen at all heights for txs not yet in block
    validateInputs(tx, inputs, frozenTXOCheckTransaction[1], "bad-txns-inputs-frozen");
    validateInputs(tx, inputs, frozenTXOCheckTransaction[2], "bad-txns-inputs-frozen");
    validateInputs(tx, inputs, frozenTXOCheckBlock[0], "bad-txns-in-belowout"); // not frozen for txs in block at height 0
    validateInputs(tx, inputs, frozenTXOCheckBlock[1], "bad-txns-inputs-frozen"); // frozen for txs in block at height 1
    validateInputs(tx, inputs, frozenTXOCheckBlock[2], "bad-txns-in-belowout"); // not frozen for txs in block at height 2

    // stop enforcing consensus and policy level freeze at height 2
    CFrozenTXODB::Instance().FreezeTXOConsensus(txTemplate.vin[0].prevout, {{1,2}}, true);
    validateInputs(tx, inputs, frozenTXOCheckTransaction[0], "bad-txns-inputs-frozen"); // frozen only at heights before stop height for txs not yet in block
    validateInputs(tx, inputs, frozenTXOCheckTransaction[1], "bad-txns-inputs-frozen");
    validateInputs(tx, inputs, frozenTXOCheckTransaction[2], "bad-txns-in-belowout");
    validateInputs(tx, inputs, frozenTXOCheckBlock[0], "bad-txns-in-belowout"); // not frozen for txs in block at height 0
    validateInputs(tx, inputs, frozenTXOCheckBlock[1], "bad-txns-inputs-frozen"); // frozen for txs in block at height 1
    validateInputs(tx, inputs, frozenTXOCheckBlock[2], "bad-txns-in-belowout"); // not frozen for txs in block at height 2
}

BOOST_AUTO_TEST_SUITE_END()

max_stars_repo_name: ZacharyPainter/bitcoin-sv

max_issues_repo_name: wjxee/bitcoin-sv

max_forks_repo_name: wjxee/bitcoin-sv


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

// Copyright (c) 2020 Barend Gehrels, Amsterdam, the Netherlands.

// This file was modified by Oracle on 2020.
// Modifications copyright (c) 2020 Oracle and/or its affiliates.
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle

// Use, modification and distribution is subject to the Boost Software License,
// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
// htt//www.boost.org/LICENSE_1_0.txt)

#include "geometry_test_common.hpp"

#include <boost/geometry/algorithms/area.hpp>
#include <boost/geometry/algorithms/detail/overlay/turn_info.hpp>
#include <boost/geometry/algorithms/detail/buffer/piece_border.hpp>
#include <boost/geometry/util/range.hpp>

#include <boost/geometry/io/wkt/wkt.hpp>
#include <boost/geometry/geometries/geometries.hpp>

#include <boost/geometry/strategies/comparable_distance_result.hpp>
#include <boost/geometry/strategies/cartesian/distance_pythagoras.hpp>
#include <boost/geometry/strategies/cartesian/distance_projected_point.hpp>

// TEMP
#include <boost/geometry/strategies/cartesian.hpp>

#if defined(TEST_WITH_SVG)
#include <boost/geometry/io/svg/svg_mapper.hpp>
#include <boost/geometry/algorithms/detail/buffer/buffer_box.hpp>
#endif

static const std::string rectangle_offsetted = "POLYGON((0 0,0 1,0 2,0 3,1 3,2 3,3 3,3 2,3 1,3 0,2 0,1 0,0 0))";
static const std::string rectangle_original = "POLYGON((1 1,1 2,2 2,2 1,1 1))";

static const std::string diamond_offsetted = "POLYGON((3 0,2 1,1 2,0 3,1 4,2 5,3 6,4 5,5 4,6 3,5 2,4 1,3 0))";
static const std::string diamond_original = "POLYGON((3 2,2 3,3 4,4 3,3 2))";

template <typename Geometry>
inline std::string as_wkt(Geometry const& geometry)
{
    std::ostringstream out;
    out << bg::wkt(geometry);
    return out.str();
}

// Return the border. It needs to pass the ring because its address
// is preserved in the border. It also gets the original.
template <typename Border, typename Ring>
Border setup_piece_border(Ring& ring, Ring& original,
                          std::string const& offsetted_wkt,
                          std::string const& original_wkt,
                          char piece)
{
    // The test setup (here fore case "rectangle"; "diamond" has same layout).
    // Points on offsetted ring are numbered, points on original like this [0].
    // a, b, c: test pieces.
    //
    //           4       5
    //   3 +-----+-------+-----+ 6
    //     |     |       |     |
    //     |     |   a   |  b  |
    //     |     |       |     |
    //   2 +----[1]-----[2]----+ 7
    //     |     |       |     |
    //     |     |  ori  |  c  |
    //     |     |       |     |
    //   1 +----[0,4]---[3]----+ 8
    //     |     |       |     |
    //     |     |       |     |
    //     |     |       |     |
    //   0 +-----+-------+-----+ 9
    //    12    11       10

    bg::read_wkt(original_wkt, original);
    bg::read_wkt(offsetted_wkt, ring);

    Border border;

    switch(piece)
    {
        case 'a' :
            border.set_offsetted(ring, 4, 6);
            border.add_original_point(original[2]);
            border.add_original_point(original[1]);
            break;
        case 'b' :
            border.set_offsetted(ring, 5, 8);
            border.add_original_point(original[2]);
            break;
        case 'c' :
            border.set_offsetted(ring, 7, 9);
            border.add_original_point(original[3]);
            border.add_original_point(original[2]);
            break;
        default :
            return border;
    }

    typedef typename bg::point_type<Ring>::type point_type;

    bg::strategies::relate::cartesian<> strategy;

    border.get_properties_of_border(false, point_type(), strategy);

    border.get_properties_of_offsetted_ring_part(strategy);

    return border;
}

template <typename Point>
void test_rectangle_properties()
{
    typedef bg::model::ring<Point> ring_type;
    typedef bg::detail::buffer::piece_border<ring_type, Point> border_type;

    border_type empty;
    BOOST_CHECK_MESSAGE(empty.ring_or_original_empty(), "piece should be empty");

    ring_type offsetted, original;
    border_type border = setup_piece_border<border_type>(offsetted, original,
            rectangle_offsetted, rectangle_original, 'a');

    // Check get_full_ring functionality (used for testing, writing SVG)
    ring_type const outline = border.get_full_ring();
    std::string const outline_wkt = as_wkt(outline);

    std::string const expected = "POLYGON((1 3,2 3,2 2,1 2,1 3))";
    BOOST_CHECK_MESSAGE(outline_wkt == expected,
                        "expected: " << expected
                        << " detected: " << outline_wkt);

    BOOST_CHECK_MESSAGE(! border.ring_or_original_empty(),
                        "piece should not be empty");

    // Check border-properties functionality (envelope, radius)
    auto const area = bg::area(border.m_envelope);
    BOOST_CHECK_MESSAGE(area > 1.0 && area < 1.01,
                        "detected: " << area);

    // Check offsetted-ring-properties functionality (convexity, monotonicity)
    BOOST_CHECK_MESSAGE(border.m_is_convex == true,
                        "piece should be convex");
    BOOST_CHECK_MESSAGE(border.m_is_monotonic_increasing == true,
                        "piece should be monotonically increasing");
    BOOST_CHECK_MESSAGE(border.m_is_monotonic_decreasing == false,
                        "piece NOT should be monotonically decreasing");
}

template <typename Point, typename Border, typename Mapper>
void test_point(std::string const& wkt, bool expected, Border const& border,
                Mapper& mapper, std::string const& color)
{
    typename Border::state_type state;
    Point point;
    bg::read_wkt(wkt, point);
    border.point_on_piece(point, false, false, state);
    bool const on_piece = state.code() == 1;
    BOOST_CHECK_MESSAGE(on_piece == expected,
                        " expected: " << expected
                        << " detected: " << on_piece
                        << " wkt: " << wkt);

#ifdef TEST_WITH_SVG
    std::string style = "fill:" + color + ";stroke:rgb(0,0,0);stroke-width:1";
    mapper.map(point, style);

    // Mark on-piece as T or F
    std::ostringstream out;
    out << (on_piece ? "T" : "F");
    mapper.text(point, out.str(), "fill:rgb(0,0,0);font-family='Arial';font-size:9px;", 10, -10);
#else
    boost::ignore_unused(mapper, color);
#endif
}

#ifdef TEST_WITH_SVG
template <typename Mapper, typename Ring, typename Border>
void start_svg(Mapper& mapper, Ring const& original, Ring const& offsetted,
               Border const& border)
{
    typedef typename bg::point_type<Ring>::type point_type;
    typedef bg::model::box<point_type> box_type;
    box_type box = border.m_envelope;
    bg::detail::buffer::buffer_box(box, 1.0, box);
    mapper.add(box);
    mapper.add(offsetted);
    mapper.map(offsetted,
               "fill-opacity:0.5;fill:rgb(255,255,128);stroke:none");
    mapper.map(original,
               "fill-opacity:0.5;fill:rgb(51,51,153);stroke:none");
    mapper.map(border.get_full_ring(),
               "fill-opacity:0.5;fill:rgb(153,204,0);"
               "stroke:rgb(153,204,0);stroke-width:1");
}
#endif

template <typename Point>
void test_rectangle_point_on_piece_a()
{
    typedef bg::model::ring<Point> ring_type;
    typedef bg::detail::buffer::piece_border<ring_type, Point> border_type;

    ring_type offsetted, original;
    border_type border = setup_piece_border<border_type>(offsetted, original,
            rectangle_offsetted, rectangle_original, 'a');

#ifdef TEST_WITH_SVG
    std::ofstream svg("/tmp/border_rectangle_a.svg");
    bg::svg_mapper<Point> mapper(svg, 500, 500);
    start_svg(mapper, original, offsetted, border);
#else
    int mapper = 0;
#endif

    // Points inside
    test_point<Point>("POINT(1.5 2.5)", true, border, mapper, "red");

    // Points outside
    test_point<Point>("POINT(0.5 2.5)", false, border, mapper, "green");
    test_point<Point>("POINT(2.5 2.5)", false, border, mapper, "green");
    test_point<Point>("POINT(1.5 1.5)", false, border, mapper, "green");
    test_point<Point>("POINT(1.5 3.5)", false, border, mapper, "green");

    // Points on the original (should be INSIDE)
    test_point<Point>("POINT(1.0 2.0)", true, border, mapper, "orange");
    test_point<Point>("POINT(1.5 2.0)", true, border, mapper, "orange");
    test_point<Point>("POINT(2.0 2.0)", true, border, mapper, "orange");

    // Points on the offsetted ring (should be OUTSIDE)
    test_point<Point>("POINT(1.0 3.0)", false, border, mapper, "blue");
    test_point<Point>("POINT(1.5 3.0)", false, border, mapper, "blue");
    test_point<Point>("POINT(2.0 3.0)", false, border, mapper, "blue");

    // Points on between original and offsetted ring (should be INSIDE)
    test_point<Point>("POINT(1.0 2.5)", true, border, mapper, "cyan");
    test_point<Point>("POINT(2.0 2.5)", true, border, mapper, "cyan");
}

template <typename Point>
void test_rectangle_point_on_piece_c()
{
    typedef bg::model::ring<Point> ring_type;
    typedef bg::detail::buffer::piece_border<ring_type, Point> border_type;

    ring_type offsetted, original;
    border_type border = setup_piece_border<border_type>(offsetted, original,
            rectangle_offsetted, rectangle_original, 'c');

#ifdef TEST_WITH_SVG
    std::ofstream svg("/tmp/border_rectangle_c.svg");
    bg::svg_mapper<Point> mapper(svg, 500, 500);
    start_svg(mapper, original, offsetted, border);
#else
    int mapper = 0;
#endif

    // Piece labeled as 'c' : from (2,1) to (3,2)

    // Points inside
    test_point<Point>("POINT(2.5 1.5)", true, border, mapper, "red");

    // Points outside
    test_point<Point>("POINT(1.5 1.5)", false, border, mapper, "green");
    test_point<Point>("POINT(3.5 1.5)", false, border, mapper, "green");
    test_point<Point>("POINT(2.5 0.5)", false, border, mapper, "green");
    test_point<Point>("POINT(2.5 2.5)", false, border, mapper, "green");

    // Points on the original (should be INSIDE)
    test_point<Point>("POINT(2.0 1.0)", true, border, mapper, "orange");
    test_point<Point>("POINT(2.0 1.5)", true, border, mapper, "orange");
    test_point<Point>("POINT(2.0 2.0)", true, border, mapper, "orange");

    // Points on the offsetted ring (should be OUTSIDE)
    test_point<Point>("POINT(3.0 1.0)", false, border, mapper, "blue");
    test_point<Point>("POINT(3.0 1.5)", false, border, mapper, "blue");
    test_point<Point>("POINT(3.0 2.0)", false, border, mapper, "blue");

    // Points on between original and offsetted ring (should be INSIDE)
    test_point<Point>("POINT(2.5 2.0)", true, border, mapper, "cyan");
    test_point<Point>("POINT(2.5 1.0)", true, border, mapper, "cyan");
}

template <typename Point>
void test_diamond_point_on_piece_a()
{
    typedef bg::model::ring<Point> ring_type;
    typedef bg::detail::buffer::piece_border<ring_type, Point> border_type;

    ring_type offsetted, original;
    border_type border = setup_piece_border<border_type>(offsetted, original,
            diamond_offsetted, diamond_original, 'a');

#ifdef TEST_WITH_SVG
    std::ofstream svg("/tmp/border_diamond_a.svg");
    bg::svg_mapper<Point> mapper(svg, 500, 500);
    start_svg(mapper, original, offsetted, border);
#else
    int mapper = 0;
#endif

    // Points inside
    test_point<Point>("POINT(2 4)", true, border, mapper, "red");

    // Points outside
    test_point<Point>("POINT(1 3)", false, border, mapper, "green");
    test_point<Point>("POINT(1 5)", false, border, mapper, "green");
    test_point<Point>("POINT(3 3)", false, border, mapper, "green");
    test_point<Point>("POINT(3 5)", false, border, mapper, "green");

    // Points on the original (should be INSIDE)
    test_point<Point>("POINT(2 3)", true, border, mapper, "orange");
    test_point<Point>("POINT(2.5 3.5)", true, border, mapper, "orange");
    test_point<Point>("POINT(3 4)", true, border, mapper, "orange");

    // Points on the offsetted ring (should be OUTSIDE)
    test_point<Point>("POINT(1 4)", false, border, mapper, "blue");
    test_point<Point>("POINT(1.5 4.5)", false, border, mapper, "blue");
    test_point<Point>("POINT(2 5)", false, border, mapper, "blue");

    // Points on between original and offsetted ring (should be INSIDE)
    test_point<Point>("POINT(1.5 3.5)", true, border, mapper, "cyan");
    test_point<Point>("POINT(2.5 4.5)", true, border, mapper, "cyan");
}

template <typename Point>
void test_diamond_point_on_piece_c()
{
    typedef bg::model::ring<Point> ring_type;
    typedef bg::detail::buffer::piece_border<ring_type, Point> border_type;

    ring_type offsetted, original;
    border_type border = setup_piece_border<border_type>(offsetted, original,
            diamond_offsetted, diamond_original, 'c');

#ifdef TEST_WITH_SVG
    std::ofstream svg("/tmp/border_diamond_c.svg");
    bg::svg_mapper<Point> mapper(svg, 500, 500);
    start_svg(mapper, original, offsetted, border);
#else
    int mapper = 0;
#endif

    // Points inside
    test_point<Point>("POINT(4 4)", true, border, mapper, "red");

    // Points outside
    test_point<Point>("POINT(3 3)", false, border, mapper, "green");
    test_point<Point>("POINT(3 5)", false, border, mapper, "green");
    test_point<Point>("POINT(5 3)", false, border, mapper, "green");
    test_point<Point>("POINT(5 5)", false, border, mapper, "green");

    // Points on the original (should be INSIDE)
    test_point<Point>("POINT(3 4)", true, border, mapper, "orange");
    test_point<Point>("POINT(3.5 3.5)", true, border, mapper, "orange");
    test_point<Point>("POINT(4 3)", true, border, mapper, "orange");

    // Points on the offsetted ring (should be OUTSIDE)
    test_point<Point>("POINT(4 5)", false, border, mapper, "blue");
    test_point<Point>("POINT(4.5 4.5)", false, border, mapper, "blue");
    test_point<Point>("POINT(5 4)", false, border, mapper, "blue");

    // Points on between original and offsetted ring (should be INSIDE)
    test_point<Point>("POINT(3.5 4.5)", true, border, mapper, "cyan");
    test_point<Point>("POINT(4.5 3.5)", true, border, mapper, "cyan");
}

int test_main(int, char* [])
{
    BoostGeometryWriteTestConfiguration();

    typedef bg::model::point<default_test_type, 2, bg::cs::cartesian> point_type;

    test_rectangle_properties<point_type>();

    test_rectangle_point_on_piece_a<point_type>();
    test_rectangle_point_on_piece_c<point_type>();

    test_diamond_point_on_piece_a<point_type>();
    test_diamond_point_on_piece_c<point_type>();

    return 0;
}


text:
//==================================================================================================
/**
  EVE - Expressive Vector Engine
  Copyright : EVE Contributors & Maintainers
  SPDX-License-Identifier: MIT
**/
//==================================================================================================
#include "test.hpp"
#include <eve/module/bessel.hpp>
#include <boost/math/special_functions/bessel.hpp>
#include <boost/math/special_functions/bessel_prime.hpp>
#include <eve/detail/diff_div.hpp>

EVE_TEST_TYPES( "Check return types of sph_bessel_j1"
            , eve::test::simd::ieee_reals
            )
<typename T>(eve::as<T>)
{
  using v_t = eve::element_type_t<T>;
  TTS_EXPR_IS(eve::sph_bessel_j1(T(0)), T);
  TTS_EXPR_IS(eve::sph_bessel_j1(v_t(0)), v_t);
};

 EVE_TEST( "Check behavior of sph_bessel_j1 on wide"
        , eve::test::simd::ieee_reals
        , eve::test::generate(eve::test::randoms(0.0, 5.5),
                              eve::test::randoms(5.5, 9.5),
                              eve::test::randoms(9.5, 60.0))
         )
   <typename T>(T const& a0, T const& a1, T const& a2)
{
  using v_t = eve::element_type_t<T>;
  auto eve__sph_bessel_j1 =  [](auto x) { return eve::sph_bessel_j1(x); };
  auto std__sph_bessel_j1 =  [](auto x)->v_t { return boost::math::sph_bessel(double(1), double(x)); };

  if constexpr( eve::platform::supports_invalids )
  {
    TTS_ULP_EQUAL(eve__sph_bessel_j1(eve::inf(eve::as<v_t>())), eve::zero(eve::as<v_t>()), 0);
    TTS_ULP_EQUAL(eve__sph_bessel_j1(eve::nan(eve::as<v_t>())), eve::nan(eve::as<v_t>()), 0);
    TTS_ULP_EQUAL(eve__sph_bessel_j1(eve::inf(eve::as< T>())),  eve::zero(eve::as< T>()), 0);
    TTS_ULP_EQUAL(eve__sph_bessel_j1(eve::nan(eve::as< T>())), eve::nan(eve::as< T>()), 0);
  }
  TTS_ULP_EQUAL(eve__sph_bessel_j1(v_t(500)), std__sph_bessel_j1(v_t(500)), 6.0);
  TTS_ULP_EQUAL(eve__sph_bessel_j1(v_t(10)), std__sph_bessel_j1(v_t(10))  , 6.0);
  TTS_ULP_EQUAL(eve__sph_bessel_j1(v_t(5)),  std__sph_bessel_j1(v_t(5))   , 6.0);
  TTS_ULP_EQUAL(eve__sph_bessel_j1(v_t(2)),  std__sph_bessel_j1(v_t(2))   , 6.0);
  TTS_ULP_EQUAL(eve__sph_bessel_j1(v_t(1.5)),std__sph_bessel_j1(v_t(1.5)) , 6.0);
  TTS_ULP_EQUAL(eve__sph_bessel_j1(v_t(0.5)),std__sph_bessel_j1(v_t(0.5)) , 6.0);
  TTS_ULP_EQUAL(eve__sph_bessel_j1(v_t(1)),  std__sph_bessel_j1(v_t(1))   , 6.0);
  TTS_ULP_EQUAL(eve__sph_bessel_j1(v_t(0)),  eve::zero(eve::as<v_t>()), 0.0);

  TTS_ULP_EQUAL(eve__sph_bessel_j1( T(500)),  T(std__sph_bessel_j1(v_t(500)) ), 6.0);
  TTS_ULP_EQUAL(eve__sph_bessel_j1( T(10)) ,  T(std__sph_bessel_j1( v_t(10)) ), 6.0);
  TTS_ULP_EQUAL(eve__sph_bessel_j1( T(5))  ,  T(std__sph_bessel_j1( v_t(5))  ), 6.0);
  TTS_ULP_EQUAL(eve__sph_bessel_j1( T(2))  ,  T(std__sph_bessel_j1( v_t(2))  ), 6.0);
  TTS_ULP_EQUAL(eve__sph_bessel_j1( T(1.5)),  T(std__sph_bessel_j1( v_t(1.5))), 6.0);
  TTS_ULP_EQUAL(eve__sph_bessel_j1( T(0.5)),  T(std__sph_bessel_j1( v_t(0.5))), 6.0);
  TTS_ULP_EQUAL(eve__sph_bessel_j1( T(1))  ,  T(std__sph_bessel_j1( v_t(1))  ), 6.0);
  TTS_ULP_EQUAL(eve__sph_bessel_j1( T(0))  , eve::zero(eve::as< T>()), 0.0);

  TTS_ULP_EQUAL(eve__sph_bessel_j1(a0), map(std__sph_bessel_j1, a0), 40.0);
  TTS_ULP_EQUAL(eve__sph_bessel_j1(a1), map(std__sph_bessel_j1, a1), 40.0);
  TTS_ULP_EQUAL(eve__sph_bessel_j1(a2), map(std__sph_bessel_j1, a2), 40.0);

};

EVE_TEST( "Check behavior of diff(sph_bessel_j1) on wide"
        , eve::test::simd::ieee_reals
        , eve::test::generate(eve::test::randoms(1.0, 10.0))
        )
  <typename T>(T a0 )
{
  auto eve__diff_bessel_j1 =  [](auto x) { return eve::diff(eve::sph_bessel_j1)(x); };
  auto df = [](auto x){return eve::detail::centered_diffdiv(eve::sph_bessel_j1, x); };

  TTS_RELATIVE_EQUAL(eve__diff_bessel_j1(a0),   df(a0), 2.0e-2);
};

max_stars_repo_name: clayne/eve

max_issues_repo_name: clayne/eve

max_forks_repo_name: clayne/eve


text:
// Copyright (C) 2008  Davis E. King (davis@dlib.net)
// License: Boost Software License   See LICENSE.txt for the full license.

#include <sstream>
#include <string>
#include <cstdlib>
#include <ctime>
#include <dlib/algs.h>

#include "tester.h"

namespace
{
    using namespace test;
    using namespace dlib;
    using namespace std;

    logger dlog("test.metaprogramming");

    void metaprogramming_test (
    )
    /*!
        ensures
            - runs tests on template metaprogramming objects and functions for compliance with the specs
    !*/
    {

        print_spinner();

        DLIB_TEST(is_signed_type<signed char>::value == true);
        DLIB_TEST(is_signed_type<signed short>::value == true);
        DLIB_TEST(is_signed_type<signed int>::value == true);
        DLIB_TEST(is_signed_type<signed long>::value == true);
        DLIB_TEST(is_unsigned_type<signed char>::value == false);
        DLIB_TEST(is_unsigned_type<signed short>::value == false);
        DLIB_TEST(is_unsigned_type<signed int>::value == false);
        DLIB_TEST(is_unsigned_type<signed long>::value == false);

        DLIB_TEST(is_unsigned_type<unsigned char>::value == true);
        DLIB_TEST(is_unsigned_type<unsigned short>::value == true);
        DLIB_TEST(is_unsigned_type<unsigned int>::value == true);
        DLIB_TEST(is_unsigned_type<unsigned long>::value == true);
        DLIB_TEST(is_signed_type<unsigned char>::value == false);
        DLIB_TEST(is_signed_type<unsigned short>::value == false);
        DLIB_TEST(is_signed_type<unsigned int>::value == false);
        DLIB_TEST(is_signed_type<unsigned long>::value == false);

        COMPILE_TIME_ASSERT(is_signed_type<signed char>::value == true);
        COMPILE_TIME_ASSERT(is_signed_type<signed short>::value == true);
        COMPILE_TIME_ASSERT(is_signed_type<signed int>::value == true);
        COMPILE_TIME_ASSERT(is_signed_type<signed long>::value == true);
        COMPILE_TIME_ASSERT(is_unsigned_type<signed char>::value == false);
        COMPILE_TIME_ASSERT(is_unsigned_type<signed short>::value == false);
        COMPILE_TIME_ASSERT(is_unsigned_type<signed int>::value == false);
        COMPILE_TIME_ASSERT(is_unsigned_type<signed long>::value == false);

        COMPILE_TIME_ASSERT(is_unsigned_type<unsigned char>::value == true);
        COMPILE_TIME_ASSERT(is_unsigned_type<unsigned short>::value == true);
        COMPILE_TIME_ASSERT(is_unsigned_type<unsigned int>::value == true);
        COMPILE_TIME_ASSERT(is_unsigned_type<unsigned long>::value == true);
        COMPILE_TIME_ASSERT(is_signed_type<unsigned char>::value == false);
        COMPILE_TIME_ASSERT(is_signed_type<unsigned short>::value == false);
        COMPILE_TIME_ASSERT(is_signed_type<unsigned int>::value == false);
        COMPILE_TIME_ASSERT(is_signed_type<unsigned long>::value == false);

    }

    class metaprogramming_tester : public tester
    {
    public:
        metaprogramming_tester (
        ) :
            tester ("test_metaprogramming",
                    "Runs tests on the metaprogramming objects and functions.")
        {}

        void perform_test (
        )
        {
            metaprogramming_test();
        }
    } a;

}


text:
////////////////////////////////////////////////////////////////////////////////
//
// Copyright (C) 2018 Intel Corporation
//
// 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.
//
// SPDX-License-Identifier: MIT
//
////////////////////////////////////////////////////////////////////////////////

#include <sys/types.h>
#include <sys/stat.h>
#include <fcntl.h>
#include <errno.h>
#include <string.h>
#include <poll.h>
#include <boost/format.hpp>

#include "EALog.h"
#include "CBCEvent.hpp"
#include "CBCEventDevice.hpp"

// A log tag for CBCEventDevice.
#define TAG "DEV"

namespace earlyapp
{
    CBCEventDevice::CBCEventDevice(const char* cbcDevice)
    {
        m_PollFd = {0,};

        if(cbcDevice)
        {
            LINF_(TAG, "Opening device " << cbcDevice);
            m_fdCBCDev = open(cbcDevice, O_RDWR | O_NOCTTY | O_NDELAY);

            if(m_fdCBCDev > 0)
            {
                // Initialization for CBC signal.
                const char initCmd = 0xff;
                write(m_fdCBCDev, &initCmd, 1);

                // Initialize Poll FD.
                m_PollFd.fd = m_fdCBCDev;
                m_PollFd.events = POLLIN;

                m_bOpenSuccess = true;
                LINF_(TAG, " : success");
            }
            else
            {
                m_bOpenSuccess = false;
                LERR_(TAG,
                      boost::str(
                          boost::format("Failed to open device %s (%s)")
                          % cbcDevice
                          % strerror(errno)));
            }
        }
    }

    CBCEventDevice::~CBCEventDevice(void)
    {
        if(m_fdCBCDev > 0)
        {
            close(m_fdCBCDev);
        }
    }

    /*
      openSuccessfully.
      Returns whether the CBC event device open successfully.
     */
    bool CBCEventDevice::openSuccessfully(void)
    {
        return m_bOpenSuccess;
    }

    /*
      readEvent.
      Polls the device node until reads something from it.
     */
    std::shared_ptr<CBCEvent> CBCEventDevice::readEvent(void)
    {
        if(m_fdCBCDev < 0)
        {
            LERR_(TAG, "CBC device not open");
            return nullptr;
        }

        // Poll the device node.
        if(poll(&m_PollFd, 1, CBCEVENT_POLL_INTERVAL) < 0)
        {
            LERR_(TAG, "Failed to read CBC device");

            char* errMsg = strerror(errno);
            if(errMsg)
            {
                LERR_(TAG, ": " << errMsg);
            }
            return nullptr;
        }

        // Event received.
        if(m_PollFd.revents & POLLIN)
        {
            unsigned char cbcSignalBuffer[CBCBUFFER_SIZE];
            ssize_t readChars;

            // Read error.
            if((readChars = read(m_PollFd.fd, &cbcSignalBuffer, CBCBUFFER_SIZE)) < 0)
            {
                char* errMsg = strerror(errno);
                LERR_(TAG, "Failed to read CBC signal buffer");
                if(errMsg)
                {
                    LERR_(TAG, ": " << errMsg);
                }
                return nullptr;
            }
            // Nothing happend.
            else if(readChars == 0)
            {
                return nullptr;
            }

            // Allocate a CBCEvent to be sent.
            CBCEvent::eCBCEvent cbcEv = CBCEvent::eGEARSTATUS_UNKNOWN;
            unsigned char cbcDataPayload = cbcSignalBuffer[CBC_DATA_INDEX];
            switch(cbcDataPayload)
            {
                 case 0x01:
                     cbcEv = CBCEvent::eGEARSTATUS_REVERSE;
                     break;
                 case 0x03:
                     cbcEv = CBCEvent::eGEARSTATUS_FORWARD;
                     break;
                 case 0x05:
                     cbcEv = CBCEvent::eAPPLICATION_EXIT;
                     break;
                 default:
                     cbcEv = CBCEvent::eGEARSTATUS_UNKNOWN;
                     LWRN_(TAG,"Wrong button, Please input correct button");
            }
            std::shared_ptr<CBCEvent> e = std::make_shared<CBCEvent>(cbcEv);

            LINF_(TAG, "Buffer: " << std::hex << cbcSignalBuffer[0] << cbcSignalBuffer[1] << cbcSignalBuffer[2] << cbcSignalBuffer[3] << cbcSignalBuffer[4] << cbcSignalBuffer[5]);
            LINF_(TAG,
                  boost::str(
                      boost::format("Got an event %s:%d")
                      %e->toString()
                      %e->toEnum()));

            return e;
        }

        return nullptr;
    }
} // namespace

max_stars_repo_name: yewwayne/earlyapp

max_issues_repo_name: yewwayne/earlyapp

max_forks_repo_name: yewwayne/earlyapp


text:
/**
 * Copyright (c) 2017 Melown Technologies SE
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *
 * *  Redistributions of source code must retain the above copyright notice,
 *    this list of conditions and the following disclaimer.
 *
 * *  Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 * POSSIBILITY OF SUCH DAMAGE.
 */
/**
 * @file geometry.hpp
 * @author Ondrej Prochazka <ondrej.prochazka@citationtech.net>
 *
 * Analytical geometry functions
 */

#ifndef MATH_GEOMETRY_HPP
#define MATH_GEOMETRY_HPP

#include "geometry_core.hpp"

#include <algorithm>
#include <array>

#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/vector_proxy.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/triangular.hpp>
#include <boost/numeric/ublas/lu.hpp>
#include <boost/numeric/ublas/io.hpp>

namespace math {

namespace ublas = boost::numeric::ublas;

typedef std::array<math::Point2d, 3> Triangle2d;
typedef std::array<math::Point3d, 3> Triangle3d;

typedef std::vector<Triangle2d> Triangles2d;
typedef std::vector<Triangle3d> Triangles3d;

typedef Points2d Polygon; // single CCW ring, not closed
typedef std::vector<Polygon> MultiPolygon; // multiple rings, holes CW

/**
 * Find the point where two lines get closest in 3D space.
 * The lines are defined parametrically in euclidian coordinates.
 * Return value is the euclidian metrics of the line distance, plus
 * the two parameter values for each of the lines.
 */

float lineDistance(
    const ublas::vector<float> & p1, const ublas::vector<float> & u1,
    const ublas::vector<float> & p2, const ublas::vector<float> & u2,
    float & r1, float & r2 );

/** normalize vector */
template <class T>
inline ublas::vector<typename T::value_type> normalize( const T & v ) {
    return v / ublas::norm_2( v );
}

/** normalize vector */
template <class T>
inline Point2_<T> normalize( const Point2_<T> & v ) {
    return v / ublas::norm_2( v );
}

/** normalize vector */
template <class T>
inline Point3_<T> normalize( const Point3_<T> & v ) {
    return v / ublas::norm_2( v );
}

/** length of the vector */
template <typename T>
inline double length( const T & v ) {
    return ublas::norm_2(v);
}

/** angle between two 3D vectors (in radians)
 * See Kahan W.: How Futile are Mindless Assessments of Roundoff in
 * Floating-Point Computation?, page 46
 * (https://people.eecs.berkeley.edu/~wkahan/Mindless.pdf)
 *  **/
template <class T>
inline T angle(const Point3_<T>& v1, const Point3_<T>& v2)
{
    return 2
           * std::atan2(length(v1 * length(v2) - length(v1) * v2),
                        length(v1 * length(v2) + length(v1) * v2));
}

/** homogeneous coordinates (from euclidian) */
template <class T>
inline ublas::vector<T, ublas::bounded_array<T, 4> >
homogeneous( const Point3_<T> & src )
{
    ublas::vector<T, ublas::bounded_array<T, 4> > dst(4);
    dst(0) = src(0); dst(1) = src(1); dst(2) = src(2);
    dst(3) = T(1);
    return dst;
}

/** homogeneous coordinates (from euclidian) */
template <class T>
inline ublas::vector<typename T::value_type> homogeneous( const T & src ) {

    ublas::vector<typename T::value_type> ret(
        ublas::unit_vector<typename T::value_type>(
            src.size() + 1, src.size() ) );
    ublas::subrange( ret, 0, src.size() ) = src;
    return ret;
}

/** euclidian coordinates (from homogenous) */
template <class T>
inline Point3_<T>
euclidian( const ublas::vector<T, ublas::bounded_array<T, 4> > & src )
{
    auto div(src(3));
    return Point3_<T>(src(0) / div, src(1) / div, src(2) / div);
}

/** euclidian coordinates (from homogenous) */
template <class T>
inline ublas::vector<typename T::value_type> euclidian( const T & src ) {
    return ublas::subrange( src, 0, src.size() - 1 ) / src( src.size() - 1 );
}

/** cross product, in euclidian 3D */

template<typename T, typename U>
inline ublas::vector<typename T::value_type> crossProduct(
    const T & u,
    const U & v ) {

    assert( u.size() == 3 && v.size() == 3 );
    ublas::vector<typename T::value_type> retval( 3 );
    retval(0) = u(1) * v(2) - v(1) * u(2);
    retval(1) = -u(0) * v(2) + v(0) * u(2);
    retval(2) = u(0) * v(1) - v(0) * u(1);
    return retval;
}

template <typename T>
inline Point3_<T> crossProduct(const Point3_<T> & u, const Point3_<T> & v )
{
    return Point3_<T>(u(1) * v(2) - v(1) * u(2)
                     , -u(0) * v(2) + v(0) * u(2)
                     , u(0) * v(1) - v(0) * u(1));
}

template<typename T>
inline Point3_<T> crossProduct(const Point3_<T> & u, const ublas::vector<T> & v )
{
    assert(v.size() == 3);
    return Point3_<T>(u(1) * v(2) - v(1) * u(2)
                     , -u(0) * v(2) + v(0) * u(2)
                     , u(0) * v(1) - v(0) * u(1));
}

template<typename T>
inline Point3_<T> crossProduct(const ublas::vector<T> & u, const Point3_<T> & v )
{
    assert(u.size() == 3);
    return Point3_<T>(u(1) * v(2) - v(1) * u(2)
                     , -u(0) * v(2) + v(0) * u(2)
                     , u(0) * v(1) - v(0) * u(1));
}

/** 2D version of cross product; result is not an vector but crossproduct's
 *  z-component.
 */
template <typename T>
inline T crossProduct(const Point2_<T> & u, const Point2_<T> & v )
{
    return u(0) * v(1) - v(0) * u(1);
}

//! return area of a triangle in 2D plane
template <typename T>
inline double triangleArea(const Point2_<T>& a, const Point2_<T>& b,
                           const Point2_<T>& c)
{
    return std::abs(crossProduct(Point2(b - a), Point2(c - a))) * 0.5;
}

//! return area of a triangle in 3D space
template <typename T>
inline double triangleArea(const Point3_<T>& a, const Point3_<T>& b,
                           const Point3_<T>& c)
{
    return norm_2(crossProduct(b - a, c - a)) * 0.5;
}

#ifdef MATH_HAS_OPENCV
    //! return area of a triangle in 2D plane
    template <typename T>
    inline double triangleArea(const cv::Point_<T>& a, const cv::Point_<T>& b,
                               const cv::Point_<T>& c)
    {
        return triangleArea(Point2_<T>(a), Point2_<T>(b), Point2_<T>(c));
    }

    //! return area of a triangle in 3D space
    template <typename T>
    inline double triangleArea(const cv::Point3_<T>& a, const cv::Point3_<T>& b,
                               const cv::Point3_<T>& c)
    {
        return triangleArea(Point3_<T>(a), Point3_<T>(b), Point3_<T>(c));
    }
#endif // MATH_HAS_OPENCV

/** Returns a positive number if the sequence of points {a, b, c} turns
 *  counter-clockwise in the XY plane (negative number otherwise).
 *  Zero (within numerical tolerance) means collinear.
 */
template <typename T>
inline T ccw(const Point2_<T> &a, const Point2_<T> &b, const Point2_<T> &c)
{
    return (b(0) - a(0))*(c(1) - a(1)) - (b(1) - a(1))*(c(0) - a(0));
}

/** Parametric line, in euclidian 2D
 */

template <typename T>
struct Line2_
{
    T p, u;

    Line2_(const T p, const T u) : p(p), u(u) {};
};

typedef Line2_<Point2> Line2;

/** Returns a point where two lines `l1` and `l2` intersect.
 */
template <typename T>
inline T lineIntersection(const Line2_<T>& l1, const Line2_<T>& l2)
{
    T pd = l2.p - l1.p; // diff of origins
    auto num = -l2.u(1) * pd(0) + l2.u(0) * pd(1);
    auto den = l2.u(0) * l1.u(1) - l1.u(0) * l2.u(1);
    return l1.p + (num / den) * l1.u;
}

/**
 * Line segment represented by start and end points
 */
template <typename T>
struct Segment2_
{
    T p1, p2;

    Segment2_(const T p1, const T p2) : p1(p1), p2(p2) {};
};

typedef Segment2_<Point2> Segment2;

/**
 * Parametric line, in euclidian 3D
 */

struct Line3 {
    Point3 p, u;

    Line3(const Point3 p = Point3(), const Point3 u = Point3() )
        : p( p ), u( u ) {}

    /** Returns line's point at given parameter (t)
     */
    Point3 point(double t) const { return p + u * t; }
};

template <typename E, typename T>
inline std::basic_ostream<E, T> & operator << (
        std::basic_ostream<E,T> & os,
        const Line3 & line ) {

    os << line.p << " + t * " << line.u;
    return os;
}

Point3 midpoint( const Line3 & line1, const Line3 & line2
               , double minAngleCos = 0.9962 );

double pointLineDistance(const Point3 &p, const Line3 &line);

/**
 * Parametric plane, in euclidian 3D
 * Legacy representation of plane by 3 points.
 */
namespace legacy {
struct Plane3 {

    Point3 p, u, v;

    Plane3(
        const ublas::vector<double> p = ublas::zero_vector<double>( 3 ),
        const ublas::vector<double> u = ublas::zero_vector<double>( 3 ),
        const ublas::vector<double> v = ublas::zero_vector<double>( 3 ) )
        : p( p ), u( u ), v( v ) {}

    template <class Matrix>
    Plane3 transform( const Matrix & trafo ) const {

        Plane3 tplane;

        tplane.p = euclidian( prod( trafo, homogeneous( p ) ) );
        tplane.u = euclidian( prod( trafo, homogeneous( p + u ) ) ) - tplane.p;
        tplane.v = euclidian( prod( trafo, homogeneous( p + v ) ) ) - tplane.p;

        return tplane;

    }
};

template <typename E, typename T>
inline std::basic_ostream<E, T> & operator << (
        std::basic_ostream<E,T> & os,
        const Plane3 & plane ) {

    os << plane.p << " + t * " << plane.u << " + s * " << plane.v;
    return os;
}
}  // namespace legacy

/** line and plane intersection */

Point3 intersection( const Line3 & line, const legacy::Plane3 & plane );

/** line and plane intersection
 *  instead of point returns 3 coefficients:
 *      * lines t-parameter
 *      * planes t-parameter
 *      * planes s-parameter
 */
Point3 intersectionParams(const Line3 &line, const legacy::Plane3 &plane);

/**
* Plane represented in general form
* Representation as (a,b,c,d) in equation: ax+by+cz+d=0
*/
struct Plane3
{
    math::Point3 n_;  // vector with (a,b,c) i.e. normal vector
    double d_;

    // Plane from parameters
    Plane3(double a, double b, double c, double d) : n_(a, b, c), d_(d) { }

    // Plane from point and normal
    Plane3(const math::Point3d& pt, const math::Point3d& n)
        : n_(n),
          d_(-ublas::inner_prod(pt, n)) {};

    // Plane from three points
    Plane3(const math::Point3d& p1,
           const math::Point3d& p2,
           const math::Point3d& p3)
        : Plane3(p1, crossProduct(p2 - p1, p3 - p1)) {};

    // Plane from legacy representation (point-normal packed as Line3)
    Plane3(const math::Line3& l) : Plane3(l.p, l.u) {};

    // Returns a plane with opposite orientation
    Plane3 opposite() const {
        Plane3 pl2(*this);
        pl2.n_ *= -1;
        pl2.d_ *= -1;
        return pl2;
    }
};

/**
 * Returns a distance of point to a plane (perpendicular)
*/
double pointPlaneDistance(const Point3 &p, const Plane3 &plane);

/**
 * Returns an orthogonal projection of a point to a plane
 */
Point3 pointPlaneProjection(const Point3 &p, const Plane3 &plane);

/**
 * Returns a point of intersection between plane and line
 */
Point3 linePlaneIntersection(const Line3 &l, const Plane3 &plane);

/**
 * Returns a line of intersection between two planes
 */
Line3 planeIntersection(const Plane3 &p1, const Plane3 &p2);

/**
 * Returns a point of intersection between three planes
 */
Point3 planeIntersection(const Plane3 &p1, const Plane3 &p2, const Plane3 &p3);

/**
 * Returns a measure of triangular polyface regularity. Value of 1.0
 * indicates an equilateral triangle, value of 0.0 indicates a triangle
 * with at least one degenerate edge.
 */
double polygonRegularity(
    const Point3 & v0, const Point3 & v1, const Point3 & v2  );

/**
 * Returns a measure of quad polyface regularity. Value of 1.0
 * indicates a square, value of 0.0 indicates a triangle
 * with at least one degenerate edge.
 */
double polygonRegularity(
    const Point3 & v0, const Point3 & v1, const Point3 & v2, const Point3 & v3 );

/**
 * Returns whether the triangle and rectagle collide.
 */
bool triangleRectangleCollision( math::Point2 triangle[3]
                               , math::Point2 ll, math::Point2 ur);

/**
 * Convert cartesian coordinates (r) to barycentric with respect to the (a,b,c)
 * triangle. See http://en.wikipedia.org/wiki/Barycentric_coordinate_system
 */
math::Point3 barycentricCoords(const math::Point2 &r, const math::Point2 &a
                               , const math::Point2 &b, const math::Point2 &c);
inline
math::Point3 barycentricCoords(const math::Point2 &r, const Triangle2d &tri) {
    return barycentricCoords(r, tri[0], tri[1], tri[2]);
}

namespace detail {
    template <typename T, typename Q> struct ExtentsTypeTraits;

    template <typename T> struct ExtentsTypeTraits<T, Point2_<T>> {
        typedef Extents2_<T> type;
    };

    template <typename T> struct ExtentsTypeTraits<T, Point3_<T>> {
        typedef Extents3_<T> type;
    };

#ifdef MATH_HAS_OPENCV
    template <typename T> struct ExtentsTypeTraits<T, cv::Point_<T>> {
        typedef Extents2_<T> type;
    };

    template <typename T> struct ExtentsTypeTraits<T, cv::Point3_<T>> {
        typedef Extents3_<T> type;
    };
#endif
} // namespace detail

template <typename Iterator>
inline auto computeExtents(Iterator begin, Iterator end)
    -> typename detail::ExtentsTypeTraits
    <typename std::iterator_traits<Iterator>::value_type::value_type
    , typename std::iterator_traits<Iterator>::value_type>::type
{
    typedef typename detail::ExtentsTypeTraits
        <typename std::iterator_traits<Iterator>::value_type::value_type
         , typename std::iterator_traits<Iterator>::value_type>::type Extents;
    typedef typename Extents::point_type point_type;

    if (begin == end) {
        return Extents(InvalidExtents{});
    }

    Extents extents(*begin++);
    std::for_each(begin, end
                  , [&extents](const point_type &p) {
                      update(extents, p);
                  });
    return extents;
}

/** Simplified version for STL container.
 */
template <typename Container>
inline auto computeExtents(const Container &c)
    -> typename detail::ExtentsTypeTraits
    <typename Container::value_type::value_type
     , typename Container::value_type>::type
{
    return computeExtents(c.begin(), c.end());
}

// prefered natural aliases for above

template <typename Iterator>
inline auto extents(Iterator begin, Iterator end)
    -> typename detail::ExtentsTypeTraits
    <typename std::iterator_traits<Iterator>::value_type::value_type
    , typename std::iterator_traits<Iterator>::value_type>::type
{
    return computeExtents<Iterator>(begin,end);
}

template <typename Container>
inline auto extents(const Container &c)
    -> typename detail::ExtentsTypeTraits
    <typename Container::value_type::value_type
     , typename Container::value_type>::type
{
    return computeExtents<Container>(c);
}

} // namespace math

#endif // MATH_GEOMETRY_HPP

max_stars_repo_name: melowntech/libmath

max_issues_repo_name: melowntech/libmath

max_forks_repo_name: melowntech/libmath


text:
/**
 * @file
 * @copyright This code is licensed under the 3-clause BSD license.\n
 *            Copyright ETH Zurich, Laboratory for Physical Chemistry, Reiher Group.\n
 *            See LICENSE.txt for details.
 */

#include "Ediis.h"
#include "EdiisCoefficientOptimizer.h"
#include <Eigen/QR>
#include <iostream>

namespace Scine {
namespace Utils {

Ediis::Ediis() {
  setSubspaceSize(6);
}

void Ediis::setUnrestricted(bool b) {
  unrestricted_ = b;
}

void Ediis::setSubspaceSize(int n) {
  bool resizeNeeded = n != subspaceSize_;

  subspaceSize_ = n;

  if (resizeNeeded)
    resizeMembers();
}

void Ediis::setNAOs(int n) {
  bool resizeNeeded = n != nAOs_;

  nAOs_ = n;

  if (resizeNeeded)
    resizeMembers();
}

void Ediis::resizeMembers() {
  fockMatrices.resize(subspaceSize_);
  densityMatrices.resize(subspaceSize_);
  energies.resize(subspaceSize_);

  B = Eigen::MatrixXd::Zero(subspaceSize_, subspaceSize_);

  restart();
}

void Ediis::restart() {
  iterationNo_ = 0;
  index_ = 0;
}

void Ediis::addMatrices(double energy, const SpinAdaptedMatrix& F, const DensityMatrix& P) {
  iterationNo_++;
  lastAdded_ = index_;

  fockMatrices[index_] = F;
  densityMatrices[index_] = P;
  energies[index_] = energy;

  updateBMatrix();

  index_ = (index_ + 1) % subspaceSize_;
}

void Ediis::updateBMatrix() {
  int activeSize = iterationNo_ > subspaceSize_ ? subspaceSize_ : iterationNo_;

  // Bii element
  B(lastAdded_, lastAdded_) = 0;

  // Bij elements
  for (int i = 0; i < activeSize; i++) {
    if (i == lastAdded_)
      continue;
    double v = getBMatrixElement(lastAdded_, i);
    B(lastAdded_, i) = v;
    B(i, lastAdded_) = B(lastAdded_, i);
  }
}

double Ediis::getBMatrixElement(int i, int j) const {
  if (unrestricted_) {
    double va = ((fockMatrices[i].alphaMatrix() - fockMatrices[j].alphaMatrix()).selfadjointView<Eigen::Lower>() *
                 (densityMatrices[i].alphaMatrix() - densityMatrices[j].alphaMatrix()))
                    .trace();
    double vb = ((fockMatrices[i].betaMatrix() - fockMatrices[j].betaMatrix()).selfadjointView<Eigen::Lower>() *
                 (densityMatrices[i].betaMatrix() - densityMatrices[j].betaMatrix()))
                    .trace();
    return (va + vb) / 2;
  }
  else {
    double v = ((fockMatrices[i].restrictedMatrix() - fockMatrices[j].restrictedMatrix()).selfadjointView<Eigen::Lower>() *
                (densityMatrices[i].restrictedMatrix() - densityMatrices[j].restrictedMatrix()))
                   .trace();
    // Divide v by two because density matrix is RHF
    return v / 2;
  }
}

SpinAdaptedMatrix Ediis::getMixedFockMatrix() {
  if (iterationNo_ > subspaceSize_)
    iterationNo_ = subspaceSize_;

  // If we have only one Fock matrix
  if (iterationNo_ < 2) {
    return fockMatrices[0];
  }
  else {
    EdiisCoefficientOptimizer opt(energies, B.block(0, 0, iterationNo_, iterationNo_));
    auto coefs = opt.getCoefficients();

    return calculateLinearCombination(coefs);
  }
}

SpinAdaptedMatrix Ediis::calculateLinearCombination(const Eigen::VectorXd& coefs) {
  if (unrestricted_) {
    Eigen::MatrixXd FAlpha = Eigen::MatrixXd::Zero(nAOs_, nAOs_);
    Eigen::MatrixXd FBeta = Eigen::MatrixXd::Zero(nAOs_, nAOs_);
    for (int i = 0; i < iterationNo_; i++) {
      FAlpha += coefs[i] * fockMatrices[i].alphaMatrix();
      FBeta += coefs[i] * fockMatrices[i].betaMatrix();
    }
    return SpinAdaptedMatrix::createUnrestricted(std::move(FAlpha), std::move(FBeta));
  }
  else {
    Eigen::MatrixXd Fsol = Eigen::MatrixXd::Zero(nAOs_, nAOs_);
    for (int i = 0; i < iterationNo_; i++)
      Fsol += coefs[i] * fockMatrices[i].restrictedMatrix();
    return SpinAdaptedMatrix::createRestricted(std::move(Fsol));
  }
}
} // namespace Utils
} // namespace Scine

max_stars_repo_name: DockBio/utilities

max_issues_repo_name: DockBio/utilities

max_forks_repo_name: DockBio/utilities