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

#ifndef CRYPTO3_ALGEBRA_CURVES_HASH_TO_CURVE_EXPAND_HPP
#define CRYPTO3_ALGEBRA_CURVES_HASH_TO_CURVE_EXPAND_HPP

#include <nil/crypto3/algebra/curves/detail/h2c/h2c_suites.hpp>
#include <nil/crypto3/algebra/curves/detail/h2c/h2c_sgn0.hpp>

#include <nil/crypto3/algebra/algorithms/strxor.hpp>

#include <nil/crypto3/algebra/curves/bls12.hpp>

#include <nil/crypto3/hash/algorithm/hash.hpp>
#include <nil/crypto3/hash/accumulators/hash.hpp>

#include <boost/assert.hpp>
#include <boost/static_assert.hpp>
#include <boost/concept/assert.hpp>

#include <array>
#include <type_traits>
#include <iterator>

namespace nil {
    namespace crypto3 {
        namespace algebra {
            namespace curves {
                namespace detail {
                    using namespace nil::crypto3::detail;
                    template<std::size_t k, typename HashType,
                             /// HashType::digest_type is required to be uint8_t[]
                             typename = typename std::enable_if<
                                 std::is_same<std::uint8_t, typename HashType::digest_type::value_type>::value>::type>
                    class expand_message_xmd {
                        // https://tools.ietf.org/html/draft-irtf-cfrg-hash-to-curve-10#section-5.4.1
                        BOOST_STATIC_ASSERT_MSG(HashType::block_bits % 8 == 0, "r_in_bytes is not a multiple of 8");
                        BOOST_STATIC_ASSERT_MSG(HashType::digest_bits % 8 == 0, "b_in_bytes is not a multiple of 8");
                        BOOST_STATIC_ASSERT_MSG(HashType::digest_bits >= 2 * k,
                                                "k-bit collision resistance is not fulfilled");

                        constexpr static const std::size_t b_in_bytes = HashType::digest_bits / 8;
                        constexpr static const std::size_t r_in_bytes = HashType::block_bits / 8;

                        constexpr static const std::array<std::uint8_t, r_in_bytes> Z_pad {0};

                    public:
                        template<typename InputMsgType, typename InputDstType, typename OutputType,
                                 typename = typename std::enable_if<
                                     std::is_same<std::uint8_t, typename InputMsgType::value_type>::value &&
                                     std::is_same<std::uint8_t, typename InputDstType::value_type>::value &&
                                     std::is_same<std::uint8_t, typename OutputType::value_type>::value>::type>
                        static inline void process(const std::size_t len_in_bytes, const InputMsgType &msg,
                                                   const InputDstType &dst, OutputType &uniform_bytes) {
                            BOOST_CONCEPT_ASSERT((boost::SinglePassRangeConcept<InputMsgType>));
                            BOOST_CONCEPT_ASSERT((boost::SinglePassRangeConcept<InputDstType>));
                            BOOST_CONCEPT_ASSERT((boost::SinglePassRangeConcept<OutputType>));
                            BOOST_CONCEPT_ASSERT((boost::WriteableRangeConcept<OutputType>));

                            // https://tools.ietf.org/html/draft-irtf-cfrg-hash-to-curve-10#section-5.4.1
                            BOOST_ASSERT(len_in_bytes < 0x10000);
                            BOOST_ASSERT(std::distance(dst.begin(), dst.end()) >= 16 &&
                                         std::distance(dst.begin(), dst.end()) <= 255);
                            BOOST_ASSERT(std::distance(uniform_bytes.begin(), uniform_bytes.end()) >= len_in_bytes);

                            const std::array<std::uint8_t, 2> l_i_b_str = {
                                static_cast<std::uint8_t>(len_in_bytes >> 8u),
                                static_cast<std::uint8_t>(len_in_bytes % 0x100)};
                            const std::size_t ell = static_cast<std::size_t>(len_in_bytes / b_in_bytes) +
                                                    static_cast<std::size_t>(len_in_bytes % b_in_bytes != 0);

                            // https://tools.ietf.org/html/draft-irtf-cfrg-hash-to-curve-10#section-5.4.1
                            BOOST_ASSERT(ell <= 255);

                            // TODO: use accumulators when they will be fixed
                            // accumulator_set<HashType> b0_acc;
                            // hash<HashType>(Z_pad, b0_acc);
                            // hash<HashType>(msg, b0_acc);
                            // hash<HashType>(l_i_b_str, b0_acc);
                            // hash<HashType>(std::array<std::uint8_t, 1> {0}, b0_acc);
                            // hash<HashType>(dst, b0_acc);
                            // hash<HashType>(std::array<std::uint8_t, 1> {static_cast<std::uint8_t>(dst.size())},
                            // b0_acc); typename HashType::digest_type b0 =
                            // accumulators::extract::hash<HashType>(b0_acc);
                            std::vector<std::uint8_t> msg_prime;
                            msg_prime.insert(msg_prime.end(), Z_pad.begin(), Z_pad.end());
                            msg_prime.insert(msg_prime.end(), msg.begin(), msg.end());
                            msg_prime.insert(msg_prime.end(), l_i_b_str.begin(), l_i_b_str.end());
                            msg_prime.insert(msg_prime.end(), static_cast<std::uint8_t>(0));
                            msg_prime.insert(msg_prime.end(), dst.begin(), dst.end());
                            msg_prime.insert(msg_prime.end(),
                                             static_cast<std::uint8_t>(std::distance(dst.begin(), dst.end())));
                            typename HashType::digest_type b0 = hash<HashType>(msg_prime);

                            // TODO: use accumulators when they will be fixed
                            // accumulator_set<HashType> bi_acc;
                            // hash<HashType>(b0, bi_acc);
                            // hash<HashType>(std::array<std::uint8_t, 1> {1}, bi_acc);
                            // hash<HashType>(dst, bi_acc);
                            // hash<HashType>(std::array<std::uint8_t, 1> {static_cast<std::uint8_t>(dst.size())},
                            // bi_acc); typename HashType::digest_type bi =
                            // accumulators::extract::hash<HashType>(bi_acc); std::copy(bi.begin(), bi.end(),
                            // uniform_bytes.begin());
                            std::vector<std::uint8_t> b_i_str;
                            b_i_str.insert(b_i_str.end(), b0.begin(), b0.end());
                            b_i_str.insert(b_i_str.end(), static_cast<std::uint8_t>(1));
                            b_i_str.insert(b_i_str.end(), dst.begin(), dst.end());
                            b_i_str.insert(b_i_str.end(),
                                           static_cast<std::uint8_t>(std::distance(dst.begin(), dst.end())));
                            typename HashType::digest_type bi = hash<HashType>(b_i_str);
                            std::copy(bi.begin(), bi.end(), uniform_bytes.begin());

                            typename HashType::digest_type xored_b;
                            for (std::size_t i = 2; i <= ell; i++) {
                                // TODO: use accumulators when they will be fixed
                                // accumulator_set<HashType> bi_acc;
                                // strxor(b0, bi, xored_b);
                                // hash<HashType>(xored_b, bi_acc);
                                // hash<HashType>(std::array<std::uint8_t, 1> {static_cast<std::uint8_t>(i)}, bi_acc);
                                // hash<HashType>(dst, bi_acc);
                                // hash<HashType>(std::array<std::uint8_t, 1> {static_cast<std::uint8_t>(dst.size())},
                                //                bi_acc);
                                // bi = accumulators::extract::hash<HashType>(bi_acc);
                                // std::copy(bi.begin(), bi.end(), uniform_bytes.begin() + (i - 1) * b_in_bytes);
                                strxor(b0, bi, xored_b);
                                std::vector<std::uint8_t> b_i_str;
                                b_i_str.insert(b_i_str.end(), xored_b.begin(), xored_b.end());
                                b_i_str.insert(b_i_str.end(), static_cast<std::uint8_t>(i));
                                b_i_str.insert(b_i_str.end(), dst.begin(), dst.end());
                                b_i_str.insert(b_i_str.end(),
                                               static_cast<std::uint8_t>(std::distance(dst.begin(), dst.end())));
                                bi = hash<HashType>(b_i_str);
                                std::copy(bi.begin(), bi.end(), uniform_bytes.begin() + (i - 1) * b_in_bytes);
                            }
                        }
                    };
                }    // namespace detail
            }        // namespace curves
        }            // namespace algebra
    }                // namespace crypto3
}    // namespace nil

#endif    // CRYPTO3_ALGEBRA_CURVES_HASH_TO_CURVE_EXPAND_HPP

max_forks_repo_name: NilFoundation/algebra


text:
//
// $Id$
//
//
// Original author: Matt Chambers <matt.chambers .@. vanderbilt.edu>
//
// Copyright 2014 Vanderbilt University - Nashville, TN 37232
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
//

#ifndef _SHIMADZUREADER_HPP_
#define _SHIMADZUREADER_HPP_

#include "pwiz/utility/misc/Export.hpp"
#include "pwiz/utility/misc/BinaryData.hpp"
#include <string>
#include <vector>
#include <set>
#include <boost/shared_ptr.hpp>
#include <boost/date_time.hpp>

namespace pwiz {
namespace vendor_api {
namespace Shimadzu {

struct PWIZ_API_DECL TimeRange { double start, end; };

PWIZ_API_DECL enum Polarity
{
    Positive = 0,
    Negative = 1,
    Undefined = 2
};

struct PWIZ_API_DECL SRMTransition
{
    short id;
    short channel;
    short event;
    short segment;
    double collisionEnergy;
    short polarity;
    double Q1;
    double Q3;
    //TimeRange acquiredTimeRange;

    bool operator< (const SRMTransition& rhs) const;
};

struct PWIZ_API_DECL Chromatogram
{
    virtual int getTotalDataPoints() const = 0;
    virtual void getXArray(pwiz::util::BinaryData<double>& x) const = 0;
    virtual void getYArray(pwiz::util::BinaryData<double>& y) const = 0;

    virtual ~Chromatogram() {}
};

typedef boost::shared_ptr<Chromatogram> ChromatogramPtr;

struct PWIZ_API_DECL SRMChromatogram : public Chromatogram
{
    virtual const SRMTransition& getTransition() const = 0;
};

typedef boost::shared_ptr<SRMChromatogram> SRMChromatogramPtr;

struct PWIZ_API_DECL SpectrumInfo
{
    double scanTime;
    int msLevel;
    double precursorMz;
    unsigned int precursorScan;
    Polarity polarity;
    int segment;
    int event;
};

struct PWIZ_API_DECL Spectrum
{
    virtual double getScanTime() const = 0;
    virtual int getMSLevel() const = 0;
    virtual Polarity getPolarity() const = 0;

    virtual double getSumY() const = 0;
    virtual double getBasePeakX() const = 0;
    virtual double getBasePeakY() const = 0;
    virtual double getMinX() const = 0;
    virtual double getMaxX() const = 0;

    virtual bool getHasIsolationInfo() const = 0;
    virtual void getIsolationInfo(double& centerMz, double& lowerLimit, double& upperLimit) const = 0;

    virtual bool getHasPrecursorInfo() const = 0;
    virtual void getPrecursorInfo(double& selectedMz, double& intensity, int& charge) const = 0;

    virtual int getTotalDataPoints(bool doCentroid = false) const = 0;
    virtual void getProfileArrays(pwiz::util::BinaryData<double>& x, pwiz::util::BinaryData<double>& y) const = 0;
    virtual void getCentroidArrays(pwiz::util::BinaryData<double>& x, pwiz::util::BinaryData<double>& y) const = 0;

    virtual ~Spectrum() {}
};

typedef boost::shared_ptr<Spectrum> SpectrumPtr;

class PWIZ_API_DECL ShimadzuReader
{
public:
    typedef boost::shared_ptr<ShimadzuReader> Ptr;
    static Ptr create(const std::string& filepath);

    //virtual std::string getVersion() const = 0;
    //virtual DeviceType getDeviceType() const = 0;
    //virtual std::string getDeviceName(DeviceType deviceType) const = 0;
    virtual boost::local_time::local_date_time getAnalysisDate(bool adjustToHostTime) const = 0;

    virtual const std::set<SRMTransition>& getTransitions() const = 0;
    virtual SRMChromatogramPtr getSRM(const SRMTransition& transition) const = 0;

    virtual ChromatogramPtr getTIC(bool ms1Only = false) const = 0;

    virtual int getScanCount() const = 0;
    virtual SpectrumPtr getSpectrum(int scanNumber, bool profileDesired) const = 0;
    virtual SpectrumInfo getSpectrumInfo(int scanNumber) const = 0;

    virtual const std::set<int>& getMSLevels() const = 0;

    virtual ~ShimadzuReader() {}
};

typedef ShimadzuReader::Ptr ShimadzuReaderPtr;

} // Shimadzu
} // vendor_api
} // pwiz

#endif // _SHIMADZUREADER_HPP_

ext: hpp

max_stars_repo_name: jeleclaire/pwiz

max_issues_repo_name: jeleclaire/pwiz

max_forks_repo_name: jeleclaire/pwiz


text:
#ifdef __GNUC__
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wdeprecated-declarations"
#endif
#include <boost/locale/encoding.hpp>
#include <boost/locale/util.hpp>
#ifdef __GNUC__
#pragma GCC diagnostic pop
#endif

#include "pkzip_io.h"

using namespace std;

using boost::locale::conv::from_utf;
using boost::locale::conv::to_utf;
using boost::locale::conv::utf_to_utf;
using boost::locale::util::get_system_locale;

static inline auto read(istream &is, void *data, size_t size)
{
    return is.read(static_cast<char *>(data), size)
        && is.gcount() == static_cast<streamsize>(size);
}
template <typename value_type>
static inline auto read(istream &is, value_type &value)
{
    return read(is, &value, sizeof value);
}

static inline auto write(ostream &os, const void *data, size_t size)
{
    os.write(static_cast<const char *>(data), size);
}
template <typename value_type>
static inline auto write(ostream &os, const value_type &value)
{
    write(os, &value, sizeof value);
}

string zz::pkzip::charset = "cp932"; // TODO: resolve from OS environment

istream & zz::pkzip::operator >> (istream &is, local_file_header &header)
{
    if (!read(is, header.signature) || !header ||
        !read(is, header.version_needed_to_extract) ||
        !read(is, header.general_purpose_bit_flag) ||
        !read(is, header.compression_method) ||
        !read(is, header.last_mod_file_time) ||
        !read(is, header.last_mod_file_date) ||
        !read(is, header.crc32) ||
        !read(is, header.compressed_size) ||
        !read(is, header.uncompressed_size) ||
        !read(is, header.file_name_length) || header.file_name_length == 0 ||
        !read(is, header.extra_field_length))
        return is;

    string file_name(header.file_name_length, '\0');
    if (!read(is, &file_name[0], header.file_name_length))
        return is;
    const auto use_utf8 = header.general_purpose_bit_flag & general_purpose_bit_flags::use_utf8;
    header.file_name = use_utf8 ? utf_to_utf<char_type>(file_name) : to_utf<char_type>(file_name, charset);
    if (header.extra_field_length) {
        header.extra_field.resize(header.extra_field_length);
        if (!read(is, &header.extra_field[0], header.extra_field_length))
            return is;
    }
    return is;
}

ostream & zz::pkzip::operator << (ostream &os, const local_file_header &header)
{
    const auto use_utf8 = header.general_purpose_bit_flag & general_purpose_bit_flags::use_utf8;
    const auto file_name = use_utf8 ? utf_to_utf<char>(header.file_name) : from_utf(header.file_name, charset);
    if (file_name.size() > numeric_limits<decltype(header.file_name_length)>::max())
        throw runtime_error("too long file name: " + file_name);
    if (header.extra_field.size() > numeric_limits<decltype(header.extra_field_length)>::max())
        throw runtime_error("too long extra field: " + file_name);

    write(os, header.signature);
    write(os, header.version_needed_to_extract);
    write(os, header.general_purpose_bit_flag);
    write(os, header.compression_method);
    write(os, header.last_mod_file_time);
    write(os, header.last_mod_file_date);
    write(os, header.crc32);
    write(os, header.compressed_size);
    write(os, header.uncompressed_size);
    write(os, static_cast<decltype(header.file_name_length)>(file_name.size()));
    write(os, static_cast<decltype(header.extra_field_length)>(header.extra_field.size()));

    write(os, file_name.data(), file_name.size());
    write(os, header.extra_field.data(), header.extra_field.size());
    return os;
}

istream & zz::pkzip::operator >> (istream &is, central_file_header &header)
{
    if (!read(is, header.signature) || !header ||
        !read(is, header.version_made_by) ||
        !read(is, header.version_needed_to_extract) ||
        !read(is, header.general_purpose_bit_flag) ||
        !read(is, header.compression_method) ||
        !read(is, header.last_mod_file_time) ||
        !read(is, header.last_mod_file_date) ||
        !read(is, header.crc32) ||
        !read(is, header.compressed_size) ||
        !read(is, header.uncompressed_size) ||
        !read(is, header.file_name_length) || header.file_name_length == 0 ||
        !read(is, header.extra_field_length) ||
        !read(is, header.file_comment_length) ||
        !read(is, header.disk_number_start) ||
        !read(is, header.internal_file_attributes) ||
        !read(is, header.external_file_attributes) ||
        !read(is, header.relative_offset_of_local_header))
        return is;

    string file_name(header.file_name_length, '\0');
    if (!is.read(&file_name[0], header.file_name_length))
        return is;
    const auto use_utf8 = header.general_purpose_bit_flag & general_purpose_bit_flags::use_utf8;
    header.file_name = use_utf8 ? utf_to_utf<char_type>(file_name) : to_utf<char_type>(file_name, charset);
    if (header.extra_field_length) {
        header.extra_field.resize(header.extra_field_length);
        if (!read(is, &header.extra_field[0], header.extra_field_length))
            return is;
    }
    if (header.file_comment_length) {
        string file_comment(header.file_comment_length, '\0');
        if (!is.read(&file_comment[0], header.file_comment_length))
            return is;
        header.file_comment = use_utf8 ? utf_to_utf<char_type>(file_comment) : to_utf<char_type>(file_comment, charset);
    }
    return is;
}

ostream & zz::pkzip::operator << (ostream &os, const central_file_header &header)
{
    const auto use_utf8 = header.general_purpose_bit_flag & general_purpose_bit_flags::use_utf8;
    const auto file_name = use_utf8 ? utf_to_utf<char>(header.file_name) : from_utf(header.file_name, charset);
    if (file_name.size() > numeric_limits<decltype(header.file_name_length)>::max())
        throw runtime_error("too long file name: " + file_name);
    if (header.extra_field.size() > numeric_limits<decltype(header.extra_field_length)>::max())
        throw runtime_error("too long extra field: " + file_name);
    const auto file_comment = use_utf8 ? utf_to_utf<char>(header.file_comment) : from_utf(header.file_comment, charset);
    if (file_comment.size() > numeric_limits<decltype(header.file_comment_length)>::max())
        throw runtime_error("too long file comment: " + file_comment);

    write(os, header.signature);
    write(os, header.version_made_by);
    write(os, header.version_needed_to_extract);
    write(os, header.general_purpose_bit_flag);
    write(os, header.compression_method);
    write(os, header.last_mod_file_time);
    write(os, header.last_mod_file_date);
    write(os, header.crc32);
    write(os, header.compressed_size);
    write(os, header.uncompressed_size);
    write(os, static_cast<decltype(header.file_name_length)>(file_name.size()));
    write(os, static_cast<decltype(header.extra_field_length)>(header.extra_field.size()));
    write(os, static_cast<decltype(header.file_comment_length)>(file_comment.size()));
    write(os, header.disk_number_start);
    write(os, header.internal_file_attributes);
    write(os, header.external_file_attributes);
    write(os, header.relative_offset_of_local_header);

    write(os, file_name.data(), file_name.size());
    write(os, header.extra_field.data(), header.extra_field.size());
    write(os, file_comment.data(), file_comment.size());
    return os;
}

istream & zz::pkzip::operator >> (istream &is, end_of_central_directory_record &record)
{
    if (!read(is, record.signature) || !record ||
        !read(is, record.number_of_this_disk) ||
        !read(is, record.number_of_the_disk_with_the_start_of_the_central_directory) ||
        !read(is, record.total_number_of_entries_in_the_central_directory_on_this_disk) ||
        !read(is, record.total_number_of_entries_in_the_central_directory) ||
        !read(is, record.size_of_the_central_directory) ||
        !read(is, record.offset_of_start_of_central_directory_with_respect_to_the_starting_disk_number) ||
        !read(is, record.zip_file_comment_length))
        return is;

    if (record.zip_file_comment_length) {
        record.zip_file_comment.resize(record.zip_file_comment_length);
        if (!is.read(&record.zip_file_comment[0], record.zip_file_comment_length))
            return is;
    }
    return is;
}

ostream & zz::pkzip::operator << (ostream &os, const end_of_central_directory_record &record)
{
    if (record.zip_file_comment.size() > numeric_limits<decltype(record.zip_file_comment_length)>::max())
        throw runtime_error("too long zip file comment: " + record.zip_file_comment);

    write(os, record.signature);
    write(os, record.number_of_this_disk);
    write(os, record.number_of_the_disk_with_the_start_of_the_central_directory);
    write(os, record.total_number_of_entries_in_the_central_directory_on_this_disk);
    write(os, record.total_number_of_entries_in_the_central_directory);
    write(os, record.size_of_the_central_directory);
    write(os, record.offset_of_start_of_central_directory_with_respect_to_the_starting_disk_number);
    write(os, static_cast<decltype(record.zip_file_comment_length)>(record.zip_file_comment.size()));

    write(os, record.zip_file_comment.data(), record.zip_file_comment.size());
    return os;
}

ext: cc


text:
#include <Eigen/Dense>
#include <Eigen/Sparse>
#include <igl/read_triangle_mesh.h>
#include <igl/write_triangle_mesh.h>
#include <igl/writeOFF.h>
#include <igl/per_vertex_normals.h>
#include <igl/cotmatrix.h>
#include <igl/cotmatrix_entries.h>
#include <igl/adjacency_list.h>
#include <igl/triangle_triangle_adjacency.h>
#include <igl/barycenter.h>
#include <igl/massmatrix.h>
#include <igl/writeOBJ.h>
#include <fstream>
#include <cmath>
#include <array>

void findRotations(const Eigen::MatrixXd& N0,
                   const Eigen::MatrixXd& N1,
                   std::vector<Eigen::Matrix3d>& rot) {

    const auto n = N0.rows();
    rot.resize(n);

    for(int i = 0; i < n; ++i) {
        Eigen::Vector3d n1 = N0.row(i);
        Eigen::Vector3d n2 = N1.row(i);
        Eigen::Vector3d v = n1.cross(n2);
        const double c = n1.dot(n2);

        if(c > -1 + 1e-8) {
            const double coeff = 1 / (1 + c);
            Eigen::Matrix3d v_x;
            v_x << 0.0, -v(2), v(1), v(2), 0.0, -v(0), -v(1), v(0), 0.0;
            rot[i] = Eigen::Matrix3d::Identity() + v_x + coeff * v_x * v_x;
        } else{
            rot[i] = -Eigen::Matrix3d::Identity();
        }
    }
}

std::vector<std::vector<int>> collectNeighbours(const std::vector<std::vector<int>>& adj,
                                                const Eigen::MatrixXd& V,
                                                const Eigen::MatrixXd& N,
                                                const double r,
                                                const double nr) {

    std::vector<int> stack;
    std::vector<int> flag(V.rows(), -1);
    std::vector<std::vector<int>> result(V.rows());
    const double normalConeThreshold = cos(nr * M_PI / 180.);

    for(int i = 0; i < V.rows(); ++i) {

        stack.push_back(i);
        flag[i] = i;

        while(!stack.empty()) {
            auto id = stack.back();
            stack.pop_back();

            result[i].push_back(id);

            for (int j : adj[id]) {
                if(flag[j] != i && (V.row(i) - V.row(j)).norm() < r && (N.row(i).dot(N.row(j))) > normalConeThreshold) {
                    stack.push_back(j);
                    flag[j] = i;
                }
            }
        }
    }

    return result;
}

void fitNormals(const std::vector<std::vector<int>>& nbh,
                const Eigen::MatrixXd& V,
                const Eigen::MatrixXd& N,
                Eigen::MatrixXd& N2,
                const double cosineThreshold,
                const double sigma = 1.) {

    const auto nv = nbh.size();
    N2.resize(nv, 3);
    double angleThreshold = cosineThreshold * M_PI / 180.;

    for(int i = 0; i < nv; ++i) {

        const auto& nbi = nbh[i];

        Eigen::MatrixXd NN(nbi.size(), 3);

        for (int k = 0; k < nbi.size(); ++k) {
            NN.row(k) = N.row(nbi[k]);
        }

        Eigen::DiagonalMatrix<double, -1> W(nbi.size());

        if(sigma < 10.) {
            for(int i = 0; i < W.diagonal().size(); ++i) {
                double dot = NN.row(0).dot(NN.row(i));
                if (dot >= 1.){
                    W.diagonal()(i) = 1;
                } else if(dot < 0) {
                    W.diagonal()(i) = 0;
                } else {
                    W.diagonal()(i) = std::exp(-std::pow(acos(dot) / angleThreshold / sigma, 2));
                }
            }
        } else {
            W.diagonal().setOnes();
        }

        Eigen::JacobiSVD<Eigen::Matrix3d> svd(NN.transpose() * W * NN, Eigen::ComputeFullV);
        Eigen::Matrix3d frame = svd.matrixV();
        N2.row(i) = (frame.leftCols(2) * frame.leftCols(2).transpose() * N.row(i).transpose()).normalized();
    }
}

void assembleRHS(const Eigen::MatrixXd& C,
                 const Eigen::MatrixXd& V,
                 const Eigen::MatrixXi& F,
                 const std::vector<Eigen::Matrix3d>& R,
                 Eigen::MatrixXd& rhs) {

    const auto nv = V.rows();
    rhs.resize(nv, 3);
    rhs.setZero();

    for(int i = 0; i < F.rows(); ++i)  {
        for(int j = 0; j < 3; ++j)  {
            int v0 = F(i, (j + 1) % 3);
            int v1 = F(i, (j + 2) % 3);

            Eigen::Vector3d b = C(i,j) * R[i] * (V.row(v0) - V.row(v1)).transpose();
            rhs.row(v0) -= b.transpose();
            rhs.row(v1) += b.transpose();
        }
    }
}

std::vector<std::vector<int>> triangleAdjacency(const Eigen::MatrixXi& F, const size_t nv) {

    std::vector<std::vector<int>> vnbhs(nv);
    const auto nf = F.rows();

    for(int i = 0; i < nf; ++i) {
        for(int j = 0; j < 3; ++j) {
            vnbhs[F(i, j)].push_back(i);
        }
    }

    std::vector<int> flags(nf, -1);
    std::vector<std::vector<int>> ret(nf);

    for(int i = 0; i < nf; ++i) {
        for(int j = 0; j < 3; ++j) {
            for(int k : vnbhs[F(i, j)]) {
                if(k != i && flags[k] != i) {
                    ret[i].push_back(k);
                    flags[k] = i;
                }
            }
        }
    }

    return ret;
}

void center(Eigen::MatrixXd& V) {
    V.rowwise() -= V.colwise().mean();;
    V /= 2. * V.rowwise().norm().maxCoeff();
}

void gaussThinning(const std::string &mesh_folder,
                   const Eigen::MatrixXd &V_in,
                   const Eigen::MatrixXi &F,
                   Eigen::MatrixXd &V,
                   const int number_iterations = 100,
                   double minConeAngle = 2.5,
                   double smooth = 1e-5,
                   double start_angle = 25,
                   double radius = 0.1,
                   double sigma = 2.) {

    double coneAngle = start_angle;
    double r = radius;
    double eps = 1e-3;

    V = V_in;
    const auto nv = V.rows();
    center(V);

    igl::writeOFF(mesh_folder + "/normalized.off", V, F);

    Eigen::SparseMatrix<double> I(nv, nv);
    I.setIdentity();

    Eigen::MatrixXi TT;
    Eigen::MatrixXd B, b, C, N, N2;
    std::vector<Eigen::Matrix3d> rot;
    Eigen::SparseMatrix<double> L, M;
    std::vector<std::vector<int>> nbhs;
    Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> chol;

    auto tt = triangleAdjacency(F, nv);
    igl::triangle_triangle_adjacency(F, TT);
    igl::cotmatrix_entries(V, F, C);
    igl::cotmatrix(V, F, L);
    igl::massmatrix(V, F, igl::MASSMATRIX_TYPE_BARYCENTRIC, M);

    if(smooth) {
        chol.compute(-L + smooth * L.transpose() * L + eps * M);
    } else {
        chol.compute(-L + eps * M);
    }

    for(int k = 0; k < number_iterations; ++k) {

        igl::per_face_normals(V, F, N);
        igl::barycenter(V, F, B);

        nbhs = collectNeighbours(tt, B, N, r, coneAngle);
        if(coneAngle > minConeAngle) coneAngle *= .95;

        fitNormals(nbhs, V, N, N2, coneAngle, sigma);
        findRotations(N, N2, rot);
        assembleRHS(C, V, F, rot, b);

        V = chol.solve(eps * M * V - b);

        if (k % std::max(1, (number_iterations / 10)) == 0) {
            std::cout << "writing " + mesh_folder + ": " << k << "\n";
            igl::writeOFF(mesh_folder + "/out" + std::to_string(k) + ".off", V, F);
        }
    }

    return;
}

void runExperiment(std::string folder, std::string inputFile, std::string outputFile, const int iters, const double minAngle, const double start_angle = 25, const double radius = 0.1, const double smooth = 1e-5) {
    Eigen::MatrixXd V_in, V_out;
    Eigen::MatrixXi F;
    igl::read_triangle_mesh(folder + "/" + inputFile, V_in, F);
    gaussThinning(folder, V_in, F, V_out, iters, minAngle, smooth, start_angle, radius);
    igl::write_triangle_mesh(folder + "/" + outputFile, V_out, F);
}

int main(int argc, const char * argv[]) {

    if(argc < 6) {
        std::cout << "Need input file, output file, output directory, number of iterations and minimum search cone. Running default experiments..." << std::endl;

        /* run default experiments here .... */
        runExperiment("./examples/architecture", "input.off", "out.obj", 100, 2.5);

        runExperiment("./examples/boat", "input.off", "out.obj", 100, 2.5);

        runExperiment("./examples/bumpy", "input.off", "out.obj", 150, 7.5);

        runExperiment("./examples/bunny", "input.off", "out.obj", 500, 2.5);

        runExperiment("./examples/bunny_high", "input.off", "out.obj", 100, 2.5);

        runExperiment("./examples/bunny_small", "input.off", "out.obj", 500, 5.0);

        runExperiment("./examples/coffee", "input.off", "out.obj", 500, 2.5);

        runExperiment("./examples/cone", "input.off", "out.obj", 100, 2.5);

        runExperiment("./examples/cone_high", "input.off", "out.obj", 100, 2.5, 25, 0.015);

        runExperiment("./examples/curved_fold", "input.off", "out.obj", 100, 2.5);

        runExperiment("./examples/cylinder", "input.off", "out.obj", 300, 7.5);

        runExperiment("./examples/dog", "input.off", "out.obj", 100, 5.0);

        runExperiment("./examples/dome", "input.off", "out.obj", 100, 7.5);

        runExperiment("./examples/dress_high", "input.off", "out.obj", 100, 7.5);

        runExperiment("./examples/drill", "input.off", "out.obj", 100, 7.5);

        runExperiment("./examples/einstein", "input.off", "out.obj", 300, 7.5, 60, 0.015);

        runExperiment("./examples/face", "input.off", "out.obj", 100, 5.0);

        runExperiment("./examples/fandisk", "input.off", "out.obj", 1000, 5.0);

        runExperiment("./examples/fertility", "input.off", "out.obj", 100, 7.5);

        runExperiment("./examples/guitar", "input.off", "out.obj", 100, 2.5);

        runExperiment("./examples/lilium", "input.off", "out.obj", 100, 5.0);

        runExperiment("./examples/mask", "input.off", "out.obj", 500, 2.5);

        runExperiment("./examples/nut", "input.off", "out.obj", 100, 2.5);

        runExperiment("./examples/swing", "input.off", "out.obj", 500, 5.0);
    } else
    {
        std::string  infile = argv[1];
        std::string  outfile = argv[2];
        std::string  folder = argv[3];

        auto numIters = std::atoi(argv[4]);
        auto minAngle = std::stold(argv[5]);

        std::cout << "Processing " << infile << " with " << numIters << " iterations and mimimum cone angle " << minAngle << ". Output directory is " << folder << std::endl;

        runExperiment(folder, infile, outfile, numIters, minAngle);
    }

    return 0;
}

ext: cpp

max_stars_repo_name: PhHerholz/GaussImageThinning

max_issues_repo_name: PhHerholz/GaussImageThinning

max_forks_repo_name: PhHerholz/GaussImageThinning


text:
//  (C) Copyright John Maddock 2005.
//  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)

#ifndef BOOST_TR1_TYPE_TRAITS_HPP_INCLUDED
#  define BOOST_TR1_TYPE_TRAITS_HPP_INCLUDED
#  include <boost/tr1/detail/config.hpp>

#ifdef BOOST_HAS_TR1_TYPE_TRAITS
#  include BOOST_TR1_HEADER(type_traits)
#else
// Boost Type Traits:
#include <boost/type_traits.hpp>

namespace std { namespace tr1{

   using ::boost::integral_constant;
   using ::boost::true_type;
   using ::boost::false_type;
   using ::boost::is_void;
   using ::boost::is_integral;
   using ::boost::is_floating_point;
   using ::boost::is_array;
   using ::boost::is_pointer;
   using ::boost::is_reference;
   using ::boost::is_member_object_pointer;
   using ::boost::is_member_function_pointer;
   using ::boost::is_enum;
   using ::boost::is_union;
   using ::boost::is_class;
   using ::boost::is_function;
   using ::boost::is_arithmetic;
   using ::boost::is_fundamental;
   using ::boost::is_object;
   using ::boost::is_scalar;
   using ::boost::is_compound;
   using ::boost::is_member_pointer;
   using ::boost::is_const;
   using ::boost::is_volatile;
   using ::boost::is_pod;
   using ::boost::is_empty;
   using ::boost::is_polymorphic;
   using ::boost::is_abstract;
   using ::boost::has_trivial_constructor;
   using ::boost::has_trivial_copy;
   using ::boost::has_trivial_assign;
   using ::boost::has_trivial_destructor;
   using ::boost::has_nothrow_constructor;
   using ::boost::has_nothrow_copy;
   using ::boost::has_nothrow_assign;
   using ::boost::has_virtual_destructor;
   using ::boost::is_signed;
   using ::boost::is_unsigned;
   using ::boost::alignment_of;
   using ::boost::rank;
   using ::boost::extent;
   using ::boost::is_same;
   using ::boost::is_base_of;
   using ::boost::is_convertible;
   using ::boost::remove_const;
   using ::boost::remove_volatile;
   using ::boost::remove_cv;
   using ::boost::add_const;
   using ::boost::add_volatile;
   using ::boost::add_cv;
   using ::boost::remove_reference;
   using ::boost::add_reference;
   using ::boost::remove_extent;
   using ::boost::remove_all_extents;
   using ::boost::remove_pointer;
   using ::boost::add_pointer;
   using ::boost::aligned_storage;

} }

#endif

#endif

max_stars_repo_name: dstrigl/mcotf

max_issues_repo_name: dstrigl/mcotf

max_forks_repo_name: dstrigl/mcotf


text:
/*
 * Copyright 2012 Open Source Robotics Foundation
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 *
*/
/* Desc: The base joint class
 * Author: Nate Koenig, Andrew Howard
 * Date: 21 May 2003
 */

#ifndef _JOINT_HH_
#define _JOINT_HH_

#include <string>

#include <boost/any.hpp>

#include "gazebo/common/Event.hh"
#include "gazebo/common/Events.hh"
#include "gazebo/math/Angle.hh"
#include "gazebo/math/Vector3.hh"
#include "gazebo/msgs/MessageTypes.hh"

#include "gazebo/physics/JointState.hh"
#include "gazebo/physics/Base.hh"
#include "gazebo/physics/JointWrench.hh"

/// \brief maximum number of axis per joint anticipated.
/// Currently, this is 2 as 3-axis joints (e.g. ball)
/// actuation, control is not there yet.
#define MAX_JOINT_AXIS 2

namespace gazebo
{
  namespace physics
  {
    /// \addtogroup gazebo_physics
    /// \{

    /// \class Joint Joint.hh physics/physics.hh
    /// \brief Base class for all joints
    class Joint : public Base
    {
      /// \enum Attribute
      /// \brief Joint attribute types.
      public: enum Attribute
              {
                /// \brief Fudge factor.
                FUDGE_FACTOR,

                /// \brief Suspension error reduction parameter.
                SUSPENSION_ERP,

                /// \brief Suspension constraint force mixing.
                SUSPENSION_CFM,

                /// \brief Stop limit error reduction parameter.
                STOP_ERP,

                /// \brief Stop limit constraint force mixing.
                STOP_CFM,

                /// \brief Error reduction parameter.
                ERP,

                /// \brief Constraint force mixing.
                CFM,

                /// \brief Maximum force.
                FMAX,

                /// \brief Velocity.
                VEL,

                /// \brief High stop angle.
                HI_STOP,

                /// \brief Low stop angle.
                LO_STOP
              };

      /// \brief Constructor
      /// \param[in] Joint parent
      public: explicit Joint(BasePtr _parent);

      /// \brief Destructor
      public: virtual ~Joint();

      /// \brief Set pose, parent and child links of a physics::Joint
      /// \param[in] _parent Parent link.
      /// \param[in] _child Child link.
      /// \param[in] _pose Pose containing Joint Anchor offset from child link.
      public: void Load(LinkPtr _parent, LinkPtr _child,
                        const math::Pose &_pose);

      /// \brief Set parent and child links of a physics::Joint and its
      /// anchor offset position.
      /// This funciton is deprecated, use
      /// Load(LinkPtr _parent, LinkPtr _child, const math::Pose &_pose)
      /// \param[in] _parent Parent link.
      /// \param[in] _child Child link.
      /// \param[in] _pos Joint Anchor offset from child link.
      public: void Load(LinkPtr _parent, LinkPtr _child,
                        const math::Vector3 &_pos) GAZEBO_DEPRECATED(1.5);

      /// \brief Load physics::Joint from a SDF sdf::Element.
      /// \param[in] _sdf SDF values to load from.
      public: virtual void Load(sdf::ElementPtr _sdf);

      /// \brief Initialize a joint.
      public: virtual void Init();

      /// \brief Update the joint.
      public: void Update();

      /// \brief Update the parameters using new sdf values.
      /// \param[in] _sdf SDF values to update from.
      public: virtual void UpdateParameters(sdf::ElementPtr _sdf);

      /// \brief Reset the joint.
      public: virtual void Reset();

      /// \brief Set the joint state.
      /// \param[in] _state Joint state
      public: void SetState(const JointState &_state);

      /// \brief Set the model this joint belongs too.
      /// \param[in] _model Pointer to a model.
      public: void SetModel(ModelPtr _model);

      /// \brief Get the link to which the joint is attached according
      /// the _index.
      /// \param[in] _index Index of the link to retreive.
      /// \return Pointer to the request link. NULL if the index was
      /// invalid.
      public: virtual LinkPtr GetJointLink(int _index) const = 0;

      /// \brief Determines of the two bodies are connected by a joint.
      /// \param[in] _one First link.
      /// \param[in] _two Second link.
      /// \return True if the two links are connected by a joint.
      public: virtual bool AreConnected(LinkPtr _one, LinkPtr _two) const = 0;

      /// \brief Attach the two bodies with this joint.
      /// \param[in] _parent Parent link.
      /// \param[in] _child Child link.
      public: virtual void Attach(LinkPtr _parent, LinkPtr _child);

      /// \brief Detach this joint from all links.
      public: virtual void Detach();

      /// \brief Set the axis of rotation.
      /// \param[in] _index Index of the axis to set.
      /// \param[in] _axis Vector in world frame of axis direction
      ///                  (must have length greater than zero).
      public: virtual void SetAxis(int _index, const math::Vector3 &_axis) = 0;

      /// \brief Set the joint damping.
      /// \param[in] _index Index of the axis to set.
      /// \param[in] _damping Damping value for the axis.
      public: virtual void SetDamping(int _index, double _damping) = 0;

      /// \brief Callback to apply damping force to joint.
      public: virtual void ApplyDamping();

      /// \brief Connect a boost::slot the the joint update signal.
      /// \param[in] _subscriber Callback for the connection.
      /// \return Connection pointer, which must be kept in scope.
      public: template<typename T>
              event::ConnectionPtr ConnectJointUpdate(T _subscriber)
              {return jointUpdate.Connect(_subscriber);}

      /// \brief Disconnect a boost::slot the the joint update signal.
      /// \param[in] _conn Connection to disconnect.
      public: void DisconnectJointUpdate(event::ConnectionPtr &_conn)
              {jointUpdate.Disconnect(_conn);}

      /// \brief Get the axis of rotation.
      /// \param[in] _index Index of the axis to get.
      /// \return Axis value for the provided index.
      public: math::Vector3 GetLocalAxis(int _index) const;

      /// \brief Get the axis of rotation in global cooridnate frame.
      /// \param[in] _index Index of the axis to get.
      /// \return Axis value for the provided index.
      public: virtual math::Vector3 GetGlobalAxis(int _index) const = 0;

      /// \brief Set the anchor point.
      /// \param[in] _index Indx of the axis.
      /// \param[in] _anchor Anchor value.
      public: virtual void SetAnchor(int _index,
                                     const math::Vector3 &_anchor) = 0;

      /// \brief Get the anchor point.
      /// \param[in] _index Index of the axis.
      /// \return Anchor value for the axis.
      public: virtual math::Vector3 GetAnchor(int _index) const = 0;

      /// \brief Set the high stop of an axis(index).
      /// \param[in] _index Index of the axis.
      /// \param[in] _angle High stop angle.
      public: virtual void SetHighStop(int _index,
                                       const math::Angle &_angle);

      /// \brief Set the low stop of an axis(index).
      /// \param[in] _index Index of the axis.
      /// \param[in] _angle Low stop angle.
      public: virtual void SetLowStop(int _index,
                                      const math::Angle &_angle);

      /// \brief Get the high stop of an axis(index).
      /// This function is replaced by GetUpperLimit(unsigned int).
      /// If you are interested in getting the value of dParamHiStop*,
      /// use GetAttribute(hi_stop, _index)
      /// \param[in] _index Index of the axis.
      /// \return Angle of the high stop value.
      public: virtual math::Angle GetHighStop(int _index) = 0;

      /// \brief Get the low stop of an axis(index).
      /// This function is replaced by GetLowerLimit(unsigned int).
      /// If you are interested in getting the value of dParamHiStop*,
      /// use GetAttribute(hi_stop, _index)
      /// \param[in] _index Index of the axis.
      /// \return Angle of the low stop value.
      public: virtual math::Angle GetLowStop(int _index) = 0;

      /// \brief Get the effort limit on axis(index).
      /// \param[in] _index Index of axis, where 0=first axis and 1=second axis
      /// \return Effort limit specified in SDF
      public: virtual double GetEffortLimit(int _index);

      /// \brief Get the velocity limit on axis(index).
      /// \param[in] _index Index of axis, where 0=first axis and 1=second axis
      /// \return Velocity limit specified in SDF
      public: virtual double GetVelocityLimit(int _index);

      /// \brief Set the velocity of an axis(index).
      /// \param[in] _index Index of the axis.
      /// \param[in] _vel Velocity.
      public: virtual void SetVelocity(int _index, double _vel) = 0;

      /// \brief Get the rotation rate of an axis(index)
      /// \param[in] _index Index of the axis.
      /// \return The rotaional velocity of the joint axis.
      public: virtual double GetVelocity(int _index) const = 0;

      /// \brief Set the force applied to this physics::Joint.
      /// Note that the unit of force should be consistent with the rest
      /// of the simulation scales.  E.g.  if you are using
      /// metric units, the unit for force is Newtons.  If using
      /// imperial units (sorry), then unit of force is lb-force
      /// not (lb-mass), etc.
      /// \param[in] _index Index of the axis.
      /// \param[in] _force Force value.
      public: virtual void SetForce(int _index, double _force);

      /// \brief @todo: not yet implemented.
      /// Get the forces applied at this Joint.
      /// Note that the unit of force should be consistent with the rest
      /// of the simulation scales.  E.g.  if you are using
      /// metric units, the unit for force is Newtons.  If using
      /// imperial units (sorry), then unit of force is lb-force
      /// not (lb-mass), etc.
      /// \param[in] _index Index of the axis.
      /// \return The force applied to an axis.
      public: virtual double GetForce(int _index) GAZEBO_DEPRECATED(1.5);

      /// \brief @todo: not yet implemented.
      /// Get the forces applied at this Joint.
      /// Note that the unit of force should be consistent with the rest
      /// of the simulation scales.  E.g.  if you are using
      /// metric units, the unit for force is Newtons.  If using
      /// imperial units (sorry), then unit of force is lb-force
      /// not (lb-mass), etc.
      /// \param[in] _index Index of the axis.
      /// \return The force applied to an axis.
      public: virtual double GetForce(unsigned int _index);

      /// \brief get internal force and torque values at a joint
      /// Note that you must set
      ///   <provide_feedback>true<provide_feedback>
      /// in the joint sdf to use this.
      /// \param[in] _index Force and torque on child link if _index = 0
      /// and on parent link of _index = 1
      /// \return The force and torque at the joint
      public: virtual JointWrench GetForceTorque(int _index)
        GAZEBO_DEPRECATED(1.5) = 0;

      /// \brief get internal force and torque values at a joint
      /// Note that you must set
      ///   <provide_feedback>true<provide_feedback>
      /// in the joint sdf to use this.
      /// \param[in] _index Force and torque on child link if _index = 0
      /// and on parent link of _index = 1
      /// \return The force and torque at the joint
      public: virtual JointWrench GetForceTorque(unsigned int _index) = 0;

      /// \brief Set the max allowed force of an axis(index).
      /// Note that the unit of force should be consistent with the rest
      /// of the simulation scales.  E.g.  if you are using
      /// metric units, the unit for force is Newtons.  If using
      /// imperial units (sorry), then unit of force is lb-force
      /// not (lb-mass), etc.
      /// \param[in] _index Index of the axis.
      /// \param[in] _force Maximum force that can be applied to the axis.
      public: virtual void SetMaxForce(int _index, double _force) = 0;

      /// \brief Get the max allowed force of an axis(index).
      /// Note that the unit of force should be consistent with the rest
      /// of the simulation scales.  E.g.  if you are using
      /// metric units, the unit for force is Newtons.  If using
      /// imperial units (sorry), then unit of force is lb-force
      /// not (lb-mass), etc.
      /// \param[in] _index Index of the axis.
      /// \return The maximum force.
      public: virtual double GetMaxForce(int _index) = 0;

      /// \brief Get the angle of rotation of an axis(index)
      /// \param[in] _index Index of the axis.
      /// \return Angle of the axis.
      public: math::Angle GetAngle(int _index) const;

      /// \brief Get the angle count.
      /// \return The number of DOF for the joint.
      public: virtual unsigned int GetAngleCount() const = 0;

      /// \brief If the Joint is static, Gazebo stores the state of
      /// this Joint as a scalar inside the Joint class, so
      /// this call will NOT move the joint dynamically for a static Model.
      /// But if this Model is not static, then it is updated dynamically,
      /// all the conencted children Link's are moved as a result of the
      /// Joint angle setting.  Dynamic Joint angle update is accomplished
      /// by calling JointController::SetJointPosition.
      /// \param[in] _index Index of the axis.
      /// \param[in] _angle Angle to set the joint to.
      public: void SetAngle(int _index, math::Angle _angle);

      /// \brief Get the forces applied to the center of mass of a physics::Link
      /// due to the existence of this Joint.
      /// Note that the unit of force should be consistent with the rest
      /// of the simulation scales.  E.g.  if you are using
      /// metric units, the unit for force is Newtons.  If using
      /// imperial units (sorry), then unit of force is lb-force
      /// not (lb-mass), etc.
      /// \param[in] index The index of the link(0 or 1).
      /// \return Force applied to the link.
      public: virtual math::Vector3 GetLinkForce(unsigned int _index) const = 0;

      /// \brief Get the torque applied to the center of mass of a physics::Link
      /// due to the existence of this Joint.
      /// Note that the unit of torque should be consistent with the rest
      /// of the simulation scales.  E.g.  if you are using
      /// metric units, the unit for force is Newtons-Meters.  If using
      /// imperial units (sorry), then unit of force is lb-force-inches
      /// not (lb-mass-inches), etc.
      /// \param[in] index The index of the link(0 or 1)
      /// \return Torque applied to the link.
      public: virtual math::Vector3 GetLinkTorque(
                  unsigned int _index) const = 0;

      /// \brief Set a non-generic parameter for the joint.
      /// replaces SetAttribute(Attribute, int, double)
      /// \param[in] _key String key.
      /// \param[in] _index Index of the axis.
      /// \param[in] _value Value of the attribute.
      public: virtual void SetAttribute(const std::string &_key, int _index,
                                        const boost::any &_value) = 0;

      /// \brief Get a non-generic parameter for the joint.
      /// \param[in] _key String key.
      /// \param[in] _index Index of the axis.
      /// \param[in] _value Value of the attribute.
      public: virtual double GetAttribute(const std::string &_key,
                                                unsigned int _index) = 0;

      /// \brief Get the child link
      /// \return Pointer to the child link.
      public: LinkPtr GetChild() const;

      /// \brief Get the parent link.
      /// \return Pointer to the parent link.
      public: LinkPtr GetParent() const;

      /// \brief Fill a joint message.
      /// \param[out] _msg Message to fill with this joint's properties.
      public: void FillMsg(msgs::Joint &_msg);

      /// \brief Accessor to inertia ratio across this joint.
      /// \param[in] _index Joint axis index.
      /// \return returns the inertia ratio across specified joint axis.
      public: double GetInertiaRatio(unsigned int _index) const;

      /// \brief:  get the joint upper limit
      /// (replaces GetLowStop and GetHighStop)
      /// \param[in] _index Index of the axis.
      /// \return Upper limit of the axis.
      public: math::Angle GetLowerLimit(unsigned int _index) const
        {
          if (_index < this->GetAngleCount())
            return this->lowerLimit[_index];

          gzwarn << "requesting lower limit of joint index out of bound\n";
          return math::Angle();
        }

      /// \brief:  get the joint lower limit
      /// (replacee GetLowStop and GetHighStop)
      /// \param[in] _index Index of the axis.
      /// \return Upper limit of the axis.
      public: math::Angle GetUpperLimit(unsigned int _index) const
        {
          if (_index < this->GetAngleCount())
            return this->upperLimit[_index];

          gzwarn << "requesting upper limit of joint index out of bound\n";
          return math::Angle();
        }

      /// \brief Get the angle of an axis helper function.
      /// \param[in] _index Index of the axis.
      /// \return Angle of the axis.
      protected: virtual math::Angle GetAngleImpl(int _index) const = 0;

      /// \brief Helper function to load a joint.
      /// \param[in] _pose Pose of the anchor.
      private: void LoadImpl(const math::Pose &_pose);

      /// \brief Helper function to load a joint.
      /// This function is deprecated, use LoadImpl(math::Pose &)
      /// \param[in] _pos Position of the anchor.
      private: void LoadImpl(const math::Vector3 &_pos) GAZEBO_DEPRECATED(1.5);

      /// \brief Computes inertiaRatio for this joint during Joint::Init
      /// The inertia ratio for each joint between [1, +inf] gives a sense
      /// of how well this model will perform in iterative LCP methods.
      private: void ComputeInertiaRatio();

      /// \brief The first link this joint connects to
      protected: LinkPtr childLink;

      /// \brief The second link this joint connects to
      protected: LinkPtr parentLink;

      /// \brief Pointer to the parent model.
      protected: ModelPtr model;

      /// \brief Anchor pose.  This is the xyz offset of the joint frame from
      /// child frame specified in the parent link frame
      protected: math::Vector3 anchorPos;

      /// \brief Anchor pose specified in SDF <joint><pose> tag.
      /// AnchorPose is the transform from child link frame to joint frame
      /// specified in the child link frame.
      /// AnchorPos is more relevant in normal usage, but sometimes,
      /// we do need this (e.g. GetForceTorque and joint visualization).
      protected: math::Pose anchorPose;

      /// \brief Anchor link.
      protected: LinkPtr anchorLink;

      /// \brief Joint update event.
      private: event::EventT<void ()> jointUpdate;

      /// \brief joint dampingCoefficient
      protected: double dampingCoefficient;

      /// \brief Angle used when the joint is paret of a static model.
      private: math::Angle staticAngle;

      /// \brief apply damping for adding viscous damping forces on updates
      protected: gazebo::event::ConnectionPtr applyDamping;

      /// \brief Save force applied by user
      /// This plus the joint feedback (joint contstraint forces) is the
      /// equivalent of simulated force torque sensor reading
      /// Allocate a 2 vector in case hinge2 joint is used.
      protected: double forceApplied[MAX_JOINT_AXIS];

      /// \brief Store Joint effort limit as specified in SDF
      protected: double effortLimit[MAX_JOINT_AXIS];

      /// \brief Store Joint velocity limit as specified in SDF
      protected: double velocityLimit[MAX_JOINT_AXIS];

      /// \brief Store Joint inertia ratio.  This is a measure of how well
      /// this model behaves using interative LCP solvers.
      protected: double inertiaRatio[MAX_JOINT_AXIS];

      /// \brief Store Joint position lower limit as specified in SDF
      protected: math::Angle lowerLimit[MAX_JOINT_AXIS];

      /// \brief Store Joint position upper limit as specified in SDF
      protected: math::Angle upperLimit[MAX_JOINT_AXIS];

      /// \brief option to use CFM damping
      protected: bool useCFMDamping;
    };
    /// \}
  }
}
#endif

ext: hh

max_stars_repo_name: thomas-moulard/gazebo-deb

max_issues_repo_name: thomas-moulard/gazebo-deb

max_forks_repo_name: thomas-moulard/gazebo-deb


text:
#include <facter/logging/logging.hpp>
#include <leatherman/util/environment.hpp>
#include <leatherman/logging/logging.hpp>
#include <leatherman/locale/locale.hpp>
#include <boost/filesystem.hpp>

// Mark string for translation (alias for leatherman::locale::format)
using leatherman::locale::_;

using namespace std;
using namespace leatherman::util;
namespace lm = leatherman::logging;

#ifdef LEATHERMAN_USE_LOCALES
static void setup_logging_internal(ostream& os, bool use_locale)
{
    // Initialize boost filesystem's locale to a UTF-8 default.
    // Logging gets setup the same way via the default 2nd argument.
    if (use_locale) {
        boost::filesystem::path::imbue(leatherman::locale::get_locale());
    }
    lm::setup_logging(os, "", PROJECT_NAME, use_locale);
}
#else
static void setup_logging_internal(ostream& os, bool)
{
    lm::setup_logging(os, "", PROJECT_NAME, false);
}
#endif

static const char* lc_vars[] = {
    "LC_CTYPE",
    "LC_NUMERIC",
    "LC_TIME",
    "LC_COLLATE",
    "LC_MONETARY",
    "LC_MESSAGES",
    "LC_PAPER",
    "LC_ADDRESS",
    "LC_TELEPHONE",
    "LC_MEASUREMENT",
    "LC_IDENTIFICATION",
    "LC_ALL"
};

namespace facter { namespace logging {
    istream& operator>>(istream& in, level& lvl)
    {
        lm::log_level lm_level;
        in >> lm_level;
        lvl = static_cast<level>(lm_level);
        return in;
    }

    ostream& operator<<(ostream& os, level lvl)
    {
        os << static_cast<lm::log_level>(lvl);
        return os;
    }

    void setup_logging(ostream& os)
    {
        try {
            setup_logging_internal(os, true);
        } catch (exception const&) {
            for (auto var : lc_vars) {
                environment::clear(var);
            }
            environment::set("LANG", "C");
            environment::set("LC_ALL", "C");
            try {
                setup_logging_internal(os, true);
                // We can't log the issue until after logging is setup. Rely on the exception for any other
                // error reporting.
                log(level::warning, "locale environment variables were bad; continuing with LANG=C LC_ALL=C");
            } catch (exception const& e) {
                // If we fail again even with a clean environment, we
                // need to signal to our consumer that things went
                // sideways.
                //
                // Since logging is busted, we raise an exception that
                // signals to the consumer that a special action must
                // be taken to alert the user.
                try {
                    setup_logging_internal(os, false);
                    log(level::warning, "Could not initialize locale, even with LC_* variables cleared. Continuing without localization support");
                } catch (exception const& e) {
                    throw locale_error(_("could not initialize logging, even with locale variables reset to LANG=C LC_ALL=C: {1}", e.what()));
                }
            }
        }
    }

    void set_level(level lvl)
    {
        lm::set_level(static_cast<lm::log_level>(lvl));
    }

    level get_level()
    {
        return static_cast<level>(lm::get_level());
    }

    void set_colorization(bool color)
    {
        lm::set_colorization(color);
    }

    bool get_colorization()
    {
        return lm::get_colorization();
    }

    bool is_enabled(level lvl)
    {
        return lm::is_enabled(static_cast<lm::log_level>(lvl));
    }

    bool error_logged()
    {
        return lm::error_has_been_logged();
    }

    void clear_logged_errors()
    {
        lm::clear_error_logged_flag();
    }

    std::string translate(std::string const& msg)
    {
#ifdef LEATHERMAN_I18N
        return leatherman::locale::translate(msg);
#else
        // When locales are disabled, use boost::format, which expects %N% style formatting
        static const boost::regex match{"\\{(\\d+)\\}"};
        static const std::string repl{"%\\1%"};
        return boost::regex_replace(leatherman::locale::translate(msg), match, repl);
#endif
    }

    void log(level lvl, string const& message)
    {
        lm::log(LOG_NAMESPACE, static_cast<lm::log_level>(lvl), 0, translate(message));
    }

    void colorize(ostream &os, level lvl)
    {
        lm::colorize(os, static_cast<lm::log_level>(lvl));
    }

}}  // namespace facter::logging

max_stars_repo_name: mihaibuzgau/facter

max_issues_repo_name: mihaibuzgau/facter

max_forks_repo_name: mihaibuzgau/facter


text:
/**
 * \file dcs/math/stats/distribution/normal.hpp
 *
 * \brief The Normal probability distribution.
 *
 * \author Marco Guazzone (marco.guazzone@gmail.com)
 *
 * <hr/>
 *
 * Copyright 2009 Marco Guazzone (marco.guazzone@gmail.com)
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

#ifndef DCS_MATH_STATS_DISTRIBUTION_NORMAL_HPP
#define DCS_MATH_STATS_DISTRIBUTION_NORMAL_HPP

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

#if !DCS_DETAIL_CONFIG_BOOST_CHECK_VERSION(103500) // 1.35
# 	error "Required Boost library version >= 1.35"
#endif

#include <boost/math/distributions/normal.hpp>
#include <boost/random/normal_distribution.hpp>
#include <boost/random/variate_generator.hpp>
#include <cstddef>
#include <dcs/math/policies/policy.hpp>
#include <iostream>
#include <vector>

namespace dcs { namespace math { namespace stats {

using ::std::size_t;

/**
 * The Normal distribution with parameter \f$\mu\f$ (the mean) and \f$\sigma\f$
 * (the standard deviation).
 *
 * \tparam RealT The type used for real numbers.
 * \tparam PolicyT The policy type.
 *
 * \author Marco Guazzone (marco.guazzone@gmail.com)
 */
template < typename RealT=double, typename PolicyT=::dcs::math::policies::policy<> >
class normal_distribution
{
	public: typedef RealT support_type;
	public: typedef RealT value_type;
	public: typedef PolicyT policy_type;

	public: explicit normal_distribution(support_type mean=0, support_type stddev=1)
		: dist_(mean,stddev)
//		: mean_(mean),
//		  stddev_(stddev)
	{
		// empty
	}

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

	/**
	 * \brief Generate a vector of random numbers distributed according to this
	 * normal distribution.
	 *
	 * \param rng A uniform random number generator.
	 * \param n The number of random numbers to generate.
	 * \return A random number distributed according to this normal
	 * distribution.
	 *
	 * A \c normal random number distribution produces random numbers * \f$x\f$
	 * distributed according to the probability density function:
	 * \f[
	 *   \Pr(x|\mu,\sigma) = \frac{1}{\sigma\sqrt{2\pi}}\exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right)
	 * \f]
	 */
	public: template <typename UniformRandomGeneratorT>
		support_type rand(UniformRandomGeneratorT& rng) const
	{
		typedef ::boost::normal_distribution<support_type> rdist_type;
		typedef ::boost::variate_generator<UniformRandomGeneratorT&, rdist_type> variate_type;

		return variate_type(rng, rdist_type(dist_.mean(), dist_.standard_deviation()))();
	}

	/**
	 * \brief Generate a vector of random numbers distributed according to this
	 * normal distribution.
	 *
	 * \param rng A uniform random number generator.
	 * \param n The number of random numbers to generate.
	 * \return A vector of random numbers distributed according to this
	 * normal distribution.
	 *
	 * A \c normal random number distribution produces random numbers * \f$x\f$
	 * distributed according to the probability density function:
	 * \f[
	 *   \Pr(x|\mu,\sigma) = \frac{1}{\sigma\sqrt{2\pi}}\exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right)
	 * \f]
	 */
	public: template <typename UniformRandomGeneratorT>
		::std::vector<support_type> rand(UniformRandomGeneratorT& rng, size_t n)
	{
		typedef ::boost::normal_distribution<support_type> rdist_type;
        typedef ::boost::variate_generator<UniformRandomGeneratorT&, rdist_type> variate_type;

		::std::vector<support_type> rnds(n);

        for ( ; n > 0; --n)
		{
			rnds.push_back(variate_type(rng, rdist_type(dist_.mean(), dist_.standard_deviation()))());
		}

		return rnds;
	}

	public: support_type mean() const
	{
		return dist_.mean();
	}

	public: support_type stddev() const
	{
		return dist_.standard_deviation();
	}

	public: support_type location() const
	{
		return dist_.location();
	}

	public: support_type scale() const
	{
		return dist_.scale();
	}

	public: support_type quantile(value_type p) const
	{
		return ::boost::math::quantile(dist_, p);
	}

	//private: support_type mean_;
	//private: support_type stddev_;
	private: ::boost::math::normal_distribution<value_type,policy_type> dist_;
};

template <
	typename CharT,
	typename CharTraitsT,
	typename RealT,
	typename PolicyT
>
::std::basic_ostream<CharT,CharTraitsT>& operator<<(::std::basic_ostream<CharT,CharTraitsT>& os, normal_distribution<RealT,PolicyT> const& dist)
{
	return os << "Normal("
			  << "mean=" <<  dist.mean()
			  << ", stddev=" <<  dist.stddev()
			  << ")";
}

}}} // Namespace dcs::math::stats

#endif // DCS_MATH_STATS_DISTRIBUTION_NORMAL_HPP

max_stars_repo_name: sguazt/dcsxx-commons

max_issues_repo_name: sguazt/dcsxx-commons

max_forks_repo_name: sguazt/dcsxx-commons


text:
//  (C) Copyright John Maddock 2006.
//  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)

#include "test_0F1.hpp"

#include <boost/multiprecision/cpp_bin_float.hpp>
#include <boost/multiprecision/cpp_dec_float.hpp>

void expected_results()
{
   //
   // Define the max and mean errors expected for
   // various compilers and platforms.
   //
   const char* largest_type;
#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
   if(boost::math::policies::digits<double, boost::math::policies::policy<> >() == boost::math::policies::digits<long double, boost::math::policies::policy<> >())
   {
      largest_type = "(long\\s+)?double|real_concept|cpp_bin_float_quad|dec_40";
   }
   else
   {
      largest_type = "long double|real_concept|cpp_bin_float_quad|dec_40";
   }
#else
   largest_type = "(long\\s+)?double|cpp_bin_float_quad|dec_40";
#endif

   add_expected_result(
      ".*",                          // compiler
      ".*",                          // stdlib
      ".*",                          // platform
      largest_type,                  // test type(s)
      "Integer.*",                   // test data group
      ".*", 300, 100);               // test function
   add_expected_result(
      ".*",                          // compiler
      ".*",                          // stdlib
      ".*",                          // platform
      largest_type,                  // test type(s)
      "Real.*",                   // test data group
      ".*", 2000, 1000);               // test function

   add_expected_result(
      ".*",                          // compiler
      ".*",                          // stdlib
      ".*",                          // platform
      largest_type,                  // test type(s)
      "Large.*",                   // test data group
      ".*", 400, 100);               // test function

   //
   // Finish off by printing out the compiler/stdlib/platform names,
   // we do this to make it easier to mark up expected error rates.
   //
   std::cout << "Tests run with " << BOOST_COMPILER << ", "
      << BOOST_STDLIB << ", " << BOOST_PLATFORM << std::endl;
}

BOOST_AUTO_TEST_CASE( test_main )
{
   expected_results();
   BOOST_MATH_CONTROL_FP;

#if !defined(TEST) || (TEST == 1)
#ifndef BOOST_MATH_BUGGY_LARGE_FLOAT_CONSTANTS
   test_spots(0.0F, "float");
#endif
   test_spots(0.0, "double");
#endif

#if !defined(TEST) || (TEST == 2)
#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
   test_spots(0.0L, "long double");
#ifndef BOOST_MATH_NO_REAL_CONCEPT_TESTS
   test_spots(boost::math::concepts::real_concept(0.1), "real_concept");
#endif
#endif
#endif
#if !defined(TEST) || (TEST == 3)
   test_spots(boost::multiprecision::cpp_bin_float_quad(), "cpp_bin_float_quad");
#endif
#if !defined(TEST) || (TEST == 4)
   typedef boost::multiprecision::number<boost::multiprecision::cpp_dec_float<40> > dec_40;
   test_spots(dec_40(), "dec_40");
#endif
}

max_stars_repo_name: oleg-alexandrov/math

max_issues_repo_name: qingkouwei/mediaones

max_forks_repo_name: qingkouwei/mediaones