text:
//
//  main.cpp
//  fontgen - x2d
//
//  Created by Alex Kremer on 03/11/2012.
//  Copyright (c) 2012 godexsoft. All rights reserved.
//

#include <iostream>
#include <fstream>
#include <boost/program_options.hpp>

#include <ft2build.h>
#include FT_FREETYPE_H

#include "glm.hpp"
#include <png.h>

#define CHAR_CR ((uint32_t)'\r')
#define CHAR_LF ((uint32_t)'\n')

namespace po = boost::program_options;
using namespace std;

struct pixel
{
    pixel()
    : red(0)
    , green(0)
    , blue(0)
    , alpha(0)
    {
    }

    pixel(uint8_t r, uint8_t g, uint8_t b, uint8_t a)
    : red(r)
    , green(g)
    , blue(b)
    , alpha(a)
    {
    }

    uint8_t red;
    uint8_t green;
    uint8_t blue;
    uint8_t alpha;
};

class bitmap
{
public:
    bitmap(uint32_t w, uint32_t h)
    : width_(w)
    , height_(h)
    , pixels_(w*h)
    {
    }

    void set_pixel(uint32_t x, uint32_t y, const pixel& p);
    pixel get_pixel(uint32_t x, uint32_t y) const;

    const uint32_t get_width() const { return width_; }
    const uint32_t get_height() const { return height_; }

private:
    vector<pixel>       pixels_;
    uint32_t            width_;
    uint32_t            height_;
};

void bitmap::set_pixel(uint32_t x, uint32_t y, const pixel& p)
{
    assert(x < width_ && y < height_);
    pixels_[y*width_ + x] = p;
}

pixel bitmap::get_pixel(uint32_t x, uint32_t y) const
{
    return pixels_[y*width_ + x];
}

struct glyph_info
{
    uint32_t index;
    uint32_t code_point;

    glyph_info(uint32_t i, uint32_t cp)
    : index(i)
    , code_point(cp)
    {
    }
};

class generator
{
public:
    generator(const string& fnt, const vector<string>& txt_lst, int sz, int pd, bool mono, bool adv, string output)
    : font_file_(fnt)
    , txt_files_(txt_lst)
    , char_size_(sz, sz)
    , padding_(pd * 2.0f)
    , monospace_(mono)
    , static_advance_(adv)
    , font_name_(output)
    {
    }

    ~generator()
    {
        if(font_face_)
        {
            FT_Done_Face(font_face_);
        }

        if(font_library_)
        {
            FT_Done_FreeType(font_library_);
        }
    }

    int run();

private:

    void load_font();
    void process_character(int c);
    void parse_text(const string& file);
    void save_image(const bitmap& bm);
    void save_config();

    map<uint32_t, glyph_info>    gen_;
    string                       font_file_;
    vector<string>               txt_files_;

    string          font_name_;
    glm::vec2       char_size_;
    glm::vec2       char_spacing_;
    float           font_scale_;
    float           padding_;

    bool            monospace_;
    bool            static_advance_;

    FT_Library      font_library_;
    FT_Face         font_face_;

    stringstream    cfg_characters_;
    stringstream    cfg_widths_;
};

int generator::run()
{
    cout << "[+] Running generator...\n";

    cout << "    Font size:          " << char_size_.y << "\n";
    cout << "    Padding:            " << padding_ << "\n";
    cout << "    Font file:          " << font_file_ << "\n";

    for(vector<string>::iterator it = txt_files_.begin(); it != txt_files_.end(); ++it)
    {
        cout << "    Text file to parse: " << *it << "\n";
    }

    // open font
    load_font();

    // parse all text files to generate a map with glyphs
    for(vector<string>::iterator it = txt_files_.begin(); it != txt_files_.end(); ++it)
    {
        parse_text(*it);
    }

    uint8_t largest_width = 0;

    if(!monospace_)
    {
        for(map<uint32_t, glyph_info>::iterator it = gen_.begin(); it != gen_.end(); ++it)
        {
            glyph_info info = it->second;

            FT_Error err = FT_Load_Glyph(font_face_, info.index, FT_LOAD_RENDER);
            assert(!err);

            // load width of character
            uint8_t char_width = 0;

            if(static_advance_)
            {
                char_width = glm::ceil(font_scale_ * font_face_->glyph->advance.x);
            }
            else
            {
                char_width = glm::max(font_face_->glyph->bitmap_left, 0) + font_face_->glyph->bitmap.width;
            }

            largest_width = glm::max(largest_width, char_width);
        }

        char_size_.x = largest_width;
        cout << "[debug] Updated font width to: " << char_size_.x << "\n";
    }

    float w = glm::floor( glm::sqrt( (float)gen_.size() ) );
    float h = glm::ceil( (float)gen_.size() / w );

    uint32_t width =  (w * (char_size_.x + padding_)) + (char_spacing_.x * glm::max( w - 1.0f, 0.0f ));
    uint32_t height = (h * (char_size_.y + padding_)) + (char_spacing_.y * glm::max( h - 1.0f, 0.0f ));

    if(!monospace_)
    {
        uint32_t x, y;
        x = y = 0;
        y += char_spacing_.y;

        for(map<uint32_t, glyph_info>::iterator it = gen_.begin(); it != gen_.end(); ++it)
        {
            glyph_info info = it->second;
            cfg_characters_ << (char)it->first;

            FT_Error err = FT_Load_Glyph(font_face_, info.index, FT_LOAD_RENDER);
            assert(!err);

            // load width of character
            uint8_t char_width = 0;

            if(static_advance_)
            {
                char_width = glm::ceil(font_scale_ * font_face_->glyph->advance.x);
            }
            else
            {
                char_width = glm::max(font_face_->glyph->bitmap_left, 0) + font_face_->glyph->bitmap.width;

                if(char_width == 0)
                {
                    char_width = glm::ceil(font_scale_ * font_face_->glyph->advance.x);
                }
            }

            char_width = padding_ + ((char_width > 0) ? char_width : char_size_.x);

            // new line?
            if(x + char_width > width)
            {
                x = 0;
                y += char_size_.y + padding_ + char_spacing_.y;
            }

            x += char_width + char_spacing_.x;
        }

        height = y + char_size_.y + padding_ + char_spacing_.y;
    }

    cout << "[debug] width and height: " << width << "x" << height << "\n";

    uint32_t baseline = glm::ceil(font_scale_ * font_face_->ascender);
    cout << "[debug] baseline: " << baseline << "\n";

    bitmap bm(width, height);

    uint32_t idx = 0;
    uint32_t x = 0;
    uint32_t y = baseline;

    for(map<uint32_t, glyph_info>::iterator it = gen_.begin(); it != gen_.end(); ++it, ++idx)
    {
        glyph_info info = it->second;
        cfg_characters_ << (char)it->first;

        FT_Error err = FT_Load_Glyph(font_face_, info.index, FT_LOAD_RENDER);
        assert(!err);

        uint32_t char_width = 0;

        if(monospace_)
        {
            char_width = char_size_.x + padding_;
        }
        else
        {
            if(static_advance_)
            {
                char_width = glm::ceil(font_scale_ * font_face_->glyph->advance.x);
            }
            else
            {
                char_width = glm::max(font_face_->glyph->bitmap_left, 0) + font_face_->glyph->bitmap.width;

                if(char_width == 0)
                {
                    char_width = glm::ceil(font_scale_ * font_face_->glyph->advance.x);
                }
            }

            char_width = padding_ + ((char_width > 0) ? char_width : char_size_.x);
            cfg_widths_ << char_width << " ";
        }

        // new line?
        if(x + char_width > width)
        {
            x = 0;
            y += char_size_.y + padding_ + char_spacing_.y;
        }

        if(font_face_->glyph->bitmap.buffer)
        {
            uint32_t adj_x = x + (0.5f * padding_) + glm::max(font_face_->glyph->bitmap_left, 0);
            uint32_t adj_y = y + (0.5f * padding_) - glm::min((uint32_t)font_face_->glyph->bitmap_top, baseline);

            uint8_t* src = (unsigned char *)(font_face_->glyph->bitmap.buffer);

            for(uint32_t i = 0; i < font_face_->glyph->bitmap.rows; ++i,
                src += glm::max(font_face_->glyph->bitmap.pitch, -font_face_->glyph->bitmap.pitch) - font_face_->glyph->bitmap.width)
            {
                for(uint32_t j = 0; j < font_face_->glyph->bitmap.width; j++, src++)
                {
                    bm.set_pixel(j + adj_x, i + adj_y, pixel(255, 255, 255, src[0]));
                }
            }
        }

        x += char_width + char_spacing_.x;
    }

    // now we have a filled bitmap. need to save it into a file
    char_size_.x += padding_;
    char_size_.y += padding_;

    // time to save image and config
    cout << "[+] Writing files...\n";
    save_image(bm);
    save_config();

    return 0;
}

void generator::save_config()
{
    stringstream ss;

    ss << "<x2d version=\"1.0\">\n";

    ss << "    <texture n = \"" << font_name_ << "_sheet\" path = \""
       << "res/fonts/" << font_name_ << ".png\" />\n\n";

    ss << "    <font    n = \"" << font_name_ << "\"\n";
    ss << "        height = \"" << char_size_.y << "\"\n";
    ss << "       spacing = \"" << char_spacing_.x << " " << char_spacing_.y << "\"\n";
    ss << "       texture = \"" << font_name_ << "_sheet\"\n";
    ss << "    >\n";

    ss << "        <string n = \"characters\">\n";
    ss << "            <![CDATA[" << cfg_characters_.str() << "]]>\n";
    ss << "        </string>\n";
    ss << "        <string n = \"widths\">\n";
    ss << "            " << cfg_widths_.str() << "\n";
    ss << "        </string>\n";

    ss << "    </font>\n";
    ss << "</x2d>\n";

    ofstream out(string(font_name_ + ".xml").c_str());
    if(!out.is_open())
    {
        throw runtime_error("Couldn't open configuration output file. aborting.");
    }

    string str(ss.str());
    out.write(str.c_str(), str.size());
}

void generator::save_image(const bitmap& bm)
{
    FILE* fp = fopen(string(font_name_ + ".png").c_str(), "wb");
    if(!fp)
    {
        throw runtime_error("Can't open png output file.");
    }

    png_structp png_ptr = NULL;
    png_infop info_ptr = NULL;
    png_byte ** row_pointers = NULL;

    int status = -1;

    int pixel_size = 4;
    int depth = 8;

    png_ptr = png_create_write_struct (PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
    if (png_ptr == NULL)
    {
        throw runtime_error("Can't open png output structure.");
    }

    info_ptr = png_create_info_struct (png_ptr);
    if (info_ptr == NULL)
    {
        throw runtime_error("Can't create png info struct.");
    }

    png_set_IHDR (png_ptr, info_ptr, bm.get_width(), bm.get_height(), depth,
        PNG_COLOR_TYPE_RGBA, PNG_INTERLACE_NONE,
        PNG_COMPRESSION_TYPE_DEFAULT, PNG_FILTER_TYPE_DEFAULT);

    row_pointers = static_cast<png_byte**>(png_malloc (png_ptr, bm.get_height() * sizeof (png_byte *)));

    for (uint32_t y = 0; y < bm.get_height(); ++y)
    {
        png_byte *row =
            static_cast<png_byte*>(png_malloc (png_ptr, sizeof (uint8_t) * bm.get_width() * pixel_size));
        row_pointers[y] = row;
        for (uint32_t x = 0; x < bm.get_width(); ++x) {
            pixel p = bm.get_pixel(x, y);
            *row++ = p.red;
            *row++ = p.green;
            *row++ = p.blue;
            *row++ = p.alpha;
        }
    }

    png_init_io (png_ptr, fp);
    png_set_rows (png_ptr, info_ptr, row_pointers);
    png_write_png (png_ptr, info_ptr, PNG_TRANSFORM_IDENTITY, NULL);

    status = 0;

    for (size_t y = 0; y < bm.get_height(); y++) {
        png_free (png_ptr, row_pointers[y]);
    }

    png_free (png_ptr, row_pointers);

    png_destroy_write_struct (&png_ptr, &info_ptr);
    fclose (fp);
}

void generator::load_font()
{
    if(!FT_Init_FreeType( &font_library_ ))
    {
        if(!FT_New_Face(font_library_, font_file_.c_str(), 0, &font_face_))
        {
            if(!FT_Select_Charmap(font_face_, ft_encoding_unicode))
            {
                char_spacing_.x = 2.0f * glm::ceil(char_size_.x * (1.0f / 32.0f));
                char_spacing_.y = 2.0f * glm::ceil(char_size_.y * (1.0f / 32.0f));

                cout << "[?] Spacing: x=" << char_spacing_.x << " / y=" << char_spacing_.y << "\n";

                font_scale_ = char_size_.y / (font_face_->bbox.yMax - font_face_->bbox.yMin);

                cout << "[?] Setting pixel sizes: " << char_size_.x - 2 << " x " << (FT_UInt)char_size_.y - 2 << "\n";

                FT_Set_Pixel_Sizes(font_face_, (FT_UInt)char_size_.x - 2, (FT_UInt)char_size_.y - 2);
                cout << "[+] Font loaded.\n";
            }
        }
    }
}

void generator::process_character(int c)
{
    uint32_t char_code_point = c;

    if(!gen_.count(char_code_point))
    {
        if(char_code_point != CHAR_CR && char_code_point != CHAR_LF)
        {
            uint32_t glyph_index = FT_Get_Char_Index(font_face_, char_code_point);
            if(glyph_index)
            {
//                cout << "[debug] Found glyph index " << glyph_index << " for character '" << (char)c << "'\n";

                // add to map
                gen_.insert( make_pair( char_code_point, glyph_info( glyph_index, char_code_point ) ) );
            }
        }
    }
}

void generator::parse_text(const string& file)
{
    cout << "[+] Parsing text file: " << file << "...\n";
    ifstream ifs(file.c_str(), ifstream::in);

    if(ifs.is_open())
    {
        while(!ifs.eof())
        {
            // TODO: most likely this wont work for anything other than ascii. need some UTF-8 magic :)
            process_character(ifs.get());
        }

        cout << "[+] Done. Added " << gen_.size() << " characters.\n";
    }
    else
    {
        cout << "[-] Could'n open text file " << file << ".\n";
    }
}

int main(int argc, const char * argv[])
{
    po::options_description desc("Allowed options");
    desc.add_options()
    ("help,h", "produce help message")
    ("output,o", po::value<std::string>()->default_value("x2d_font"), "output name")
    ("size,s", po::value<int>()->default_value(32), "height of characters")
    ("padding,p", po::value<int>()->default_value(0), "character padding")
    ("font,f", po::value<string>(), "input font file (.ttf)")
    ("text,t", po::value<vector<string> >()->composing(), "input text file")
    ("monospace,m", "monospaced font")
    ("advance,a", "use glyph advance values for non-monospace fonts")
    ;

    po::variables_map vm;
    po::store(po::parse_command_line(argc, argv, desc), vm);
    po::notify(vm);

    string fnt;
    vector<string> txt_files;
    int size;
    int padding;

    if (vm.count("help"))
    {
        cout << desc << "\n";
        return 1;
    }

    // check font file
    if (vm.count("font"))
    {
        fnt = vm["font"].as<string>();
    }
    else
    {
        cout << "[-] Input font is not passed.\n";
        return 127;
    }

    // set size and padding
    size = vm["size"].as<int>();
    padding = vm["padding"].as<int>();

    // check input text files
    if (vm.count("text"))
    {
        txt_files = vm["text"].as<vector<string> >();
    }
    else
    {
        cout << "[-] Text files not specified. abort.\n";
        return 127;
    }

    // create a generator
    generator gen(fnt, txt_files, size, padding, vm.count("monospace")>0,
                  vm.count("advance")>0, vm["output"].as<string>());
    return gen.run();
}


text:
/*
    Copyright (C) 2008 Wei Dong <wdong@princeton.edu>. All Rights Reserved.

    This file is part of LSHKIT.

    LSHKIT is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    LSHKIT is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with LSHKIT.  If not, see <http://www.gnu.org/licenses/>.
*/

#include <boost/program_options.hpp>
#include <boost/progress.hpp>
#include <boost/format.hpp>
#include <boost/timer.hpp>
#include <lshkit.h>

/**
  * \file mplsh-run.cpp
  * \brief Example of using MPLSH.
  *
  * This program is an example of using MPLSH index.
  *
  * The program reconstruct the LSH index by default.  You can give the
  * --index option to make the program save the LSH index.  The next
  * time you run the program with the same --index option, the program
  * will try to load the previously saved index.  When a saved index is
  * used, you need to make sure that the dataset and other parameters match
  * the previous run.  However, the benchmark file, Q and K can be different.
  *
\verbatim
Allowed options:
  -h [ --help ]                   produce help message.
  -W [ -- ] arg (=1)
  -M [ -- ] arg (=1)
  -T [ -- ] arg (=1)              # probes
  -L [ -- ] arg (=1)              # hash tables
  -Q [ -- ] arg (=100)            # queries
  -K [ -- ] arg (=50)             # nearest neighbor to retrieve
  -R [ -- ] arg (=3.40282347e+38) R-NN distance range
  --recall arg                    desired recall
  -D [ --data ] arg               data file
  -B [ --benchmark ] arg          benchmark file
  --index arg                     index file
  -H [ -- ] arg (=1017881)        hash table size, use the default value.
\endverbatim
  */

using namespace std;
using namespace lshkit;
namespace po = boost::program_options;

/*
    You must provide an access class to query the MPLSH.
    MPLSH only saves keys (pointers to the real feature vectors) in the
    hash tables and it will rely on the accessor class to retrieve
    the feature vector.

    An accessor must provide three methods:

        bool mark (unsigned key);

        mark that key has been accessed.  If key has already been marked, return false,
        otherwise return true.  MPLSH will use this to avoid scanning the data more than
        one time per query.

        void reset ();

        to clear all the marks.

        const float *operator () (unsigned key);

        given a key, return the pointer to a feature vector.

    The MatrixAccessor class is used to access feature vectors stored in a Matrix.
*/

/* This class has been merged to include/lshkit/matrix.h */
/*
class MatrixAccessor
{
    const Matrix<float> &matrix_;
    boost::dynamic_bitset<> flags_;
public:
    typedef unsigned Key;
    MatrixAccessor(const Matrix<float> &matrix)
        : matrix_(matrix), flags_(matrix.getSize()) {}

    MatrixAccessor(const MatrixAccessor &accessor)
        : matrix_(accessor.matrix_), flags_(matrix_.getSize()) {}

    void reset ()
    {
        flags_.reset();
    }

    bool mark (unsigned key)
    {
        if (flags_[key]) return false;
        flags_.set(key);
        return true;
    }

    const float *operator () (unsigned key)
    {
        return matrix_[key];
    }
};
*/

int main (int argc, char *argv[])
{
    string data_file;
    string benchmark;
    string index_file;

    float W, R, desired_recall = 1.0;
    unsigned M, L, H;
    unsigned Q, K, T;
    bool do_recall = false;
    bool do_benchmark = true;
    bool use_index = false; // load the index from a file

    boost::timer timer;

    po::options_description desc("Allowed options");
    desc.add_options()
        ("help,h", "produce help message.")
        (",W", po::value<float>(&W)->default_value(1.0), "")
        (",M", po::value<unsigned>(&M)->default_value(1), "")
        (",T", po::value<unsigned>(&T)->default_value(1), "# probes")
        (",L", po::value<unsigned>(&L)->default_value(1), "# hash tables")
        (",Q", po::value<unsigned>(&Q)->default_value(100), "# queries")
        (",K", po::value<unsigned>(&K)->default_value(0), "# nearest neighbor to retrieve")
        ("radius,R", po::value<float>(&R)->default_value(numeric_limits<float>::max()), "R-NN distance range (L2)")
        ("recall", po::value<float>(&desired_recall), "desired recall")
        ("data,D", po::value<string>(&data_file), "data file")
        ("benchmark,B", po::value<string>(&benchmark), "benchmark file")
        ("index", po::value<string>(&index_file), "index file")
        (",H", po::value<unsigned>(&H)->default_value(1017881), "hash table size, use the default value.")
        ;

    po::variables_map vm;
    po::store(po::parse_command_line(argc, argv, desc), vm);
    po::notify(vm);

    if (vm.count("help") || (vm.count("data") < 1))
    {
        cout << desc;
        return 0;
    }

    if (vm.count("radius") >= 1) {
        R *= R; // we use L2sqr in the program.
    }

    if (vm.count("recall") >= 1)
    {
        do_recall = true;
        if (K == 0) {
            cerr << "Automatic probing does not support R-NN query." << endl;
        }
    }

    if ((Q == 0) || (vm.count("benchmark") == 0)) {
        do_benchmark = false;
    }

    if (vm.count("index") == 1) {
        use_index = true;
    }

    cout << "LOADING DATA..." << endl;
    timer.restart();
    FloatMatrix data(data_file);
    cout << boost::format("LOAD TIME: %1%s.") % timer.elapsed() << endl;

    typedef MultiProbeLshIndex<unsigned> Index;

    FloatMatrix::Accessor accessor(data);
    Index index;

    // try loading index
    bool index_loaded = false;

    if (use_index) {
        ifstream is(index_file.c_str(), ios_base::binary);
        if (is) {
            is.exceptions(ios_base::eofbit | ios_base::failbit | ios_base::badbit);
            cout << "LOADING INDEX..." << endl;
            timer.restart();
            index.load(is);
            BOOST_VERIFY(is);
            cout << boost::format("LOAD TIME: %1%s.") % timer.elapsed() << endl;
            index_loaded = true;
        }
    }

    if (!index_loaded) {
        // We define a short name for the MPLSH index.
        Index::Parameter param;

        // Setup the parameters.  Note that L is not provided here.
        param.W = W;
        param.range = H; // See H in the program parameters.  You can just use the default value.
        param.repeat = M;
        param.dim = data.getDim();
        DefaultRng rng;

        index.init(param, rng, L);
        // The accessor.

        // Initialize the index structure.  Note L is passed here.
        cout << "CONSTRUCTING INDEX..." << endl;

        timer.restart();
        {
            boost::progress_display progress(data.getSize());
            for (int i = 0; i < data.getSize(); ++i)
            {
                // Insert an item to the hash table.
                // Note that only the key is passed in here.
                // MPLSH will get the feature from the accessor.
                index.insert(i, data[i]);
                ++progress;
            }
        }
        cout << boost::format("CONSTRUCTION TIME: %1%s.") % timer.elapsed() << endl;

        if (use_index) {
            timer.restart();
            cout << "SAVING INDEX..." << endl;
            {
                ofstream os(index_file.c_str(), ios_base::binary);
                os.exceptions(ios_base::eofbit | ios_base::failbit | ios_base::badbit);
                index.save(os);
            }
            cout << boost::format("SAVING TIME: %1%s") % timer.elapsed() << endl;
        }
    }

    if (do_benchmark) {

        Benchmark<> bench;
        cout << "LOADING BENCHMARK..." << endl;
        bench.load(benchmark);
        bench.resize(Q, K);
        cout << "DONE." << endl;

        for (unsigned i = 0; i < Q; ++i)
        {
            for (unsigned j = 0; j < K; ++j)
            {
                assert(bench.getAnswer(i)[j].key < data.getSize());
            }
        }

        cout << "RUNNING QUERIES..." << endl;

        Stat recall;
        Stat cost;
        metric::l2sqr<float> l2sqr(data.getDim());
        TopkScanner<FloatMatrix::Accessor, metric::l2sqr<float> > query(accessor, l2sqr, K, R);
        vector<Topk<unsigned> > topks(Q);

        timer.restart();
        if (do_recall)
            // Specify the required recall
            // and let MPLSH to guess how many bins to probe.
        {
            boost::progress_display progress(Q);
            for (unsigned i = 0; i < Q; ++i)
            {
                // Query for one point.
                query.reset(data[bench.getQuery(i)]);
                index.query_recall(data[bench.getQuery(i)], desired_recall, query);
                cost << double(query.cnt())/double(data.getSize());
                topks[i].swap(query.topk());
                ++progress;
            }
        }
        else
            // specify how many bins to probe.
        {
            boost::progress_display progress(Q);
            for (unsigned i = 0; i < Q; ++i)
            {
                query.reset(data[bench.getQuery(i)]);
                index.query(data[bench.getQuery(i)], T, query);
                cost << double(query.cnt())/double(data.getSize());
                topks[i].swap(query.topk());
                ++progress;
            }
        }

        for (unsigned i = 0; i < Q; ++i) {
            recall << bench.getAnswer(i).recall(topks[i]);
        }

        cout << boost::format("QUERY TIME: %1%s.") % timer.elapsed() << endl;

        cout << "[RECALL] " << recall.getAvg() << " +/- " << recall.getStd() << endl;
        cout << "[COST] " << cost.getAvg() << " +/- " << cost.getStd() << endl;

    }

    return 0;
}

ext: cpp

max_stars_repo_name: huonw/nmslib

max_issues_repo_name: huonw/nmslib

max_forks_repo_name: huonw/nmslib


text:
// Boost.Bimap
//
// Copyright (c) 2006-2007 Matias Capeletto
//
// 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)

//  VC++ 8.0 warns on usage of certain Standard Library and API functions that
//  can be cause buffer overruns or other possible security issues if misused.
//  See http://msdn.microsoft.com/msdnmag/issues/05/05/SafeCandC/default.aspx
//  But the wording of the warning is misleading and unsettling, there are no
//  portable alternative functions, and VC++ 8.0's own libraries use the
//  functions in question. So turn off the warnings.
#define _CRT_SECURE_NO_DEPRECATE
#define _SCL_SECURE_NO_DEPRECATE

// Boost.Bimap Example
//-----------------------------------------------------------------------------

#include <boost/config.hpp>

#include <iostream>
#include <string>

#include <boost/bimap/bimap.hpp>
#include <boost/bimap/unordered_set_of.hpp>
#include <boost/bimap/multiset_of.hpp>

#include <boost/optional.hpp>
#include <boost/none.hpp>
#include <boost/foreach.hpp>
#include <boost/assign/list_inserter.hpp>

using namespace boost::bimaps;
using namespace boost;
using namespace std;

int main()
{
    {

    typedef bimap<

        string,
        multiset_of< optional<string> >

    > bm_type;

    bm_type bm;

    assign::insert( bm )

        ( "John" , string("lazarus" ) )
        ( "Peter", string("vinicius") )
        ( "Simon", string("vinicius") )
        ( "Brian", none               )
    ;

    cout << "John is working in "
         << bm.left.at( "John" ).get_value_or( "no project" )
         << endl;

    cout << "Project vinicius is being developed by " << endl;
    BOOST_FOREACH( bm_type::right_reference rp,
                   bm.right.equal_range( std::string("vinicius") ) )
    {
        cout << rp.second << endl;
    }

    cout << "This workers need a project " << endl;
    BOOST_FOREACH( bm_type::right_reference rp,
                   bm.right.equal_range(none) )
    {
        cout << rp.second << endl;
    }

}

    //[ code_population_bimap

    typedef bimap<

        unordered_set_of< std::string >,
        multiset_of< long, std::greater<long> >

    > population_bimap;

    typedef population_bimap::value_type population;

    population_bimap pop;
    pop.insert( population("China",          1321000000) );
    pop.insert( population("India",          1129000000) );
    pop.insert( population("United States",   301950000) );
    pop.insert( population("Indonesia",       234950000) );
    pop.insert( population("Brazil",          186500000) );
    pop.insert( population("Pakistan",        163630000) );

    std::cout << "Countries by their population:" << std::endl;

    // First requirement
    /*<< The right map view works like a
         `std::multimap< long, std::string, std::greater<long> >`,
         We can iterate over it to print the results in the required order. >>*/
    for( population_bimap::right_const_iterator
            i = pop.right.begin(), iend = pop.right.end();
            i != iend ; ++i )
    {
        std::cout << i->second << " with " << i->first << std::endl;
    }

    // Second requirement
    /*<< The left map view works like a `std::unordered_map< std::string, long >`,
         given the name of the country we can use it to search for the population
         in constant time >>*/
    std::cout << "Population of China: " << pop.left.at("China") << std::endl;
    //]

    return 0;
}

max_stars_repo_name: zyiacas/boost-doc-zh

max_issues_repo_name: sdfict/boost-doc-zh

max_forks_repo_name: sdfict/boost-doc-zh


text:
/**
 * \file OutWavFilter.cpp
 */

#include "OutWavFilter.h"

#include <cstring>
#include <stdexcept>

#include <boost/lexical_cast.hpp>

#include <ATK/Core/Utilities.h>

namespace
{
  template<typename DataType1, typename DataType2>
  void convert(std::vector<std::vector<DataType1> >& outputs, const std::vector<char>& inputs)
  {
    int nbChannels = outputs.size();
    int64_t size = outputs[0].size();
    for(int j = 0; j < nbChannels; ++j)
    {
      ATK::ConversionUtilities<DataType2, DataType1>::convert_array(reinterpret_cast<const DataType2*>(inputs.data() + j * sizeof(DataType2)), outputs[j].data(), size, nbChannels);
    }
  }
}

namespace ATK
{
  template<typename DataType>
  OutWavFilter<DataType>::OutWavFilter(const std::string& filename)
  :TypedBaseFilter<DataType>(0, 0)
  {
    wavstream.open(filename.c_str(), std::ios_base::binary);
    if(!wavstream.good())
    {
      throw std::runtime_error("Could not WAV file " + filename);
    }
  }

  template<typename DataType>
  void OutWavFilter<DataType>::process_impl(int64_t size) const
  {
    std::size_t nb_inputs = converted_inputs.size();
    std::vector<DataType> buffer(nb_inputs * size);

    for(int64_t i = 0; i < size; ++i)
    {
      for(std::size_t j = 0; j < nb_inputs; ++j)
      {
        buffer[j + i * nb_inputs] = converted_inputs[j][i];
      }
    }

    wavstream.write(reinterpret_cast<const char*>(buffer.data()), buffer.size() * sizeof(DataType));
    write_header();
  }

  template<typename DataType>
  void OutWavFilter<DataType>::set_nb_input_ports(int nb_ports)
  {
    Parent::set_nb_input_ports(nb_ports);
    setup();
  }

  template<typename DataType>
  void OutWavFilter<DataType>::setup()
  {
    write_header();
  }

  template<typename DataType>
  void OutWavFilter<DataType>::write_header() const
  {
    WavHeader header;
    WavFormat format;
    WavData data;

    std::strncpy(header.FileTypeBlocID, "RIFF", 4);
    std::strncpy(header.FileFormatID, "WAVE", 4);
    std::strncpy(format.FormatBlocID, "fmt ", 4);
    std::strncpy(data.DataBlocID, "data", 4);
    format.AudioFormat = WavTraits<DataType>::get_wav_type();
    format.Frequence = input_sampling_rate;
    format.BitsPerSample = sizeof(DataType)* 8;
    format.BytePerBloc = format.NbChannels * format.BitsPerSample / 8;
    format.BytePerSec = format.BytePerBloc * input_sampling_rate;
    format.NbChannels = nb_input_ports;

    std::size_t total_size = wavstream.tellp();
    std::size_t bloc_size = sizeof(WavFormat);
    std::size_t data_size = total_size - sizeof(WavFormat) - sizeof(WavHeader) - sizeof(WavData);
    wavstream.seekp(0);
    header.FileSize = static_cast<std::int32_t>(total_size - 8);
    format.BlocSize = static_cast<std::int32_t>(bloc_size - 8);
    data.DataSize = static_cast<std::int32_t>(data_size);
    wavstream.write(reinterpret_cast<char*>(&header), sizeof(WavHeader));
    wavstream.write(reinterpret_cast<char*>(&format), sizeof(WavFormat));
    wavstream.write(reinterpret_cast<char*>(&data), sizeof(WavData));

    wavstream.seekp(0, std::ios_base::end);
  }

  template class OutWavFilter<std::int16_t>;
  template class OutWavFilter<float>;
  template class OutWavFilter<double>;
}


text:
#include "latex.hpp"

#include <boost/filesystem.hpp>
#include <cstdlib>
#include <fstream>
#include <iostream>
#include <libplatform/libplatform.h>
#include <wkhtmltox/image.h>

std::string Latex::_find_katex_path()
{
	for (std::string dir = "./", end = "./../../"; dir != end; dir += "../")
	{
		for (const auto& file : boost::filesystem::directory_iterator(dir))
		{
			auto path = file.path();

			if (path.stem() == "katex" && boost::filesystem::is_directory(path))
			{
				return path.string();
			}
		}
	}

	// We don't need the CSS necessarily, but do need the JavaScript.
	throw ExistentialException("Could not find KaTeX directory!");
}

std::string Latex::_katex_path;

Latex::V8 Latex::_v8;

Latex::Latex(WarningBehavior behavior)
: Latex("katex/katex.min.css", behavior)
{ }

Latex::Latex(const std::string& stylesheet, WarningBehavior behavior)
: _stylesheet(stylesheet)
, _warning_behaviour(behavior)
, _isolate(_new_isolate())
{
	if (_katex_path.empty())
	{
		_katex_path = _find_katex_path();
	}

	v8::HandleScope handle_scope(_isolate);

	v8::Isolate::Scope isolate_scope(_isolate);

	auto context = v8::Context::New(_isolate);

	v8::Context::Scope context_scope(context);

	_load_katex(context);

	_persistent_context = v8::UniquePersistent<v8::Context>(_isolate, context);

	wkhtmltoimage_init(false);
}

Latex::Latex(const Latex& other)
: Latex(other._stylesheet,
		other._warning_behaviour)
{
	_additional_css = other._additional_css;
}

Latex::Latex(Latex&& other) noexcept
: Latex()
{
	swap(other);
}

Latex& Latex::operator=(Latex other)
{
	swap(other);

	return *this;
}

void Latex::swap(Latex &other) noexcept
{
	// Enable ADL
	using std::swap;

	swap(_allocator, other._allocator);

	swap(_isolate, other._isolate);

	swap(_persistent_context, other._persistent_context);

	swap(_stylesheet, other._stylesheet);

	swap(_additional_css, other._additional_css);

	swap(_warning_behaviour, other._warning_behaviour);
}

void swap(Latex& first, Latex& second) noexcept
{
	first.swap(second);
}

Latex::~Latex()
{
	wkhtmltoimage_deinit();
}

std::string Latex::to_html(const std::string& latex) const
{
	static const std::string arguments = "{'displayMode': true}";

	v8::Isolate::Scope isolate_scope(_isolate);

	// Stack-allocated handle-scope (takes care of handles such
	// that object are garbage-collected after the scope ends)
	v8::HandleScope handle_scope(_isolate);

	// Get a local context handle from the persistent handle.
	auto context = v8::Local<v8::Context>::New(_isolate,
											   _persistent_context);

	v8::Context::Scope context_scope(context);

	std::string source = "katex.renderToString('";

	source += _escape(latex) + "', ";
	source += arguments + ");";

	auto value = _run(source, context);

	std::string html = *static_cast<v8::String::Utf8Value>(value);

	return "
\n" + html + "
\n";
}

std::string Latex::to_complete_html(const std::string &latex) const
{
	auto snippet = to_html(latex);

	std::string html = "<!DOCTYPE html>\n\n";

	html += "\n<meta charset='utf-8'/>\n";
	html += "<link rel='stylesheet' type='text/css' ";
 	html += "href='" + _stylesheet + "'>\n";

	if (! _additional_css.empty())
	{
		html += "\n";
	}

	html += "\n\n";
	html += snippet;
	html += "\n";

	return html;
}

void Latex::to_image(const std::string &latex,
				  const std::string &filepath,
				  ImageFormat format) const
{
	std::ofstream temp("temp.html");

	if (! temp) throw FileException("Could not open temporary file!");

	temp << to_complete_html(latex);

	temp.close();

	auto converter = _new_converter(filepath, format);

	if (! wkhtmltoimage_convert(converter))
	{
		throw ConversionException("Could not convert to png!");
	}

	wkhtmltoimage_destroy_converter(converter);

	boost::filesystem::remove("temp.html");
}

void Latex::to_png(const std::string &latex,
				const std::string &filepath) const
{
	to_image(latex, filepath, ImageFormat::PNG);
}

void Latex::to_jpg(const std::string &latex,
				const std::string &filepath) const
{
	to_image(latex, filepath, ImageFormat::JPG);
}

void Latex::to_svg(const std::string &latex,
				const std::string &filepath) const
{
	to_image(latex, filepath, ImageFormat::SVG);
}

void Latex::add_css(const std::string& css)
{
	_additional_css += css;
}

const std::string& Latex::additional_css() const
{
	return _additional_css;
}

void Latex::clear_additional_css()
{
	_additional_css.clear();
}

const std::string& Latex::stylesheet() const
{
	return _stylesheet;
}

void Latex::stylesheet(const std::string& stylesheet)
{
	_stylesheet = stylesheet;
}

const Latex::WarningBehavior& Latex::warning_behavior() const
{
	return _warning_behaviour;
}

void Latex::warning_behavior(WarningBehavior behavior)
{
	_warning_behaviour = behavior;
}

v8::Isolate* Latex::_new_isolate() const
{
	v8::Isolate::CreateParams parameters;

	parameters.array_buffer_allocator = &_allocator;

	// Isolated JavaScript Virtual Environment
	return v8::Isolate::New(parameters);
}

void Latex::_load_katex(const v8::Local<v8::Context>& context) const
{
	std::ifstream file(_katex_path + "/katex.min.js");

	std::string source;
	std::string input;

	while (std::getline(file, input))
	{
		source += input + " ";
	}

	_run(source, context);
}

v8::Local<v8::Value> Latex::_run(const std::string& source,
								 const v8::Local<v8::Context>& context) const
{
	v8::EscapableHandleScope handle_scope(_isolate);

	auto unchecked = v8::String::NewFromUtf8(_isolate,
											 source.c_str(),
											 v8::NewStringType::kNormal);

	auto checked = unchecked.ToLocalChecked();

	// Compile the source code.
	auto script = v8::Script::Compile(context, checked).ToLocalChecked();

	// V8 engine's try-catch mechanism
	v8::TryCatch try_catch(_isolate);

	auto result = script->Run(context);

	if (result.IsEmpty())
	{
		// Grab last exception
		auto exception = try_catch.Exception();

		std::string what = *static_cast<v8::String::Utf8Value>(exception);

		// Remove the 'ParseError' (redundant)
		throw ParseException(what.substr(12));
	}

	// Allows us to return local-scope objects to the outside scope
	return handle_scope.Escape(result.ToLocalChecked());
}

std::string Latex::_escape(std::string source) const
{
	for (auto i = source.begin(); i != source.end(); ++i)
	{
		if (*i == '\\')
		{
			i = source.insert(i, '\\');

			++i;
		}
	}

	return source;
}

void _throw(wkhtmltoimage_converter*, const char* message)
{
	throw Latex::ConversionException(message);
}

void _log(wkhtmltoimage_converter*, const char* message)
{
	std::clog << message << std::endl;
}

wkhtmltoimage_converter*
Latex::_new_converter(const std::string& filepath, ImageFormat format) const
{
	auto settings = _new_converter_settings(filepath, format);

	auto converter = wkhtmltoimage_create_converter(settings, nullptr);

	wkhtmltoimage_set_error_callback(converter, _throw);

	if (_warning_behaviour == WarningBehavior::Strict)
	{
		wkhtmltoimage_set_error_callback(converter, _throw);
	}

	else if (_warning_behaviour == WarningBehavior::Log)
	{
		wkhtmltoimage_set_warning_callback(converter, _log);
	}

	return converter;
}

wkhtmltoimage_global_settings*
Latex::_new_converter_settings(const std::string& filepath,
					 Latex::ImageFormat format) const
{
	auto settings = wkhtmltoimage_create_global_settings();

	wkhtmltoimage_set_global_setting(settings,
									 "transparent",
									 "false");

	wkhtmltoimage_set_global_setting(settings,
									 "in",
									 "temp.html");

	wkhtmltoimage_set_global_setting(settings,
									 "out",
									 filepath.c_str());
	const char* fmt{ "png" };

	switch (format)
	{
		case ImageFormat::PNG: fmt = "png"; break;

		case ImageFormat::JPG: fmt = "jpg"; break;

		case ImageFormat::SVG: fmt = "svg"; break;
	}

	wkhtmltoimage_set_global_setting(settings,
									 "fmt",
									 fmt);

	wkhtmltoimage_set_global_setting(settings,
									 "screenWidth",
									 "0");

	wkhtmltoimage_set_global_setting(settings,
									 "quality",
									 "100");

	return settings;
}

Latex::V8::V8()
: platform(v8::platform::CreateDefaultPlatform())
{
	v8::V8::InitializeICU();

	v8::V8::InitializePlatform(platform.get());

	v8::V8::Initialize();
}

Latex::V8::~V8()
{
	v8::V8::Dispose();

	v8::V8::ShutdownPlatform();
}

void* Latex::Allocator::Allocate(size_t length)
{
	auto data = AllocateUninitialized(length);

	return data ? std::memset(data, 0, length) : data;
}

void Latex::Allocator::Free(void *data, size_t)
{
	free(data);
}

void* Latex::Allocator::AllocateUninitialized(size_t length)
{
	return malloc(length);
}

max_stars_repo_name: ggarciatisa/latexpp

max_issues_repo_name: ggarciatisa/latexpp

max_forks_repo_name: ggarciatisa/latexpp


text:
/*
 * Copyright 2015 David A. Boyuka II
 *
 * 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.
 */

/*
 * cblq-setops.hpp
 *
 *  Created on: Mar 25, 2014
 *      Author: David A. Boyuka II
 */
#ifndef CBLQ_SETOPS_HPP_
#define CBLQ_SETOPS_HPP_

#include <cstdlib>
#include <vector>
#include <boost/smart_ptr.hpp>

#include "pique/region/region-encoding.hpp"
#include "pique/region/cblq/cblq.hpp"
#include "pique/setops/setops.hpp"

struct CBLQSetOperationsConfig {
	CBLQSetOperationsConfig(bool compact_after_setop, bool no_compact_during_nary = true) :
		compact_after_setop(compact_after_setop), no_compact_during_nary(no_compact_during_nary)
	{}

	bool compact_after_setop;
	bool no_compact_during_nary;
};

template<int ndim>
class CBLQSetOperations : public SetOperations< CBLQRegionEncoding<ndim> > {
public:
	static constexpr int NDIM = ndim;

	typedef CBLQSetOperationsConfig SetOperationsConfig;
	CBLQSetOperations(CBLQSetOperationsConfig conf) : conf(conf), suppress_compact(false) {}
	virtual ~CBLQSetOperations() {}

protected:
    using cblq_word_t = typename CBLQRegionEncoding<ndim>::cblq_word_t;
    using typename SetOperations< CBLQRegionEncoding<ndim> >::RegionEncodingCPtrCIter;

private:
    virtual boost::shared_ptr< CBLQRegionEncoding<ndim> > unary_set_op_impl(boost::shared_ptr< const CBLQRegionEncoding<ndim> > region, UnarySetOperation op) const;
    virtual boost::shared_ptr< CBLQRegionEncoding<ndim> > binary_set_op_impl(boost::shared_ptr< const CBLQRegionEncoding<ndim> > left, boost::shared_ptr< const CBLQRegionEncoding<ndim> > right, NArySetOperation op) const;

    virtual boost::shared_ptr< CBLQRegionEncoding<ndim> > nary_set_op_impl(RegionEncodingCPtrCIter region_it, RegionEncodingCPtrCIter region_end_it, NArySetOperation op) const;

protected:
    const CBLQSetOperationsConfig conf;
    mutable bool suppress_compact;

    template<int ndim2> friend class CBLQSetOperationsFast; // Allows delegation upward to this class
};

template<int ndim>
class CBLQSetOperations2 : public CBLQSetOperations<ndim> {
public:
	CBLQSetOperations2(CBLQSetOperationsConfig conf) : CBLQSetOperations<ndim>(conf) {}
	virtual ~CBLQSetOperations2() {}

private:
    virtual boost::shared_ptr< CBLQRegionEncoding<ndim> > unary_set_op_impl(boost::shared_ptr< const CBLQRegionEncoding<ndim> > region, UnarySetOperation op) const;
    virtual boost::shared_ptr< CBLQRegionEncoding<ndim> > binary_set_op_impl(boost::shared_ptr< const CBLQRegionEncoding<ndim> > left, boost::shared_ptr< const CBLQRegionEncoding<ndim> > right, NArySetOperation op) const;

protected:
    using typename CBLQSetOperations<ndim>::RegionEncodingCPtrCIter;
    using cblq_word_t = typename CBLQRegionEncoding<ndim>::cblq_word_t;
};

template<int ndim>
class CBLQSetOperationsFast : public CBLQSetOperations<ndim> {
public:
	CBLQSetOperationsFast(CBLQSetOperationsConfig conf) : CBLQSetOperations<ndim>(conf) {}
	virtual ~CBLQSetOperationsFast() {}

private:
    virtual boost::shared_ptr< CBLQRegionEncoding<ndim> > unary_set_op_impl(boost::shared_ptr< const CBLQRegionEncoding<ndim> > region, UnarySetOperation op) const;
    virtual boost::shared_ptr< CBLQRegionEncoding<ndim> > binary_set_op_impl(boost::shared_ptr< const CBLQRegionEncoding<ndim> > left, boost::shared_ptr< const CBLQRegionEncoding<ndim> > right, NArySetOperation op) const;

    //virtual boost::shared_ptr< const CBLQRegionEncoding<ndim> > nary_set_op(RegionEncodingCPtrCIter region_it, RegionEncodingCPtrCIter region_end_it, NArySetOperation op) const;

protected:
    using typename CBLQSetOperations<ndim>::RegionEncodingCPtrCIter;
    using cblq_word_t = typename CBLQRegionEncoding<ndim>::cblq_word_t;
};

// Code-by-code nary set ops
template<int ndim>
class CBLQSetOperationsNAry1 : public CBLQSetOperations<ndim> {
public:
	CBLQSetOperationsNAry1(CBLQSetOperationsConfig conf) : CBLQSetOperations<ndim>(conf) {}
	virtual ~CBLQSetOperationsNAry1() {}

private:
	using typename CBLQSetOperations<ndim>::RegionEncodingCPtrCIter;
    virtual boost::shared_ptr< CBLQRegionEncoding<ndim> > nary_set_op_impl(RegionEncodingCPtrCIter region_it, RegionEncodingCPtrCIter region_end_it, NArySetOperation op) const;

private:
    typedef typename CBLQSetOperations<ndim>::cblq_word_t cblq_word_t;
};

// Word-by-word nary set ops with dense suffix support
template<int ndim>
class CBLQSetOperationsNAry2Dense : public CBLQSetOperations<ndim> {
public:
	CBLQSetOperationsNAry2Dense(CBLQSetOperationsConfig conf) : CBLQSetOperations<ndim>(conf) {}
	virtual ~CBLQSetOperationsNAry2Dense() {}

private:
	using typename CBLQSetOperations<ndim>::RegionEncodingCPtrCIter;
    virtual boost::shared_ptr< CBLQRegionEncoding<ndim> > nary_set_op_impl(RegionEncodingCPtrCIter region_it, RegionEncodingCPtrCIter region_end_it, NArySetOperation op) const;

private:
    typedef typename CBLQSetOperations<ndim>::cblq_word_t cblq_word_t;
    typedef typename CBLQRegionEncoding<ndim>::cblq_semiword_t cblq_semiword_t;
};

// Level-by-level nary set ops with dense suffix support
template<int ndim>
class CBLQSetOperationsNAry3Dense : public CBLQSetOperations<ndim> {
public:
	CBLQSetOperationsNAry3Dense(CBLQSetOperationsConfig conf) : CBLQSetOperations<ndim>(conf) {}
	virtual ~CBLQSetOperationsNAry3Dense() {}

private:
	using typename CBLQSetOperations<ndim>::RegionEncodingCPtrCIter;
    virtual boost::shared_ptr< CBLQRegionEncoding<ndim> > nary_set_op_impl(RegionEncodingCPtrCIter region_it, RegionEncodingCPtrCIter region_end_it, NArySetOperation op) const;

private:
    typedef typename CBLQSetOperations<ndim>::cblq_word_t cblq_word_t;
    typedef typename CBLQRegionEncoding<ndim>::cblq_semiword_t cblq_semiword_t;
};

// Level-by-level nary set ops with dense suffix support, optimized for union
template<int ndim>
class CBLQSetOperationsNAry3Fast : public CBLQSetOperationsFast<ndim> {
private:
    typedef typename CBLQRegionEncoding<ndim>::cblq_word_t cblq_word_t;
    typedef typename CBLQRegionEncoding<ndim>::cblq_semiword_t cblq_semiword_t;
    typedef typename CBLQSemiwords<ndim>::block_t semiword_block_t;

public:
	CBLQSetOperationsNAry3Fast(CBLQSetOperationsConfig conf) : CBLQSetOperationsFast<ndim>(conf) {}
	virtual ~CBLQSetOperationsNAry3Fast() {}

private:
	//virtual boost::shared_ptr< CBLQRegionEncoding<ndim> > binary_set_op(boost::shared_ptr< const CBLQRegionEncoding<ndim> > left, boost::shared_ptr< const CBLQRegionEncoding<ndim> > right, NArySetOperation op) const;

	using typename CBLQSetOperations<ndim>::RegionEncodingCPtrCIter;
    virtual boost::shared_ptr< CBLQRegionEncoding<ndim> > nary_set_op_impl(RegionEncodingCPtrCIter region_it, RegionEncodingCPtrCIter region_end_it, NArySetOperation op) const;

    template<NArySetOperation op>
    inline boost::shared_ptr< CBLQRegionEncoding<ndim> > nary_set_op_impl(RegionEncodingCPtrCIter cblq_it, RegionEncodingCPtrCIter cblq_end_it) const;
};

#endif /* CBLQ_SETOPS_HPP_ */

ext: hpp

max_stars_repo_name: daboyuka/PIQUE

max_issues_repo_name: daboyuka/PIQUE

max_forks_repo_name: daboyuka/PIQUE


text:
/*
 * Copyright 2009-2011 The VOTCA Development Team (http://www.votca.org)
 *
 * 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.
 *
 */

#include <boost/algorithm/string/trim.hpp>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <votca/tools/rangeparser.h>
#include <votca/tools/table.h>
#include <boost/lexical_cast.hpp>
#include <votca/tools/tokenizer.h>
#include <votca/tools/getline.h>
#include <iostream>
#include <votca/csg/imcio.h>

namespace votca { namespace csg {

typedef ub::matrix<double> group_matrix;
using namespace std;

void imcio_write_dS(const string &file, ub::vector<double> &r, ub::vector<double> &dS, std::list<int> *list)
{
    // write the dS
    ofstream out_dS;
    out_dS.open(file.c_str());
    out_dS << setprecision(8);
    if (!out_dS)
        throw runtime_error(string("error, cannot open file ") + file);

    if(list == NULL) {
        for (size_t i = 0; i < dS.size(); ++i) {
            out_dS << r[i] << " " << dS[i] << endl;
        }
    }
    else {
        for(std::list<int>::iterator i = list->begin(); i!=list->end(); ++i) {
            out_dS << r[*i] << " " << dS[*i] << endl;
        }
    }

    out_dS.close();
    cout << "written " << file << endl;
}

void imcio_write_matrix(const string &file, ub::matrix<double> &gmc, std::list<int> *list)
{
    ofstream out_A;
    out_A.open(file.c_str());
    out_A << setprecision(8);

    if (!out_A)
        throw runtime_error(string("error, cannot open file ") + file);

    if(list == NULL) {
        for (group_matrix::size_type i = 0; i < gmc.size1(); ++i) {
            for (group_matrix::size_type j = 0; j < gmc.size2(); ++j) {
                out_A << gmc(min(i,j), max(i,j)) << " ";
            }
            out_A << endl;
        }
    }
    else {
        for(std::list<int>::iterator i = list->begin(); i!=list->end(); ++i) {
            for(std::list<int>::iterator j = list->begin(); j!=list->end(); ++j) {
                out_A << gmc(*i, *j) << " ";
            }
            out_A << endl;
        }
    }
out_A.close();
    cout << "written " << file << endl;
}

void imcio_write_index(const string &file, vector<string> &names, vector<RangeParser> &ranges)
{
    // write the index

    ofstream out_idx;
    out_idx.open(file.c_str());

    if (!out_idx)
        throw runtime_error(string("error, cannot open file ") + file);

    for (size_t i = 0; i < names.size(); ++i)
        out_idx << names[i] << " " << ranges[i] << endl;

    out_idx.close();
    cout << "written " << file << endl;
}

void imcio_read_dS(const string &filename, ub::vector<double> &r, ub::vector<double> &dS)
{
    Table tbl;
    tbl.Load(filename);

    r.resize(tbl.size());
    dS.resize(tbl.size());

    for(int i=0; i<tbl.size(); ++i) {
        r(i) = tbl.x(i);
        dS(i) = tbl.y(i);
    }
}

void imcio_read_matrix(const string &filename, ub::matrix<double> &gmc)
{
    ifstream in;
    in.open(filename.c_str());

    bool is_initialized = false;
    if(!in)
        throw runtime_error(string("error, cannot open file ") + filename);

    size_t line_count =0;
    string line;
    // read till the first data line
    while(getline(in, line)) {
        // remove comments and xmgrace stuff
        line = line.substr(0, line.find("#"));
        line = line.substr(0, line.find("@"));

        // tokenize string and put it to vector
        Tokenizer tok(line, " \t");
        vector<string> tokens;
        tok.ToVector(tokens);

        // skip empty lines
        if(tokens.size()==0) continue;
        if(!is_initialized)
            gmc.resize(tokens.size(),tokens.size());
        is_initialized=true;

        if(gmc.size1()!=tokens.size())
            throw runtime_error(string("error loading ")
                    + filename + ": size mismatchm, number of columns differ");
        for(size_t i=0; i<tokens.size(); ++i)
            gmc(line_count,i) = boost::lexical_cast<double>(tokens[i]);
        ++line_count;
    }
    if(line_count != gmc.size1())
            throw runtime_error(string("error loading ")
                    + filename + ": size mismatch, not enough lines");
    in.close();
}

void imcio_read_index(const string &filename, vector<string> &names, vector<RangeParser> &ranges)
{
    ifstream in;
    in.open(filename.c_str());
    if(!in)
        throw runtime_error(string("error, cannot open file ") + filename);

    names.clear();
    ranges.clear();
    string line;

    // read till the first data line
    while(getline(in, line)) {
        // remove comments and xmgrace stuff
        line = line.substr(0, line.find("#"));
        line = line.substr(0, line.find("@"));

        boost::trim(line);

        size_t found;
        found = line.find(" ");
        if(found==string::npos)
            throw runtime_error(string("wrong format in ") + filename);

        string name = line.substr(0, found);

        string range = line.substr(found);

        RangeParser rp;
        rp.Parse(range);
        names.push_back(name);
        ranges.push_back(rp);

    }
    in.close();
}

}}

ext: cc


text:
//
//  ShaderReloader.hpp
//  OpenGLTest
//
//  Created by Samuel Hall on 27/07/2016.
//
//

#ifndef ShaderReloader_hpp
#define ShaderReloader_hpp

#include <stdio.h>

#include <vector>
#include <thread>
#include <ctime>
#include <boost/filesystem.hpp>

class ShaderReloader
{
private:

    std::thread m_monitoring_thread;
    bool m_is_running;
    std::condition_variable m_thread_conditional;
    std::mutex m_thread_mutex;
    std::timed_mutex m_deinit_mutex;

    struct File {
        boost::filesystem::path path;
        std::time_t last_write_time;
        std::function<void()> callback;
    };

    std::vector<std::shared_ptr<File>> m_watched_files;

    void initialize();
    void deinitialize();

    void add_files_to_watch(std::function<void ()> callback);

public:

    ShaderReloader() {

        initialize();
    };

    ~ShaderReloader() {

        deinitialize();
    };

    template<typename ... Args>
	    void add_files_to_watch(std::function<void ()> callback, std::string file_path, Args ... args);
};

template<typename ... Args>
inline void ShaderReloader::add_files_to_watch(std::function<void ()> callback, std::string file_path, Args ... args)
{
    std::shared_ptr<File> file(new File());

    file->path = boost::filesystem::path(file_path);
    file->last_write_time = boost::filesystem::last_write_time(file->path);
    file->callback = callback;

    m_watched_files.push_back(file);

    add_files_to_watch(callback, args...);
}

#endif /* ShaderReloader_hpp */

max_stars_repo_name: BrassLion/opengl-opencl-framework

max_issues_repo_name: BrassLion/opengl-opencl-framework

max_forks_repo_name: BrassLion/opengl-opencl-framework


text:
// Project: libv.parse, File: src/libv/parse/bool.cpp, Author: Császár Mátyás [Vader]

// hpp
#include <libv/parse/bool.hpp>
// ext
#include <boost/spirit/home/x3.hpp>
#include <boost/spirit/home/x3/char/unicode.hpp>
// libv
#include <libv/utility/concat.hpp>
// std
#include <algorithm>
#include <stdexcept>

namespace libv {
namespace parse {

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

bool parse_bool_or(const std::string_view str, const bool fallback) {
	const auto result = parse_bool_optional(str);
	return result ? *result : fallback;
}

bool parse_bool_or_throw(const std::string_view str) {
	const auto result = parse_bool_optional(str);
	if (result)
		return *result;
	else
		throw std::invalid_argument(libv::concat("\"", str, "\" is not a valid boolean"));
}

std::optional<bool> parse_bool_optional(const std::string_view str) {
	namespace x3 = boost::spirit::x3;

	if (std::all_of(str.begin(), str.end(), static_cast<int(*)(int)>(std::isspace)))
		return false;

	const auto value = [](bool value) {
		return [value](auto& ctx) {
			_val(ctx) = value;
		};
	};

	const auto value_from_int = [](auto& ctx) {
		if (_attr(ctx) == 0)
			_val(ctx) = false;
		else if (_attr(ctx) == 1)
			_val(ctx) = true;
		else if (_attr(ctx) == -1)
			_val(ctx) = true;
		else
			_pass(ctx) = false;
	};

	const auto bool_rule = x3::rule<class bool_rule_, bool>{} =
			x3::no_case[x3::lit("true")][value(true)] |
			x3::no_case[x3::lit("false")][value(false)] |
			x3::int_[value_from_int];

	bool result = false;

	auto it = str.begin();
	// NOTE: x3::unicode::space is required for certain invalid UTF8 inputs as x3::space skipper would assert
	auto success = x3::phrase_parse(it, str.end(), bool_rule, x3::unicode::space, result);
	success = success && it == str.end();

	return success ? std::optional<bool>{result} : std::nullopt;
}

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

} // namespace parse
} // namespace libv

max_stars_repo_name: cpplibv/libv

max_issues_repo_name: cpplibv/libv

max_forks_repo_name: cpplibv/libv


text:
#ifndef SUPERGENIUS_API_REQUEST_UNSUBSCRIBE_STORAGE
#define SUPERGENIUS_API_REQUEST_UNSUBSCRIBE_STORAGE

#include <jsonrpc-lean/request.h>

#include <boost/optional.hpp>

#include "api/service/state/state_api.hpp"
#include "base/buffer.hpp"
#include "outcome/outcome.hpp"
#include "primitives/block_id.hpp"

namespace sgns::api::state::request {

  class UnsubscribeStorage final {
   public:
    UnsubscribeStorage(const UnsubscribeStorage &) = delete;
    UnsubscribeStorage &operator=(const UnsubscribeStorage &) = delete;

    UnsubscribeStorage(UnsubscribeStorage &&) = default;
    UnsubscribeStorage &operator=(UnsubscribeStorage &&) = default;

    explicit UnsubscribeStorage(std::shared_ptr<StateApi> api)
        : api_(std::move(api)) { };
    ~UnsubscribeStorage() = default;

    outcome::result<void> init(const jsonrpc::Request::Parameters &params);
    outcome::result<void> execute();

   private:
    std::shared_ptr<StateApi> api_;
    std::vector<uint32_t> subscriber_id_;
  };

}  // namespace sgns::api::state::request

#endif  // SUPERGENIUS_API_REQUEST_UNSUBSCRIBE_STORAGE

max_stars_repo_name: GeniusVentures/SuperGenius

max_issues_repo_name: GeniusVentures/SuperGenius

max_forks_repo_name: GeniusVentures/SuperGenius


text:
//==================================================================================================
/**
  Copyright 2016 NumScale SAS

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

#if defined(DOXYGEN_ONLY)
namespace boost { namespace simd
{
  /*!
    @ingroup group-swar
    Function object implementing split_low capabilities

    SIMD register type-based split_low

    @par Header <boost/simd/function/split_low.hpp>

    @c split_low extract the lower half of a SIMD register and convert it
    to the appropriate SIMD register type of corresponding cardinal.

    @see split_high, split, slice
  **/
  upgrade_t<Value> split_low(Value const& x);
} }
#endif

#include <boost/simd/function/simd/split_low.hpp>

#endif

max_stars_repo_name: SylvainCorlay/pythran

max_issues_repo_name: SylvainCorlay/pythran

max_forks_repo_name: SylvainCorlay/pythran


text:
#include <memo/model/faith/Faith.hh>

#include <boost/uuid/random_generator.hpp>

#include <elle/log.hh>
#include <elle/serialization/binary.hh>

#include <elle/cryptography/hash.hh>

#include <memo/model/Address.hh>
#include <memo/model/MissingBlock.hh>
#include <memo/model/blocks/MutableBlock.hh>
#include <memo/silo/MissingKey.hh>

ELLE_LOG_COMPONENT("memo.model.faith.Faith");

using Serializer = elle::serialization::Binary;

namespace memo
{
  namespace model
  {
    namespace faith
    {
      Faith::Faith(std::unique_ptr<silo::Silo> storage,
                   boost::optional<elle::Version> version)
        : Model(std::move(version))
        , _storage(std::move(storage))
      {}

      std::unique_ptr<blocks::Block>
      Faith::_fetch(Address address, boost::optional<int>) const
      {
        ELLE_TRACE_SCOPE("%s: fetch block at %x", *this, address);
        try
        {
          auto data = this->_storage->get(address);
          return elle::serialization::deserialize<
            Serializer, std::unique_ptr<blocks::Block>>(data, false);
        }
        catch (memo::silo::MissingKey const&)
        {
          throw MissingBlock(address);
        }
      }

      void
      Faith::_insert(std::unique_ptr<blocks::Block> block,
                     std::unique_ptr<ConflictResolver> resolver)
      {
        ELLE_TRACE_SCOPE("%s: insert %f", *this, *block);
        this->_storage->set(
          block->address(),
          elle::serialization::serialize<Serializer>(block, false),
          true,
          false);
      }

      void
      Faith::_update(std::unique_ptr<blocks::Block> block,
                     std::unique_ptr<ConflictResolver> resolver)
      {
        ELLE_TRACE_SCOPE("%s: update %f", *this, *block);
        this->_storage->set(
          block->address(),
          elle::serialization::serialize<Serializer>(block, false),
          false,
          true);
      }

      void
      Faith::_remove(Address address, blocks::RemoveSignature)
      {
        ELLE_TRACE_SCOPE("%s: remove block at %x", *this, address);
        try
        {
          this->_storage->erase(address);
        }
        catch (memo::silo::MissingKey const&)
        {
          throw MissingBlock(address);
        }
      }

      struct FaithModelConfig:
        public ModelConfig
      {
      public:
        FaithModelConfig(elle::serialization::SerializerIn& input,
                         elle::Version version)
          : ModelConfig(nullptr, std::move(version))
        {
          this->serialize(input);
        }

        void
        serialize(elle::serialization::Serializer& s) override
        {
          ModelConfig::serialize(s);
        }

        virtual
        std::unique_ptr<memo::model::Model>
        make(bool,
             boost::filesystem::path const&) override
        {
          return std::make_unique<memo::model::faith::Faith>(
            this->storage->make(), this->version);
        }
      };

      static const elle::serialization::Hierarchy<ModelConfig>::
      Register<FaithModelConfig> _register_FaithModelConfig("faith");
    }
  }
}

max_stars_repo_name: infinit/memo

max_issues_repo_name: infinit/memo

max_forks_repo_name: infinit/memo


text:
#if HAVE_CONFIG_H
# include <config.h>
#endif

#include <boost/bind.hpp>
using boost::cref;

#include "Cosmology.h"
#include "PowerSpectrum.h"
#include "Quadrature.h"
#include "BSPT.h"
#include "BSPTN.h"
#include "SPT.h"
#include "RegPT.h"

#include "Spline.h"
#include "SpecialFunctions.h"
#include "LinearPS.h"
#include "CMBc.h"

#include <cerrno>
#include <cstdio>
#include <cstdlib>
#include <sys/stat.h>
#include <math.h>       /* pow */
#include <boost/math/tools/roots.hpp>
#include <random>

#include <gsl/gsl_math.h>
#include <gsl/gsl_deriv.h>

// bottom limit of comoving distance integral
const double xim = 10.;

CMBc::CMBc(const Cosmology& C_, const PowerSpectrum& P_L_, real epsrel_)
: C(C_), P_L(P_L_)
{
    epsrel = epsrel_;
}

//////////////////////////////////////////
 /* Functions used in initialisation */
//////////////////////////////////////////

// Hubble function
double HAc(double omega0, double z){
	double omegaL= 1.-omega0;
	return  1./sqrt(omega0*pow(1.+z,3)+omegaL);}

// Functions for GM or SC formula
// neff(k)
  static double myfuncnw(double k,void * params){
    (void)(params);
    return log(neffnw(k));
  }

  static double myfuncehu(double k,void * params){
    (void)(params);
    return log(neffehu(k));
  }

  static double myfunclin(double k,void * params){
    (void)(params);
    return log(mylinearps(k));
  }

static double myneff(double k){
  gsl_function F;
  double result, abserr;
    F.params = 0;
    F.function = &myfuncnw;
    gsl_deriv_central(&F, k, 1e-6, &result, &abserr);
return k*result;
}

// Find K_nl
  double myknls(const PowerSpectrum& P_L)
  {
      std::mt19937 gen(2);
       const double lower_bound =  1e-4;
       const double upper_bound =  1000.;
      std::uniform_real_distribution<> dis(lower_bound, upper_bound);

      double pos_pt = dis(gen);
      double neg_pt = dis(gen);

      while ( (pow2(D_spt/dnorm_spt)*P_L(pos_pt)*pow3(pos_pt)/2./M_PI/M_PI  - 1.) < 0.){
          pos_pt = dis(gen);}

      while ( (pow2(D_spt/dnorm_spt)*P_L(neg_pt)*pow3(neg_pt)/2./M_PI/M_PI- 1.) > 0.){
          neg_pt = dis(gen);}

       const double about_zero_mag = 1e-5;
      for (;;)
      {
          double mid_pt = (pos_pt + neg_pt)/2.0;
          double f_mid_pt =  pow2(D_spt/dnorm_spt)*P_L(mid_pt)*pow3(mid_pt)/2./M_PI/M_PI -1.;

          if (fabs(f_mid_pt)  < about_zero_mag){
              return mid_pt;
                }

          if (f_mid_pt >= 0.){
              pos_pt = mid_pt;}
          else{
              neg_pt = mid_pt;}

      }
  }

// regpt damp term for post born prediction
static real regpt_exp(const PowerSpectrum& P_l, double q){
    return  P_l(q)/(6.*pow2(M_PI));
  }

// Calculate sigma_8 at some redshift
static double sigma8int(const PowerSpectrum& P_L, double k){
  return pow2(3.*((sin(k*8.) - k*8.*cos(k*8.))/pow2(k*8.))/k/8.) * pow2(k) * pow2(D_spt/dnorm_spt) * P_L(k)/2./M_PI/M_PI;
}

static double mysigma8(const PowerSpectrum& P_L){
  return sqrt(Integrate(bind(sigma8int,cref(P_L),_1),0.0001,30.,1e-4));
}

// Spline function z->comoving distance
Spline zofc;
// Spline linear growth factor + F2_dgp
Spline Dlin, Dlingr, F2DGP;
// Splines for HALOFIT
Spline phpars0,phpars1,phpars2,phpars3,phpars4,phpars5,phpars6,phpars7,phpars8,phpars9,phpars10,phpars11,phpars12;
// Splines for dgp 1-loop matter power Spectrum
Spline fdgps, cdgps, idgps, kldgps, p22spl,p13spl, damp_termcmb;
// Comoving distance to z=1000 (~CMB)
double xis;
// comoving distance to z=40 as in TN's computation
double xis40;

// Splines for GM formula
Spline myneffspl, mys8spl, myknlspl;

// Initialize z(comoving_dist) and xi_star = comoving distance to CMB
// a = 1 linear
// a = 2 1-loop (+ initialises regpt components)
// a = 3 Gil-Marin or Scoccimaro [+ initialises all halofit params, n_eff(k), k_nl(z) and sigma_8(z)]
 void CMBc::comdist_init(int a, double omega0, double mg1){
  SPT spt(C, P_L, 1e-3);
  RegPT rpt(C, P_L, 1e-3);
  BSPT bspt(C, P_L, 1e-3);
  IOW iow;

  // initialise damping factor for regpt
  if(a==2){
  rpt.sigmad_init();
    }

   // initialise EHU PS and Fonseca NW
    if(a==3){
      bspt.mypspline(2);
    }

  double h0inv = 2997.92;
  const double zmin = 1e-4;
  const double zmax = 1100.;
  const double kmin = 1e-4;
  const double kmax = 100.;

  vector<double> ztable, ktable, ctable, dtable, dgrtable, f2table;
  // Halofit tables
  vector<double> phtable0,phtable1,phtable2,phtable3,phtable4,phtable5,phtable6,phtable7,phtable8,phtable9,phtable10,phtable11,phtable12;
  // 1-LOOP SPT nDGP GROWTH TABLES
  vector<double> ftab,ctab,itab,kltab;
  // 1-Loop k-components
  vector<double> p13tab,p22tab, dampt;
  // GM formula tables
  vector<double> nefftab, knltab, sig8tab;

  xis =  h0inv*Integrate(bind(HAc,omega0,_1),0.,zmax,1e-4);
  xis40 = h0inv*Integrate(bind(HAc,omega0,_1),0.,40.,1e-4);

  const int n1 = 1000;
  double cval,zval,scalef,kval,sigd;

  for(int i = 0; i<n1; i++ ){
    zval = zmin*exp(i*log(zmax/zmin)/(n1-1.));
    kval = kmin*exp(i*log(kmax/kmin)/(n1-1.));

    scalef = 1./(1.+zval);
    iow.inite(scalef,omega0,mg1,1.,1.);

    cval = h0inv*Integrate(bind(HAc,omega0,_1),zmin,zval,1e-4);

    ztable.push_back(zval);
    ktable.push_back(kval);
    ctable.push_back(cval);

    dtable.push_back(D_spt/dnorm_spt);
    dgrtable.push_back(Dl_spt/dnorm_spt);
    f2table.push_back(F_spt/pow2(dnorm_spt));

    if (a==2) {
        p22tab.push_back(pow4(dnorm_spt/Dl_spt)*spt.P22_dd(kval));
        p13tab.push_back(pow4(dnorm_spt/Dl_spt)*spt.P13_dd(kval));
        sigd = pow2(kval)*Integrate<ExpSub>(bind(regpt_exp, cref(P_L), _1), kmin, kval/2., 1e-4);
        dampt.push_back(sigd);
    }

    if (a ==3 ) {
      spt.phinit(scalef,omega0);
      phtable0.push_back(phpars[0]);
      phtable1.push_back(phpars[1]);
      phtable2.push_back(phpars[2]);
      phtable3.push_back(phpars[3]);
      phtable4.push_back(phpars[4]);
      phtable5.push_back(phpars[5]);
      phtable6.push_back(phpars[6]);
      phtable7.push_back(phpars[7]);
      phtable8.push_back(phpars[8]);
      phtable9.push_back(phpars[9]);
      phtable10.push_back(phpars[10]);
      phtable11.push_back(phpars[11]);
      phtable12.push_back(phpars[12]);

        if (i==0) {
          knltab.push_back(myknls(cref(P_L)));
          }
        else if( knltab[i-1] > 200. ){
           knltab.push_back(knltab[i-1]);
        }
        else{
           knltab.push_back(myknls(cref(P_L)));
              }
             printf("%d %e \n", i, knltab[i]);

      sig8tab.push_back(mysigma8(cref(P_L)));
      nefftab.push_back(myneff(kval));

    }
      }

      zofc = LinearSpline(ctable,ztable);
      Dlin = LinearSpline(ztable,dtable);
      Dlingr = LinearSpline(ztable,dgrtable);
      F2DGP = LinearSpline(ztable,f2table);

      if (a==2) {
        p22spl=LinearSpline(ktable,p22tab);
        p13spl=LinearSpline(ktable,p13tab);
        damp_termcmb = LinearSpline(ktable,dampt);
            }
      if (a==3) {
        phpars0 = LinearSpline(ztable,phtable0);
        phpars1 = LinearSpline(ztable,phtable1);
        phpars2 = LinearSpline(ztable,phtable2);
        phpars3 = LinearSpline(ztable,phtable3);
        phpars4 = LinearSpline(ztable,phtable4);
        phpars5 = LinearSpline(ztable,phtable5);
        phpars6 = LinearSpline(ztable,phtable6);
        phpars7 = LinearSpline(ztable,phtable7);
        phpars8 = LinearSpline(ztable,phtable8);
        phpars9 = LinearSpline(ztable,phtable9);
        phpars10 = LinearSpline(ztable,phtable10);
        phpars11 = LinearSpline(ztable,phtable11);
        phpars12 = LinearSpline(ztable,phtable12);

        myknlspl = LinearSpline(ztable,knltab);
        mys8spl = LinearSpline(ztable,sig8tab);
        myneffspl = LinearSpline(ktable,nefftab);

      }
}

// B_LSS integrands
// tree analytic
static double btlssa_integrand(const Cosmology& C, const PowerSpectrum& P_L, int a, double params[], double l1, double l2, double x, double xi ){
   BSPT bspt(C, P_L, 1e-4);
   IOW iow;
   double bis, kernel, F2A, F2B, F2C, growth;
   double h0 = 1./2997.92;
   double lensingk = (xis-xi)/xi/xis;
   double myz = zofc(xi);
   double scalef = 1./(1.+myz);
   double prefactor = pow3(3.*params[0]*pow2(h0)/(2.*scalef) * lensingk);
   double k1 = l1/xi;
   double k2 = l2/xi;

   D_spt = Dlin(myz) * dnorm_spt;
   Dl_spt = Dlin(myz) * dnorm_spt; //Dlingr(myz) * dnorm_spt;
   F_spt = F2DGP(myz) * dnorm_spt;
   bis = bspt.Btree(a, k1, k2, x);

   return prefactor * bis * pow2(xi);
}

// 1-loop analytic
static double bllssa_integrand(const Cosmology& C,const PowerSpectrum& P_L, int a, double params[], double l1, double l2, double x, double xi ){
   BSPT bspt(C, P_L, 1e-2);
   IOW iow;
   double bis, kernel;
   double h0 = 1./2997.92;
   double lensingk = (xis-xi)/xi/xis;
   double myz = zofc(xi);
   double scalef = 1./(1.+myz);
   double prefactor = pow3(3.*params[0]*pow2(h0)/(2.*scalef) * lensingk);
   switch (a) {
     case 1:
      Dl_spt = Dlin(myz) * dnorm_spt;
      break;
      case 2:
   iow.inite2(scalef,params[0],params[1],params[2],params[3]);
      break;
 }
   double k1 = l1/xi;
   double k2 = l2/xi;
   bis = bspt.Bloop(a, k1, k2, x);
   return prefactor * bis * pow2(xi);
}

// Gil-Marin formula
static double bgmlssa_integrand(const Cosmology& C,const PowerSpectrum& P_L,  double params[], double l1, double l2, double x, double xi ){
   SPT spt(C, P_L, 1e-3);
   double vars[6],F2A,F2B,F2C,p1,p2,p3;

   double h0 = 1./2997.92;
   double lensingk = (xis-xi)/xi/xis;
   double myz = zofc(xi);
   double scalef = 1./(1.+myz);
   double prefactor = pow3(3.*params[0]*pow2(h0)/(2.*scalef) * lensingk);

   double k1 = l1/xi;
   double k2 = l2/xi;
   double k3 = sqrt(pow2(k1) + pow2(k2) + 2.*k2*k1*x);

    D_spt =  Dlin(myz) * dnorm_spt;
    F_spt =  F2DGP(myz) * pow2(dnorm_spt);
    vars[0] = mys8spl(myz);
    vars[1] = myknlspl(myz);
    vars[2] = 1.;
    vars[3] = (1.- F_spt * 7./2./pow2(D_spt));

     vars[4] = myneffspl(k1);
     vars[5] = myneffspl(k2);
     F2A = F2fit(vars,k1,k2,x);

     vars[4] = myneffspl(k1);
     vars[5] = myneffspl(k3);
     F2B = F2fit(vars,k1,k3,-(k2*x+k1)/k3);

     vars[4] = myneffspl(k2);
     vars[5] = myneffspl(k3);
     F2C = F2fit(vars,k2,k3,-(k1*x+k2)/k3);

    phpars[0] = phpars0(myz);
    phpars[1] = phpars1(myz);
    phpars[2] = phpars2(myz);
    phpars[3] = phpars3(myz);
    phpars[4] = phpars4(myz);
    phpars[5] = phpars5(myz);
    phpars[6] = phpars6(myz);
    phpars[7] = phpars7(myz);
    phpars[8] = phpars8(myz);
    phpars[9] = phpars9(myz);
    phpars[10] = phpars10(myz);
    phpars[11] = phpars11(myz);
    phpars[12] = phpars12(myz);

     if (k1>100.) {
      p1 = spt.PHALO(100.);
     }
     else{
     p1 = spt.PHALO(k1);
        }
        if (k2>100.) {
          p2 = spt.PHALO(100.);
        }
        else{
        p2 = spt.PHALO(k2);
           }

           if (k3>100.) {
              p3 = spt.PHALO(100.);
           }
           else{
           p3 = spt.PHALO(k3);
              }

    double bis =2.*(p1*p2*F2A + p3*p1*F2B + p2*p3*F2C);

    return prefactor * bis * pow2(xi);
}

// SC fitting formula
static double bsclssa_integrand(const Cosmology& C,const PowerSpectrum& P_L,  double params[], double l1, double l2, double x, double xi ){
   SPT spt(C, P_L, 1e-3);
   double vars[6],F2A,F2B,F2C,p1,p2,p3;

   double h0 = 1./2997.92;
   double lensingk = (xis-xi)/xi/xis;
   double myz = zofc(xi);
   double scalef = 1./(1.+myz);
   double prefactor = pow3(3.*params[0]*pow2(h0)/(2.*scalef) * lensingk);

   double k1 = l1/xi;
   double k2 = l2/xi;
   double k3 = sqrt(pow2(k1) + pow2(k2) + 2.*k2*k1*x);

    D_spt =  Dlin(myz) * dnorm_spt;
    F_spt =  F2DGP(myz) * pow2(dnorm_spt);
    vars[0] = mys8spl(myz);
    vars[1] = myknlspl(myz);
    vars[2] = 1.;
    vars[3] = (1.- F_spt * 7./2./pow2(D_spt));

     vars[4] = myneffspl(k1);
     vars[5] = myneffspl(k2);
     F2A = F2fitsc(vars,k1,k2,x);

     vars[4] = myneffspl(k1);
     vars[5] = myneffspl(k3);
     F2B = F2fitsc(vars,k1,k3,-(k2*x+k1)/k3);

     vars[4] = myneffspl(k2);
     vars[5] = myneffspl(k3);
     F2C = F2fitsc(vars,k2,k3,-(k1*x+k2)/k3);

    phpars[0] = phpars0(myz);
    phpars[1] = phpars1(myz);
    phpars[2] = phpars2(myz);
    phpars[3] = phpars3(myz);
    phpars[4] = phpars4(myz);
    phpars[5] = phpars5(myz);
    phpars[6] = phpars6(myz);
    phpars[7] = phpars7(myz);
    phpars[8] = phpars8(myz);
    phpars[9] = phpars9(myz);
    phpars[10] = phpars10(myz);
    phpars[11] = phpars11(myz);
    phpars[12] = phpars12(myz);

     if (k1>100.) {
      p1 = spt.PHALO(100.);
     }
     else{
     p1 = spt.PHALO(k1);
        }
        if (k2>100.) {
          p2 = spt.PHALO(100.);
        }
        else{
        p2 = spt.PHALO(k2);
           }

           if (k3>100.) {
              p3 = spt.PHALO(100.);
           }
           else{
           p3 = spt.PHALO(k3);
              }

    double bis =2.*(p1*p2*F2A + p3*p1*F2B + p2*p3*F2C);

    return prefactor * bis * pow2(xi);
}

// Numerical tree level
static double btnlssa_integrand(const Cosmology& C, const PowerSpectrum& P_L, int a, double params[], double l1, double l2, double x, double xi ){
   BSPT bspt(C, P_L, 1e-3);
   IOW iow;
   double bis, kernel, F2A, F2B, F2C, growth;
   double h0 = 1./2997.92;
   double lensingk = (xis-xi)/xi/xis;
   double myz = zofc(xi);
   double scalef = 1./(1.+myz);
   double prefactor = pow3(3.*params[0]*pow2(h0)/(2.*scalef) * lensingk);
   double k1 = l1/xi;
   double k2 = l2/xi;
   double k3 = sqrt(pow2(k1) + pow2(k2) + 2.*k2*k1*x);

   double vars[6];
   vars[0] = params[0];
   vars[1] = params[1];
   vars[2] = params[2];
   vars[3] = params[3];
   vars[4] = scalef;
   vars[5] = 1.*a; // doubling as choice of dgp (1), ide (2) , fr(3)

   bis = bspt.Btreen(vars, k1, k2, x);
   return prefactor * bis * pow2(xi);
}

// numerical 1-loop level
static double blnlssa_integrand(const Cosmology& C, const PowerSpectrum& P_L, int a, double params[], double l1, double l2, double x, double xi ){
   BSPT bspt(C, P_L, 1e-4);
   IOW iow;
   double bis, kernel, F2A, F2B, F2C, growth;
   double h0 = 1./2997.92;
   double lensingk = (xis-xi)/xi/xis;
   double myz = zofc(xi);
   double scalef = 1./(1.+myz);
   double prefactor = pow3(3.*params[0]*pow2(h0)/(2.*scalef) * lensingk);
   double k1 = l1/xi;
   double k2 = l2/xi;
   double vars[9];
   vars[0] = params[0];
   vars[1] = params[1];
   vars[2] = params[2];
   vars[3] = params[3];
   vars[4] = scalef;
   vars[5] = 1.*a; // doubling as choice of dgp (1), ide (2) , fr(3)
   vars[6] = params[4];
   vars[7] = params[5];
   vars[8] = params[6];
   bis = bspt.Bloopn(vars, k1, k2, x);
   return prefactor * bis * pow2(xi);
}

// POST BORN STUFF

// C_l(xi1, xi2) integrands
// linear
static double cl_integrandt(const Cosmology& C, const PowerSpectrum& P_L, double omega0, double l2, double xi1, double xi2){
  double h0 = 1./2997.92;
  double k2 = l2/xi2;
  double lenker1 = (xi1 - xi2)/xi2/xi1 ;
  double lenker2 = (xis - xi2)/xi2/xis;
  double weylpot2 =  pow2(Dlin(zofc(xi2))) * P_L(k2) * pow2(3. * omega0 * pow2(h0) * (1.+zofc(xi2)) / 2. / pow2(k2));

  return  lenker1 * lenker2 / pow2(xi2) * weylpot2 ;
}

// 1-loop (regpt)
static double cl_integrandl(const Cosmology& C, const PowerSpectrum& P_L, double omega0, double l2, double xi1, double xi2){
  double D0, D0sqr, D2, expon, p13, pl,  ploop,ploopr,h0,k2,weylpot2,lenker1,lenker2;
   D0 =  pow2(Dlin(zofc(xi2)));
   h0 = 1./2997.92;
   k2 = l2/xi2;

  expon = D0 * damp_termcmb(k2) ;
  if(expon >= 80.){
   return 0.;
  }
  else{

   lenker1 = (xi1 - xi2)/xi2/xi1;
   lenker2 = (xis - xi2)/xi2/xis;

  D0sqr = pow2(D0);
  pl = D0 * P_L(k2);
  p13 =  D0sqr*p13spl(k2);
  ploop = pl+ D0sqr*p22spl(k2) + p13;
  D2 = pow2(p13/pl/2.)*D0;

  ploopr = exp(-expon)*(ploop + expon/2.*(2.*pl + p13) + pow2(expon)/4.*pl + P_L(k2)*D2);

  weylpot2 =   ploopr * pow2(3. * omega0 * pow2(h0) * (1.+zofc(xi2)) / 2. / pow2(k2));

  return  lenker1 * lenker2 / pow2(xi2) * weylpot2 ;
}
}

// 1-loop (spt)
static double cl_integrandspt(const Cosmology& C, const PowerSpectrum& P_L, double omega0, double l2, double xi1, double xi2){
  double D0, D0sqr, D2, expon, p13, pl,  ploop,ploopr,h0,k2,weylpot2,lenker1,lenker2;
   D0 =  pow2(Dlin(zofc(xi2)));
   h0 = 1./2997.92;
   k2 = l2/xi2;
   lenker1 = (xi1 - xi2)/xi2/xi1;
   lenker2 = (xis - xi2)/xi2/xis;

  D0sqr = pow2(D0);
  pl = D0 * P_L(k2);
  p13 =  D0sqr*p13spl(k2);
  ploop = pl+ D0sqr*p22spl(k2) + p13;

  weylpot2 =   ploop * pow2(3. * omega0 * pow2(h0) * (1.+zofc(xi2)) / 2. / pow2(k2));

  return  lenker1 * lenker2 / pow2(xi2) * weylpot2 ;
}

// halofit
static double cl_integrandh(const Cosmology& C, const PowerSpectrum& P_L, double omega0, double l2, double xi1, double xi2){
  SPT spt(C, P_L, 1e-3);

  double h0 = 1./2997.92;
  double k2 = l2/xi2;
  double lenker1 = (xi1 - xi2)/xi2/xi1 ;
  double lenker2 = (xis - xi2)/xi2/xis;
  double myz = zofc(xi2);

    D_spt =  Dlin(myz) * dnorm_spt;
    phpars[0] = phpars0(myz);
    phpars[1] = phpars1(myz);
    phpars[2] = phpars2(myz);
    phpars[3] = phpars3(myz);
    phpars[4] = phpars4(myz);
    phpars[5] = phpars5(myz);
    phpars[6] = phpars6(myz);
    phpars[7] = phpars7(myz);
    phpars[8] = phpars8(myz);
    phpars[9] = phpars9(myz);
    phpars[10] = phpars10(myz);
    phpars[11] = phpars11(myz);
    phpars[12] = phpars12(myz);

  double weylpot2 =   spt.PHALO(k2) * pow2(3. * omega0 * pow2(h0) * (1.+zofc(xi2)) / 2. / pow2(k2));

  return  lenker1 * lenker2 / pow2(xi2) * weylpot2 ;
}

// M(l1,l2) function integrands
// tree level
static double mfunc_integrandt(const Cosmology& C, const PowerSpectrum& P_L, double omega0, double l1, double l2, double xi1){
  double h0 = 1./2997.92;
  double mycl = pow4(l2)*Integrate(bind(cl_integrandt,cref(C),cref(P_L), omega0, l2, xi1, _1), xim,xi1,1e-4);  // C_l2 (xi1, xis)
  double k1 = l1/xi1;
  double weylpot1 =  pow2(Dlin(zofc(xi1))) * P_L(k1) * pow2(3. * omega0 * pow2(h0) * (1.+zofc(xi1)) / 2. / pow2(k1));

  return pow2((xis - xi1)/pow2(xi1)/xis) * weylpot1 * mycl;
}

// 1-loop level (RegPT implementation)
static double mfunc_integrandl(const Cosmology& C, const PowerSpectrum& P_L, double omega0, double l1, double l2, double xi1){
  double D0,D0sqr, D2, expon, p13, pl,  ploop,ploopr,h0,mycl,k1,weylpot1;
  h0 = 1./2997.92;
  D0 =  pow2(Dlin(zofc(xi1)));
  k1 = l1/xi1;
  expon = D0 * damp_termcmb(k1) ;
  if(expon >= 80. ){
   return 0.;
  }
  else{
     mycl = pow4(l2) * Integrate(bind(cl_integrandl, cref(C), cref(P_L), omega0, l2, xi1, _1), xim,xi1,1e-2);   // C_l2 (xi1, xis)
     D0sqr = pow2(D0);
     pl = D0 * P_L(k1);
     p13 =  D0sqr*p13spl(k1);
     ploop = pl+ D0sqr*p22spl(k1) + p13;
     D2 = pow2(p13/pl/2.)*D0;

     ploopr = exp(-expon)*(ploop + expon/2.*(2.*pl + p13) + pow2(expon)/4.*pl + P_L(k1)*D2);

     weylpot1 =  ploopr * pow2(3. * omega0 * pow2(h0) * (1.+zofc(xi1)) / 2. / pow2(k1));

  return pow2( (xis - xi1)/pow2(xi1)/xis ) * weylpot1 * mycl;
}
}

// pure spt
static double mfunc_integrandspt(const Cosmology& C, const PowerSpectrum& P_L, double omega0, double l1, double l2, double xi1){
  double D0,D0sqr, D2, expon, p13, pl,  ploop,ploopr,h0,mycl,k1,weylpot1;
  h0 = 1./2997.92;
  D0 =  pow2(Dlin(zofc(xi1)));
  k1 = l1/xi1;

     mycl = pow4(l2) * Integrate(bind(cl_integrandspt, cref(C), cref(P_L), omega0, l2, xi1, _1), xim,xi1,1e-2);   // C_l2 (xi1, xis)
     D0sqr = pow2(D0);
     pl = D0 * P_L(k1);
     p13 =  D0sqr*p13spl(k1);
     ploop = pl+ D0sqr*p22spl(k1) + p13;

     weylpot1 =  ploop * pow2(3. * omega0 * pow2(h0) * (1.+zofc(xi1)) / 2. / pow2(k1));

  return pow2( (xis - xi1)/pow2(xi1)/xis ) * weylpot1 * mycl;
}

// halofit
static double mfunc_integrandh(const Cosmology& C, const PowerSpectrum& P_L, double omega0, double l1, double l2, double xi1){
  SPT spt(C, P_L, 1e-2);

  double h0 = 1./2997.92;
  double mycl = pow4(l2)*Integrate(bind(cl_integrandh,cref(C), cref(P_L), omega0, l2, xi1, _1), xim,xi1,1e-4);   // C_l2 (xi1, xis)
  double k1 = l1/xi1;
  double myz = zofc(xi1);

    D_spt =  Dlin(myz) * dnorm_spt;
    phpars[0] = phpars0(myz);
    phpars[1] = phpars1(myz);
    phpars[2] = phpars2(myz);
    phpars[3] = phpars3(myz);
    phpars[4] = phpars4(myz);
    phpars[5] = phpars5(myz);
    phpars[6] = phpars6(myz);
    phpars[7] = phpars7(myz);
    phpars[8] = phpars8(myz);
    phpars[9] = phpars9(myz);
    phpars[10] = phpars10(myz);
    phpars[11] = phpars11(myz);
    phpars[12] = phpars12(myz);

  double weylpot1 =  spt.PHALO(k1) * pow2(3. * omega0 * pow2(h0) * (1.+zofc(xi1)) / 2. / pow2(k1));

  return pow2( (xis - xi1)/pow2(xi1)/xis ) * weylpot1 * mycl;
}

// M(l1,l2) as in Eq.4.5 of 1605.05662
static double Mfunc(int a, const Cosmology& C, const PowerSpectrum& P_L, double omega0, double l1, double l2){
    switch (a) {
      case 1:
      return pow4(l1) * Integrate(bind(mfunc_integrandt,cref(C),cref(P_L), omega0, l1, l2, _1), xim,xis40,1e-3); // tree level
      break;
      case 2:
      return pow4(l1) * Integrate(bind(mfunc_integrandl,cref(C),cref(P_L), omega0, l1, l2, _1), xim,xis40,1e-3); // 1-loop (regpt)
      break;
      case 3:
      return pow4(l1) * Integrate(bind(mfunc_integrandh,cref(C),cref(P_L), omega0, l1, l2, _1), xim,xis,1e-3); // halofit
      break;
      case 4:
      return pow4(l1) * Integrate(bind(mfunc_integrandspt,cref(C),cref(P_L), omega0, l1, l2, _1), xim,xis40,1e-3); // 1-loop (spt)
      break;
}}

// B_LSS contribution
// a chooses LCDM (1) or nDGP/MG analytic (2) (tree or 1-loop only) or 3(IDE) or 4 [f(R)] (only numerical)
// b chooses tree analytic (1) or 1-loop analytic (2) or Gil-Marin (3) or Scoccimaro (4) or tree numerical (5) or 1-loop numerical (6)
// params[0] = omega0, params[1]=mg1, params[2]=mg2, params[3] =mg3

// LSS bispectrum
double CMBc::BLSS(int a, int b,  double params[], double l1, double l2, double x) const{
switch (b) {
  case 1:
   return   Integrate(bind(btlssa_integrand,cref(C),cref(P_L), a, params, l1, l2, x, _1),xim,xis,epsrel); // tree level EdS
   break;
  case 2:
    return  Integrate(bind(bllssa_integrand,cref(C),cref(P_L),a, params, l1, l2, x, _1),xim,xis40,epsrel); // 1-loop EdS
    break;
  case 3:
    return  Integrate(bind(bgmlssa_integrand,cref(C),cref(P_L), params, l1, l2, x, _1),xim,xis,epsrel); // Gil-Marin fit
    break;
  case 4:
    return  Integrate(bind(bsclssa_integrand,cref(C),cref(P_L), params, l1, l2, x, _1),xim,xis40,epsrel); // Scoccimaro fit
    break;
  case 5:
  if (a==1) {
    return  Integrate(bind(btlssa_integrand,cref(C),cref(P_L), a, params, l1, l2, x, _1),xim,xis40,epsrel); // For GR we use EdS
    }
  else{
    return  Integrate(bind(btnlssa_integrand,cref(C),cref(P_L),a-1, params, l1, l2, x, _1),xim,xis40,epsrel); // tree level numerical
  }
  break;
  case 6:
  if (a==1) {
    return  Integrate(bind(bllssa_integrand,cref(C),cref(P_L), a, params, l1, l2, x, _1),xim,xis40,epsrel);  // For GR we use EdS
    }
  else{
    return  Integrate(bind(blnlssa_integrand,cref(C),cref(P_L),a-1, params, l1, l2, x, _1),xim,xis40,epsrel); // 1-loop level numerical
  }
    break;
  }
}

// BLSS lensing kernel
double CMBc::BLSS_kernel(int a, int b, double params[], double l1, double l2, double x, double xi) const{
  switch (b) {
    case 1:
     return   btlssa_integrand(cref(C),cref(P_L), a, params, l1, l2, x, xi); // tree level EdS
     break;
    case 2:
      return   bllssa_integrand(cref(C),cref(P_L), a, params, l1, l2, x, xi); // 1-loop EdS
      break;
    case 3:
      return   bgmlssa_integrand(cref(C),cref(P_L), params, l1, l2, x, xi); // Gil-Marin fit
      break;
    }
}

// post born correction
double CMBc::BPB(int a, double omega0, double l1, double l2, double x) const{
   double l1sqr = pow2(l1);
   double l2sqr = pow2(l2);
   double l3 = sqrt(l1sqr + l2sqr + 2.*l1*l2*x);
   double l3sqr = pow2(l3);
   double l1l3 = -(l1sqr + l1*l2*x);
   double l2l3 = -(l2sqr + l1*l2*x);
   return   2.*x/l1/l2 * (l1l3 * Mfunc(a, cref(C), cref(P_L), omega0, l1,l2) + l2l3 * Mfunc(a,cref(C), cref(P_L), omega0, l2,l1))
           +2.*l2l3/l2sqr/l3sqr * (l2*l1*x * Mfunc(a,cref(C), cref(P_L), omega0, l2,l3) + l1l3 * Mfunc(a,cref(C), cref(P_L), omega0, l3,l2))
           +2.*l1l3/l3sqr/l1sqr * (l2l3 * Mfunc(a,cref(C), cref(P_L), omega0, l3,l1) + l2*l1*x * Mfunc(a,cref(C), cref(P_L), omega0, l1,l3));
    }

// CMB lensing spectrum
// nDGP or LCDM are chosen via comdist_init(int a, double omega0, double mg1) function (Only linear growth dependence - 1-loop uses UsA + EdS)
double CMBc::CMBla(int a, int b, double params[],double l1, double l2, double x) const{
  double l1sqr = pow2(l1);
  double l2sqr = pow2(l2);
  double l3 = sqrt(l1sqr + l2sqr + 2.*l1*l2*x);
  int lw1,lw2,lw3;
  lw1 = (int)l1;
  lw2 = (int)l2;
  lw3 = (int)round(l3);
  return  wigner3j(lw1, lw2, lw3) * sqrt((2.*l1+1.)*(2.*l2+1.)*(2.*l3+1.)/4./M_PI) * (BLSS(a,b, params, l1,l2,x) + BPB(a, params[0], l1,l2,x));
}

max_stars_repo_name: PedroCarrilho/ReACT

max_issues_repo_name: PedroCarrilho/ReACT

max_forks_repo_name: PedroCarrilho/ReACT


text:
//==============================================================================
//         Copyright 2014        LASMEA UMR 6602 CNRS/Univ. Clermont II
//         Copyright 2014        LRI    UMR 8623 CNRS/Univ Paris Sud XI
//         Copyright 2014        NumScale SAS
//
//          Distributed under the Boost Software License, Version 1.0.
//                 See accompanying file LICENSE.txt or copy at
//                     http://www.boost.org/LICENSE_1_0.txt
//==============================================================================
#define NT2_UNIT_MODULE "nt2 arithmetic toolbox - fast_toint/simd Mode"

//////////////////////////////////////////////////////////////////////////////
// unit test behavior of arithmetic components in simd mode
//////////////////////////////////////////////////////////////////////////////
///
#include <nt2/arithmetic/include/functions/fast_toint.hpp>
#include <boost/simd/sdk/simd/native.hpp>
#include <nt2/sdk/functor/meta/call.hpp>
#include <nt2/sdk/meta/as_integer.hpp>
#include <nt2/sdk/unit/tests.hpp>
#include <nt2/sdk/unit/tests/relation.hpp>
#include <nt2/sdk/unit/module.hpp>
#include <boost/simd/sdk/simd/io.hpp>
#include <nt2/include/constants/inf.hpp>
#include <nt2/include/constants/minf.hpp>
#include <nt2/include/constants/nan.hpp>
#include <nt2/include/constants/zero.hpp>
#include <nt2/include/constants/one.hpp>
#include <nt2/include/constants/mone.hpp>
#include <boost/simd/include/functions/splat.hpp>

NT2_TEST_CASE_TPL ( fast_toint_real__1_0,  NT2_SIMD_REAL_TYPES)
{
  using nt2::fast_toint;
  using nt2::tag::fast_toint_;
  using boost::simd::native;
  typedef NT2_SIMD_DEFAULT_EXTENSION  ext_t;
  typedef native<T,ext_t>           vT;
  typedef typename nt2::meta::call<fast_toint_(vT)>::type r_t;

  // specific values tests
  NT2_TEST_EQUAL(fast_toint(nt2::Mone<vT>()), nt2::Mone<r_t>());
  NT2_TEST_EQUAL(fast_toint(nt2::One<vT>()), nt2::One<r_t>());
  NT2_TEST_EQUAL(fast_toint(nt2::Zero<vT>()), nt2::Zero<r_t>());

  NT2_TEST_EQUAL(fast_toint(boost::simd::splat<vT>(T(4194304))), boost::simd::splat<r_t>(4194304));
  NT2_TEST_EQUAL(fast_toint(boost::simd::splat<vT>(T(16777216))), boost::simd::splat<r_t>(16777216));

  NT2_TEST_EQUAL(fast_toint(boost::simd::splat<vT>(T(-4194304))), boost::simd::splat<r_t>(-4194304));
  NT2_TEST_EQUAL(fast_toint(boost::simd::splat<vT>(T(-8388608))), boost::simd::splat<r_t>(-8388608));
  NT2_TEST_EQUAL(fast_toint(boost::simd::splat<vT>(T(-16777216))), boost::simd::splat<r_t>(-16777216));
}


text:
/*

Copyright (c) 2005-2019, University of Oxford.
All rights reserved.

University of Oxford means the Chancellor, Masters and Scholars of the
University of Oxford, having an administrative office at Wellington
Square, Oxford OX1 2JD, UK.

This file is part of Chaste.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
 * Redistributions of source code must retain the above copyright notice,
   this list of conditions and the following disclaimer.
 * Redistributions in binary form must reproduce the above copyright notice,
   this list of conditions and the following disclaimer in the documentation
   and/or other materials provided with the distribution.
 * Neither the name of the University of Oxford nor the names of its
   contributors may be used to endorse or promote products derived from this
   software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

*/

#ifndef CHASTECUBOID_HPP_
#define CHASTECUBOID_HPP_
#include "ChasteSerialization.hpp"
#include <boost/serialization/base_object.hpp>
#include "AbstractChasteRegion.hpp"
#include "ChastePoint.hpp"

/**
 * This class defines a 3D cuboid and provides a method to check
 * if a given point is contained in the volume.
 */

template <unsigned SPACE_DIM>
class ChasteCuboid : public AbstractChasteRegion<SPACE_DIM>
{
    /** Needed for serialization. */
    friend class boost::serialization::access;
    /**
     * Archive the member variables.
     *
     * @param archive
     * @param version
     */
    template<class Archive>
    void serialize(Archive & archive, const unsigned int version)
    {
        archive & boost::serialization::base_object<AbstractChasteRegion<SPACE_DIM> >(*this);
    }

private:
    /** Lower vertex of the cuboid. */
    ChastePoint<SPACE_DIM> mLowerCorner;

    /** Upper vertex of the cuboid.  The space-diagonal opposite corner of mLowerCorner. */
    ChastePoint<SPACE_DIM> mUpperCorner;

public:
    /**
     * The cuboid is defined by any of its two space-diagonal opposite corners.
     *
     * @param rLowerPoint Lower vertex of the cuboid.
     * @param rUpperPoint Upper vertex of the cuboid.
     */
    ChasteCuboid(ChastePoint<SPACE_DIM>& rLowerPoint, ChastePoint<SPACE_DIM>& rUpperPoint);

    /**
     * @return true if a given point is contained in the cuboid.
     *
     * @param rPointToCheck Point to be checked to be contained in the cuboid.
     */
    bool DoesContain(const ChastePoint<SPACE_DIM>& rPointToCheck) const;

    /** @return the upper vertex of the cuboid */
    const ChastePoint<SPACE_DIM>& rGetUpperCorner() const;

    /** @return the lower vertex of the cuboid */
    const ChastePoint<SPACE_DIM>& rGetLowerCorner() const;

    /**
     * @param rDimension dimension
     * @return the width in a particular dimension */
    double GetWidth(unsigned rDimension) const;

    /**
     * @return the longest axis of the cuboid (<SPACE_DIM).
     * In the case where there's a tie, then any of the longest
     * axes are returned.
     */
     unsigned GetLongestAxis() const;
};

// Declare identifier for the serializer
#include "SerializationExportWrapper.hpp"
EXPORT_TEMPLATE_CLASS_SAME_DIMS(ChasteCuboid)

namespace boost
{
namespace serialization
{

template<class Archive, unsigned SPACE_DIM>
inline void save_construct_data(
    Archive & ar, const ChasteCuboid<SPACE_DIM> * t, const unsigned int file_version)
{
    const ChastePoint<SPACE_DIM>* p_upper_corner =  &(t->rGetUpperCorner());
    const ChastePoint<SPACE_DIM>* p_lower_corner =  &(t->rGetLowerCorner());
    ar & p_upper_corner;
    ar & p_lower_corner;
}

/**
 * Allow us to not need a default constructor, by specifying how Boost should
 * instantiate an instance (using existing constructor)
 */
template<class Archive, unsigned SPACE_DIM>
inline void load_construct_data(
    Archive & ar, ChasteCuboid<SPACE_DIM> * t, const unsigned int file_version)
{
    ChastePoint<SPACE_DIM>* p_upper_corner;
    ChastePoint<SPACE_DIM>* p_lower_corner;

    ar & p_upper_corner;
    ar & p_lower_corner;

    ::new(t)ChasteCuboid<SPACE_DIM>((*p_lower_corner), (*p_upper_corner));

    delete p_upper_corner;
    delete p_lower_corner;
}
}
} // namespace ...

#endif /*CHASTECUBOID_HPP_*/