text:
// Boost.Geometry Index
// Additional tests

// Copyright (c) 2011-2015 Adam Wulkiewicz, Lodz, Poland.

// Use, modification and distribution is 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)

#define BOOST_GEOMETRY_INDEX_DETAIL_EXPERIMENTAL

#include <iostream>

#include <boost/chrono.hpp>
#include <boost/foreach.hpp>
#include <boost/random.hpp>
#include <boost/range/algorithm/copy.hpp>

#include <boost/geometry.hpp>
#include <boost/geometry/index/rtree.hpp>
#include <boost/geometry/geometries/linestring.hpp>
#include <boost/geometry/geometries/segment.hpp>

#include <boost/geometry/index/detail/rtree/utilities/are_levels_ok.hpp>
#include <boost/geometry/index/detail/rtree/utilities/are_boxes_ok.hpp>

namespace bg = boost::geometry;
namespace bgi = bg::index;

typedef bg::model::point<double, 2, bg::cs::cartesian> P;
typedef bg::model::box<P> B;
typedef bg::model::linestring<P> LS;
typedef bg::model::segment<P> S;
//typedef P V;
typedef B V;
//typedef S V;
//#define SEGMENT_INDEXABLE

template <typename V>
struct generate_value {};

template <>
struct generate_value<B>
{
    static inline B apply(float x, float y) { return B(P(x - 0.5f, y - 0.5f), P(x + 0.5f, y + 0.5f)); }
};

template <>
struct generate_value<S>
{
    static inline S apply(float x, float y) { return S(P(x - 0.5f, y - 0.5f), P(x + 0.5f, y + 0.5f)); }
};

template <>
struct generate_value<P>
{
    static inline P apply(float x, float y) { return P(x, y); }
};

//#include <boost/geometry/extensions/nsphere/nsphere.hpp>
//typedef bg::model::nsphere<P, double> NS;
//typedef NS V;
//
//template <>
//struct generate_value<NS>
//{
//    static inline NS apply(float x, float y) { return NS(P(x, y), 0.5); }
//};

template <typename I1, typename I2, typename O>
void mycopy(I1 first, I2 last, O o)
{
    for ( ; first != last ; ++o, ++first )
        *o = *first;
}

//#define BOOST_GEOMETRY_INDEX_BENCHMARK_DEBUG

int main()
{
    typedef boost::chrono::thread_clock clock_t;
    typedef boost::chrono::duration<float> dur_t;

#ifndef BOOST_GEOMETRY_INDEX_BENCHMARK_DEBUG
    size_t values_count = 1000000;
    size_t queries_count = 100000;
    size_t nearest_queries_count = 20000;
    unsigned neighbours_count = 10;
    size_t path_queries_count = 2000;
    size_t path_queries_count2 = 20000;
    unsigned path_values_count = 10;
#else
    size_t values_count = 1000;
    size_t queries_count = 1;
    size_t nearest_queries_count = 1;
    unsigned neighbours_count = 10;
    size_t path_queries_count = 1;
    size_t path_queries_count2 = 1;
    unsigned path_values_count = 10;
#endif

    float max_val = static_cast<float>(values_count / 2);
    std::vector< std::pair<float, float> > coords;
    std::vector<V> values;

    //randomize values
    {
        boost::mt19937 rng;
        //rng.seed(static_cast<unsigned int>(std::time(0)));
        boost::uniform_real<float> range(-max_val, max_val);
        boost::variate_generator<boost::mt19937&, boost::uniform_real<float> > rnd(rng, range);

        coords.reserve(values_count);

        std::cout << "randomizing data\n";
        for ( size_t i = 0 ; i < values_count ; ++i )
        {
            float x = rnd();
            float y = rnd();
            coords.push_back(std::make_pair(x, y));
            values.push_back(generate_value<V>::apply(x, y));
        }
        std::cout << "randomized\n";
    }

    typedef bgi::rtree<V, bgi::linear<16, 4> > RT;
    //typedef bgi::rtree<V, bgi::quadratic<16, 4> > RT;
    //typedef bgi::rtree<V, bgi::rstar<16, 4> > RT;

    std::cout << "sizeof rtree: " << sizeof(RT) << std::endl;

    for (;;)
    {
        std::vector<V> result;
        result.reserve(100);
        B result_one;

        // packing test
        {
            clock_t::time_point start = clock_t::now();

            RT t(values.begin(), values.end());

            dur_t time = clock_t::now() - start;
            std::cout << time << " - pack " << values_count /*<< '\n'*/;

            std::cout << (bgi::detail::rtree::utilities::are_levels_ok(t) ? " ok" : " NOK")
                      << (bgi::detail::rtree::utilities::are_boxes_ok(t) ? " ok\n" : "NOK\n");

            {
                clock_t::time_point start = clock_t::now();
                size_t temp = 0;
                for (size_t i = 0 ; i < queries_count ; ++i )
                {
                    float x = coords[i].first;
                    float y = coords[i].second;
                    result.clear();
                    t.query(bgi::intersects(B(P(x - 10, y - 10), P(x + 10, y + 10))), std::back_inserter(result));
                    temp += result.size();
                }
                dur_t time = clock_t::now() - start;
                std::cout << time << " - query(B) " << queries_count << " found " << temp << '\n';
            }
        }

        RT t;

        // inserting test
        {
            clock_t::time_point start = clock_t::now();
            t.insert(values);
            dur_t time = clock_t::now() - start;
            std::cout << time << " - insert " << values_count /*<< '\n'*/;

            std::cout << (bgi::detail::rtree::utilities::are_levels_ok(t) ? " ok" : " NOK")
                      << (bgi::detail::rtree::utilities::are_boxes_ok(t) ? " ok\n" : "NOK\n");
        }

        {
            clock_t::time_point start = clock_t::now();
            size_t temp = 0;
            for (size_t i = 0 ; i < queries_count ; ++i )
            {
                float x = coords[i].first;
                float y = coords[i].second;
                result.clear();
                t.query(bgi::intersects(B(P(x - 10, y - 10), P(x + 10, y + 10))), std::back_inserter(result));
                temp += result.size();
            }
            dur_t time = clock_t::now() - start;
            std::cout << time << " - query(B) " << queries_count << " found " << temp << '\n';
        }

        {
            clock_t::time_point start = clock_t::now();
            size_t temp = 0;
            for (size_t i = 0 ; i < queries_count ; ++i )
            {
                float x = coords[i].first;
                float y = coords[i].second;
                result.clear();
                boost::copy(t | bgi::adaptors::queried(bgi::intersects(B(P(x - 10, y - 10), P(x + 10, y + 10)))),
                    std::back_inserter(result));
                temp += result.size();
            }
            dur_t time = clock_t::now() - start;
            std::cout << time << " - range queried(B) " << queries_count << " found " << temp << '\n';
        }

#ifdef BOOST_GEOMETRY_INDEX_DETAIL_EXPERIMENTAL
        {
            clock_t::time_point start = clock_t::now();
            size_t temp = 0;
            for (size_t i = 0 ; i < queries_count ; ++i )
            {
                float x = coords[i].first;
                float y = coords[i].second;
                result.clear();
                std::copy(
                    t.qbegin_(bgi::intersects(B(P(x - 10, y - 10), P(x + 10, y + 10)))),
                    t.qend_(bgi::intersects(B(P(x - 10, y - 10), P(x + 10, y + 10)))),
                           std::back_inserter(result));
                temp += result.size();
            }
            dur_t time = clock_t::now() - start;
            std::cout << time << " - qbegin(B) qend(B) " << queries_count << " found " << temp << '\n';
        }
        {
            clock_t::time_point start = clock_t::now();
            size_t temp = 0;
            for (size_t i = 0 ; i < queries_count ; ++i )
            {
                float x = coords[i].first;
                float y = coords[i].second;
                result.clear();
                mycopy(
                    t.qbegin_(bgi::intersects(B(P(x - 10, y - 10), P(x + 10, y + 10)))),
                    t.qend_(),
                    std::back_inserter(result));
                temp += result.size();
            }
            dur_t time = clock_t::now() - start;
            std::cout << time << " - qbegin(B) qend() " << queries_count << " found " << temp << '\n';
        }
        {
            clock_t::time_point start = clock_t::now();
            size_t temp = 0;
            for (size_t i = 0 ; i < queries_count ; ++i )
            {
                float x = coords[i].first;
                float y = coords[i].second;
                result.clear();
                boost::copy(
                    std::make_pair(
                        t.qbegin_(bgi::intersects(B(P(x - 10, y - 10), P(x + 10, y + 10)))),
                        t.qend_(bgi::intersects(B(P(x - 10, y - 10), P(x + 10, y + 10))))
                    ), std::back_inserter(result));
                temp += result.size();
            }
            dur_t time = clock_t::now() - start;
            std::cout << time << " - range qbegin(B) qend(B)" << queries_count << " found " << temp << '\n';
        }
#endif // BOOST_GEOMETRY_INDEX_DETAIL_EXPERIMENTAL

        {
            clock_t::time_point start = clock_t::now();
            size_t temp = 0;
            for (size_t i = 0 ; i < queries_count ; ++i )
            {
                float x = coords[i].first;
                float y = coords[i].second;
                result.clear();
                RT::const_query_iterator first = t.qbegin(bgi::intersects(B(P(x - 10, y - 10), P(x + 10, y + 10))));
                RT::const_query_iterator last = t.qend();
                std::copy(first, last, std::back_inserter(result));
                temp += result.size();
            }
            dur_t time = clock_t::now() - start;
            std::cout << time << " - type-erased qbegin(B) qend() " << queries_count << " found " << temp << '\n';
        }
        {
            clock_t::time_point start = clock_t::now();
            size_t temp = 0;
            for (size_t i = 0 ; i < queries_count ; ++i )
            {
                float x = coords[i].first;
                float y = coords[i].second;
                result.clear();
                RT::const_query_iterator first = t.qbegin(bgi::intersects(B(P(x - 10, y - 10), P(x + 10, y + 10))));
                RT::const_query_iterator last = t.qend();
                boost::copy(std::make_pair(first, last), std::back_inserter(result));
                temp += result.size();
            }
            dur_t time = clock_t::now() - start;
            std::cout << time << " - range type-erased qbegin(B) qend() " << queries_count << " found " << temp << '\n';
        }

#ifndef SEGMENT_INDEXABLE
        {
            clock_t::time_point start = clock_t::now();
            size_t temp = 0;
            for (size_t i = 0 ; i < queries_count / 2 ; ++i )
            {
                float x1 = coords[i].first;
                float y1 = coords[i].second;
                float x2 = coords[i+1].first;
                float y2 = coords[i+1].second;
                float x3 = coords[i+2].first;
                float y3 = coords[i+2].second;
                result.clear();
                t.query(
                    bgi::intersects(B(P(x1 - 10, y1 - 10), P(x1 + 10, y1 + 10)))
                    &&
                    !bgi::within(B(P(x2 - 10, y2 - 10), P(x2 + 10, y2 + 10)))
                    &&
                    !bgi::covered_by(B(P(x3 - 10, y3 - 10), P(x3 + 10, y3 + 10)))
                    ,
                    std::back_inserter(result)
                    );
                temp += result.size();
            }
            dur_t time = clock_t::now() - start;
            std::cout << time << " - query(i && !w && !c) " << queries_count << " found " << temp << '\n';
        }
#endif

        result.clear();

        {
            clock_t::time_point start = clock_t::now();
            size_t temp = 0;
            for (size_t i = 0 ; i < nearest_queries_count ; ++i )
            {
                float x = coords[i].first + 100;
                float y = coords[i].second + 100;
                result.clear();
                temp += t.query(bgi::nearest(P(x, y), neighbours_count), std::back_inserter(result));
            }
            dur_t time = clock_t::now() - start;
            std::cout << time << " - query(nearest(P, " << neighbours_count << ")) " << nearest_queries_count << " found " << temp << '\n';
        }

#ifdef BOOST_GEOMETRY_INDEX_DETAIL_EXPERIMENTAL
        {
            clock_t::time_point start = clock_t::now();
            size_t temp = 0;
            for (size_t i = 0 ; i < nearest_queries_count ; ++i )
            {
                float x = coords[i].first + 100;
                float y = coords[i].second + 100;
                result.clear();
                std::copy(
                    t.qbegin_(bgi::nearest(P(x, y), neighbours_count)),
                    t.qend_(bgi::nearest(P(x, y), neighbours_count)),
                    std::back_inserter(result));
                temp += result.size();
            }
            dur_t time = clock_t::now() - start;
            std::cout << time << " - qbegin(nearest(P, " << neighbours_count << ")) qend(n) " << nearest_queries_count << " found " << temp << '\n';
        }
        {
            clock_t::time_point start = clock_t::now();
            size_t temp = 0;
            for (size_t i = 0 ; i < nearest_queries_count ; ++i )
            {
                float x = coords[i].first + 100;
                float y = coords[i].second + 100;
                result.clear();
                mycopy(
                    t.qbegin_(bgi::nearest(P(x, y), neighbours_count)),
                    t.qend_(),
                    std::back_inserter(result));
                temp += result.size();
            }
            dur_t time = clock_t::now() - start;
            std::cout << time << " - qbegin(nearest(P, " << neighbours_count << ")) qend() " << nearest_queries_count << " found " << temp << '\n';
        }
#endif // BOOST_GEOMETRY_INDEX_DETAIL_EXPERIMENTAL

        {
            clock_t::time_point start = clock_t::now();
            size_t temp = 0;
            for (size_t i = 0 ; i < nearest_queries_count ; ++i )
            {
                float x = coords[i].first;
                float y = coords[i].second;
                result.clear();
                RT::const_query_iterator first = t.qbegin(bgi::nearest(P(x, y), neighbours_count));
                RT::const_query_iterator last = t.qend();
                std::copy(first, last, std::back_inserter(result));
                temp += result.size();
            }
            dur_t time = clock_t::now() - start;
            std::cout << time << " - type-erased qbegin(nearest(P, " << neighbours_count << ")) qend() " << nearest_queries_count << " found " << temp << '\n';
        }

#ifdef BOOST_GEOMETRY_INDEX_DETAIL_EXPERIMENTAL
#ifndef SEGMENT_INDEXABLE

        {
            LS ls;
            ls.resize(6);

            clock_t::time_point start = clock_t::now();
            size_t temp = 0;
            for (size_t i = 0 ; i < path_queries_count ; ++i )
            {
                float x = coords[i].first;
                float y = coords[i].second;
                for ( int i = 0 ; i < 3 ; ++i )
                {
                    float foo = i*max_val/300;
                    ls[2*i] = P(x, y+foo);
                    ls[2*i+1] = P(x+max_val/100, y+foo);
                }
                result.clear();
                t.query(bgi::path(ls, path_values_count), std::back_inserter(result));
                temp += result.size();
            }
            dur_t time = clock_t::now() - start;
            std::cout << time << " - query(path(LS6, " << path_values_count << ")) " << path_queries_count << " found " << temp << '\n';
        }

        {
            LS ls;
            ls.resize(2);

            clock_t::time_point start = clock_t::now();
            size_t temp = 0;
            for (size_t i = 0 ; i < path_queries_count2 ; ++i )
            {
                float x = coords[i].first;
                float y = coords[i].second;
                ls[0] = P(x, y);
                ls[1] = P(x+max_val/100, y+max_val/100);
                result.clear();
                t.query(bgi::path(ls, path_values_count), std::back_inserter(result));
                temp += result.size();
            }
            dur_t time = clock_t::now() - start;
            std::cout << time << " - query(path(LS2, " << path_values_count << ")) " << path_queries_count2 << " found " << temp << '\n';
        }

        {
            clock_t::time_point start = clock_t::now();
            size_t temp = 0;
            for (size_t i = 0 ; i < path_queries_count2 ; ++i )
            {
                float x = coords[i].first;
                float y = coords[i].second;
                S seg(P(x, y), P(x+max_val/100, y+max_val/100));
                result.clear();
                t.query(bgi::path(seg, path_values_count), std::back_inserter(result));
                temp += result.size();
            }
            dur_t time = clock_t::now() - start;
            std::cout << time << " - query(path(S, " << path_values_count << ")) " << path_queries_count2 << " found " << temp << '\n';
        }
#endif
#endif
        {
            clock_t::time_point start = clock_t::now();
            for (size_t i = 0 ; i < values_count / 10 ; ++i )
            {
                float x = coords[i].first;
                float y = coords[i].second;

                t.remove(generate_value<V>::apply(x, y));
            }
            dur_t time = clock_t::now() - start;
            std::cout << time << " - remove " << values_count / 10 << '\n';

            std::cout << (bgi::detail::rtree::utilities::are_boxes_ok(t) ? " boxes ok\n" : "boxes NOT ok\n");
        }

        std::cout << "------------------------------------------------\n";
    }

    return 0;
}


text:
/*=============================================================================
    Copyright (c) 2001-2011 Joel de Guzman

    Distributed under the Boost Software License, Version 1.0. (See accompanying
    file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

    Autogenerated by MultiStageTable.py (Unicode multi-stage
    table builder) (c) Peter Kankowski, 2008
==============================================================================*/
#if !defined(BOOST_SPIRIT_UNICODE_QUERY_FEBRUARY_2_2010)
#define BOOST_SPIRIT_UNICODE_QUERY_FEBRUARY_2_2010

#include <boost/cstdint.hpp>

# include "category_table.hpp"
# include "script_table.hpp"
# include "lowercase_table.hpp"
# include "uppercase_table.hpp"

namespace boost { namespace spirit { namespace ucd
{
    // This header provides Basic (Level 1) Unicode Support
    // See http://unicode.org/reports/tr18/ for details

    struct properties
    {
        // bit pattern: xxMMMCCC
        // MMM: major_category
        // CCC: category

        enum major_category
        {
            letter,
            mark,
            number,
            separator,
            other,
            punctuation,
            symbol
        };

        enum category
        {
            uppercase_letter = 0,   // [Lu] an uppercase letter
            lowercase_letter,       // [Ll] a lowercase letter
            titlecase_letter,       // [Lt] a digraphic character, with first part uppercase
            modifier_letter,        // [Lm] a modifier letter
            other_letter,           // [Lo] other letters, including syllables and ideographs

            nonspacing_mark = 8,    // [Mn] a nonspacing combining mark (zero advance width)
            enclosing_mark,         // [Me] an enclosing combining mark
            spacing_mark,           // [Mc] a spacing combining mark (positive advance width)

            decimal_number = 16,    // [Nd] a decimal digit
            letter_number,          // [Nl] a letterlike numeric character
            other_number,           // [No] a numeric character of other type

            space_separator = 24,   // [Zs] a space character (of various non-zero widths)
            line_separator,         // [Zl] U+2028 LINE SEPARATOR only
            paragraph_separator,    // [Zp] U+2029 PARAGRAPH SEPARATOR only

            control = 32,           // [Cc] a C0 or C1 control code
            format,                 // [Cf] a format control character
            private_use,            // [Co] a private-use character
            surrogate,              // [Cs] a surrogate code point
            unassigned,             // [Cn] a reserved unassigned code point or a noncharacter

            dash_punctuation = 40,  // [Pd] a dash or hyphen punctuation mark
            open_punctuation,       // [Ps] an opening punctuation mark (of a pair)
            close_punctuation,      // [Pe] a closing punctuation mark (of a pair)
            connector_punctuation,  // [Pc] a connecting punctuation mark, like a tie
            other_punctuation,      // [Po] a punctuation mark of other type
            initial_punctuation,    // [Pi] an initial quotation mark
            final_punctuation,      // [Pf] a final quotation mark

            math_symbol = 48,       // [Sm] a symbol of primarily mathematical use
            currency_symbol,        // [Sc] a currency sign
            modifier_symbol,        // [Sk] a non-letterlike modifier symbol
            other_symbol            // [So] a symbol of other type
        };

        enum derived_properties
        {
            alphabetic = 64,
            uppercase = 128,
            lowercase = 256,
            white_space = 512,
            hex_digit = 1024,
            noncharacter_code_point = 2048,
            default_ignorable_code_point = 4096
        };

        enum script
        {
            arabic = 0,
            imperial_aramaic = 1,
            armenian = 2,
            avestan = 3,
            balinese = 4,
            bamum = 5,
            bengali = 6,
            bopomofo = 7,
            braille = 8,
            buginese = 9,
            buhid = 10,
            canadian_aboriginal = 11,
            carian = 12,
            cham = 13,
            cherokee = 14,
            coptic = 15,
            cypriot = 16,
            cyrillic = 17,
            devanagari = 18,
            deseret = 19,
            egyptian_hieroglyphs = 20,
            ethiopic = 21,
            georgian = 22,
            glagolitic = 23,
            gothic = 24,
            greek = 25,
            gujarati = 26,
            gurmukhi = 27,
            hangul = 28,
            han = 29,
            hanunoo = 30,
            hebrew = 31,
            hiragana = 32,
            katakana_or_hiragana = 33,
            old_italic = 34,
            javanese = 35,
            kayah_li = 36,
            katakana = 37,
            kharoshthi = 38,
            khmer = 39,
            kannada = 40,
            kaithi = 41,
            tai_tham = 42,
            lao = 43,
            latin = 44,
            lepcha = 45,
            limbu = 46,
            linear_b = 47,
            lisu = 48,
            lycian = 49,
            lydian = 50,
            malayalam = 51,
            mongolian = 52,
            meetei_mayek = 53,
            myanmar = 54,
            nko = 55,
            ogham = 56,
            ol_chiki = 57,
            old_turkic = 58,
            oriya = 59,
            osmanya = 60,
            phags_pa = 61,
            inscriptional_pahlavi = 62,
            phoenician = 63,
            inscriptional_parthian = 64,
            rejang = 65,
            runic = 66,
            samaritan = 67,
            old_south_arabian = 68,
            saurashtra = 69,
            shavian = 70,
            sinhala = 71,
            sundanese = 72,
            syloti_nagri = 73,
            syriac = 74,
            tagbanwa = 75,
            tai_le = 76,
            new_tai_lue = 77,
            tamil = 78,
            tai_viet = 79,
            telugu = 80,
            tifinagh = 81,
            tagalog = 82,
            thaana = 83,
            thai = 84,
            tibetan = 85,
            ugaritic = 86,
            vai = 87,
            old_persian = 88,
            cuneiform = 89,
            yi = 90,
            inherited = 91,
            common = 92,
            unknown = 93
        };
    };

    inline properties::category get_category(::boost::uint32_t ch)
    {
        return static_cast<properties::category>(detail::category_lookup(ch) & 0x3F);
    }

    inline properties::major_category get_major_category(::boost::uint32_t ch)
    {
        return static_cast<properties::major_category>(get_category(ch) >> 3);
    }

    inline bool is_punctuation(::boost::uint32_t ch)
    {
        return get_major_category(ch) == properties::punctuation;
    }

    inline bool is_decimal_number(::boost::uint32_t ch)
    {
        return get_category(ch) == properties::decimal_number;
    }

    inline bool is_hex_digit(::boost::uint32_t ch)
    {
        return (detail::category_lookup(ch) & properties::hex_digit) != 0;
    }

    inline bool is_control(::boost::uint32_t ch)
    {
        return get_category(ch) == properties::control;
    }

    inline bool is_alphabetic(::boost::uint32_t ch)
    {
        return (detail::category_lookup(ch) & properties::alphabetic) != 0;
    }

    inline bool is_alphanumeric(::boost::uint32_t ch)
    {
        return is_decimal_number(ch) || is_alphabetic(ch);
    }

    inline bool is_uppercase(::boost::uint32_t ch)
    {
        return (detail::category_lookup(ch) & properties::uppercase) != 0;
    }

    inline bool is_lowercase(::boost::uint32_t ch)
    {
        return (detail::category_lookup(ch) & properties::lowercase) != 0;
    }

    inline bool is_white_space(::boost::uint32_t ch)
    {
        return (detail::category_lookup(ch) & properties::white_space) != 0;
    }

    inline bool is_blank(::boost::uint32_t ch)
    {
        switch (ch)
        {
            case '\n': case '\v': case '\f': case '\r':
                return false;
            default:
                return is_white_space(ch)
                && !(   get_category(ch) == properties::line_separator
                    ||  get_category(ch) == properties::paragraph_separator
                    );
        }
    }

    inline bool is_graph(::boost::uint32_t ch)
    {
        return !(   is_white_space(ch)
                ||  get_category(ch) == properties::control
                ||  get_category(ch) == properties::surrogate
                ||  get_category(ch) == properties::unassigned
                );
    }

    inline bool is_print(::boost::uint32_t ch)
    {
        return (is_graph(ch) || is_blank(ch)) && !is_control(ch);
    }

    inline bool is_noncharacter_code_point(::boost::uint32_t ch)
    {
        return (detail::category_lookup(ch) & properties::noncharacter_code_point) != 0;
    }

    inline bool is_default_ignorable_code_point(::boost::uint32_t ch)
    {
        return (detail::category_lookup(ch) & properties::default_ignorable_code_point) != 0;
    }

    inline properties::script get_script(::boost::uint32_t ch)
    {
        return static_cast<properties::script>(detail::script_lookup(ch) & 0x7F);
    }

    inline ::boost::uint32_t to_lowercase(::boost::uint32_t ch)
    {
        // The table returns 0 to signal that this code maps to itself
        ::boost::uint32_t r = detail::lowercase_lookup(ch);
        return (r == 0)? ch : r;
    }

    inline ::boost::uint32_t to_uppercase(::boost::uint32_t ch)
    {
        // The table returns 0 to signal that this code maps to itself
        ::boost::uint32_t r = detail::uppercase_lookup(ch);
        return (r == 0)? ch : r;
    }
}}}

#endif

ext: hpp


text:
//  io_ex1.cpp  ----------------------------------------------------------//

//  Copyright 2010 Howard Hinnant
//  Copyright 2010 Vicente J. Botet Escriba

//  Distributed under the Boost Software License, Version 1.0.
//  See http://www.boost.org/LICENSE_1_0.txt

/*
This code was adapted by Vicente J. Botet Escriba from Hinnant's html documentation.
Many thanks to Howard for making his code available under the Boost license.

*/

#include <boost/chrono/chrono_io.hpp>
#include <sstream>
#include <iostream>
#include <boost/assert.hpp>

int main()
{
    using namespace boost::chrono;
    using std::cout;

    high_resolution_clock::time_point t0 = high_resolution_clock::now();
    std::stringstream io;
    io << t0;
    BOOST_ASSERT(!io.fail());
    cout << io.str() << '\n';
    BOOST_ASSERT(!io.fail());
    high_resolution_clock::time_point t1;
    io >> t1;
    BOOST_ASSERT(!io.fail());
    cout << io.str() << '\n';
    cout << t0 << '\n';
    cout << t1 << '\n';
    high_resolution_clock::time_point t = high_resolution_clock::now();
    cout << t << '\n';

    cout << "That took " << t - t0 << '\n';
    cout << "That took " << t - t1 << '\n';

    return 0;
}

//~ 50908679121461 nanoseconds since boot
//~ That took 649630 nanoseconds

ext: cpp

max_stars_repo_name: randolphwong/mcsema

max_issues_repo_name: randolphwong/mcsema

max_forks_repo_name: randolphwong/mcsema


text:
/**
 * Genome
 * @author astaroth
 * @date 10/8/20
 *
 * @brief
 */

#ifndef FA_GENOME_HPP
#define FA_GENOME_HPP

#include <string>
#include <array>
#include <utility>
#include <map>

#include <boost/algorithm/string.hpp>

#include "../utils/searcher.hpp"

class Genome {
protected:
    /// constructor
    std::string name, sequence;
    std::size_t threshold;
    bool is_rna;
    std::string start_codon;
    std::array<std::string, 3> end_codons;
    /// !constructor
    std::map< std::size_t, std::pair<std::size_t, std::string> > fragments;  // proteins + islands - aligned version
    std::size_t start_pad {0}, end_pad {0};  // useful when a query starts/ends before/after the reference
    /// flags
    bool proteins_aligned {false}, islands_aligned {false}, is_island {false};
    /// everyone needs a friend
    friend class FileManager;
    friend class Aligner;
    friend class Reference;
    friend class Query;
public:
    /// constructors
    Genome(std::pair<std::string, std::string> ns, const std::size_t &protein_threshold)
    : name(std::move(ns.first)), sequence(std::move(ns.second)), threshold(protein_threshold)
    {
        // make the sequence all uppercase to avoid issues
        boost::to_upper(sequence);
        sequence.shrink_to_fit();
        is_rna = sequence.find('U') != std::string::npos;
        set_codons();  // may be reused if we convert a RNA into a DNA and vice versa
    }
    /// member methods
    // setters
    void set_codons();
    // flags
    void set_proteins_aligned() { proteins_aligned = true; }
    void set_islands_aligned() { islands_aligned = true; }
    void set_whole_island() { is_island = true; }
    // analyzers
    std::map< std::size_t, std::pair<std::size_t, std::string> > find_orfs();
    /// abstract methods
    virtual void search_proteins() = 0;
    virtual void search_islands() = 0;
};

#endif //FA_GENOME_HPP

max_stars_repo_name: ovaltriangle/fingerprint-aligner

max_issues_repo_name: ovaltriangle/fingerprint-aligner

max_forks_repo_name: ovaltriangle/fingerprint-aligner


text:
// Author: Daisuke Kanaizumi
// Affiliation: Department of Applied Mathematics, Waseda University
// Date: Dec 28th, 2017

// verification program for q-Airy functions

#ifndef QAIRY_HPP
#define QAIRY_HPP

#include <kv/interval.hpp>
#include <kv/rdouble.hpp>
#include <kv/constants.hpp>
#include <kv/Pochhammer.hpp>
#include <kv/qBessel.hpp>
#include <cmath>
#include <limits>
#include <kv/convert.hpp> // this was included to use complex numbers
#include <kv/complex.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/io.hpp>
namespace ub = boost::numeric::ublas;
namespace kv {

  // calculate basic hypergeometric series 1phi1

  // reference
  // Fredrik Johansson, Computing hypergeometric functions rigorously, arXiv, 2016

template <class T> interval<T> _1phi_1(const interval<T>& a, const interval<T>& c,const interval<T>& q, const interval<T> & z) {
  int N;
    N=1000;
    interval<T> mid,rad,res,ratio,first;
    mid=1.;
   while(abs(c)>pow(1/q,N)){
     N=N+500;
     //throw std::domain_error("value of N not large enough");
    }   if (q>=1){
     throw std::domain_error("value of q must be under 1");
   }
   if (q<=0){
     throw std::domain_error("q must be positive");
   }
    for(int n=1;n<=N-1;n++){
      mid=mid+qPochhammer(interval<T>(a),interval<T>(q),int (n))*pow(-z,n)*pow(q,n*(n-1)/2)
	/qPochhammer(interval<T>(c),interval<T>(q),int (n))/qPochhammer(interval<T>(q),interval<T>(q),int (n));
    }
    first=abs(qPochhammer(interval<T>(a),interval<T>(q),int (N))*pow(-z,N)*pow(q,N*(N-1)/2)
	      /qPochhammer(interval<T>(c),interval<T>(q),int (N))/qPochhammer(interval<T>(q),interval<T>(q),int (N)));
    ratio=abs(z)*abs(pow(q,N))*(1+abs(pow(q,N)*(c-a)/(1-c*pow(q,N))))*abs(1/(1-pow(q,N+1)));
    if(ratio<1){
      rad=(first/(1-ratio)).upper();
      res=mid+rad*interval<T>(-1.,1.);
      // Alternative implementation
      // Reference
      // Koekoek-Swarttouw,The Askey-scheme of hypergeometric orthogonal polynomials and its q-analogue, arXiv (1996)
      // formula 0.6.12 and 0.6.13
        if((abs(res)).upper()==std::numeric_limits<T>::infinity()&&abs(a)<1&&abs(a)>0){
      res=infinite_qPochhammer(interval<T>(a),interval<T>(q))*infinite_qPochhammer(interval<T>(z),interval<T>(q))
	*Heine(interval<T>(c/a),interval<T>(0.),interval<T>(z),interval<T>(q),interval<T>(a))/infinite_qPochhammer(interval<T>(c),interval<T>(q));
	}
	if((abs(res)).upper()==std::numeric_limits<T>::infinity()&&abs(a*z/c)<1&&abs(a)>0){
	  res=infinite_qPochhammer(interval<T>(a*z/c),interval<T>(q))*Heine(interval<T>(c/a),interval<T>(0.),interval<T>(c),interval<T>(q),interval<T>(a*z/c));
	}

	return res;
    }
    else{
      throw std::domain_error("ratio is more than 1");
    }

}template <class T> interval<T> _1phi_1_bs(const interval<T>& a, const interval<T>& c,const interval<T>& q, const interval<T> & z) {
  int N,nn;
  nn=10;
  N=std::pow(2,nn);
  interval<T> mid,rad,res,ratio,first;
  mid=0.;

    while(abs(c)>pow(1/q,N)){
     nn=nn+5;
     N=std::pow(2,nn);
     //throw std::domain_error("value of N not large enough");
    }
    ub::vector<interval<T> >A(N),B(N),C(N);
    for(int i=0;i<=N/2-1;i++){
      A(i)=pow(-z,i)*pow(q,i*(i-1)*0.5)/qPochhammer(interval<T>(q),interval<T>(q),int (i));
    }
    if (q>=1){
     throw std::domain_error("value of q must be under 1");
    }
    if (q<=0){
      throw std::domain_error("q must be positive");
    }
    for(int n=1;n<=N/2-1;n++){
      B(0)=qPochhammer(interval<T>(a),interval<T>(q),int (n));
      C(0)=qPochhammer(interval<T>(c),interval<T>(q),int (n));
      for(int j=1;j<=N/2-1;j++){
	B(j)=1.;
      }
      for(int k=1;k<=N/2-1;k++){
	C(k)=1.;
      }
      mid=mid+B(0)/C(0)*(A(2*n)*C(2*n+1)+B(2*n)*A(2*n+1));

    }
    first=abs(qPochhammer(interval<T>(a),interval<T>(q),int (N))*pow(-z,N)*pow(q,N*(N-1)/2)
	      /qPochhammer(interval<T>(c),interval<T>(q),int (N))/qPochhammer(interval<T>(q),interval<T>(q),int (N)));
    ratio=abs(z)*abs(pow(q,N))*(1+abs(pow(q,N)*(c-a)/(1-c*pow(q,N))))*abs(1/(1-pow(q,N+1)));
    if(ratio<1){
      rad=(first/(1-ratio)).upper();
      res=mid+rad*interval<T>(-1.,1.);
      // Alternative implementation
      // Reference
      // Koekoek-Swarttouw,The Askey-scheme of hypergeometric orthogonal polynomials and its q-analogue, arXiv (1996)
      // formula 0.6.12 and 0.6.13
        if((abs(res)).upper()==std::numeric_limits<T>::infinity()&&abs(a)<1&&abs(a)>0){
      res=infinite_qPochhammer(interval<T>(a),interval<T>(q))*infinite_qPochhammer(interval<T>(z),interval<T>(q))
	*Heine(interval<T>(c/a),interval<T>(0.),interval<T>(z),interval<T>(q),interval<T>(a))/infinite_qPochhammer(interval<T>(c),interval<T>(q));
	}
	if((abs(res)).upper()==std::numeric_limits<T>::infinity()&&abs(a*z/c)<1&&abs(a)>0){
	  res=infinite_qPochhammer(interval<T>(a*z/c),interval<T>(q))*Heine(interval<T>(c/a),interval<T>(0.),interval<T>(c),interval<T>(q),interval<T>(a*z/c));
	}
	if((abs(res)).upper()==std::numeric_limits<T>::infinity()&&a<1&&a>=0&&c>=0&&c<1){
  // Alternative implementation
  // Reference: Plancherel-Rotach asymptotics for certain basic hypergeometric series
  // Zhang, 2008, Advances in Mathematics
	  int nn,mu;
	  nn=5;
	  mu=std::floor(nn/2);
	  interval<T>r,amid;

	  T arad;
	  arad=(32*Euler(interval<T>(q))*infinite_qPochhammer(interval<T>(-q*pow(q,-nn)/abs(z)),interval<T>(q))
		*infinite_qPochhammer(interval<T>(-q*pow(q,nn)*abs(z)),interval<T>(q))
		/infinite_qPochhammer(interval<T>(a),interval<T>(q))
		*(pow(q,mu+1)/(1-q)+pow(q,mu*mu/2+mu/2)/pow(abs(z*pow(q,nn)),mu))).upper();
	  amid=Euler(interval<T>(q))*infinite_qPochhammer(interval<T>(q*pow(q,-nn)/abs(z)),interval<T>(q))
	    *infinite_qPochhammer(interval<T>(q*pow(q,nn)*abs(z)),interval<T>(q));
	  r=interval<T>(amid.lower()-arad,amid.upper()+arad);
	  res=Euler(interval<T>(q))*infinite_qPochhammer(interval<T>(c),interval<T>(q))*pow(-z*pow(q,nn),nn)*r
	    /infinite_qPochhammer(interval<T>(a),interval<T>(q))/pow(q,nn*(nn+1)/2);
	}
	return res;
    }
    else{
      throw std::domain_error("ratio is more than 1");
    }

}
  template <class T> complex<interval<T> >_1phi_1(const complex<interval<T> >& a, const complex<interval<T> >& c,const interval<T>& q, const complex<interval<T> >& z) {
     int N;
     N=1000;
     while(abs(c)>pow(1/q,N)){
       N=N+500;
       //throw std::domain_error("value of N not large enough");
     }
     if (q>=1){
       throw std::domain_error("value of q must be under 1");
     }
     if (q<=0){
       throw std::domain_error("q must be positive");
     }
     complex<interval<T> > mid,res;
     interval<T>ratio,first;
     T rad;
     mid=1.;
     for(int n=1;n<=N-1;n++){
       mid=mid+qPochhammer(complex<interval<T> >(a),interval<T>(q) ,int (n))*pow(-z,n)*pow(q,n*(n-1)/2)
	 /qPochhammer(complex<interval<T> >(c),interval<T>(q),int (n))/qPochhammer(interval<T>(q),interval<T>(q),int (n));
     }
     first=abs(qPochhammer(complex<interval<T> >(a),interval<T>(q),int (N))*pow(-z,N)*pow(q,N*(N-1)/2)
	       /qPochhammer(complex<interval<T> >(c),interval<T>(q),int (N))/qPochhammer(interval<T>(q),interval<T>(q),int (N)));
     ratio=abs(z)*pow(q,N)*(1+abs(pow(q,N)*(c-a)/(1-c*pow(q,N))))*abs(1/(1-pow(q,N+1)));
     if(abs(ratio)<1){
       rad=(first/(1-ratio)).upper();
       res=complex_nbd(mid,rad);
       // Alternative implementation
       // Reference
       // Koekoek-Swarttouw,The Askey-scheme of hypergeometric orthogonal polynomials and its q-analogue, arXiv (1996)
       // formula 0.6.12 and 0.6.13
       if((abs(res)).upper()==std::numeric_limits<T>::infinity()&&abs(a)<1&&abs(a)>0){
	 res=infinite_qPochhammer(complex<interval<T> >(a),interval<T> (q))*infinite_qPochhammer(complex<interval<T> >(z),interval<T>(q))
	   *Heine(complex<interval<T> >(c/a),complex<interval<T> >(0.),complex<interval<T> >(z),interval<T>(q),complex<interval<T> >(a))/infinite_qPochhammer(complex<interval<T> >(c),interval<T>(q));
       }
       if((abs(res)).upper()==std::numeric_limits<T>::infinity()&&abs(a*z/c)<1&&abs(a)>0){
	 res=infinite_qPochhammer(complex<interval<T> >(a*z/c),interval<T>(q))*Heine(complex<interval<T> >(c/a),complex<interval<T> >(0.),complex<interval<T> >(c),interval<T>(q),complex<interval<T> >(a*z/c));
       }
      return res;
     }
     else{
       throw std::domain_error("ratio is more than 1");
     }
  }
  // El-Guindy, A., & Ismail, M. E. (2016). On Certain Generalizations of Rogers-Ramanujan Type Identities. arXiv preprint arXiv:1602.00316.
  template <class T> interval<T> _1phi_1_ver2(const interval<T> & a, const interval<T> & c,const interval<T>& q, const interval<T> & z) {
    interval<T>  res;
    res=_1phi_1(interval<T> (a),interval<T> (c), interval<T>(q),interval<T> (z));
    if((abs(res)).upper()==std::numeric_limits<T>::infinity()){
      res=infinite_qPochhammer(interval<T> (z),interval<T>(q))
	/infinite_qPochhammer(interval<T> (c),interval<T>(q))
	*_1phi_1(interval<T> (a*z/c),interval<T> (z),interval<T>(q),interval<T> (c));

      }
    return res;
  }
  template <class T> complex<interval<T> >_1phi_1_ver2(const complex<interval<T> >& a, const complex<interval<T> >& c,const interval<T>& q, const complex<interval<T> >& z) {
    complex<interval<T> > res;
    res=_1phi_1(complex<interval<T> >(a),complex<interval<T> >(c), interval<T>(q),complex<interval<T> >(z));
    if((abs(res)).upper()==std::numeric_limits<T>::infinity()){
      res=infinite_qPochhammer(complex<interval<T> >(z),interval<T>(q))
	/infinite_qPochhammer(complex<interval<T> >(c),interval<T>(q))
	*_1phi_1(complex<interval<T> >(a*z/c),complex<interval<T> >(z),interval<T>(q),complex<interval<T> >(c));
    }
    return res;
  }
  // calculate basic hypergeometric series 0phi1

  // reference
  // Fredrik Johansson, Computing hypergeometric functions rigorously, arXiv, 2016

  template <class T> interval<T> _0phi_1(const interval<T>& c,const interval<T>& q, const interval<T> & z) {
    int N;
    N=1000;
    interval<T> mid,res,ratio,first;
    T rad;
    while(abs(c)>pow(1/q,N)){
      N=N+500;
      //throw std::domain_error("value of N not large enough");
    }
    if (q>=1){
      throw std::domain_error("value of q must be under 1");
    }
    if (q<=0){
      throw std::domain_error("q must be positive");
    }
    mid=1.;
    for(int n=1;n<=N-1;n++){
      mid=mid+pow(z,n)*pow(q,n*(n-1))
	/qPochhammer(interval<T>(c),interval<T>(q),int (n))/qPochhammer(interval<T>(q),interval<T>(q),int (n));
    }
    first=abs(pow(z,N)*pow(q,N*(N-1))
	      /qPochhammer(interval<T>(c),interval<T>(q),int (N))/qPochhammer(interval<T>(q),interval<T>(q),int (N)));
    ratio=abs(z)*pow(q,2*N)*abs(1/(1-c*pow(q,N)))*abs(1/(1-pow(q,N+1)));
    if(ratio<1){
      rad=(first/(1-ratio)).upper();
      res=mid+rad*interval<T>(-1.,1.);
      if((abs(res)).upper()==std::numeric_limits<T>::infinity()&&abs(c)>0){
	// Alternative implementation
	// Reference
	// Koekoek-Swarttouw,The Askey-scheme of hypergeometric orthogonal polynomials and its q-analogue, arXiv (1996)
	// formula 0.6.8 and 0.6.9
	// H. T Koelink, Hansen-Lommel orthogonality Relations for Jackson`s q-Bessel functions, formula 3.2
	// Journal of Mathematical Analysis and Applications 175, 425-437 (1993)
	res=_1phi_1(interval<T>(z/c),interval<T>(0.),interval<T>(q),interval<T>(c))
	  /infinite_qPochhammer(interval<T>(c),interval<T>(q));
      }

      return res;
    }
    else{
      throw std::domain_error("ratio is more than 1");
    }
 }
  template <class T> complex<interval<T> >_0phi_1(const complex<interval<T> >& c,const interval<T>& q, const complex<interval<T> >& z) {
    int N;
    N=1000;
    complex<interval<T> > mid,res;
    interval<T>ratio,first;
    T rad;
    mid=1.;
  while(abs(c)>pow(1/q,N)){
     N=N+500;
     //throw std::domain_error("value of N not large enough");
    }
   if (q>=1){
     throw std::domain_error("value of q must be under 1");
   }
   if (q<=0){
     throw std::domain_error("q must be positive");
   }
    for(int n=1;n<=N-1;n++){
      mid=mid+pow(z,n)*pow(q,n*(n-1))
	/qPochhammer(complex<interval<T> >(c),interval<T>(q),int (n))/qPochhammer(interval<T>(q),interval<T>(q),int (n));
    }
    first=abs(pow(z,N)*pow(q,N*(N-1))
	      /qPochhammer(complex<interval<T> >(c),interval<T>(q),int (N))/qPochhammer(interval<T>(q),interval<T>(q),int (N)));
    ratio=abs(z)*pow(q,2*N)*abs(1/(1-c*pow(q,N)))*abs(1/(1-pow(q,N+1)));
    if(abs(ratio)<1){
      rad=(first/(1-ratio)).upper();
      res=complex_nbd(mid,rad);
      if((abs(res)).upper()==std::numeric_limits<T>::infinity()&&abs(c)>0){
	// Alternative implementation
	// Reference
	// Koekoek-Swarttouw,The Askey-scheme of hypergeometric orthogonal polynomials and its q-analogue, arXiv (1996)
	// formula 0.6.8 and 0.6.9
	// H. T Koelink, Hansen-Lommel orthogonality Relations for Jackson`s q-Bessel functions, formula 3.2
	// Journal of Mathematical Analysis and Applications 175, 425-437 (1993)
	res=_1phi_1(complex<interval<T> >(z/c),complex<interval<T> >(0.),interval<T>(q),complex<interval<T> >(c))
	  /infinite_qPochhammer(complex<interval<T> >(c),interval<T>(q));
 }

      return res;
    }
    else{
      throw std::domain_error("ratio is more than 1");
    }
 }
 template <class T> complex<interval<T> >HKW_qAiry(const interval<T>& q, const complex<interval<T> >& z) {
   // calculate Hamamoto-Kajiwara-Witte q-Airy function
   // reference
   // Hamamoto, Kajiwara, Witte, Hypergeometric solutions to the q-Painlev\'e equation of type $(A_1+A^{\prime}_1)~{(1)}$
   // International Mathematics Research Notices, 2006
   complex<interval<T> >res,qq;
   qq=q;
   res=_1phi_1(complex<interval<T> >(0),complex<interval<T> >(-qq),interval<T>(q),complex<interval<T> >(-z));
   if((abs(res)).upper()==std::numeric_limits<T>::infinity()){
     // implementation to avoid numerical divergence
     // reference
     // T. H. Koornwinder and R. F. Swarttouw, On $q$-analogues of the Fourier and Hankel transforms
     // Transactions of the American Mathematical Society, 333 (1992), 445-461
     // symmetry relation is used
     res=infinite_qPochhammer(complex<interval<T> >(-z),interval<T> (q))*_1phi_1(complex<interval<T> >(0),complex<interval<T> >(-z),interval<T>(q),complex<interval<T> >(-qq))/infinite_qPochhammer(interval<T>(-q),interval<T>(q));
   }
   return res;
 }
  template <class T> complex<interval<T> >KMNOY_qAiry(const interval<T>& q, const complex<interval<T> >& z) {
    // calculate Kajiwara-Masuda-Noumi-Ohta-Yamada q-Airy function
    // reference
    // Kajiwara, Masuda, Noumi, Ohta, Yamada, Hypergeometric Solutions to the q-Painlev\'e Equations
    // International Mathematics Research Notices, 2004
    complex<interval<T> >res,qq;
    qq=q;
    res=_1phi_1(complex<interval<T> >(0),complex<interval<T> >(-qq),interval<T>(q),complex<interval<T> >(-z*pow(q,0.5)));
   if((abs(res)).upper()==std::numeric_limits<T>::infinity()){
     // implementation to avoid numerical divergence
     // reference
     // T. H. Koornwinder and R. F. Swarttouw, On $q$-analogues of the Fourier and Hankel transforms
     // Transactions of the American Mathematical Society, 333 (1992), 445-461
     // symmetry relation is used
     res=infinite_qPochhammer(complex<interval<T> >(-z*pow(qq,0.5)),interval<T> (q))
       *_1phi_1(complex<interval<T> >(0),complex<interval<T> >(-z*pow(qq,0.5)),interval<T>(q),complex<interval<T> >(-qq))
       /infinite_qPochhammer(interval<T>(-q),interval<T>(q));
   }
    return res;
 }
template <class T> interval<T> Ramanujan_qAiry(const interval<T>& q, const interval<T> & z) {
  // calculate Ramanujan`s q-Airy function
  interval<T> res;
  if(z>0){
    // reference: ENCYCLOPEDIA OF MATHEMATICS AND ITS APPLICATOINS 71, SPECIAL FUNCTIONS
    // G.E.ANDREWS, RICHARD ASKEY, RANJAN ROY
    // CAMBRIDGE, 1999, Page 551, Exercise 39
    int K=1000;
    int j=1;
    interval<T>a,b,c,series;
    a=1.;
    b=1.;
    c=1.;
    // Leibniz criterion
    for (int k=1; k<=K; k++){
    j = -1*j;
    b=b*(1-pow(q,2*k));
    c=c*(1-z*pow(q,2*k));
    a = a+j*pow(q,k*k)*pow(z,k)/(b*c);
    }
    b=b*(1-pow(q,2*K+2));
    c=c*(1-z*pow(q,2*K+2));
    T rad=abs(pow(q,(K+1)*(K+1))*pow(z,K+1)/(b*c)).upper();
    series=a+rad*interval<T>(-1.,1.);
    res=series*infinite_qPochhammer(interval<T>(z*q*q),interval<T>(q*q));
  }
  else{
   res=_0phi_1(interval<T> (0),interval<T>(q),interval<T> (-q*z));
  }

  if((abs(res)).upper()==std::numeric_limits<T>::infinity()){
    // implementation to avoid numerical divergence
    // reference
    // Mourad E. H. Ismail, Changgui Zhang, Zeros of entire functions and a problem of Ramanujan, Theorem 3.1
    // Advances in Mathematics 209 (2007) 363–380

    // In this paper, the proof of theorem 3.1 is partially abbreviated.
    // By using the following identities, the proof will be completed.
    // Gasper and Rahman, Basic Hypergeometric Series, Cambridge University Press, 1990, Appendix I, (I.2), (I.5), (I.17), (I.25)
    // Weisstein, Eric W. "q-Pochhammer Symbol." From MathWorld--A Wolfram Web Resource. http://mathworld.wolfram.com/q-PochhammerSymbol.html, formula (10)

    res=(infinite_qPochhammer(interval<T>(q*z),interval<T>(q*q))*infinite_qPochhammer(interval<T>(q/z),interval<T>(q*q))
	 *_1phi_1(interval<T>(0.),interval<T>(q),interval<T>(q*q),interval<T>(pow(q,2)/z))-
	 q*infinite_qPochhammer(interval<T>(q*q*z),interval<T>(q*q))*infinite_qPochhammer(interval<T>(1/z),interval<T>(q*q))/(1-q)
	 *_1phi_1(interval<T>(0.),interval<T>(q*q*q),interval<T>(q*q),interval<T>(pow(q,3)/z)))
      /infinite_qPochhammer(interval<T>(q),interval<T>(q*q));
  }
   return res;
 }
  template <class T> complex<interval<T> >Ramanujan_qAiry(const interval<T>& q, const complex<interval<T> >& z) {
    // calculate Ramanujan`s q-Airy function
    complex<interval<T> >res,qq,i;
    qq=q;
    res=_0phi_1(complex<interval<T> >(0),interval<T>(q),complex<interval<T> >(-qq*z));
    if((abs(res)).upper()==std::numeric_limits<T>::infinity()){
      // implementation to avoid numerical divergence
      // reference
      // Mourad E. H. Ismail, Ruiming Zhang - Proceedings of the American Mathematical Society, 2017
      // Integral and series representations of $q$-polynomials and functions: Part II Schur polynomials and the Rogers-Ramanujan identities
     i=complex<interval<T> >::i();
     res=infinite_qPochhammer(complex<interval<T> >(-i*sqrt(z)*pow(q,0.25)),interval<T>(sqrt(q)))
       *_1phi_1(complex<interval<T> >(0.),complex<interval<T> >(-i*sqrt(z)*pow(q,0.25)),interval<T>(sqrt(q)),complex<interval<T> >(i*sqrt(z)*pow(q,0.25)));
    }
    if((abs(res)).upper()==std::numeric_limits<T>::infinity()){
      // implementation to avoid numerical divergence
      // reference
      // Mourad E. H. Ismail, Changgui Zhang, Zeros of entire functions and a problem of Ramanujan, Theorem 3.1
      // Advances in Mathematics 209 (2007) 363–380

      // In this paper, the proof of theorem 3.1 is partially abbreviated.
      // By using the following identities, the proof will be completed.
      // Gasper and Rahman, Basic Hypergeometric Series, Cambridge University Press, 1990, Appendix I, (I.2), (I.5), (I.17), (I.25)
      // Weisstein, Eric W. "q-Pochhammer Symbol." From MathWorld--A Wolfram Web Resource. http://mathworld.wolfram.com/q-PochhammerSymbol.html, formula (10)
      res=(infinite_qPochhammer(complex<interval<T> >(qq*z),interval<T>(q*q))*infinite_qPochhammer(complex<interval<T> >(qq/z),interval<T>(q*q))
	   *_1phi_1(complex<interval<T> >(0.),complex<interval<T> >(qq),interval<T>(q*q),complex<interval<T> >(pow(qq,2)/z))-
	   q*infinite_qPochhammer(complex<interval<T> >(qq*qq*z),interval<T>(q*q))*infinite_qPochhammer(complex<interval<T> >(1/z),interval<T>(q*q))/(1-q)
	   *_1phi_1(complex<interval<T> >(0.),complex<interval<T> >(qq*qq*qq),interval<T>(q*q),complex<interval<T> >(pow(qq,3)/z)))
	/infinite_qPochhammer(interval<T>(q),interval<T>(q*q));
       }
    return res;
  }

  template <class T> complex<interval<T> >Ramanujan_qAiry_Morita(const interval<T>& q, const complex<interval<T> >& z) {
    // Morita, T. (2011). A connection formula between the Ramanujan function and the $ q $-Airy function. arXiv preprint arXiv:1104.0755.
    complex<interval<T> >res,i,x;
    i=complex<interval<T> >::i();
    x=i*pow(q,0.75)/sqrt(z);
    res=(infinite_qPochhammer(complex<interval<T> >(-x/sqrt(q)),interval<T>(sqrt(q)))*infinite_qPochhammer(complex<interval<T> >(-q/x),interval<T>(sqrt(q)))*HKW_qAiry(interval<T>(sqrt(q)),complex<interval<T> >(-x))
	  +infinite_qPochhammer(complex<interval<T> >(x/sqrt(q)),interval<T>(sqrt(q)))*infinite_qPochhammer(complex<interval<T> >(q/x),interval<T>(sqrt(q)))*HKW_qAiry(interval<T>(sqrt(q)),complex<interval<T> >(x)))
      /infinite_qPochhammer(interval<T>(-1.),interval<T>(sqrt(q)));
    return res;
  }
  template <class T> complex<interval<T> >Ramanujan_qAiry_mqB(const interval<T>& q, const complex<interval<T> >& z) {
    complex<interval<T> >res;
    interval<T>pi;
    pi=constants<interval<T> >::pi();
    res=sqrt(q*z)/sqrt(1+q+q*q)/pi
      *(modified_qBesselI2(complex<interval<T> >(2/(1+q+q*q)*pow(q*z,1.5)*(1-pow(q,1/3.))),complex<interval<T> >(-1/3.,0.),interval<T>(pow(q,3/2.)))-
	modified_qBesselI2(complex<interval<T> >(2/(1+q+q*q)*pow(q*z,1.5)*(1-pow(q,1/3.))),complex<interval<T> >(1/3.,0),interval<T>(pow(q,3/2.))));
    return res;
    // to be updated
  }
  /*
  template <class T> complex<interval<T> >Ramanujan_qAiry_ae(const interval<T>& q, const complex<interval<T> >& Z) {
    // calculate Ramanujan`s q-Airy function with asymptotic formula
    // reference
    // Ismail, M. E., & Zhang, R. (2016).
    // Integral and Series Representations of $ q $-Polynomials and Functions: Part I.
    // arXiv preprint arXiv:1604.08441.
    complex<interval<T> >res,z,nbd;
    int n;
    n=10;

    z=Z*pow(q,2*n);

    res=pow(q,-n*n)*pow(-z,n)/Euler(interval<T>(q))*Euler(interval<T>(q*q))
      *infinite_qPochhammer(complex<interval<T> >(q*z),interval<T>(q*q))
      *infinite_qPochhammer(complex<interval<T> >(q/z),interval<T>(q*q));
    return res;
  }
  */
  template <class T> interval<T> IKMMS_qAiry(const interval<T>& q, const interval<T>& tau){
    // calculate Isojima-Konno-Mimura-Murata-Satsuma q-Airy function
    // reference
    // Isojima, Konno, Mimura, Murata, Satsuma, Ultradiscrete Painlev\`e II equation and a special function solution
    // JOURNAL OF PHYSICS A: MATHEMATICAL AND THEORETICAL, 2011
   if (q>=1){
     throw std::domain_error("value of q must be under 1");
   }
   if (q<=0){
     throw std::domain_error("q must be positive");
   }
   if (tau<=0){
     throw std::domain_error("tau must be positive");
   }
   int N;
   N=1000;
   interval<T> mid,res,ratio,first,pi,qq,pt;
   T rad;
   qq=1.;
   pt=1.;
   pi=constants<interval<T> >::pi();
   mid=sqrt(2)*sin(pi*(2*log(tau)/log(q)+1)*0.25);

  for(int n=1;n<=N-1;n++){
    qq*=q*q;
    pt*=tau;
    mid=mid+sqrt(2)*sin(pi*(2*n+2*log(tau)/log(q)+1)*0.25)*pow(q,n*(n+1)/2)*tau/(1-qq);
  }
  qq*=q*q;
  pt*=tau;
  first=abs(sqrt(2)*sin(pi*(2*N+2*log(tau)/log(q)+1)*0.25)*pow(q,N*(N+1)/2)*tau/(1-qq));
  ratio=abs(tau)*pow(q,N+1)/(1-pow(q,2*N+2));
  // ratio of majorant
  if(ratio<1){
    rad=(first/(1-ratio)).upper();
    res=mid+rad*interval<T>(-1.,1.);
    return res;
  }
  else{
    throw std::domain_error("ratio is more than 1");
  }
  }
template <class T> interval<T> q_Bi(const interval<T>& q, const interval<T>& tau){
    // calculate q-Bi function
    // reference
    // Isojima, Konno, Mimura, Murata, Satsuma, Ultradiscrete Painlev\`e II equation and a special function solution
    // JOURNAL OF PHYSICS A: MATHEMATICAL AND THEORETICAL, 2011
   if (q>=1){
     throw std::domain_error("value of q must be under 1");
   }
   if (q<=0){
     throw std::domain_error("q must be positive");
   }
   if (tau<=0){
     throw std::domain_error("tau must be positive");
   }

   int N;
   N=1000;
   interval<T> mid,res,ratio,first,pi,qq,pt;
   T rad;
   qq=1.;
   pt=1.;
   pi=constants<interval<T> >::pi();
   mid=sqrt(2)*sin(pi*(2*log(tau)/log(q)+3)*0.25);
   for(int n=1;n<=N-1;n++){
    qq*=q*q;
    pt*=tau;
    mid=mid+sqrt(2)*sin(pi*(2*n+2*log(tau)/log(q)+3)*0.25)*pow(q,n*(n+1)/2)*tau/(1-qq);
   }
   qq*=q*q;
   pt*=tau;
   first=abs(sqrt(2)*sin(pi*(2*N+2*log(tau)/log(q)+3)*0.25)*pow(q,N*(N+1)/2)*tau/(1-qq));
   ratio=abs(tau)*pow(q,N+1)/(1-pow(q,2*N+2));
  // ratio of majorant
   if(ratio<1){
     rad=(first/(1-ratio)).upper();
     res=mid+rad*interval<T>(-1.,1.);
     return res;
   }
   else{
     throw std::domain_error("ratio is more than 1");
   }
}
}
#endif

max_stars_repo_name: Daisuke-Kanaizumi/q-special-functions

max_issues_repo_name: Daisuke-Kanaizumi/q-special-functions

max_forks_repo_name: Daisuke-Kanaizumi/q-special-functions


text:
//
// Licensed to Green Energy Corp (www.greenenergycorp.com) under one
// or more contributor license agreements. See the NOTICE file
// distributed with this work for additional information
// regarding copyright ownership.  Green Enery Corp licenses this file
// to you 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.
//
#include <boost/test/unit_test.hpp>

#include <opendnp3/APL/Log.h>
#include <opendnp3/APL/Exception.h>
#include <tinyxml/tinyxml.h>
#include <XMLBindings/APLXML_Base.h>
#include <APLXML/tinybinding.h>
#include <APLXML/SingleNodeSaver.h>

using namespace apl;

BOOST_AUTO_TEST_SUITE(XmlTests)

BOOST_AUTO_TEST_CASE(ParseInt)
{
	TiXmlElement testNode("Test");
	BOOST_CHECK_EQUAL(IXMLDataBound::FromString_int(&testNode, "100"), 100);
	BOOST_CHECK_EQUAL(IXMLDataBound::FromString_int(&testNode, "123456789"), 123456789);
	BOOST_CHECK_EQUAL(IXMLDataBound::FromString_int(&testNode, "-100"), -100);

	//dont do implit casting to int
	BOOST_CHECK_THROW(IXMLDataBound::FromString_int(&testNode, "-1.00"), apl::Exception);
	//dont read a string as int
	BOOST_CHECK_THROW(IXMLDataBound::FromString_int(&testNode, "abcasdasd"), apl::Exception);
	//empty string where an int should fail
	BOOST_CHECK_THROW(IXMLDataBound::FromString_int(&testNode, ""), apl::Exception);

	BOOST_CHECK_EQUAL(IXMLDataBound::FromString_int(&testNode, "2147483647"), 2147483647);
	BOOST_CHECK_THROW(IXMLDataBound::FromString_int(&testNode, "2147483648"), apl::Exception);
}
BOOST_AUTO_TEST_CASE(WriteInt)
{
	BOOST_CHECK_EQUAL(IXMLDataBound::ToString_int(100), "100");
	BOOST_CHECK_EQUAL(IXMLDataBound::ToString_int(-1), "-1");
}
BOOST_AUTO_TEST_CASE(ParseBool)
{
	TiXmlElement testNode("Test");
	BOOST_CHECK_EQUAL(IXMLDataBound::FromString_bool(&testNode, "true"), true);
	//can only use "true" and "false" any other casing is an error.
	BOOST_CHECK_THROW(IXMLDataBound::FromString_bool(&testNode, "TRUE"), apl::Exception);
	BOOST_CHECK_THROW(IXMLDataBound::FromString_bool(&testNode, "True"), apl::Exception);
	BOOST_CHECK_THROW(IXMLDataBound::FromString_bool(&testNode, "tRue"), apl::Exception);

	BOOST_CHECK_EQUAL(IXMLDataBound::FromString_bool(&testNode, "false"), false);
	BOOST_CHECK_THROW(IXMLDataBound::FromString_bool(&testNode, "FALSE"), apl::Exception);
	BOOST_CHECK_THROW(IXMLDataBound::FromString_bool(&testNode, "False"), apl::Exception);
	BOOST_CHECK_THROW(IXMLDataBound::FromString_bool(&testNode, "faLse"), apl::Exception);

	BOOST_CHECK_EQUAL(IXMLDataBound::FromString_bool(&testNode, "0"), false);
	BOOST_CHECK_EQUAL(IXMLDataBound::FromString_bool(&testNode, "1"), true);
	BOOST_CHECK_THROW(IXMLDataBound::FromString_bool(&testNode, ""), apl::Exception);
}
BOOST_AUTO_TEST_CASE(WriteBool)
{
	BOOST_CHECK_EQUAL(IXMLDataBound::ToString_bool(true), "true");
	BOOST_CHECK_EQUAL(IXMLDataBound::ToString_bool(false), "false");
}
BOOST_AUTO_TEST_CASE(ParseDouble)
{
	TiXmlElement testNode("Test");
	BOOST_CHECK_EQUAL(IXMLDataBound::FromString_double(&testNode, "100"), 100);
	BOOST_CHECK_EQUAL(IXMLDataBound::FromString_double(&testNode, "3.1456"), 3.1456);
	BOOST_CHECK_EQUAL(IXMLDataBound::FromString_double(&testNode, "-100."), -100);

	BOOST_CHECK_EQUAL(IXMLDataBound::FromString_double(&testNode, "-100.0e-7"), -1.0e-5);
	BOOST_CHECK_EQUAL(IXMLDataBound::FromString_double(&testNode, "-100.0E-7"), -1.0e-5);
	BOOST_CHECK_EQUAL(IXMLDataBound::FromString_double(&testNode, "2147483648"), 2147483648UL);
	BOOST_CHECK_EQUAL(IXMLDataBound::FromString_double(&testNode, "-100.0E4"), -1.0e6);
	BOOST_CHECK_EQUAL(IXMLDataBound::FromString_double(&testNode, "-100.0E+4"), -1.0e6);
}
BOOST_AUTO_TEST_CASE(WriteDouble)
{
	BOOST_CHECK_EQUAL(IXMLDataBound::ToString_double(55), "55");
}
BOOST_AUTO_TEST_CASE(ParseString)
{
	TiXmlElement testNode("Test");
	BOOST_CHECK_EQUAL(IXMLDataBound::FromString_string(&testNode, "qwertyuiop"), "qwertyuiop");

	//consider an empty string to be a valid string
	BOOST_CHECK_EQUAL(IXMLDataBound::FromString_string(&testNode, ""), "");
}
BOOST_AUTO_TEST_CASE(WriteString)
{
	BOOST_CHECK_EQUAL(IXMLDataBound::ToString_string("qwertyuiop"), "qwertyuiop");
	BOOST_CHECK_EQUAL(IXMLDataBound::ToString_string(""), "");
}

BOOST_AUTO_TEST_SUITE_END()

BOOST_AUTO_TEST_SUITE(SingleNodeSuite)

BOOST_AUTO_TEST_CASE(SaveAndLoad)
{
	EventLog log;
	SingleNodeSaver<APLXML_Base::TCPServer_t> s("TestSingle.tiny.xml", log.GetLogger(LEV_DEBUG, "SingleNode"));
	APLXML_Base::TCPServer_t node;

	node.Name = "HelloWorld";
	node.Port = 1000;

	BOOST_CHECK(s.SaveNode(node));

	APLXML_Base::TCPServer_t node2;
	BOOST_CHECK(s.LoadNode(node2));

	BOOST_CHECK_EQUAL(node2.Name, "HelloWorld");
	BOOST_CHECK_EQUAL(node2.Port, 1000);
}
BOOST_AUTO_TEST_CASE(LoadFails)
{
	EventLog log;
	SingleNodeSaver<APLXML_Base::Serial_t> s("TestSingleFailure.tiny.xml", log.GetLogger(LEV_DEBUG, "SingleNode"));
	APLXML_Base::Serial_t node;

	BOOST_CHECK_EQUAL(s.LoadNode(node), false);
}

BOOST_AUTO_TEST_SUITE_END()

max_stars_repo_name: rahul-sachdev/DNPSimulator

max_issues_repo_name: rahul-sachdev/DNPSimulator


text:
#include "Utils.h"

#include <cctype>
#include <iostream>

#include <boost/icl/type_traits/is_numeric.hpp>

#include "Log.h"

bool is_number(const std::string& s)
{
	return !s.empty() && std::find_if(s.begin(), s.end(), [](char c) { return !std::isdigit(c); }) == s.end();
}

bool copyDir(
	boost::filesystem::path const & source,
	boost::filesystem::path const & destination
	)
{
	namespace fs = boost::filesystem;
	try
	{
		// Check whether the function call is valid
		if (
			!fs::exists(source) ||
			!fs::is_directory(source)
			)
		{
			LOG(LogLevel::Error) << "Source directory " << source.string()
				<< " does not exist or is not a directory.";
			return false;
		}
		if (fs::exists(destination))
		{
			LOG(LogLevel::Error) << "Destination directory " << destination.string()
				<< " already exists.";
			return false;
		}
		// Create the destination directory
		if (!fs::create_directory(destination))
		{
			LOG(LogLevel::Error) << "Unable to create destination directory"
				<< destination.string();
			return false;
		}
	}
	catch (fs::filesystem_error const & e)
	{
		LOG(LogLevel::Error) << e.what();
		return false;
	}
	// Iterate through the source directory
	for (
		fs::directory_iterator file(source);
		file != fs::directory_iterator(); ++file
		)
	{
		try
		{
			fs::path current(file->path());
			if (fs::is_directory(current))
			{
				// Found directory: Recursion
				if (
					!copyDir(
					current,
					destination / current.filename()
					)
					)
				{
					return false;
				}
			}
			else
			{
				// Found file: Copy
				fs::copy_file(
					current,
					destination / current.filename()
					);
			}
		}
		catch (fs::filesystem_error const & e)
		{
			LOG(LogLevel::Error) << e.what();
		}
	}
	return true;
}


text:
/**
 * @file initcondlv_main.cc
 * @brief NPDE homework InitCondLV code
 * @author lfilippo, tille, jgacon, dcasati
 * @copyright Developed at ETH Zurich
 */

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

#include "initcondlv.h"

int main() {
  /* Compute initial conditions u0, v0 such that solution has period T=5
   * using inital approximation [3,2]^T */
  /* SAM_LISTING_BEGIN_2 */
  Eigen::Vector2d y(3, 2);  // Initial guess
  double T = 5;             // Period

#if SOLUTION
  // Compute Phi and W from initial guess
  std::pair<Eigen::Vector2d, Eigen::Matrix2d> PaW =
      InitCondLV::PhiAndW(y(0), y(1), T);
  Eigen::Vector2d F = PaW.first - y;  // Value of F at initial guess
  Eigen::Matrix2d DF;                 // Declare Jacobian of F

  // Until stopping condition
  while (F.norm() > 1e-5) {  // Residual based termination
    // Calculate Jacobian
    DF = PaW.second - Eigen::MatrixXd::Identity(2, 2);
    // Use Newton iteration
    y = y - DF.lu().solve(F);
    // Test current guess
    PaW = InitCondLV::PhiAndW(y(0), y(1), T);
    // Get value of F at the current guess
    F = PaW.first - y;
  }
#else
  //====================
  // Your code goes here
  //====================
#endif

  // Calculate the evolution of the Lotka-Volterra ODE up to
  // time 100, using the initial value y, found by the Newton iterations above.
  // If the evolution is indeed 5-periodic, then the solution at time 100
  // should have the value y. Print a warning message if this is not the case.
#if SOLUTION
  PaW = InitCondLV::PhiAndW(y(0), y(1), 100);
  if ((y - PaW.first).norm() > 1e-5) {
    std::cout << "Warning: Solution not periodic, y(100) != y(0)" << std::endl;
  }
#else
  //====================
  // Your code goes here
  //====================
#endif

  // Testing obtained solution against reference "exact" solution
  Eigen::Vector2d ref(3.1098751029156, 2.08097564048345);
  double tol = 1.0e-8;
  double error = (ref - y).lpNorm<Eigen::Infinity>();
  if (std::abs(error) > tol) {
    std::cout
        << " Error w.r.t. to reference solution is GREATER than tol=1.0e-8."
        << std::endl;
  } else {
    std::cout
        << " Error w.r.t. to reference solution is SMALLER than tol=1.0e-8."
        << std::endl;
  }

  /* SAM_LISTING_END_2 */
  return 0;
}

ext: cc