text:
//-----------------------------------------------------------------------------
// boost variant/recursive_variant.hpp header file
// See http://www.boost.org for updates, documentation, and revision history.
//-----------------------------------------------------------------------------
//
// Copyright (c) 2003 Eric Friedman
// Copyright (c) 2013-2020 Antony Polukhin
//
// 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)

#ifndef BOOST_VARIANT_RECURSIVE_VARIANT_HPP
#define BOOST_VARIANT_RECURSIVE_VARIANT_HPP

#include <boost/variant/variant_fwd.hpp>
#include <boost/variant/detail/enable_recursive.hpp>
#include <boost/variant/detail/substitute_fwd.hpp>
#include <boost/variant/detail/make_variant_list.hpp>
#include <boost/variant/detail/over_sequence.hpp>

#include <boost/mpl/aux_/lambda_arity_param.hpp>

#include <boost/mpl/equal.hpp>
#include <boost/mpl/eval_if.hpp>
#include <boost/mpl/identity.hpp>
#include <boost/mpl/if.hpp>
#include <boost/mpl/protect.hpp>
#include <boost/mpl/transform.hpp>
#include <boost/type_traits/is_same.hpp>
#include <boost/preprocessor/cat.hpp>
#include <boost/preprocessor/repeat.hpp>

#include <boost/mpl/bool.hpp>
#include <boost/mpl/is_sequence.hpp>
#include <boost/variant/variant.hpp>

namespace boost {

namespace detail { namespace variant {

///////////////////////////////////////////////////////////////////////////////
// (detail) metafunction specialization substitute
//
// Handles embedded variant types when substituting for recursive_variant_.
//

#if !defined(BOOST_VARIANT_DETAIL_NO_SUBSTITUTE)

template <
      BOOST_VARIANT_ENUM_PARAMS(typename T)
    , typename RecursiveVariant
      BOOST_MPL_AUX_LAMBDA_ARITY_PARAM(typename Arity)
    >
struct substitute<
      ::boost::variant<
          recursive_flag< T0 >
        , BOOST_VARIANT_ENUM_SHIFTED_PARAMS(T)
        >
    , RecursiveVariant
    , ::boost::recursive_variant_
      BOOST_MPL_AUX_LAMBDA_ARITY_PARAM(Arity)
    >
{
    typedef ::boost::variant<
          recursive_flag< T0 >
        , BOOST_VARIANT_ENUM_SHIFTED_PARAMS(T)
        > type;
};

template <
      BOOST_VARIANT_ENUM_PARAMS(typename T)
    , typename RecursiveVariant
      BOOST_MPL_AUX_LAMBDA_ARITY_PARAM(typename Arity)
    >
struct substitute<
      ::boost::variant<
          ::boost::detail::variant::over_sequence< T0 >
        , BOOST_VARIANT_ENUM_SHIFTED_PARAMS(T)
        >
    , RecursiveVariant
    , ::boost::recursive_variant_
      BOOST_MPL_AUX_LAMBDA_ARITY_PARAM(Arity)
    >
{
private:

    typedef T0 initial_types;

    typedef typename mpl::transform<
          initial_types
        , mpl::protect< quoted_enable_recursive<RecursiveVariant,mpl::true_> >
        >::type types;

public:

    typedef typename mpl::if_<
          mpl::equal<initial_types, types, ::boost::is_same<mpl::_1, mpl::_2> >
        , ::boost::variant<
              ::boost::detail::variant::over_sequence< T0 >
            , BOOST_VARIANT_ENUM_SHIFTED_PARAMS(T)
            >
        , ::boost::variant< over_sequence<types> >
        >::type type;
};

template <
      BOOST_VARIANT_ENUM_PARAMS(typename T)
    , typename RecursiveVariant
      BOOST_MPL_AUX_LAMBDA_ARITY_PARAM(typename Arity)
    >
struct substitute<
      ::boost::variant< BOOST_VARIANT_ENUM_PARAMS(T) >
    , RecursiveVariant
    , ::boost::recursive_variant_
      BOOST_MPL_AUX_LAMBDA_ARITY_PARAM(Arity)
    >
{
#if !defined(BOOST_VARIANT_DO_NOT_USE_VARIADIC_TEMPLATES)

    typedef ::boost::variant<
        typename enable_recursive<
              T0
            , RecursiveVariant
            , mpl::true_
        >::type,
        typename enable_recursive<
              TN
            , RecursiveVariant
            , mpl::true_
        >::type...
    > type;

#else // defined(BOOST_VARIANT_DO_NOT_USE_VARIADIC_TEMPLATES)

private: // helpers, for metafunction result (below)

    #define BOOST_VARIANT_AUX_ENABLE_RECURSIVE_TYPEDEFS(z,N,_)  \
        typedef typename enable_recursive<   \
              BOOST_PP_CAT(T,N)              \
            , RecursiveVariant               \
            , mpl::true_                     \
            >::type BOOST_PP_CAT(wknd_T,N);  \
        /**/

    BOOST_PP_REPEAT(
          BOOST_VARIANT_LIMIT_TYPES
        , BOOST_VARIANT_AUX_ENABLE_RECURSIVE_TYPEDEFS
        , _
        )

    #undef BOOST_VARIANT_AUX_ENABLE_RECURSIVE_TYPEDEFS

public: // metafunction result

    typedef ::boost::variant< BOOST_VARIANT_ENUM_PARAMS(wknd_T) > type;
#endif // BOOST_VARIANT_DO_NOT_USE_VARIADIC_TEMPLATES workaround
};

#else // defined(BOOST_VARIANT_DETAIL_NO_SUBSTITUTE)

//
// no specializations: embedded variants unsupported on these compilers!
//

#endif // !defined(BOOST_VARIANT_DETAIL_NO_SUBSTITUTE)

}} // namespace detail::variant

///////////////////////////////////////////////////////////////////////////////
// metafunction make_recursive_variant
//
// See docs and boost/variant/variant_fwd.hpp for more information.
//
template < BOOST_VARIANT_ENUM_PARAMS(typename T) >
struct make_recursive_variant
{
public: // metafunction result

    typedef boost::variant<
          detail::variant::recursive_flag< T0 >
        , BOOST_VARIANT_ENUM_SHIFTED_PARAMS(T)
        > type;

};

///////////////////////////////////////////////////////////////////////////////
// metafunction make_recursive_variant_over
//
// See docs and boost/variant/variant_fwd.hpp for more information.
//
template <typename Types>
struct make_recursive_variant_over
{
private: // precondition assertions

    BOOST_STATIC_ASSERT(( ::boost::mpl::is_sequence<Types>::value ));

public: // metafunction result

    typedef typename make_recursive_variant<
          detail::variant::over_sequence< Types >
        >::type type;

};

} // namespace boost

#endif // BOOST_VARIANT_RECURSIVE_VARIANT_HPP

max_stars_repo_name: mamil/demo

max_issues_repo_name: mamil/demo

max_forks_repo_name: mamil/demo


text:
/*
 * Copyright 2020 Yamana Laboratory, Waseda University
 * Supported by JST CREST Grant Number JPMJCR1503, Japan.
 *
 * 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.
 */

#define _POSIX_SOURCE // for mkdir
#include <sys/stat.h> // for mkdir

#include <unordered_map>
#include <string>
#include <fstream>
#include <boost/algorithm/string.hpp>
#include <map>

#include "FHE.h"
#include "EncryptedArray.h"

#include <stdsc/stdsc_exception.hpp>
#include <stdsc/stdsc_log.hpp>

#include <sses_share/sses_utility.hpp>
#include <sses_share/sses_computation_param.hpp>
#include <sses_share/sses_mask.hpp>
#include <sses_share/sses_types.hpp>

#include <sses_server/sses_server_db.hpp>

#define ENABLE_LOCAL_DEBUG
#ifdef ENABLE_LOCAL_DEBUG
#include <sses_share/sses_fhe_debug.hpp>
#endif

namespace sses_server
{

static constexpr char* DEFAULT_DBBASIC_FILENAME = (char*)"dbbasics.bin";
static constexpr char* DEFAULT_MEDINV_FILENAME = (char*)"med.inv";
static constexpr char* DEFAULT_SIDEINV_FILENAME = (char*)"side.inv";
static constexpr char* DEFAULT_ENCDATA_DIRNAME = (char*)"encdata";
static constexpr char* DEFAULT_AUXDATA_DIRNAME = (char*)"auxdata";

namespace csvcolumns
{
    constexpr int ID_IDX = 0;
    constexpr int MEDICINE_ID_IDX = 1;
    constexpr int SYMPTOM_ID_IDX = 2;
    constexpr int AGE_IDX = 8;
    constexpr int GENDER_IDX = 9;
}

struct Record
{
    int maskValue;
    std::set<size_t> medicineIds;
    std::set<size_t> symptomIds;
};

DBBasicFile::DBBasicFile()
    : dbstatus(false),
      totalRecordsNum(0),
      totalMedicinesNum(0),
      totalSymptomsNum(0)
{}

DBBasicFile::DBBasicFile(std::string& filepath)
    : dbstatus(false),
      totalRecordsNum(0),
      totalMedicinesNum(0),
      totalSymptomsNum(0)
{
    try {
        load_from_file(filepath);
    } catch (const stdsc::AbstractException& e) {
        STDSC_LOG_WARN(e.what());
    }
}

void DBBasicFile::write_to_file(const std::string& filepath,
                                const size_t totalRecordsNum,
                                const size_t totalMedicinesNum,
                                const size_t totalSymptomsNum)
{
    std::ofstream ofs(filepath, std::ios::binary);
    if (!ofs.is_open())
    {
        std::ostringstream oss;
        oss << "failed to open. (" << filepath << ")";
        STDSC_THROW_FILE(oss.str());
    }
    ofs << 1 << std::endl;
    ofs << totalRecordsNum << std::endl;
    ofs << totalMedicinesNum << std::endl;
    ofs << totalSymptomsNum << std::endl;
}

void DBBasicFile::load_from_file(const std::string& filepath)
{
    std::ifstream ifs(filepath, std::ios::binary);
    if (!ifs.is_open())
    {
        std::ostringstream oss;
        oss << "failed to open. (" << filepath << ")";
        STDSC_THROW_FILE(oss.str());
    }
    ifs >> dbstatus;
    ifs >> totalRecordsNum;
    ifs >> totalMedicinesNum;
    ifs >> totalSymptomsNum;
}

bool DBBasicFile::status() const
{
    return dbstatus;
}

std::string DBBasicFile::to_string() const
{
    std::ostringstream oss;
    oss << totalRecordsNum << ", ";
    oss << totalMedicinesNum << ", ";
    oss << totalSymptomsNum;
    return oss.str();
}

struct DatasetInfo
{
    explicit DatasetInfo(std::string& dir)
        : dir_(dir)
    {}
    virtual ~DatasetInfo() = default;

    bool is_enable() const
    {
        bool ret = false;
        if (sses_share::utility::dir_exist(dir_)) {

            auto filepath = dbbasic_filepath();
            DBBasicFile dbbasic(filepath);
            if (dbbasic.status()) {
                ret = true;
            }

        }
        return ret;
    }

    const std::string dbbasic_filepath() const
    {
        std::ostringstream oss;
        oss << dir_ << "/";
        oss << DEFAULT_DBBASIC_FILENAME;
        return oss.str();
    }

    const std::string medinv_filepath() const
    {
        auto dir = auxdata_dirpath() + "/";
        return dir + DEFAULT_MEDINV_FILENAME;
    }

    const std::string sideinv_filepath() const
    {
        auto dir = auxdata_dirpath() + "/";
        return dir + DEFAULT_SIDEINV_FILENAME;
    }

    const std::string encdata_dirpath() const
    {
        std::ostringstream oss;
        oss << dir_ << "/";
        oss << DEFAULT_ENCDATA_DIRNAME;
        return oss.str();
    }

    const std::string auxdata_dirpath() const
    {
        std::ostringstream oss;
        oss << dir_ << "/";
        oss << DEFAULT_AUXDATA_DIRNAME;
        return oss.str();
    }

    std::string dir_;
};

struct DB::Impl
{
    static constexpr char* LIST_FILENAME = (char*)"list.txt";

    Impl(const std::string& db_basedir)
        : db_basedir_(db_basedir)
    {
        {
            std::ostringstream oss;
            oss << db_basedir << "/" << LIST_FILENAME;
            list_filepath_ = oss.str();
        }

        load_listfile(list_filepath_);

        std::ostringstream oss;
        for (auto m : map_) {
            auto key_id = m.first;
            auto dsinfo = m.second;
            oss << "{key_id:" << key_id;
            oss << ", dir:" << dsinfo.dir_.c_str();
            oss << ", dbstatus:" << dsinfo.is_enable();
            oss << "} ";
        }
        STDSC_LOG_INFO("Initialized DB : %s", oss.str().c_str());
    }

    bool is_enable(const int32_t key_id) const
    {
        bool ret = false;
        if (map_.count(key_id) > 0) {
            auto& dsinfo = map_.at(key_id);
            ret = dsinfo.is_enable();
        }
        return ret;
    }

    void setup(const int32_t key_id,
               const std::string& db_src_filepath,
               const FHEcontext& context,
               const FHEPubKey& pubkey)
    {
        if (is_enable(key_id)) {
            return;
        }

        std::string top_dir = db_basedir_ + "/db_" + std::to_string(key_id);

        if (sses_share::utility::dir_exist(top_dir)) {
            sses_share::utility::remove_dir(top_dir);
            STDSC_LOG_INFO("Remove directory : %s", top_dir.c_str());
        }

        STDSC_LOG_INFO("Creating new DB. [In:%s, Out: %s]",
                       db_src_filepath.c_str(),
                       top_dir.c_str());

        STDSC_THROW_FILE_IF_CHECK(mkdir(top_dir.c_str(), S_IRWXU) == 0,
                                  "Err: failed to create DB directory");
        auto encdata_dir = top_dir + "/encdata";
        STDSC_THROW_FILE_IF_CHECK(mkdir(encdata_dir.c_str(), S_IRWXU) == 0,
                                  "Err: failed to create DB directory");
        auto auxdata_dir = top_dir + "/auxdata";
        STDSC_THROW_FILE_IF_CHECK(mkdir(auxdata_dir.c_str(), S_IRWXU) == 0,
                                  "Err: failed to create DB directory");

        std::ifstream csvIfs(db_src_filepath);
        std::string line;
        std::map<size_t, Record> records;
        std::set<size_t> totalMedicines;
        std::set<size_t> totalSymptoms;

        std::vector<std::string> record;
        getline(csvIfs, line); // header
        while (getline(csvIfs, line))
        {
            boost::algorithm::split(record, line, boost::is_any_of(","));
            size_t id = stoi(record[csvcolumns::ID_IDX]);
            size_t medicineId = stoi(record[csvcolumns::MEDICINE_ID_IDX]);
            size_t symptomId = stoi(record[csvcolumns::SYMPTOM_ID_IDX]);
            totalMedicines.insert(medicineId);
            totalSymptoms.insert(symptomId);
            if (records.count(id) == 0)
            {
                int age = stoi(record[csvcolumns::AGE_IDX]);
                assert(age >= 0 && age <= 122);
                int genderData = stoi(record[csvcolumns::GENDER_IDX]);
                assert(genderData == 1 || genderData == 2 || genderData == 3);
                auto gender = sses_share::GENDER_INT_MAP.at(genderData);
                int mask = age + static_cast<int>(gender) * 128 + 5;
                std::set<size_t> medicineIds;
                medicineIds.insert(medicineId);
                std::set<size_t> symptomIds;
                symptomIds.insert(symptomId);
                records.insert(
                    std::make_pair(id, Record{mask, medicineIds, symptomIds}));
            }
            else
            {
                records.at(id).medicineIds.insert(medicineId);
                records.at(id).symptomIds.insert(symptomId);
            }
        }
        csvIfs.close();

        auto totalRecordsNum = records.size();
        auto totalMedicinesNum = totalMedicines.size();
        auto totalSymptomsNum = totalSymptoms.size();

        auto dbbasicfilepath = top_dir + "/" + std::string(DEFAULT_DBBASIC_FILENAME);
        auto medinvfilepath = auxdata_dir + "/" + std::string(DEFAULT_MEDINV_FILENAME);
        auto sideinvfilepath = auxdata_dir + "/" + std::string(DEFAULT_SIDEINV_FILENAME);

        DBBasicFile dbbasicfile;
        dbbasicfile.write_to_file(dbbasicfilepath,
                                  totalRecordsNum,
                                  totalMedicinesNum,
                                  totalSymptomsNum);
        STDSC_LOG_INFO("Created DBBasic file. [%s]", dbbasicfilepath.c_str());

        std::ofstream fmedinv(medinvfilepath, std::ios::binary);
        fmedinv << totalMedicinesNum << std::endl;
        for (const auto& medicineId : totalMedicines) {
            fmedinv << medicineId << ":" << 0 << std::endl;
        }
        fmedinv.close();

        std::ofstream fsideinv(sideinvfilepath, std::ios::binary);
        fsideinv << totalSymptomsNum << std::endl;
        for (const auto& symptomId : totalSymptoms) {
            fsideinv << symptomId << ":" << 0 << std::endl;
        }
        fsideinv.close();

        NTL::ZZX G = context.alMod.getFactorsOverZZ()[0];
        EncryptedArray ea(context, G);

        STDSC_LOG_INFO("Start generating DB data.");

        for (const auto& recordPair : records)
        {
            size_t recordId = recordPair.first;
            Record record = recordPair.second;
            int mask = record.maskValue;

             std::string encfilepath = encdata_dir + "/" + std::to_string(recordId) + ".bin";
            std::string auxfilepath = auxdata_dir + "/" + std::to_string(recordId) + ".bin";

            STDSC_LOG_INFO("  recID:%lu, numMed:%lu, numSide:%lu",
                           recordId,
                           record.medicineIds.size(),
                           record.symptomIds.size());

            Ctxt encmask(pubkey);
            pubkey.Encrypt(encmask, NTL::to_ZZX(mask));
            {
                std::ofstream ofs(encfilepath, std::ios::binary);
                ofs << encmask;
                STDSC_LOG_INFO("    | %s", encfilepath.c_str());
            }

            std::vector<size_t> meds, sides;
            for (const auto v : record.medicineIds) {
                meds.push_back(v);
            }
            for (const auto v : record.symptomIds) {
                sides.push_back(v);
            }

            {
                std::ofstream ofs(auxfilepath, std::ios::binary);
                ofs << "Medicine: [" << meds[0];
                for (size_t i = 1; i < meds.size(); ++i) {
                    ofs << ", " << meds[i];
                }
                ofs << "]" << std::endl;
                ofs << "Side Effect: [" << sides[0];
                for (size_t i = 1; i < sides.size(); ++i) {
                    ofs << ", " << sides[i];
                }
                ofs << "]" << std::endl;
                STDSC_LOG_INFO("    | %s", auxfilepath.c_str());
            }

            update(recordId,
                   record.medicineIds.size(), meds,
                   record.symptomIds.size(), sides,
                   auxfilepath,
                   medinvfilepath,
                   sideinvfilepath);
            STDSC_LOG_INFO("    | %s", medinvfilepath.c_str());
            STDSC_LOG_INFO("    | %s", sideinvfilepath.c_str());
        }

        STDSC_LOG_INFO("Finish generating DB data.");

        map_.emplace(key_id, DatasetInfo(top_dir));
        save_listfile(list_filepath_);
        STDSC_LOG_INFO("Updated List file. [%s]", list_filepath_.c_str());
    }

    std::string dbbasic_filepath(const int32_t key_id) const
    {
        return map_.at(key_id).dbbasic_filepath();
    }

    std::string medinv_filepath(const int32_t key_id) const
    {
        return map_.at(key_id).medinv_filepath();
    }

    std::string sideinv_filepath(const int32_t key_id) const
    {
        return map_.at(key_id).sideinv_filepath();
    }

    std::string encdata_dirpath(const int32_t key_id) const
    {
        return map_.at(key_id).encdata_dirpath();
    }

    std::string auxdata_dirpath(const int32_t key_id) const
    {
        return map_.at(key_id).auxdata_dirpath();
    }

private:

    void load_listfile(const std::string& filepath)
    {
        std::ifstream ifs(filepath);
        if (ifs.is_open()) {
            size_t n;
            ifs >> n;

            for (size_t i=0; i<n; ++i) {
                std::string line;
                ifs >> line;

                std::vector<std::string> elems;
                boost::algorithm::split(elems, line, boost::is_any_of(","));

                map_.emplace(std::stoi(elems[0]), DatasetInfo(elems[1]));
            }
        }
    }

    void save_listfile(const std::string& filepath)
    {
        std::ofstream ofs(filepath);
        if (ofs.is_open()) {
            ofs << map_.size() << std::endl;

            for (const auto& v : map_) {
                ofs << v.first;
                ofs << ",";
                ofs << v.second.dir_;
                ofs << std::endl;
            }
        }
    }

    void update(const size_t recordId,
                const size_t numMed,
                const std::vector<size_t>& meds,
                const size_t numSide,
                const std::vector<size_t>& sides,
                const std::string& auxfilepath,
                const std::string& medinvfilepath,
                const std::string& sideinvfilepath)
    {
        std::map<int, std::vector<int>> medIndex;
        std::map<int, std::vector<int>> sideIndex;

        std::ofstream faux(auxfilepath.c_str(), std::ios::binary);
        faux << "Medicine: [" << meds[0];
        for (size_t i = 1; i < numMed; ++i) {
            faux << ", " << meds[i];
        }
        faux << "]" << std::endl;
        faux << "Side Effect: [" << sides[0];
        for (size_t i = 1; i < numSide; ++i) {
            faux << ", " << sides[i];
        }
        faux << "]" << std::endl;
        faux.close();

        std::ifstream findexmed(medinvfilepath, std::ios::binary);
        std::string line;
        std::vector<std::string> medInfo;
        int numMedDict;
        findexmed >> numMedDict;
        for (int i = 0; i < numMedDict; ++i)
        {
            findexmed >> line;
            boost::algorithm::split(medInfo, line, boost::is_any_of(":"));
            int medId = std::stoi(medInfo[0]);
            int numindex = std::stoi(medInfo[1]);
            std::vector<int> tempindex;
            while (numindex--)
            {
                int temprec;
                findexmed >> temprec;
                tempindex.push_back(temprec);
            }
            medIndex.insert(std::make_pair(medId, tempindex));
        }
        findexmed.close();

        std::ifstream findexside(sideinvfilepath, std::ios::binary);
        std::vector<std::string> sideInfo;
        int numSideDict;
        findexside >> numSideDict;
        for (int i = 0; i < numSideDict; ++i)
        {
            findexside >> line;
            boost::algorithm::split(sideInfo, line, boost::is_any_of(":"));
            int sideId = std::stoi(sideInfo[0]);
            int numindex = std::stoi(sideInfo[1]);
            std::vector<int> tempindex;
            while (numindex--)
            {
                int temprec;
                findexside >> temprec;
                tempindex.push_back(temprec);
            }
            sideIndex.insert(make_pair(sideId, tempindex));
        }
        findexside.close();

        for (size_t i = 0; i < numMed; ++i) {
            medIndex[meds[i]].push_back(recordId);
        }
        for (size_t i = 0; i < numSide; ++i) {
            sideIndex[sides[i]].push_back(recordId);
        }
        STDSC_LOG_TRACE("Updated rev index");

        std::ofstream findexmedw(medinvfilepath, std::ios::binary);
        findexmedw << numMedDict << std::endl;

        for (const auto& medIdxPair : medIndex)
        {
            int medId = medIdxPair.first;
            std::vector<int> recordIds = medIdxPair.second;
            findexmedw << medId << ":" << recordIds.size() << std::endl;
            for (const int& recordId : recordIds)
            {
                findexmedw << recordId << std::endl;
            }
        }
        findexmedw.close();

        std::ofstream findexsidew(sideinvfilepath, std::ios::binary);
        findexsidew << numSideDict << std::endl;
        for (const auto& sideIdxPair : sideIndex)
        {
            int sideId = sideIdxPair.first;
            std::vector<int> recordIds = sideIdxPair.second;
            findexsidew << sideId << ":" << recordIds.size() << std::endl;
            for (const int& recordId : recordIds)
            {
                findexsidew << recordId << std::endl;
            }
        }
        findexsidew.close();

        STDSC_LOG_TRACE("New record id: %lu", recordId);
        STDSC_LOG_TRACE("DB successfully updated.");
    }

private:
    std::string db_basedir_;
    std::string list_filepath_;
    std::unordered_map<int32_t, DatasetInfo> map_;
};

DB::DB(const std::string& db_basedir)
  : pimpl_(new Impl(db_basedir))
{}

bool DB::is_enable(const int32_t key_id) const
{
    return pimpl_->is_enable(key_id);
}

void DB::setup(const int32_t key_id,
               const std::string& db_src_filepath,
               const FHEcontext& context,
               const FHEPubKey& pubkey)
{
    pimpl_->setup(key_id, db_src_filepath, context, pubkey);
}

std::string DB::dbbasic_filepath(const int32_t key_id) const
{
    return pimpl_->dbbasic_filepath(key_id);
}

std::string DB::medinv_filepath(const int32_t key_id) const
{
        return pimpl_->medinv_filepath(key_id);
}

std::string DB::sideinv_filepath(const int32_t key_id) const
{
        return pimpl_->sideinv_filepath(key_id);
}

std::string DB::encdata_dirpath(const int32_t key_id) const
{
    return pimpl_->encdata_dirpath(key_id);
}

std::string DB::auxdata_dirpath(const int32_t key_id) const
{
    return pimpl_->auxdata_dirpath(key_id);
}

} /* namespace sses_server */

ext: cpp

max_stars_repo_name: yamanalab/SecureSideEffectSearch

max_issues_repo_name: yamanalab/SecureSideEffectSearch

max_forks_repo_name: yamanalab/SecureSideEffectSearch


text:
/*=============================================================================
    Copyright (c) 1999-2003 Jaakko Jarvi
    Copyright (c) 2001-2011 Joel de Guzman

    Distributed under the Boost Software License, Version 1.0. (See accompanying
    file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
==============================================================================*/
#if !defined(FUSION_GREATER_05052005_0432)
#define FUSION_GREATER_05052005_0432

#include <boost/fusion/sequence/intrinsic/begin.hpp>
#include <boost/fusion/sequence/intrinsic/end.hpp>
#include <boost/fusion/sequence/intrinsic/size.hpp>
#include <boost/fusion/sequence/comparison/enable_comparison.hpp>

#if defined(FUSION_DIRECT_OPERATOR_USAGE)
#include <boost/fusion/sequence/comparison/detail/greater.hpp>
#else
#include <boost/fusion/sequence/comparison/less.hpp>
#endif

namespace boost { namespace fusion
{
    template <typename Seq1, typename Seq2>
    inline bool
    greater(Seq1 const& a, Seq2 const& b)
    {
#if defined(FUSION_DIRECT_OPERATOR_USAGE)
        return detail::sequence_greater<Seq1 const, Seq2 const>::
            call(fusion::begin(a), fusion::begin(b));
#else
        return (b < a);
#endif
    }

    namespace operators
    {
        template <typename Seq1, typename Seq2>
        inline typename
            enable_if<
                traits::enable_comparison<Seq1, Seq2>
              , bool
            >::type
        operator>(Seq1 const& a, Seq2 const& b)
        {
            return fusion::greater(a, b);
        }
    }
    using operators::operator>;
}}

#endif

ext: hpp


text:
#ifndef STAN_MATH_PRIM_SCAL_PROB_CAUCHY_RNG_HPP
#define STAN_MATH_PRIM_SCAL_PROB_CAUCHY_RNG_HPP

#include <boost/random/cauchy_distribution.hpp>
#include <boost/random/variate_generator.hpp>
#include <stan/math/prim/scal/err/check_consistent_sizes.hpp>
#include <stan/math/prim/scal/err/check_finite.hpp>
#include <stan/math/prim/scal/err/check_not_nan.hpp>
#include <stan/math/prim/scal/err/check_positive_finite.hpp>
#include <stan/math/prim/scal/fun/constants.hpp>
#include <stan/math/prim/scal/fun/log1p.hpp>
#include <stan/math/prim/scal/fun/square.hpp>
#include <stan/math/prim/scal/fun/value_of.hpp>
#include <stan/math/prim/scal/meta/include_summand.hpp>

namespace stan {
  namespace math {

    /**
     * Return a pseudorandom Cauchy variate for the given location and scale
     * using the specified random number generator.
     *
     * @tparam RNG type of random number generator
     * @param mu location parameter
     * @param sigma positive scale parameter
     * @param rng random number generator
     * @return Cauchy random variate
     * @throw std::domain_error if mu is infinite or sigma is nonpositive
     */
    template <class RNG>
    inline double
    cauchy_rng(double mu,
               double sigma,
               RNG& rng) {
      using boost::variate_generator;
      using boost::random::cauchy_distribution;

      static const char* function("cauchy_rng");

      check_finite(function, "Location parameter", mu);
      check_positive_finite(function, "Scale parameter", sigma);

      variate_generator<RNG&, cauchy_distribution<> >
        cauchy_rng(rng, cauchy_distribution<>(mu, sigma));
      return cauchy_rng();
    }

  }
}
#endif

max_stars_repo_name: yizhang-cae/torsten

max_issues_repo_name: yizhang-cae/torsten

max_forks_repo_name: yizhang-cae/torsten


text:
//
// Copyright (c) 2002--2010
// Toon Knapen, Karl Meerbergen, Kresimir Fresl,
// Thomas Klimpel and Rutger ter Borg
//
// 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)
//
// THIS FILE IS AUTOMATICALLY GENERATED
// PLEASE DO NOT EDIT!
//

#ifndef BOOST_NUMERIC_BINDINGS_LAPACK_DRIVER_GEJSV_HPP
#define BOOST_NUMERIC_BINDINGS_LAPACK_DRIVER_GEJSV_HPP

#include <boost/assert.hpp>
#include <boost/numeric/bindings/begin.hpp>
#include <boost/numeric/bindings/detail/array.hpp>
#include <boost/numeric/bindings/is_column_major.hpp>
#include <boost/numeric/bindings/is_mutable.hpp>
#include <boost/numeric/bindings/lapack/workspace.hpp>
#include <boost/numeric/bindings/remove_imaginary.hpp>
#include <boost/numeric/bindings/size.hpp>
#include <boost/numeric/bindings/stride.hpp>
#include <boost/numeric/bindings/value_type.hpp>
#include <boost/static_assert.hpp>
#include <boost/type_traits/is_same.hpp>
#include <boost/type_traits/remove_const.hpp>

//
// The LAPACK-backend for gejsv is the netlib-compatible backend.
//
#include <boost/numeric/bindings/lapack/detail/lapack.h>
#include <boost/numeric/bindings/lapack/detail/lapack_option.hpp>

namespace boost {
namespace numeric {
namespace bindings {
namespace lapack {

//
// The detail namespace contains value-type-overloaded functions that
// dispatch to the appropriate back-end LAPACK-routine.
//
namespace detail {

//
// Overloaded function for dispatching to
// * netlib-compatible LAPACK backend (the default), and
// * float value-type.
//
inline std::ptrdiff_t gejsv( const char joba, const char jobu, const char jobv,
        const char jobr, const char jobt, const char jobp,
        const fortran_int_t m, const fortran_int_t n, float* a,
        const fortran_int_t lda, float* sva, float* u,
        const fortran_int_t ldu, float* v, const fortran_int_t ldv,
        float* work, const fortran_int_t lwork, fortran_int_t* iwork ) {
    fortran_int_t info(0);
    LAPACK_SGEJSV( &joba, &jobu, &jobv, &jobr, &jobt, &jobp, &m, &n, a, &lda,
            sva, u, &ldu, v, &ldv, work, &lwork, iwork, &info );
    return info;
}

//
// Overloaded function for dispatching to
// * netlib-compatible LAPACK backend (the default), and
// * double value-type.
//
inline std::ptrdiff_t gejsv( const char joba, const char jobu, const char jobv,
        const char jobr, const char jobt, const char jobp,
        const fortran_int_t m, const fortran_int_t n, double* a,
        const fortran_int_t lda, double* sva, double* u,
        const fortran_int_t ldu, double* v, const fortran_int_t ldv,
        double* work, const fortran_int_t lwork, fortran_int_t* iwork ) {
    fortran_int_t info(0);
    LAPACK_DGEJSV( &joba, &jobu, &jobv, &jobr, &jobt, &jobp, &m, &n, a, &lda,
            sva, u, &ldu, v, &ldv, work, &lwork, iwork, &info );
    return info;
}

} // namespace detail

//
// Value-type based template class. Use this class if you need a type
// for dispatching to gejsv.
//
template< typename Value >
struct gejsv_impl {

    typedef Value value_type;
    typedef typename remove_imaginary< Value >::type real_type;

    //
    // Static member function for user-defined workspaces, that
    // * Deduces the required arguments for dispatching to LAPACK, and
    // * Asserts that most arguments make sense.
    //
    template< typename MatrixA, typename VectorSVA, typename MatrixU,
            typename MatrixV, typename WORK, typename IWORK >
    static std::ptrdiff_t invoke( const char joba, const char jobu,
            const char jobv, const char jobr, const char jobt,
            const char jobp, MatrixA& a, VectorSVA& sva, MatrixU& u,
            MatrixV& v, detail::workspace2< WORK, IWORK > work ) {
        namespace bindings = ::boost::numeric::bindings;
        BOOST_STATIC_ASSERT( (bindings::is_column_major< MatrixA >::value) );
        BOOST_STATIC_ASSERT( (bindings::is_column_major< MatrixU >::value) );
        BOOST_STATIC_ASSERT( (bindings::is_column_major< MatrixV >::value) );
        BOOST_STATIC_ASSERT( (boost::is_same< typename remove_const<
                typename bindings::value_type< MatrixA >::type >::type,
                typename remove_const< typename bindings::value_type<
                VectorSVA >::type >::type >::value) );
        BOOST_STATIC_ASSERT( (boost::is_same< typename remove_const<
                typename bindings::value_type< MatrixA >::type >::type,
                typename remove_const< typename bindings::value_type<
                MatrixU >::type >::type >::value) );
        BOOST_STATIC_ASSERT( (boost::is_same< typename remove_const<
                typename bindings::value_type< MatrixA >::type >::type,
                typename remove_const< typename bindings::value_type<
                MatrixV >::type >::type >::value) );
        BOOST_STATIC_ASSERT( (bindings::is_mutable< MatrixA >::value) );
        BOOST_STATIC_ASSERT( (bindings::is_mutable< VectorSVA >::value) );
        BOOST_STATIC_ASSERT( (bindings::is_mutable< MatrixU >::value) );
        BOOST_STATIC_ASSERT( (bindings::is_mutable< MatrixV >::value) );
        BOOST_ASSERT( bindings::size(work.select(fortran_int_t())) >=
                min_size_iwork( bindings::size_row(a),
                bindings::size_column(a) ));
        BOOST_ASSERT( bindings::size(work.select(real_type())) >=
                min_size_work( joba, jobu, jobv, bindings::size_row(a),
                bindings::size_column(a) ));
        BOOST_ASSERT( bindings::size_minor(a) == 1 ||
                bindings::stride_minor(a) == 1 );
        BOOST_ASSERT( bindings::size_minor(u) == 1 ||
                bindings::stride_minor(u) == 1 );
        BOOST_ASSERT( bindings::size_minor(v) == 1 ||
                bindings::stride_minor(v) == 1 );
        BOOST_ASSERT( bindings::size_row(a) >= 0 );
        BOOST_ASSERT( bindings::stride_major(a) >= std::max< std::ptrdiff_t >(1,
                bindings::size_row(a)) );
        return detail::gejsv( joba, jobu, jobv, jobr, jobt, jobp,
                bindings::size_row(a), bindings::size_column(a),
                bindings::begin_value(a), bindings::stride_major(a),
                bindings::begin_value(sva), bindings::begin_value(u),
                bindings::stride_major(u), bindings::begin_value(v),
                bindings::stride_major(v),
                bindings::begin_value(work.select(real_type())),
                bindings::size(work.select(real_type())),
                bindings::begin_value(work.select(fortran_int_t())) );
    }

    //
    // Static member function that
    // * Figures out the minimal workspace requirements, and passes
    //   the results to the user-defined workspace overload of the
    //   invoke static member function
    // * Enables the unblocked algorithm (BLAS level 2)
    //
    template< typename MatrixA, typename VectorSVA, typename MatrixU,
            typename MatrixV >
    static std::ptrdiff_t invoke( const char joba, const char jobu,
            const char jobv, const char jobr, const char jobt,
            const char jobp, MatrixA& a, VectorSVA& sva, MatrixU& u,
            MatrixV& v, minimal_workspace ) {
        namespace bindings = ::boost::numeric::bindings;
        bindings::detail::array< real_type > tmp_work( min_size_work( joba,
                jobu, jobv, bindings::size_row(a),
                bindings::size_column(a) ) );
        bindings::detail::array< fortran_int_t > tmp_iwork(
                min_size_iwork( bindings::size_row(a),
                bindings::size_column(a) ) );
        return invoke( joba, jobu, jobv, jobr, jobt, jobp, a, sva, u, v,
                workspace( tmp_work, tmp_iwork ) );
    }

    //
    // Static member function that
    // * Figures out the optimal workspace requirements, and passes
    //   the results to the user-defined workspace overload of the
    //   invoke static member
    // * Enables the blocked algorithm (BLAS level 3)
    //
    template< typename MatrixA, typename VectorSVA, typename MatrixU,
            typename MatrixV >
    static std::ptrdiff_t invoke( const char joba, const char jobu,
            const char jobv, const char jobr, const char jobt,
            const char jobp, MatrixA& a, VectorSVA& sva, MatrixU& u,
            MatrixV& v, optimal_workspace ) {
        namespace bindings = ::boost::numeric::bindings;
        return invoke( joba, jobu, jobv, jobr, jobt, jobp, a, sva, u, v,
                minimal_workspace() );
    }

    //
    // Static member function that returns the minimum size of
    // workspace-array work.
    //
    static std::ptrdiff_t min_size_work( const char joba, const char jobu,
            const char jobv, const std::ptrdiff_t m, const std::ptrdiff_t n ) {
        if ( jobu == 'N' && jobv == 'N' ) {
            if ( joba != 'E' && joba != 'G' )
                return std::max< std::ptrdiff_t >( std::max<
                        std::ptrdiff_t >( 2*m+n, 4*n+1), 7 );
            else
                return std::max< std::ptrdiff_t >( std::max<
                        std::ptrdiff_t >( 2*m+n, n*n+4*n), 7 );
        } else if ( jobu == 'N' || jobu == 'W' || jobv == 'N' ||
                jobv == 'W' ) {
                return std::max< std::ptrdiff_t >( 2*n+m, 7);
        } else {
            if ( jobv != 'J' )
                return 6*n+2*n*n;
            else
                return std::max< std::ptrdiff_t >( m+3*n+n*n, 7);
        }
    }

    //
    // Static member function that returns the minimum size of
    // workspace-array iwork.
    //
    static std::ptrdiff_t min_size_iwork( const std::ptrdiff_t m,
            const std::ptrdiff_t n ) {
        return m+3*n;
    }
};

//
// Functions for direct use. These functions are overloaded for temporaries,
// so that wrapped types can still be passed and used for write-access. In
// addition, if applicable, they are overloaded for user-defined workspaces.
// Calls to these functions are passed to the gejsv_impl classes. In the
// documentation, most overloads are collapsed to avoid a large number of
// prototypes which are very similar.
//

//
// Overloaded function for gejsv. Its overload differs for
// * User-defined workspace
//
template< typename MatrixA, typename VectorSVA, typename MatrixU,
        typename MatrixV, typename Workspace >
inline typename boost::enable_if< detail::is_workspace< Workspace >,
        std::ptrdiff_t >::type
gejsv( const char joba, const char jobu, const char jobv,
        const char jobr, const char jobt, const char jobp, MatrixA& a,
        VectorSVA& sva, MatrixU& u, MatrixV& v, Workspace work ) {
    return gejsv_impl< typename bindings::value_type<
            MatrixA >::type >::invoke( joba, jobu, jobv, jobr, jobt, jobp, a,
            sva, u, v, work );
}

//
// Overloaded function for gejsv. Its overload differs for
// * Default workspace-type (optimal)
//
template< typename MatrixA, typename VectorSVA, typename MatrixU,
        typename MatrixV >
inline typename boost::disable_if< detail::is_workspace< MatrixV >,
        std::ptrdiff_t >::type
gejsv( const char joba, const char jobu, const char jobv,
        const char jobr, const char jobt, const char jobp, MatrixA& a,
        VectorSVA& sva, MatrixU& u, MatrixV& v ) {
    return gejsv_impl< typename bindings::value_type<
            MatrixA >::type >::invoke( joba, jobu, jobv, jobr, jobt, jobp, a,
            sva, u, v, optimal_workspace() );
}

} // namespace lapack
} // namespace bindings
} // namespace numeric
} // namespace boost

#endif

max_stars_repo_name: ljktest/siconos

max_issues_repo_name: ljktest/siconos

max_forks_repo_name: ljktest/siconos


text:
#include "Surface.h"
#include "Ray.h"
#include "Light.h"
#include "SurfaceList.h"
#include "Intersection.h"
#include "Scene.h"
#include "KD.h"

#include <algorithm>

#include <Eigen/Dense>

void calculate_Ld_Ls(const Ray& ray, const Ray& lightray, float I,
                     const Eigen::Vector3f& n, const Material& mat,
                     Color& Ld, Color& Ls) {
    auto v = -ray.d.normalized();
    auto l = lightray.d.normalized();
    auto h = (v+l).normalized();

    Ld = mat.kd * I * std::max(0.0f, n.dot(l));
    Ls = mat.ks * I * std::pow(std::max(0.0f, n.dot(h)), mat.sp);
}

Color Surface::shade(const Ray& ray, const Eigen::Vector3f& point, const Eigen::Vector3f& n,
                     const Light& light, const Scene& scene, int depth, bool shadows, bool kd) {
    Color La, Ld, Ls, Lm;

    La = mat.ka * light.intensity;

    float I = light.intensity;

    auto corrected_point = point + n*0.01f;
    auto lightray = Ray(corrected_point, light.pos - corrected_point);

    // Apply diffuse and specular shading, with shadows, if applicable
    if (shadows) {
        std::unique_ptr<Intersection> lighthit;
        if (kd)
            lighthit = std::unique_ptr<Intersection>(kd_intersect(ray, kd_tree, 0, mat).intersect(lightray));
        else
            lighthit = std::unique_ptr<Intersection>(scene.intersect(lightray));

        if (!lighthit)
            calculate_Ld_Ls(ray, lightray, I, n, mat, Ld, Ls);
    } else {
        calculate_Ld_Ls(ray, lightray, I, n, mat, Ld, Ls);
    }

    // Apply recursive raytracing, if applicable
    if (mat.a != 0.0f && depth != 0) {
        auto v = -ray.d.normalized();
        auto r = -v - 2*((-v).dot(n)*n);

        auto reflected_ray = Ray(corrected_point, r);
        auto hit = scene.intersect(reflected_ray);
        if (hit)
            Lm = scene.shade(reflected_ray, *hit, depth-1, shadows, kd);
    }

    return (La + Ld + Ls)*(1-mat.a) + (Lm)*(mat.a);
}

max_stars_repo_name: fmenozzi/raytracer

max_issues_repo_name: fmenozzi/raytracer

max_forks_repo_name: fmenozzi/raytracer


text:
#include "TcpConnecting.hpp"

#include <faucet/tcp/TcpSocket.hpp>

#include <boost/lexical_cast.hpp>
#include <boost/bind.hpp>
#include <boost/thread/locks.hpp>

#include <faucet/resolve.hpp>

using namespace boost::asio::ip;

TcpConnecting::TcpConnecting(TcpSocket &socket) :
	ConnectionState(socket),
	resolver(Asio::getIoService()),
	abortRequested(false),
	noDelay_(TcpSocket::DEFAULT_TCP_NODELAY) {
}

void TcpConnecting::enter(const char *host, uint16_t port) {
    fct_async_resolve<tcp>(host, port, resolver, boost::bind(&TcpConnecting::handleResolve,
			this, socket->shared_from_this(), boost::asio::placeholders::error,
			boost::asio::placeholders::iterator));
}

void TcpConnecting::abort() {
	resolver.cancel();
	abortRequested = true;
}

bool TcpConnecting::setNoDelay(bool noDelay) {
    noDelay_ = noDelay;
    return true;
}

void TcpConnecting::handleResolve(std::shared_ptr<TcpSocket> socket,
		const boost::system::error_code &error,
		tcp::resolver::iterator endpointIterator) {
	boost::lock_guard<boost::recursive_mutex> guard(getCommonMutex());

	if (abortRequested) {
		return;
	}

	if (error) {
		enterErrorState(error.message());
		return;
	}

	V4FirstIterator<tcp> endpoints(endpointIterator);
	boost::system::error_code ec;
	if(startConnectionAttempt(socket, endpoints, ec)) {
		enterErrorState(ec.message());
		return;
	}
}

boost::system::error_code TcpConnecting::startConnectionAttempt(std::shared_ptr<TcpSocket> socket,
		V4FirstIterator<tcp> endpoints, boost::system::error_code &ec) {
	if (getSocket().close(ec)) {
		return ec;
	}
	tcp::endpoint endpoint = endpoints.next();
	getSocket().async_connect(endpoint, boost::bind(
			&TcpConnecting::handleConnect, this, socket,
			boost::asio::placeholders::error, endpoints));
	return ec;
}

void TcpConnecting::handleConnect(std::shared_ptr<TcpSocket> socket,
		const boost::system::error_code &error,
		V4FirstIterator<tcp> endpoints) {
	boost::lock_guard<boost::recursive_mutex> guard(getCommonMutex());

	if (abortRequested) {
		return;
	}

	if (!error) {
		enterConnectedState(noDelay_);
	} else if (endpoints.hasNext()) {
		boost::system::error_code ec;
		if(startConnectionAttempt(socket, endpoints, ec)) {
			enterErrorState(ec.message());
			return;
		}
	} else {
		enterErrorState(error.message());
		return;
	}
}


text:
#include <iostream>
#include <algorithm>
#include <iterator>
#include <cmath>
#include <fstream>
#include <valarray>
#include <sstream>

#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/math/constants/constants.hpp>

#include "vlasovpp/field.h"
#include "vlasovpp/weno.h"
#include "vlasovpp/fft.h"
#include "vlasovpp/array_view.h"
#include "vlasovpp/poisson.h"

#ifndef SIGMA
#define SIGMA (2.*std::sqrt(2.))
#endif
#define E_MAX (0.6)

/*
template <unsigned int i>
struct phi
{
  static std::complex<double>
  operator () ( std::complex<double> const & z ) {
    static std::valarray<std::complex<double>> coeff(i);
    coeff[0] = 1.;

    for ( unsigned int k=1 ; k<coeff.size() ; ++k ) {
      coeff[k] = coeff[k-1] * z / (double(k));
    }

    return (std::exp(z) - std::accumulate( std::begin(coeff) , std::end(coeff) , std::complex<double>(0.,0.) ))/(std::pow(z,i));
  }
};
*/
template <unsigned int i>
std::complex<double>
phi ( std::complex<double> const & _z )
{
  std::valarray<std::complex<double>> coeff(i);
  coeff[0] = 1.;

  std::complex<double> z = _z;
  if ( _z == 0. ) { z = std::complex<double>(1.,0.); }

  for ( unsigned int k=1 ; k<coeff.size() ; ++k ) {
    coeff[k] = coeff[k-1] * z / (double(k));
  }
  //std::copy(std::begin(coeff),std::end(coeff),std::ostream_iterator<std::complex<double>>(std::cout," . "));
  //std::cout << std::endl;

  if ( z != 0. ) {
    return (std::exp(z) - std::accumulate( std::begin(coeff) , std::end(coeff) , std::complex<double>(0.,0.) ))/(std::pow(z,i));
  }
  return coeff[i-1];
}

namespace o2 {
  template < typename _T , std::size_t NumDimsV >
  auto
  trp_v ( field<_T,NumDimsV> const & u , ublas::vector<_T> const& E )
  {
    field<_T,NumDimsV> trp(tools::array_view<const std::size_t>(u.shape(),NumDimsV+1));

    { auto k=0, km1=trp.size(0)-1;
      for ( auto i=0 ; i<trp.size(1) ; ++i ) {
        trp[k][i] = ( E(i)*(u[k+1][i]-u[km1][i])/(2.*u.step.dv) );
      }
    }
    for ( auto k=1 ; k<trp.size(0)-1 ; ++k ) {
      for ( auto i=0 ; i<trp.size(1) ; ++i ) {
        trp[k][i] = ( E(i)*(u[k+1][i]-u[k-1][i])/(2.*u.step.dv) );
      }
    }
    { auto k=trp.size(0)-1, kp1=0;
      for ( auto i=0 ; i<trp.size(1) ; ++i ) {
        trp[k][i] = ( E(i)*(u[kp1][i]-u[k-1][i])/(2.*u.step.dv) );
      }
    }

    return trp;
  }
}

namespace math = boost::math::constants;
const std::complex<double> & I = std::complex<double>(0.,1.);

#define SQ(X) ((X)*(X))
#define Xi(i) (i*f.step.dx+f.range.x_min)
#define Vk(k) (k*f.step.dv+f.range.v_min)

int main(int,char**)
{
	std::size_t Nx = 135, Nv = 256 , Nb_iter=10;
	field<double,1> f(boost::extents[Nv][Nx]);

	f.range.v_min = -8.; f.range.v_max = 8.;
	f.step.dv = (f.range.v_max-f.range.v_min)/Nv;

  const double Kx = 0.3;
  f.range.x_min = 0.; f.range.x_max = 2./Kx*math::pi<double>();
	//f.range.x_min = 0.; f.range.x_max = 20.*math::pi<double>();
  //f.range.x_min = 0.; f.range.x_max = 4.0*math::pi<double>();
  //f.range.x_min = -8.; f.range.x_max = 8.;
	f.step.dx = (f.range.x_max-f.range.x_min)/Nx;

  // SIGMA is the CFL number 0.45 is E_max in our test case
	const double dt = 0.05;//SIGMA*f.step.dv/E_MAX; //1.606*f.step.dv/0.6; //0.005; //1.606*f.step.dv/0.6; //0.5*6.*math::pi<double>()/(Nv*f.range.v_max);

  //field<double,1> f_sol = f;
  //field<double,1> f_ini = f;

	ublas::vector<double> v (Nv,0.);
  ublas::vector<double> E (Nx,0.),rho(Nx);
  for ( std::size_t k=0 ; k<Nv ; ++k ) { v[k] = Vk(k); }
  //for ( std::size_t i=0 ; i<Nx ; ++i ) { E[i] = -Xi(i); }

	const double l = f.range.x_max-f.range.x_min;
	ublas::vector<double> kx(Nx);
	//for ( auto i=0 ; i<Nx/2+1 ; ++i )   { kx[i] = 2.*math::pi<double>()*i/l; }
	//for ( auto i=0 ; i<((Nx/2)) ; ++i ) { kx[i+Nx/2+1] = -kx[Nx/2-i]; }
  for ( auto i=0 ; i<Nx/2 ; ++i ) { kx[i]    = 2.*math::pi<double>()*i/l; }
  for ( int i=-Nx/2 ; i<0 ; ++i ) { kx[Nx+i] = 2.*math::pi<double>()*i/l; }

  double ui = 3.4;
  double alpha = 0.1;
  double Tc = 0.0001;

  //double np = 0.9 , nb = 0.2 , ui = 4.5;
  for (field<double,2>::size_type k=0 ; k<f.size(0) ; ++k ) {
    for (field<double,2>::size_type i=0 ; i<f.size(1) ; ++i ) {
      //f[k][i] = ( std::exp(-0.5*SQ(Vk(k)))*np/std::sqrt(2.*math::pi<double>()) + nb/std::sqrt(2.*math::pi<double>())*std::exp(-0.5*SQ(Vk(k)-ui)/0.25) )*(1.+0.04*std::cos(0.5*Xi(i)));

      // Bump on Tail
      f[k][i] = ( std::exp(-0.5*SQ(Vk(k)))*np/std::sqrt(2.*math::pi<double>()) + nb/std::sqrt(2.*math::pi<double>())*std::exp(-0.5*SQ(Vk(k)-ui)/0.25) )*(1.+0.04*std::cos(0.3*Xi(i)));
      // Landau dumpping test
      //f[k][i] = (1./std::sqrt(2*math::pi<double>()))*std::exp(-0.5*SQ(Vk(k)))*(1.+0.001*std::cos(0.5*Xi(i)));

      //f[k][i] = std::exp( -SQ(Xi(i)-6) );
      //f[k][i] = std::exp(-SQ(Xi(i)-3)/0.5 - SQ(Vk(k))/2.);
      //f_sol[k][i] = std::exp(-SQ(Xi(i)-3)/0.5 - SQ(Vk(k))/2.);
      //f[k][i] = SQ(v(k))*std::exp(-0.5*SQ(Vk(k)))/std::sqrt(2.*math::pi<double>())*(1.+0.01*std::cos(0.5*Xi(i)));
      //f_ini[k][i] = f[k][i];
      //f_sol[k][i] =  std::exp(-SQ(Xi(i)-3.-v[k]*Nb_iter*dt)/0.5 - SQ(Vk(k)-E[i]*Nb_iter*dt)/2.);

      // triple bump
      //f[k][i] = ( 0.5*alpha/std::sqrt(2.*math::pi<double>())*std::exp(-0.5*SQ(Vk(k)-ui)) + 0.5*alpha/std::sqrt(2.*math::pi<double>())*std::exp(-0.5*SQ(Vk(k)+ui)) )*(1.+0.04*std::cos(Kx*Xi(i))) + ((1-alpha)/(std::sqrt(2.*math::pi<double>()*Tc))*std::exp(-0.5*SQ(Vk(k))/Tc));
    }
  }
  f.write("vphl/kin/init.dat");

  poisson<double> poisson_solver(Nx,l);
  rho = f.density();
  E = poisson_solver(rho);

  //double Tf = 60.;//2*math::pi<double>();
  double Tf = 40.;
  int i_t=0;

  std::cout << "Nx: " << Nx << "\n";
  std::cout << "Nv: " << Nv << "\n";
  std::cout << "v_min: " << f.range.v_min << "\n";
  std::cout << "v_max: " << f.range.v_max << "\n";
  std::cout << "x_min: " << f.range.x_min << "\n";
  std::cout << "x_max: " << f.range.x_max << "\n";
  std::cout << "dt: " << dt << "\n";
  std::cout << "dx: " << f.step.dx << "\n";
  std::cout << "dv: " << f.step.dv << "\n";
  std::cout << "Tf: " << Tf << "\n";
  std::cout << "f_0: " << "\"bot\"" << "\n";
  std::cout << std::endl;

  std::ofstream info("info.yaml");

  info << "Nx: " << Nx << "\n";
  info << "Nv: " << Nv << "\n";
  info << "v_min: " << f.range.v_min << "\n";
  info << "v_max: " << f.range.v_max << "\n";
  info << "x_min: " << f.range.x_min << "\n";
  info << "x_max: " << f.range.x_max << "\n";
  info << "dt: " << dt << "\n";
  info << "dx: " << f.step.dx << "\n";
  info << "dv: " << f.step.dv << "\n";
  info << "Tf: " << Tf << "\n";
  info << "f_0: " << "\"bot\"" << "\n";
  info << std::endl;
  info.close();

  ublas::vector<double> ee(int(std::ceil(Tf/dt))+1,0.);
  ublas::vector<double> Emax(int(std::ceil(Tf/dt))+1,0.);
  ublas::vector<double> H(int(std::ceil(Tf/dt))+1,0.);
  //ublas::vector<double> ee(Nb_iter);
  //ublas::vector<double> Emax(Nb_iter);
  //ublas::vector<double> H(Nb_iter);

  field<double,1> f1(tools::array_view<const std::size_t>(f.shape(),2)),f2(tools::array_view<const std::size_t>(f.shape(),2)),f3(tools::array_view<const std::size_t>(f.shape(),2)),f4(tools::array_view<const std::size_t>(f.shape(),2)),f5(tools::array_view<const std::size_t>(f.shape(),2));
  fft::spectrum_ hf(Nx),hf1(Nx),hf2(Nx),hf3(Nx),hf4(Nx),hf5(Nx),hEdvf(Nx),hEdvf1(Nx),hEdvf2(Nx),hEdvf3(Nx),hEdvf4(Nx);

  rho = f.density();
  E = poisson_solver(rho);
  //for ( int i=0 ; i<E.size() ; ++i ) { E[i] = 1.; }
  //for ( int i=0 ; i<E.size() ; ++i ) { std::cout << E[i] << " "; }
  //  std::cout << std::endl;

  //ee(i_t) = 0.;
  //for ( auto i=0 ; i<Nx ; ++i ) {
  //  ee(i_t) += SQ(E(i))*f.step.dx;
  //}
  //ee(i_t) = std::sqrt(ee(i_t));
  H(i_t) = energy(f,E);

#define L (-v(k)*I*kx[i])
//#define L (0.)
  while (  i_t*dt < Tf ) {
  //while (  i_t < Nb_iter ) {
  	//if (t%32==0) { std::cout<<"\r"<<t<<" "<<std::flush ; }
    std::cout<<" ["<<std::setw(5)<<i_t<<"] "<<i_t*dt<<"\r"<<std::flush;
    Emax(i_t) = std::abs(*std::max_element( E.begin() , E.end() , [](double a,double b){return (std::abs(a) < std::abs(b));} ));
    ee(i_t) = 0.;
    for ( auto i=0 ; i<Nx ; ++i ) { ee(i_t) += SQ(E(i))*f.step.dx; }
    ee(i_t) = std::sqrt(ee(i_t));

    /**
    // exprk(2,2) =============================================================
    #define SCHEME "expRK22"
    // SIGMA = 0.551 (10^-2)
    E = poisson_solver(f.density());
    field<double,1> Edvf = o2::trp_v(f,E);
    for ( auto k=0 ; k<f.size(0) ; ++k ) {
      hf.fft(&(f[k][0]));
      hEdvf.fft(&(Edvf[k][0]));

      for ( auto i=0 ; i<Nx ; ++i ) {
        hf1[i] = std::exp(-L*dt)*hf[i] + dt*phi<1>(dt*L)*hEdvf[i];
      }
      hf1.ifft(&(f1[k][0]));
    }

    E = poisson_solver(f1.density());
    field<double,1> Edvf1 = o2::trp_v(f1,E);
    for ( auto k=0 ; k<f.size(0) ; ++k ) {
      hf.fft(&f[k][0]);
      hf1.fft(&f1[k][0]);
      hEdvf.fft(&(Edvf[k][0]));
      hEdvf1.fft(&(Edvf1[k][0]));

      for ( auto i=0 ; i<Nx ; ++i ) {
        hf[i] = std::exp(-L*dt)*hf[i] + dt*( (phi<1>(dt*L)-phi<2>(dt*L))*hEdvf[i] + phi<2>(dt*L)*hEdvf1[i] );
      }

      hf.ifft(&(f[k][0]));
    }

    **/
    /**
    // Cox-Matthews ===========================================================
    #define SCHEME "CM"
    // SIGMA = 0.450 (10^-2)
    E = poisson_solver( f.density() );
    field<double,1> Edvf = o2::trp_v(f,E);
    for ( auto k=0 ; k<f.size(0) ; ++k ) {
      hf.fft(&(f[k][0]));
      hEdvf.fft(&(Edvf[k][0]));

      for ( auto i=0 ; i<Nx ; ++i ) {
        hf1[i] = std::exp(-0.5*L*dt)*hf[i] + 0.5*dt*phi<1>(-0.5*L*dt)*hEdvf[i];
      }

      hf1.ifft(&(f1[k][0]));
    }

    E = poisson_solver( f1.density() );
    field<double,1> Edvf1 = o2::trp_v(f1,E);
    for ( auto k=0 ; k<f.size(0) ; ++k ) {
      hf.fft(&(f[k][0]));
      hf1.fft(&(f1[k][0]));
      hEdvf.fft(&(Edvf[k][0]));
      hEdvf1.fft(&(Edvf1[k][0]));

      for ( auto i=0 ; i<Nx ; ++i ) {
        hf2[i] = std::exp(-0.5*L*dt)*hf[i] + 0.5*dt*phi<1>(-0.5*L*dt)*hEdvf1[i];
      }

      hf2.ifft(&(f2[k][0]));
    }

    E = poisson_solver( f2.density() );
    field<double,1> Edvf2 = o2::trp_v(f2,E);
    for ( auto k=0 ; k<f.size(0) ; ++k ) {
      hf.fft(&(f[k][0]));
      hf1.fft(&(f1[k][0]));
      hf2.fft(&(f2[k][0]));
      hEdvf.fft(&(Edvf[k][0]));
      hEdvf1.fft(&(Edvf1[k][0]));
      hEdvf2.fft(&(Edvf2[k][0]));

      for ( auto i=0 ; i<Nx ; ++i ) {
        hf3[i] = std::exp(-L*dt)*hf[i] + 0.5*dt*phi<1>(-0.5*L*dt)*(std::exp(-0.5*L*dt)-1.)*hEdvf[i] + dt*phi<1>(-0.5*L*dt)*hEdvf2[i];
      }

      hf3.ifft(&(f3[k][0]));
    }

    E = poisson_solver( f3.density() );
    field<double,1> Edvf3  = o2::trp_v(f3,E);
    for ( auto k=0 ; k<f.size(0) ; ++k ) {
      hf.fft(&(f[k][0]));
      hf1.fft(&(f1[k][0]));
      hf2.fft(&(f2[k][0]));
      hf3.fft(&(f3[k][0]));
      hEdvf.fft(&(Edvf[k][0]));
      hEdvf1.fft(&(Edvf1[k][0]));
      hEdvf2.fft(&(Edvf2[k][0]));
      hEdvf3.fft(&(Edvf3[k][0]));

      for ( auto i=0 ; i<Nx ; ++i ) {
        hf4[i] = std::exp(-L*dt)*hf[i] + dt*( (phi<1>(-L*dt)-3.*phi<2>(-L*dt)+4.*phi<3>(-L*dt))*hEdvf[i]
                                             + (2.*phi<2>(-L*dt)-4.*phi<3>(-L*dt))*(hEdvf1[i]+hEdvf2[i])
                                             + (-phi<2>(-L*dt)+4.*phi<3>(-L*dt))*hEdvf3[i] );
      }

      hf4.ifft(&(f[k][0]));
    }

    **/
    /**
    // Krogstad ===============================================================
    #define SCHEME "K"
    // SIGMA = 0.200 (10^-2)
    E = poisson_solver( f.density() );
    field<double,1> Edvf = o2::trp_v(f,E);
    for ( auto k=0 ; k<f.size(0) ; ++k ) {
      hf.fft(&(f[k][0]));
      hEdvf.fft(&(Edvf[k][0]));

      for ( auto i=0 ; i<Nx ; ++i ) {
        hf1[i] = std::exp(-0.5*L*dt)*hf[i] + 0.5*dt*phi<1>(-0.5*L*dt)*hEdvf[i];
      }

      hf1.ifft(&(f1[k][0]));
    }

    E = poisson_solver( f1.density() );
    field<double,1> Edvf1 = o2::trp_v(f1,E);
    for ( auto k=0 ; k<f.size(0) ; ++k ) {
      hf.fft(&(f[k][0]));
      hf1.fft(&(f1[k][0]));
      hEdvf.fft(&(Edvf[k][0]));
      hEdvf1.fft(&(Edvf1[k][0]));

      for ( auto i=0 ; i<Nx ; ++i ) {
        hf2[i] = std::exp(-0.5*L*dt)*hf[i] + dt*(0.5*phi<1>(-0.5*L*dt)-phi<2>(-0.5*L*dt))*hEdvf[i] + dt*phi<2>(-0.5*L*dt)*hEdvf1[i];
      }

      hf2.ifft(&(f2[k][0]));
    }

    E = poisson_solver( f2.density() );
    field<double,1> Edvf2 = o2::trp_v(f2,E);
    for ( auto k=0 ; k<f.size(0) ; ++k ) {
      hf.fft(&(f[k][0]));
      hf1.fft(&(f1[k][0]));
      hf2.fft(&(f2[k][0]));
      hEdvf.fft(&(Edvf[k][0]));
      hEdvf1.fft(&(Edvf1[k][0]));
      hEdvf2.fft(&(Edvf2[k][0]));

      for ( auto i=0 ; i<Nx ; ++i ) {
        hf3[i] = std::exp(-L*dt)*hf[i] + dt*(phi<1>(-L*dt)-2.*phi<2>(-0.5*L*dt))*hEdvf[i] + 2.*dt*phi<2>(-L*dt)*hEdvf2[i];
      }

      hf3.ifft(&(f3[k][0]));
    }

    E = poisson_solver( f3.density() );
    field<double,1> Edvf3  = o2::trp_v(f3,E);
    for ( auto k=0 ; k<f.size(0) ; ++k ) {
      hf.fft(&(f[k][0]));
      hf1.fft(&(f1[k][0]));
      hf2.fft(&(f2[k][0]));
      hf3.fft(&(f3[k][0]));
      hEdvf.fft(&(Edvf[k][0]));
      hEdvf1.fft(&(Edvf1[k][0]));
      hEdvf2.fft(&(Edvf2[k][0]));
      hEdvf3.fft(&(Edvf3[k][0]));

      for ( auto i=0 ; i<Nx ; ++i ) {
        hf4[i] = std::exp(-L*dt)*hf[i] + dt*( (phi<1>(-L*dt)-3.*phi<2>(-L*dt)+4.*phi<3>(-L*dt))*hEdvf[i]
                                             + (2.*phi<2>(-L*dt)-4.*phi<3>(-L*dt))*(hEdvf1[i]+hEdvf2[i])
                                             + (-phi<2>(-L*dt)+4.*phi<3>(-L*dt))*hEdvf3[i] );
      }

      hf4.ifft(&(f[k][0]));
    }

    **/
    /**
    // Hochbruck-Ostermann ====================================================
    #define SCHEME "HO"
    // SIGMA = 0.501 (10^-2)
    E = poisson_solver( f.density() );
    field<double,1> Edvf = o2::trp_v(f,E);
    for ( auto k=0 ; k<f.size(0) ; ++k ) {
      hf.fft(&(f[k][0]));
      hEdvf.fft(&(Edvf[k][0]));

      for ( auto i=0 ; i<Nx ; ++i ) {
        hf1[i] = std::exp(-0.5*L*dt)*hf[i] + 0.5*dt*phi<1>(-0.5*L*dt)*hEdvf[i];
      }

      hf1.ifft(&(f1[k][0]));
    }

    E = poisson_solver( f1.density() );
    field<double,1> Edvf1 = o2::trp_v(f1,E);
    for ( auto k=0 ; k<f.size(0) ; ++k ) {
      hf.fft(&(f[k][0]));
      hf1.fft(&(f1[k][0]));
      hEdvf.fft(&(Edvf[k][0]));
      hEdvf1.fft(&(Edvf1[k][0]));

      for ( auto i=0 ; i<Nx ; ++i ) {
        hf2[i] = std::exp(-0.5*L*dt)*hf[i] + dt*(0.5*phi<1>(-0.5*L*dt)-phi<2>(-0.5*L*dt))*hEdvf[i] + dt*phi<2>(-0.5*L*dt)*hEdvf1[i];
      }

      hf2.ifft(&(f2[k][0]));
    }

    E = poisson_solver( f2.density() );
    field<double,1> Edvf2 = o2::trp_v(f2,E);
    for ( auto k=0 ; k<f.size(0) ; ++k ) {
      hf.fft(&(f[k][0]));
      hf1.fft(&(f1[k][0]));
      hf2.fft(&(f2[k][0]));
      hEdvf.fft(&(Edvf[k][0]));
      hEdvf1.fft(&(Edvf1[k][0]));
      hEdvf2.fft(&(Edvf2[k][0]));

      for ( auto i=0 ; i<Nx ; ++i ) {
        hf3[i] = std::exp(-L*dt)*hf[i] + dt*(phi<1>(-L*dt)-2.*phi<2>(-L*dt))*hEdvf[i] + dt*phi<2>(-L*dt)*hEdvf1[i] + dt*phi<2>(-L*dt)*hEdvf2[i];
      }

      hf3.ifft(&(f3[k][0]));
    }

    E = poisson_solver( f3.density() );
    field<double,1> Edvf3  = o2::trp_v(f3,E);
    for ( auto k=0 ; k<f.size(0) ; ++k ) {
      hf.fft(&(f[k][0]));
      hf1.fft(&(f1[k][0]));
      hf2.fft(&(f2[k][0]));
      hf3.fft(&(f3[k][0]));
      hEdvf.fft(&(Edvf[k][0]));
      hEdvf1.fft(&(Edvf1[k][0]));
      hEdvf2.fft(&(Edvf2[k][0]));
      hEdvf3.fft(&(Edvf3[k][0]));

#define a52 (0.5*phi<2>(-0.5*L*dt)-phi<3>(-L*dt)+0.25*phi<2>(-L*dt)-0.5*phi<3>(-0.5*L*dt))
#define a54 (0.25*phi<2>(-0.5*L*dt)-a52)
      for ( auto i=0 ; i<Nx ; ++i ) {
        hf4[i] = std::exp(-0.5*L*dt)*hf[i] + dt*(0.5*phi<1>(-0.5*L*dt)-2.*a52-a54)*hEdvf[i] + dt*a52*(hEdvf1[i]+hEdvf2[i]) + dt*(0.25*phi<2>(-0.5*L*dt)-a52)*hEdvf3[i];
      }
#undef a54
#undef a52

      hf4.ifft(&(f4[k][0]));
    }

    E = poisson_solver( f4.density() );
    field<double,1> Edvf4  = o2::trp_v(f4,E);
    for ( auto k=0 ; k<f.size(0) ; ++k ) {
      hf.fft(&(f[k][0]));
      hf1.fft(&(f1[k][0]));
      hf2.fft(&(f2[k][0]));
      hf3.fft(&(f3[k][0]));
      hf4.fft(&(f4[k][0]));
      hEdvf.fft(&(Edvf[k][0]));
      hEdvf1.fft(&(Edvf1[k][0]));
      hEdvf2.fft(&(Edvf2[k][0]));
      hEdvf3.fft(&(Edvf3[k][0]));
      hEdvf4.fft(&(Edvf4[k][0]));

      for ( auto i=0 ; i<Nx ; ++i ) {
        hf5[i] = std::exp(-L*dt)*hf[i] + dt*( (phi<1>(-L*dt)-3.*phi<2>(-L*dt)+4.*phi<3>(-L*dt))*hEdvf[i]
                                             + (-phi<2>(-L*dt)+4.*phi<3>(-L*dt))*hEdvf3[i]
                                             + (4.*phi<2>(-L*dt)-8.*phi<3>(-L*dt))*hEdvf4[i] );
      }

      hf5.ifft(&(f[k][0]));
    }
    **/
    /**
    // RK(3,2) best ===========================================================
    #define SCHEME "RK32"
    // SIGMA = 2. (y_max)
    // SIGMA = 1.344 (WENO)
    E = poisson_solver(f.density());
    field<double,1> Edvf = o2::trp_v(f,E);
    for ( auto k=0 ; k<f.size(0) ; ++k ) {
      hf.fft(&(f[k][0]));
      hEdvf.fft(&(Edvf[k][0]));

      for ( auto i=0 ; i<Nx ; ++i ) {
        hf1[i] = std::exp(0.5*L*dt)*( hf[i]-0.5*dt*hEdvf[i] );
      }
      hf1.ifft(&(f1[k][0]));
    }

    E = poisson_solver(f1.density());
    Edvf = o2::trp_v(f1,E);
    for ( auto k=0 ; k<f.size(0) ; ++k ) {
      hf.fft(&(f[k][0]));
      //hf1.fft(&(f1[k][0]));
      hEdvf.fft(&(Edvf[k][0]));

      for ( auto i=0 ; i<Nx ; ++i ) {
        hf2[i] = std::exp(0.5*L*dt)*hf[i] - 0.5*dt*hEdvf[i];
      }
      hf2.ifft(&(f2[k][0]));
    }

    E = poisson_solver(f2.density());
    Edvf = o2::trp_v(f2,E);
    for ( auto k=0 ; k<f.size(0) ; ++k ) {
      hf.fft(&(f[k][0]));
      //hf1.fft(&(f1[k][0]));
      //hf2.fft(&(f2[k][0]));
      hEdvf.fft(&(Edvf[k][0]));

      for ( auto i=0 ; i<Nx ; ++i ) {
        hf[i] = std::exp(L*dt)*hf[i] - dt*std::exp(0.5*L*dt)*hEdvf[i];
      }
      hf.ifft(&(f[k][0]));
    }
    **/
    /**/
    // RK(3,3) ================================================================
    // RK(3,3) eq19 ===========================================================
    #define SCHEME "RK33"
    //#define SCHEME "RK33_eq19"
    // SIGMA = std::sqrt(3) (y_max)
    // SIGMA = 1.433 (WENO)
    E = poisson_solver(f.density());

    //std::cout << "\n" << L << " " << *std::max_element(E.begin(),E.end()) << " " << *std::min_element(E.begin(),E.end()) << std::endl;

    field<double,1> Edvf = weno::trp_v(f,E);
    for ( auto k=0 ; k<f.size(0) ; ++k ) {
      hf.fft(&(f[k][0]));
      hEdvf.fft(&(Edvf[k][0]));

      for ( auto i=0 ; i<Nx ; ++i ) {
        hf1[i] = std::exp(L*dt)*( hf[i]-dt*hEdvf[i] );
        //hf1[i] = 0.5*std::exp((2./3.)*L*dt)*hf[i] + 0.5*std::exp((2./3.)*dt*L)*( hf[i] - (4./3.)*dt*hEdvf[i] );
      }
      hf1.ifft(&(f1[k][0]));
    }

    E = poisson_solver(f1.density());
    Edvf = weno::trp_v(f1,E);
    for ( auto k=0 ; k<f.size(0) ; ++k ) {
      hf.fft(&(f[k][0]));
      hf1.fft(&(f1[k][0]));
      hEdvf.fft(&(Edvf[k][0]));

      for ( auto i=0 ; i<Nx ; ++i ) {
        hf2[i] = 0.75*std::exp(0.5*L*dt)*hf[i] + 0.25*std::exp(-0.5*L*dt)*( hf1[i]-dt*hEdvf[i] );
        //hf2[i] = (2./3.)*std::exp((2./3.)*dt*L)*hf[i] + (1./3.)*( hf1[i] - (4./3.)*dt*hEdvf[i] );
      }
      hf2.ifft(&(f2[k][0]));
    }

    E = poisson_solver(f2.density());
    Edvf = weno::trp_v(f2,E);
    for ( auto k=0 ; k<f.size(0) ; ++k ) {
      hf.fft(&(f[k][0]));
      hf1.fft(&(f1[k][0]));
      hf2.fft(&(f2[k][0]));
      hEdvf.fft(&(Edvf[k][0]));

      for ( auto i=0 ; i<Nx ; ++i ) {
        hf[i] = (1./3.)*std::exp(L*dt)*hf[i] + (2./3.)*std::exp(0.5*L*dt)*( hf2[i]-dt*hEdvf[i] );
        //hf[i] = (59./128.)*std::exp(L*dt)*hf[i] + (15./128.)*std::exp(L*dt)*( 2.*hf1[i]*std::exp(-(2./3.)*L*dt) - hf[i] ) + (27./64.)*std::exp((1./3.)*dt*L)*( hf2[i] - (4./3.)*dt*hEdvf[i] );
      }
      hf.ifft(&(f[k][0]));
    }
    /**/
    /**
    // RK(4,4) ================================================================
    // RK(4,4) 3/8 rule =======================================================
    #define SCHEME "RK44"
    //#define SCHEME "RK44_38"
    // SIGMA = 2.*std::sqrt(2.) (y_max)
    // SIGMA = 1.731 (WENO)
    E = poisson_solver( f.density() );
    field<double,1> Edvf = weno::trp_v(f,E);
    for ( auto k=0 ; k<f.size(0) ; ++k ) {
      hf.fft(&(f[k][0]));
      hEdvf.fft(&(Edvf[k][0]));

      for ( auto i=0 ; i<Nx ; ++i ) {
        hf1[i] = std::exp(0.5*L*dt)*( hf[i] - 0.5*dt*hEdvf[i] );
        //hf1[i] = std::exp((1./3.)*L*dt)*( hf[i] - (1./3.)*dt*hEdvf[i] );
      }

      hf1.ifft(&(f1[k][0]));
    }

    E = poisson_solver( f1.density() );
    Edvf = weno::trp_v(f1,E);
    for ( auto k=0 ; k<f.size(0) ; ++k ) {
      hf.fft(&(f[k][0]));
      hf1.fft(&(f1[k][0]));
      hEdvf.fft(&(Edvf[k][0]));

      for ( auto i=0 ; i<Nx ; ++i ) {
        hf2[i] = std::exp(0.5*L*dt)*hf[i] - 0.5*dt*hEdvf[i] ;
        //hf2[i] = 2.*std::exp((2./3.)*L*dt)*hf[i] - std::exp((1./3.)*L*dt)*hf1[i] - dt*std::exp((1./3.)*L*dt)*hEdvf[i];
      }

      hf2.ifft(&(f2[k][0]));
    }

    E = poisson_solver( f2.density() );
    Edvf = weno::trp_v(f2,E);
    for ( auto k=0 ; k<f.size(0) ; ++k ) {
      hf.fft(&(f[k][0]));
      hf1.fft(&(f1[k][0]));
      hf2.fft(&(f2[k][0]));
      hEdvf.fft(&(Edvf[k][0]));

      for ( auto i=0 ; i<Nx ; ++i ) {
        hf3[i] = std::exp(L*dt)*hf[i] - dt*std::exp(0.5*L*dt)*hEdvf[i];
        //hf3[i] = 2.*std::exp((2./3.)*L*dt)*hf1[i] - std::exp((1./3.)*L*dt)*hf2[i] - dt*std::exp((1./3.)*L*dt)*hEdvf[i];
      }

      hf3.ifft(&(f3[k][0]));
    }

    E = poisson_solver( f3.density() );
    Edvf  = weno::trp_v(f3,E);
    for ( auto k=0 ; k<f.size(0) ; ++k ) {
      hf.fft(&(f[k][0]));
      hf1.fft(&(f1[k][0]));
      hf2.fft(&(f2[k][0]));
      hf3.fft(&(f3[k][0]));
      hEdvf.fft(&(Edvf[k][0]));

      for ( auto i=0 ; i<Nx ; ++i ) {
        hf[i] = -(1./3.)*std::exp(L*dt)*hf[i] + (1./3.)*std::exp(0.5*L*dt)*hf1[i] + (2./3.)*std::exp(0.5*L*dt)*hf2[i] + (1./3.)*hf3[i] - (1./6.)*dt*hEdvf[i];
        //hf[i] = -(1./8.)*std::exp(L*dt)*hf[i] + 0.75*std::exp((1./3.)*L*dt)*hf2[i] + (3./8.)*hf3[i] - (1./8.)*dt*hEdvf[i];
      }

      hf.ifft(&(f[k][0]));
    }
    **/
    /**
    // RK(5,3) ================================================================
    #define SCHEME "RK53"
    E = poisson_solver(f.density());
    field<double,1> Edvf = weno::trp_v(f,E);
    for ( auto k=0 ; k<f.size(0) ; ++k ) {
      hf.fft(&(f[k][0]));
      hEdvf.fft(&(Edvf[k][0]));

      for ( auto i=0 ; i<Nx ; ++i ) {
        hf1[i] = std::exp((1./7.)*L*dt)*hf[i] - (1./7.)*dt*std::exp((1./7.)*L*dt)*hEdvf[i];
      }
      hf1.ifft(&(f1[k][0]));
    }

    E = poisson_solver(f1.density());
    Edvf = weno::trp_v(f1,E);
    for ( auto k=0 ; k<f.size(0) ; ++k ) {
      hf.fft(&(f[k][0]));
      hEdvf.fft(&(Edvf[k][0]));

      for ( auto i=0 ; i<Nx ; ++i ) {
        hf2[i] = std::exp((3./16.)*L*dt)*hf[i] - (3./16.)*dt*std::exp((5./112.)*L*dt)*hEdvf[i];
      }
      hf2.ifft(&(f2[k][0]));
    }

    E = poisson_solver(f2.density());
    Edvf = weno::trp_v(f2,E);
    for ( auto k=0 ; k<f.size(0) ; ++k ) {
      hf.fft(&(f[k][0]));
      hEdvf.fft(&(Edvf[k][0]));

      for ( auto i=0 ; i<Nx ; ++i ) {
        hf3[i] = std::exp((1./3.)*L*dt)*hf[i] - (1./3.)*dt*std::exp((7./48.)*L*dt)*hEdvf[i];
      }
      hf3.ifft(&(f3[k][0]));
    }

    E = poisson_solver(f3.density());
    Edvf = weno::trp_v(f3,E);
    for ( auto k=0 ; k<f.size(0) ; ++k ) {
      hf.fft(&(f[k][0]));
      hEdvf.fft(&(Edvf[k][0]));

      for ( auto i=0 ; i<Nx ; ++i ) {
        hf4[i] = std::exp((2./3.)*L*dt)*hf[i] - (2./3.)*dt*std::exp((1./3.)*L*dt)*hEdvf[i];
      }
      hf4.ifft(&(f4[k][0]));
    }

    E = poisson_solver(f4.density());
    Edvf = weno::trp_v(f4,E);
    for ( auto k=0 ; k<f.size(0) ; ++k ) {
      hf.fft(&(f[k][0]));
      hf1.fft(&(f1[k][0]));
      hEdvf.fft(&(Edvf[k][0]));

      for ( auto i=0 ; i<Nx ; ++i ) {
        hf[i] = -0.75*std::exp(L*dt)*hf[i] + 1.75*std::exp((6./7.)*L*dt)*hf1[i] - 0.75*dt*std::exp((1./3.)*L*dt)*hEdvf[i];
      }
      hf.ifft(&(f[k][0]));
    }
    **/
    /**
    // DP5  ===================================================================
    #define SCHEME "DP5"
    E = poisson_solver(f.density());
    field<double,1> Edvf = weno::trp_v(f,E);
    for ( auto k=0 ; k<f.size(0) ; ++k ) {
      hf.fft(&(f[k][0]));
      hEdvf.fft(&(Edvf[k][0]));

      for ( auto i=0 ; i<Nx ; ++i ) {
        hf1[i] = std::exp((1./5.)*L*dt)*hf[i] - (1./5.)*dt*std::exp((1./5.)*L*dt)*hEdvf[i];
      }
      hf1.ifft(&(f1[k][0]));
    }

    E = poisson_solver(f1.density());
    Edvf = weno::trp_v(f1,E);
    for ( auto k=0 ; k<f.size(0) ; ++k ) {
      hf.fft(&(f[k][0]));
      hf1.fft(&(f1[k][0]));
      hEdvf.fft(&(Edvf[k][0]));

      for ( auto i=0 ; i<Nx ; ++i ) {
        hf2[i] = (5./8.)*std::exp((3./10.)*L*dt)*hf[i] + (3./8.)*std::exp((1./10.)*L*dt)*hf1[i] - (9./40.)*dt*std::exp((1./10.)*L*dt)*hEdvf[i];
      }
      hf2.ifft(&(f2[k][0]));
    }

    E = poisson_solver(f2.density());
    Edvf = weno::trp_v(f2,E);
    for ( auto k=0 ; k<f.size(0) ; ++k ) {
      hf.fft(&(f[k][0]));
      hf1.fft(&(f1[k][0]));
      hf2.fft(&(f2[k][0]));
      hEdvf.fft(&(Edvf[k][0]));

      for ( auto i=0 ; i<Nx ; ++i ) {
        hf3[i] = (175./27.)*std::exp((4./5.)*L*dt)*hf[i] + (100./9.)*std::exp((3./5.)*L*dt)*hf1[i] - (448./27.)*std::exp(0.5*L*dt)*hf2[i] - (32./9.)*dt*std::exp(0.5*L*dt)*hEdvf[i];
      }
      hf3.ifft(&(f3[k][0]));
    }

    E = poisson_solver(f3.density());
    Edvf = weno::trp_v(f3,E);
    for ( auto k=0 ; k<f.size(0) ; ++k ) {
      hf.fft(&(f[k][0]));
      hf1.fft(&(f1[k][0]));
      hf2.fft(&(f2[k][0]));
      hf3.fft(&(f3[k][0]));
      hEdvf.fft(&(Edvf[k][0]));

      for ( auto i=0 ; i<Nx ; ++i ) {
        hf4[i] = (3551./6561.)*std::exp((8./9.)*L*dt)*hf[i] + (7420./2187.)*std::exp((31./45.)*L*dt)*hf1[i] - (37376./6561.)*std::exp((53./90.)*L*dt)*hf2[i] + (2014./729.)*std::exp((4./45.)*L*dt)*hf3[i] + (212./729.)*dt*std::exp((4./45.)*L*dt)*hEdvf[i];
      }
      hf4.ifft(&(f4[k][0]));
    }

    E = poisson_solver(f4.density());
    Edvf = weno::trp_v(f4,E);
    for ( auto k=0 ; k<f.size(0) ; ++k ) {
      hf.fft(&(f[k][0]));
      hf1.fft(&(f1[k][0]));
      hf2.fft(&(f2[k][0]));
      hf3.fft(&(f3[k][0]));
      hf4.fft(&(f4[k][0]));
      hEdvf.fft(&(Edvf[k][0]));

      for ( auto i=0 ; i<Nx ; ++i ) {
        hf5[i] = (313397./335808.)*std::exp(L*dt)*hf[i] + (424025./55968.)*std::exp((4./5.)*L*dt)*hf1[i] - (61400./5247.)*std::exp((7./10.)*L*dt)*hf2[i] + (96075./18656.)*std::exp((1./5.)*L*dt)*hf3[i] - (35721./37312.)*std::exp((1./9.)*L*dt)*hf4[i] + (5103./18656.)*dt*std::exp((1./9.)*L*dt)*hEdvf[i];
      }
      hf5.ifft(&(f5[k][0]));
    }

    E = poisson_solver(f5.density());
    Edvf = weno::trp_v(f5,E);
    for ( auto k=0 ; k<f.size(0) ; ++k ) {
      hf.fft(&(f[k][0]));
      hf1.fft(&(f1[k][0]));
      hf2.fft(&(f2[k][0]));
      hf3.fft(&(f3[k][0]));
      hf4.fft(&(f4[k][0]));
      hf5.fft(&(f5[k][0]));
      hEdvf.fft(&(Edvf[k][0]));

      for ( auto i=0 ; i<Nx ; ++i ) {
        hf[i] = -(563./3456.)*std::exp(L*dt)*hf[i] - (575./252.)*std::exp((4./5.)*L*dt)*hf1[i] + (31400./10017.)*std::exp((7./10.)*L*dt)*hf2[i] + (325./1344.)*std::exp((1./5.)*L*dt)*hf3[i] - (7533./6784.)*std::exp((1./9.)*L*dt)*hf4[i] + (33./28.)*hf5[i] - (11./84.)*dt*hEdvf[i];
      }
      hf.ifft(&(f[k][0]));
    }

    **/
    /**
    // RK(8,6) ================================================================
    #define SCHEME"RK86"
    E = poisson_solver(f.density());
    field<double,1> Edvf = weno::trp_v(f,E);
    for ( auto k=0 ; k<f.size(0) ; ++k ) {
      hf.fft(&(f[k][0]));
      hEdvf.fft(&(Edvf[k][0]));

      for ( auto i=0 ; i<Nx ; ++i ) {
        hf1[i] = std::exp((1./9.)*L*dt)*hf[i] - (1./9.)*dt*std::exp((1./9.)*dt*L)*hEdvf[i];
      }
      hf1.ifft(&(f1[k][0]));
    }

    E = poisson_solver(f1.density());
    Edvf = weno::trp_v(f1,E);
    for ( auto k=0 ; k<f.size(0) ; ++k ) {
      hf.fft(&(f[k][0]));
      hf1.fft(&(f1[k][0]));
      hEdvf.fft(&(Edvf[k][0]));

      for ( auto i=0 ; i<Nx ; ++i ) {
        hf2[i] = (5./8.)*std::exp((1./6.)*L*dt)*hf[i] + (3./8.)*std::exp((1./18.)*L*dt)*hf1[i] - (1./8.)*dt*std::exp((1./18.)*L*dt)*hEdvf[i];
      }
      hf2.ifft(&(f2[k][0]));
    }

    E = poisson_solver(f2.density());
    Edvf = weno::trp_v(f2,E);
    for ( auto k=0 ; k<f.size(0) ; ++k ) {
      hf.fft(&(f[k][0]));
      hf1.fft(&(f1[k][0]));
      hf2.fft(&(f2[k][0]));
      hEdvf.fft(&(Edvf[k][0]));

      for ( auto i=0 ; i<Nx ; ++i ) {
        hf3[i] = 2.*std::exp((1./3.)*L*dt)*hf[i] + 3.*std::exp((2./9.)*L*dt)*hf1[i] -4.*std::exp((1./6.)*L*dt)*hf2[i] - (2./3.)*dt*std::exp((1./6.)*L*dt)*hEdvf[i];
      }
      hf3.ifft(&(f3[k][0]));
    }

    E = poisson_solver(f3.density());
    Edvf = weno::trp_v(f3,E);
    for ( auto k=0 ; k<f.size(0) ; ++k ) {
      hf.fft(&(f[k][0]));
      hf1.fft(&(f1[k][0]));
      hf2.fft(&(f2[k][0]));
      hf3.fft(&(f3[k][0]));
      hEdvf.fft(&(Edvf[k][0]));

      for ( auto i=0 ; i<Nx ; ++i ) {
        hf4[i] = (305./1268.)*std::exp(0.5*L*dt)*hf[i] + (2817./1268.)*std::exp((7./18.)*L*dt)*hf1[i] - (927./317.)*std::exp((1./3.)*L*dt)*hf2[i] + (927./634.)*std::exp((1./6.)*L*dt)*hf3[i] - (321./1268.)*dt*std::exp((1./6.)*L*dt)*hEdvf[i];
      }
      hf4.ifft(&(f4[k][0]));
    }

    E = poisson_solver(f4.density());
    Edvf = weno::trp_v(f4,E);
    for ( auto k=0 ; k<f.size(0) ; ++k ) {
      hf.fft(&(f[k][0]));
      hf1.fft(&(f1[k][0]));
      hf2.fft(&(f2[k][0]));
      hf3.fft(&(f3[k][0]));
      hf4.fft(&(f4[k][0]));
      hEdvf.fft(&(Edvf[k][0]));

      for ( auto i=0 ; i<Nx ; ++i ) {
        hf5[i] = (3191./321.)*std::exp((2./3.)*L*dt)*hf[i] - (2436./107.)*std::exp((5./9.)*L*dt)*hf1[i] - (2404./107.)*std::exp(0.5*L*dt)*hf2[i] + (12330./107.)*std::exp((1./3.)*L*dt)*hf3[i] - (25340./321.)*std::exp((1./6.)*L*dt)*hf4[i] - 8.*dt*std::exp((1./6.)*L*dt)*hEdvf[i];
      }
      hf5.ifft(&(f5[k][0]));
    }

    E = poisson_solver(f5.density());
    Edvf = weno::trp_v(f5,E);
    for ( auto k=0 ; k<f.size(0) ; ++k ) {
      hf.fft(&(f[k][0]));
      hf1.fft(&(f1[k][0]));
      hf2.fft(&(f2[k][0]));
      hf3.fft(&(f3[k][0]));
      hf4.fft(&(f4[k][0]));
      hf5.fft(&(f5[k][0]));
      hEdvf.fft(&(Edvf[k][0]));

      for ( auto i=0 ; i<Nx ; ++i ) {
        hf6[i] = -(15130159./6286464.)*std::exp((5./6.)*L*dt)*hf[i] + (2014319./349248.)*std::exp((13./18.)*L*dt)*hf1[i] + (1194095./523872.)*std::exp((2./3.)*L*dt)*hf2[i] - (1471057./116416.)*std::exp(0.5*L*dt)*hf3[i] + (12601453./1571616.)*std::exp((1./3.)*L*dt)*hf4[i] - (433./19584.)*std::exp((1./6.)*L*dt)*hf5[i] - (33./1088.)*dt*std::exp((1./6.)*L*dt)*hEdvf[i];
      }
      hf6.ifft(&(f6[k][0]));
    }

    E = poisson_solver(f6.density());
    Edvf = weno::trp_v(f6,E);
    for ( auto k=0 ; k<f.size(0) ; ++k ) {
      hf.fft(&(f[k][0]));
      hf1.fft(&(f1[k][0]));
      hf2.fft(&(f2[k][0]));
      hf3.fft(&(f3[k][0]));
      hf4.fft(&(f4[k][0]));
      hf5.fft(&(f5[k][0]));
      hf6.fft(&(f6[k][0]));
      hEdvf.fft(&(Edvf[k][0]));

      for ( auto i=0 ; i<Nx ; ++i ) {
        hf7[i] = (805187./78966.)*std::exp(L*dt)*hf[i] - (2263766./48257.)*std::exp((8./9.)*L*dt)*hf1[i] + (2745422./144771.)*std::exp((5./6.)*L*dt)*hf2[i] + (2271108./48257.)*std::exp((2./3.)*L*dt)*hf3[i] - (13115270./434313.)*std::exp(0.5*L*dt)*hf4[i] - (227./2706.)*std::exp((1./3.)*L*dt)*hf5[i] + (888./451.)*std::exp((1./6.)*L*dt)*hf6[i] - (36./41.)*dt*std::exp((1./6.)*L*dt)*hEdvf[i];
      }
      hf7.ifft(&(f7[k][0]));
    }

    E = poisson_solver(f7.density());
    Edvf = weno::trp_v(f7,E);
    for ( auto k=0 ; k<f.size(0) ; ++k ) {
      hf.fft(&(f[k][0]));
      hf1.fft(&(f1[k][0]));
      hf2.fft(&(f2[k][0]));
      hf3.fft(&(f3[k][0]));
      hf4.fft(&(f4[k][0]));
      hf5.fft(&(f5[k][0]));
      hf6.fft(&(f6[k][0]));
      hf7.fft(&(f7[k][0]));
      hEdvf.fft(&(Edvf[k][0]));

      for ( auto i=0 ; i<Nx ; ++i ) {
        hf[i] = -(193999./179760.)*std::exp(L*dt)*hf[i] + (2487363./329560.)*std::exp((8./9.)*L*dt)*hf1[i] - (847909./164780.)*std::exp((5./6.)*L*dt)*hf2[i] - (1600251./329560.)*std::exp((2./3.)*L*dt)*hf3[i] + (362713./98868.)*std::exp(0.5*L*dt)*hf4[i] + (109./1232.)*std::exp((1./3.)*L*dt)*hf5[i] + (186./385.)*std::exp((1./6.)*L*dt)*hf6[i] + (41./140.)*hf7[i] - (41./840.)*dt*hEdvf[i];
      }
      hf.ifft(&(f[k][0]));
    }
    **/

    // end of time loop
    ++i_t;
    rho = f.density();
    E = poisson_solver(rho);
    //ee(i_t) = 0.;
    //for ( auto i=0 ; i<Nx ; ++i ) {
    //  ee(i_t) += SQ(E(i))*f.step.dx;
    //}
    //ee(i_t) = std::sqrt(ee(i_t));
    H(i_t) = energy(f,E);

//#define FOLDER "lukas/vp/"
#define FOLDER "vphl/kin/"
#define SPACE_SCHEME "weno"
/*
    if ( i_t == int(15./dt) ) {
      std::stringstream ss; ss << FOLDER << "vp_" << SCHEME << "_" << SPACE_SCHEME << "_15.dat";
      f.write(ss.str());
      std::cout << std::endl;
    }
    if ( i_t == int(20./dt) ) {
      std::stringstream ss; ss << FOLDER << "vp_" << SCHEME << "_" << SPACE_SCHEME << "_20.dat";
      f.write(ss.str());
      std::cout << std::endl;
    }
    if ( i_t == int(25./dt) ) {
      std::stringstream ss; ss << FOLDER << "vp_" << SCHEME << "_" << SPACE_SCHEME << "_25.dat";
      f.write(ss.str());
      std::cout << std::endl;
    }
    if ( i_t == int(30./dt) ) {
      std::stringstream ss; ss << FOLDER << "vp_" << SCHEME << "_" << SPACE_SCHEME << "_30.dat";
      f.write(ss.str());
      std::cout << std::endl;
    }
    if ( i_t == int(35./dt) ) {
      std::stringstream ss; ss << FOLDER << "vp_" << SCHEME << "_" << SPACE_SCHEME << "_35.dat";
      f.write(ss.str());
      std::cout << std::endl;
    }
*/
	} // while (  i_t*dt < Tf )
#undef L
  std::cout<<" ["<<std::setw(5)<<i_t<<"] "<<i_t*dt<<"   \r"<<std::endl;

  std::stringstream ss; ss << FOLDER << "vp_" << SCHEME << "_" << SPACE_SCHEME << "_dt" << dt << ".dat";
  f.write(ss.str());

  rho = f.density();
  auto dx_y = [&,count=0](auto const& y) mutable { std::stringstream ss; ss<<Xi(count++)<<" "<<y; return ss.str(); };
  ss.str(std::string());
  ss << FOLDER << "rho_" << SCHEME << "_" << SPACE_SCHEME << "_" << Tf << ".dat";
  std::ofstream of;
  of.open(ss.str()); ss.str(std::string());
  std::transform( rho.begin() , rho.end() , std::ostream_iterator<std::string>(of,"\n") , dx_y );
  of.close();

  ss.str(std::string());
  //std::ofstream of;
  std::size_t count = 0;
  auto dt_y = [&,count=0](auto const& y) mutable { std::stringstream ss; ss<<(count++)*dt<<" "<<y; return ss.str(); };
  ss << FOLDER << "ee_" << SCHEME << "_" << SPACE_SCHEME << "_10.dat";
  of.open(ss.str()); ss.str(std::string());
  for ( auto i=0; i<ee.size() ; ++i ) {
    of << i*dt <<" " << ee[i] << "\n";
  }
  of.close();
  ss << FOLDER << "H_" << SCHEME << "_" << SPACE_SCHEME << ".dat";
  of.open(ss.str()); ss.str(std::string());
  for ( auto i=0; i<H.size() ; ++i ) {
    of << i*dt <<" " << (H[i]-H[0])/std::abs(H[0]) << "\n";
  }

  //double h = std::abs(*std::max_element( H.begin() , H.end() , [&](double a,double b){return ( std::abs((a-H[0])/std::abs(H[0])) < std::abs((b-H[0])/std::abs(H[0])) );} ));
  //std::cout << dt << " " << std::abs((h-H[0])/std::abs(H[0])) << "\n";

  of.close();
  ss << FOLDER << "Emax_" << SCHEME << "_" << SPACE_SCHEME << "_10.dat";
  of.open(ss.str()); ss.str(std::string());
  std::transform( Emax.begin() , Emax.end() , std::ostream_iterator<std::string>(of,"\n") , dt_y );
  of.close();

	return 0;
}

ext: cc

max_stars_repo_name: kivvix/vlasovpp

max_issues_repo_name: kivvix/vlasovpp

max_forks_repo_name: kivvix/vlasovpp


text:
/**
 * @file loss_functions_test.cpp
 * @author Dakshit Agrawal
 * @author Sourabh Varshney
 * @author Atharva Khandait
 * @author Saksham Rastogi
 *
 * Tests for loss functions in mlpack::methods::ann:loss_functions.
 *
 * mlpack is free software; you may redistribute it and/or modify it under the
 * terms of the 3-clause BSD license.  You should have received a copy of the
 * 3-clause BSD license along with mlpack.  If not, see
 * http://www.opensource.org/licenses/BSD-3-Clause for more information.
 */
#include <mlpack/core.hpp>

#include <mlpack/methods/ann/layer/layer.hpp>
#include <mlpack/methods/ann/loss_functions/kl_divergence.hpp>
#include <mlpack/methods/ann/loss_functions/earth_mover_distance.hpp>
#include <mlpack/methods/ann/loss_functions/mean_squared_error.hpp>
#include <mlpack/methods/ann/loss_functions/sigmoid_cross_entropy_error.hpp>
#include <mlpack/methods/ann/loss_functions/cross_entropy_error.hpp>
#include <mlpack/methods/ann/loss_functions/reconstruction_loss.hpp>
#include <mlpack/methods/ann/loss_functions/mean_squared_logarithmic_error.hpp>
#include <mlpack/methods/ann/loss_functions/mean_bias_error.hpp>
#include <mlpack/methods/ann/loss_functions/dice_loss.hpp>
#include <mlpack/methods/ann/loss_functions/log_cosh_loss.hpp>
#include <mlpack/methods/ann/init_rules/nguyen_widrow_init.hpp>
#include <mlpack/methods/ann/ffn.hpp>

#include <boost/test/unit_test.hpp>
#include "test_tools.hpp"
#include "ann_test_tools.hpp"

using namespace mlpack;
using namespace mlpack::ann;

BOOST_AUTO_TEST_SUITE(LossFunctionsTest);

/**
 * Simple KL Divergence test.  The loss should be zero if input = target.
 */
BOOST_AUTO_TEST_CASE(SimpleKLDivergenceTest)
{
  arma::mat input, target, output;
  double loss;
  KLDivergence<> module(true);

  // Test the Forward function.  Loss should be 0 if input = target.
  input = arma::ones(10, 1);
  target = arma::ones(10, 1);
  loss = module.Forward(std::move(input), std::move(target));
  BOOST_REQUIRE_SMALL(loss, 0.00001);
}

/*
 * Simple test for the mean squared logarithmic error function.
 */
BOOST_AUTO_TEST_CASE(SimpleMeanSquaredLogarithmicErrorTest)
{
  arma::mat input, output, target;
  MeanSquaredLogarithmicError<> module;

  // Test the Forward function on a user generator input and compare it against
  // the manually calculated result.
  input = arma::zeros(1, 8);
  target = arma::zeros(1, 8);
  double error = module.Forward(std::move(input), std::move(target));
  BOOST_REQUIRE_SMALL(error, 0.00001);

  // Test the Backward function.
  module.Backward(std::move(input), std::move(target), std::move(output));
  // The output should be equal to 0.
  CheckMatrices(input, output);
  BOOST_REQUIRE_EQUAL(output.n_rows, input.n_rows);
  BOOST_REQUIRE_EQUAL(output.n_cols, input.n_cols);

  // Test the error function on a single input.
  input = arma::mat("2");
  target = arma::mat("3");
  error = module.Forward(std::move(input), std::move(target));
  BOOST_REQUIRE_CLOSE(error, 0.082760974810151655, 0.001);

  // Test the Backward function on a single input.
  module.Backward(std::move(input), std::move(target), std::move(output));
  BOOST_REQUIRE_CLOSE(arma::accu(output), -0.1917880483011872, 0.001);
  BOOST_REQUIRE_EQUAL(output.n_elem, 1);
}

/**
 * Test to check KL Divergence loss function when we take mean.
 */
BOOST_AUTO_TEST_CASE(KLDivergenceMeanTest)
{
  arma::mat input, target, output;
  double loss;
  KLDivergence<> module(true);

  // Test the Forward function.
  input = arma::mat("1 1 1 1 1 1 1 1 1 1");
  target = arma::exp(arma::mat("2 1 1 1 1 1 1 1 1 1"));

  loss = module.Forward(std::move(input), std::move(target));
  BOOST_REQUIRE_CLOSE_FRACTION(loss, -1.1 , 0.00001);

  // Test the Backward function.
  module.Backward(std::move(input), std::move(target), std::move(output));
  BOOST_REQUIRE_CLOSE_FRACTION(arma::as_scalar(output), -0.1, 0.00001);
}

/**
 * Test to check KL Divergence loss function when we do not take mean.
 */
BOOST_AUTO_TEST_CASE(KLDivergenceNoMeanTest)
{
  arma::mat input, target, output;
  double loss;
  KLDivergence<> module(false);

  // Test the Forward function.
  input = arma::mat("1 1 1 1 1 1 1 1 1 1");
  target = arma::exp(arma::mat("2 1 1 1 1 1 1 1 1 1"));

  loss = module.Forward(std::move(input), std::move(target));
  BOOST_REQUIRE_CLOSE_FRACTION(loss, -11, 0.00001);

  // Test the Backward function.
  module.Backward(std::move(input), std::move(target), std::move(output));
  BOOST_REQUIRE_CLOSE_FRACTION(arma::as_scalar(output), -1, 0.00001);
}

/*
 * Simple test for the mean squared error performance function.
 */
BOOST_AUTO_TEST_CASE(SimpleMeanSquaredErrorTest)
{
  arma::mat input, output, target;
  MeanSquaredError<> module;

  // Test the Forward function on a user generator input and compare it against
  // the manually calculated result.
  input = arma::mat("1.0 0.0 1.0 0.0 -1.0 0.0 -1.0 0.0");
  target = arma::zeros(1, 8);
  double error = module.Forward(std::move(input), std::move(target));
  BOOST_REQUIRE_EQUAL(error, 0.5);

  // Test the Backward function.
  module.Backward(std::move(input), std::move(target), std::move(output));
  // We subtract a zero vector, so according to the used backward formula:
  // output = 2 * (input - target) / target.n_cols,
  // output * nofColumns / 2 should be equal to input.
  CheckMatrices(input, output * output.n_cols / 2);
  BOOST_REQUIRE_EQUAL(output.n_rows, input.n_rows);
  BOOST_REQUIRE_EQUAL(output.n_cols, input.n_cols);

  // Test the error function on a single input.
  input = arma::mat("2");
  target = arma::mat("3");
  error = module.Forward(std::move(input), std::move(target));
  BOOST_REQUIRE_EQUAL(error, 1.0);

  // Test the Backward function on a single input.
  module.Backward(std::move(input), std::move(target), std::move(output));
  // Test whether the output is negative.
  BOOST_REQUIRE_EQUAL(arma::accu(output), -2);
  BOOST_REQUIRE_EQUAL(output.n_elem, 1);
}

/*
 * Simple test for the cross-entropy error performance function.
 */
BOOST_AUTO_TEST_CASE(SimpleCrossEntropyErrorTest)
{
  arma::mat input1, input2, output, target1, target2;
  CrossEntropyError<> module(1e-6);

  // Test the Forward function on a user generator input and compare it against
  // the manually calculated result.
  input1 = arma::mat("0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5");
  target1 = arma::zeros(1, 8);
  double error1 = module.Forward(std::move(input1), std::move(target1));
  BOOST_REQUIRE_SMALL(error1 - 8 * std::log(2), 2e-5);

  input2 = arma::mat("0 1 1 0 1 0 0 1");
  target2 = arma::mat("0 1 1 0 1 0 0 1");
  double error2 = module.Forward(std::move(input2), std::move(target2));
  BOOST_REQUIRE_SMALL(error2, 1e-5);

  // Test the Backward function.
  module.Backward(std::move(input1), std::move(target1), std::move(output));
  for (double el : output)
  {
    // For the 0.5 constant vector we should get 1 / (1 - 0.5) = 2 everywhere.
    BOOST_REQUIRE_SMALL(el - 2, 5e-6);
  }
  BOOST_REQUIRE_EQUAL(output.n_rows, input1.n_rows);
  BOOST_REQUIRE_EQUAL(output.n_cols, input1.n_cols);

  module.Backward(std::move(input2), std::move(target2), std::move(output));
  for (size_t i = 0; i < 8; ++i)
  {
    double el = output.at(0, i);
    if (input2.at(i) == 0)
      BOOST_REQUIRE_SMALL(el - 1, 2e-6);
    else
      BOOST_REQUIRE_SMALL(el + 1, 2e-6);
  }
  BOOST_REQUIRE_EQUAL(output.n_rows, input2.n_rows);
  BOOST_REQUIRE_EQUAL(output.n_cols, input2.n_cols);
}

/**
 * Simple test for the Sigmoid Cross Entropy performance function.
 */
BOOST_AUTO_TEST_CASE(SimpleSigmoidCrossEntropyErrorTest)
{
  arma::mat input1, input2, input3, output, target1,
            target2, target3, expectedOutput;
  SigmoidCrossEntropyError<> module;

  // Test the Forward function on a user generator input and compare it against
  // the calculated result.
  input1 = arma::mat("0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5");
  target1 = arma::zeros(1, 8);
  double error1 = module.Forward(std::move(input1), std::move(target1));
  double expected = 0.97407699;
  // Value computed using tensorflow.
  BOOST_REQUIRE_SMALL(error1 / input1.n_elem - expected, 1e-7);

  input2 = arma::mat("1 2 3 4 5");
  target2 = arma::mat("0 0 1 0 1");
  double error2 = module.Forward(std::move(input2), std::move(target2));
  expected = 1.5027283;
  BOOST_REQUIRE_SMALL(error2 / input2.n_elem - expected, 1e-6);

  input3 = arma::mat("0 -1 -1 0 -1 0 0 -1");
  target3 = arma::mat("0 -1 -1 0 -1 0 0 -1");
  double error3 = module.Forward(std::move(input3), std::move(target3));
  expected = 0.00320443;
  BOOST_REQUIRE_SMALL(error3 / input3.n_elem - expected, 1e-6);

  // Test the Backward function.
  module.Backward(std::move(input1), std::move(target1), std::move(output));
  expected = 0.62245929;
  for (size_t i = 0; i < output.n_elem; i++)
    BOOST_REQUIRE_SMALL(output(i) - expected, 1e-5);
  BOOST_REQUIRE_EQUAL(output.n_rows, input1.n_rows);
  BOOST_REQUIRE_EQUAL(output.n_cols, input1.n_cols);

  expectedOutput = arma::mat(
      "0.7310586 0.88079709 -0.04742587 0.98201376 -0.00669285");
  module.Backward(std::move(input2), std::move(target2), std::move(output));
  for (size_t i = 0; i < output.n_elem; i++)
    BOOST_REQUIRE_SMALL(output(i) - expectedOutput(i), 1e-5);
  BOOST_REQUIRE_EQUAL(output.n_rows, input2.n_rows);
  BOOST_REQUIRE_EQUAL(output.n_cols, input2.n_cols);

  module.Backward(std::move(input3), std::move(target3), std::move(output));
  expectedOutput = arma::mat("0.5 1.2689414");
  for (size_t i = 0; i < 8; ++i)
  {
    double el = output.at(0, i);
    if (std::abs(input3.at(i) - 0.0) < 1e-5)
      BOOST_REQUIRE_SMALL(el - expectedOutput[0], 2e-6);
    else
      BOOST_REQUIRE_SMALL(el - expectedOutput[1], 2e-6);
  }
  BOOST_REQUIRE_EQUAL(output.n_rows, input3.n_rows);
  BOOST_REQUIRE_EQUAL(output.n_cols, input3.n_cols);
}

/**
 * Simple test for the Earth Mover Distance Layer.
 */
BOOST_AUTO_TEST_CASE(SimpleEarthMoverDistanceLayerTest)
{
  arma::mat input1, input2, output, target1, target2, expectedOutput;
  EarthMoverDistance<> module;

  // Test the Forward function on a user generator input and compare it against
  // the manually calculated result.
  input1 = arma::mat("0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5");
  target1 = arma::zeros(1, 8);
  double error1 = module.Forward(std::move(input1), std::move(target1));
  double expected = 0.0;
  BOOST_REQUIRE_SMALL(error1 / input1.n_elem - expected, 1e-7);

  input2 = arma::mat("1 2 3 4 5");
  target2 = arma::mat("1 0 1 0 1");
  double error2 = module.Forward(std::move(input2), std::move(target2));
  expected = -1.8;
  BOOST_REQUIRE_SMALL(error2 / input2.n_elem - expected, 1e-6);

  // Test the Backward function.
  module.Backward(std::move(input1), std::move(target1), std::move(output));
  expected = 0.0;
  for (size_t i = 0; i < output.n_elem; i++)
    BOOST_REQUIRE_SMALL(output(i) - expected, 1e-5);
  BOOST_REQUIRE_EQUAL(output.n_rows, input1.n_rows);
  BOOST_REQUIRE_EQUAL(output.n_cols, input1.n_cols);

  expectedOutput = arma::mat("-1 0 -1 0 -1");
  module.Backward(std::move(input2), std::move(target2), std::move(output));
  for (size_t i = 0; i < output.n_elem; i++)
    BOOST_REQUIRE_SMALL(output(i) - expectedOutput(i), 1e-5);
  BOOST_REQUIRE_EQUAL(output.n_rows, input2.n_rows);
  BOOST_REQUIRE_EQUAL(output.n_cols, input2.n_cols);
}

/*
 * Mean Squared Error numerical gradient test.
 */
BOOST_AUTO_TEST_CASE(GradientMeanSquaredErrorTest)
{
  // Linear function gradient instantiation.
  struct GradientFunction
  {
    GradientFunction()
    {
      input = arma::randu(10, 1);
      target = arma::randu(2, 1);

      model = new FFN<MeanSquaredError<>, NguyenWidrowInitialization>();
      model->Predictors() = input;
      model->Responses() = target;
      model->Add<IdentityLayer<> >();
      model->Add<Linear<> >(10, 2);
      model->Add<SigmoidLayer<> >();
    }

    ~GradientFunction()
    {
      delete model;
    }

    double Gradient(arma::mat& gradient) const
    {
      arma::mat output;
      double error = model->Evaluate(model->Parameters(), 0, 1);
      model->Gradient(model->Parameters(), 0, gradient, 1);
      return error;
    }

    arma::mat& Parameters() { return model->Parameters(); }

    FFN<MeanSquaredError<>, NguyenWidrowInitialization>* model;
    arma::mat input, target;
  } function;

  BOOST_REQUIRE_LE(CheckGradient(function), 1e-4);
}

/*
 * Reconstruction Loss numerical gradient test.
 */
BOOST_AUTO_TEST_CASE(GradientReconstructionLossTest)
{
  // Linear function gradient instantiation.
  struct GradientFunction
  {
    GradientFunction()
    {
      input = arma::randu(10, 1);
      target = arma::randu(2, 1);

      model = new FFN<ReconstructionLoss<>, NguyenWidrowInitialization>();
      model->Predictors() = input;
      model->Responses() = target;
      model->Add<IdentityLayer<> >();
      model->Add<Linear<> >(10, 2);
      model->Add<SigmoidLayer<> >();
    }

    ~GradientFunction()
    {
      delete model;
    }

    double Gradient(arma::mat& gradient) const
    {
      arma::mat output;
      double error = model->Evaluate(model->Parameters(), 0, 1);
      model->Gradient(model->Parameters(), 0, gradient, 1);
      return error;
    }

    arma::mat& Parameters() { return model->Parameters(); }

    FFN<ReconstructionLoss<>, NguyenWidrowInitialization>* model;
    arma::mat input, target;
  } function;

  BOOST_REQUIRE_LE(CheckGradient(function), 1e-4);
}

/*
 * Simple test for the dice loss function.
 */
BOOST_AUTO_TEST_CASE(DiceLossTest)
{
  arma::mat input1, input2, target, output;
  double loss;
  DiceLoss<> module;

  // Test the Forward function. Loss should be 0 if input = target.
  input1 = arma::ones(10, 1);
  target = arma::ones(10, 1);
  loss = module.Forward(std::move(input1), std::move(target));
  BOOST_REQUIRE_SMALL(loss, 0.00001);

  // Test the Forward function. Loss should be 0.185185185.
  input2 = arma::ones(10, 1) * 0.5;
  loss = module.Forward(std::move(input2), std::move(target));
  BOOST_REQUIRE_CLOSE(loss, 0.185185185, 0.00001);

  // Test the Backward function for input = target.
  module.Backward(std::move(input1), std::move(target), std::move(output));
  for (double el : output)
  {
    // For input = target we should get 0.0 everywhere.
    BOOST_REQUIRE_CLOSE(el, 0.0, 0.00001);
  }
  BOOST_REQUIRE_EQUAL(output.n_rows, input1.n_rows);
  BOOST_REQUIRE_EQUAL(output.n_cols, input1.n_cols);

  // Test the Backward function.
  module.Backward(std::move(input2), std::move(target), std::move(output));
  for (double el : output)
  {
    // For the 0.5 constant vector we should get -0.0877914951989026 everywhere.
    BOOST_REQUIRE_CLOSE(el, -0.0877914951989026, 0.00001);
  }
  BOOST_REQUIRE_EQUAL(output.n_rows, input2.n_rows);
  BOOST_REQUIRE_EQUAL(output.n_cols, input2.n_cols);
}

/*
 * Simple test for the mean bias error performance function.
 */
BOOST_AUTO_TEST_CASE(SimpleMeanBiasErrorTest)
{
  arma::mat input, output, target;
  MeanBiasError<> module;

  // Test the Forward function on a user generator input and compare it against
  // the manually calculated result.
  input = arma::mat("1.0 0.0 1.0 -1.0 -1.0 0.0 -1.0 0.0");
  target = arma::zeros(1, 8);
  double error = module.Forward(std::move(input), std::move(target));
  BOOST_REQUIRE_EQUAL(error, 0.125);

  // Test the Backward function.
  module.Backward(std::move(input), std::move(target), std::move(output));
  // We should get a vector with -1 everywhere.
  for (double el : output)
  {
    BOOST_REQUIRE_EQUAL(el, -1);
  }
  BOOST_REQUIRE_EQUAL(output.n_rows, input.n_rows);
  BOOST_REQUIRE_EQUAL(output.n_cols, input.n_cols);

  // Test the error function on a single input.
  input = arma::mat("2");
  target = arma::mat("3");
  error = module.Forward(std::move(input), std::move(target));
  BOOST_REQUIRE_EQUAL(error, 1.0);

  // Test the Backward function on a single input.
  module.Backward(std::move(input), std::move(target), std::move(output));
  // Test whether the output is negative.
  BOOST_REQUIRE_EQUAL(arma::accu(output), -1);
  BOOST_REQUIRE_EQUAL(output.n_elem, 1);
}

/**
 * Simple test for the Log-Hyperbolic-Cosine loss function.
 */
BOOST_AUTO_TEST_CASE(LogCoshLossTest)
{
  arma::mat input, target, output;
  double loss;
  LogCoshLoss<> module(2);

  // Test the Forward function. Loss should be 0 if input = target.
  input = arma::ones(10, 1);
  target = arma::ones(10, 1);
  loss = module.Forward(std::move(input), std::move(target));
  BOOST_REQUIRE_EQUAL(loss, 0);

  // Test the Backward function for input = target.
  module.Backward(std::move(input), std::move(target), std::move(output));
  for (double el : output)
  {
    // For input = target we should get 0.0 everywhere.
    BOOST_REQUIRE_CLOSE(el, 0.0, 1e-5);
  }

  BOOST_REQUIRE_EQUAL(output.n_rows, input.n_rows);
  BOOST_REQUIRE_EQUAL(output.n_cols, input.n_cols);

  // Test the Forward function. Loss should be 0.546621.
  input = arma::mat("1 2 3 4 5");
  target = arma::mat("1 2.4 3.4 4.2 5.5");
  loss = module.Forward(std::move(input), std::move(target));
  BOOST_REQUIRE_CLOSE(loss, 0.546621, 1e-3);

  // Test the Backward function.
  module.Backward(std::move(input), std::move(target), std::move(output));
  BOOST_REQUIRE_CLOSE(arma::accu(output), 2.46962, 1e-3);
  BOOST_REQUIRE_EQUAL(output.n_rows, input.n_rows);
  BOOST_REQUIRE_EQUAL(output.n_cols, input.n_cols);
}
BOOST_AUTO_TEST_SUITE_END();

max_stars_repo_name: mhmohona/mlpack

max_issues_repo_name: mhmohona/mlpack

max_forks_repo_name: mhmohona/mlpack