text:
/*
 * Copyright (c) 2013 Juniper Networks, Inc. All rights reserved.
 */

#include "base/os.h"
#include "testing/gunit.h"

#include <boost/scoped_array.hpp>
#include <base/logging.h>

#include <pugixml/pugixml.hpp>
#include <io/event_manager.h>
#include "http/http_request.h"
#include "http/http_session.h"
#include "http/http_server.h"
#include "http/client/http_client.h"
#include "http/client/http_curl.h"
#include <cmn/agent_cmn.h>
#include <oper/operdb_init.h>
#include <oper/interface_common.h>
#include <pkt/pkt_init.h>
#include <services/services_init.h>
#include <test/test_cmn_util.h>
#include <services/services_sandesh.h>
#include <services/metadata_proxy.h>
#include "pkt/test/test_pkt_util.h"

#define MAX_WAIT_COUNT 5000
#define BUF_SIZE 8192
#define vm1_ip "10.1.1.3"
#define METADATA_CHECK(condition)                                              \
                    count = 0;                                                 \
                    do {                                                       \
                      client->WaitForIdle();                                   \
                      stats = Agent::GetInstance()->services()->               \
                              metadataproxy()->metadatastats();                \
                      if (++count == MAX_WAIT_COUNT)                           \
                          assert(0);                                           \
                      usleep(1000);                                            \
                    } while (condition);                                       \

IpamInfo ipam_info[] = {
    {"10.1.1.0", 24, "10.1.1.1"},
};

pthread_t curl_thread;
// send multiple requests using curl
void *SendMultiGetRequest(void *arg) {
    uint16_t port = Agent::GetInstance()->metadata_server_port();
    stringstream str;
    str << "curl http://127.0.0.1:";
    str << port;
    str << "/{req1,req2}";
    string cmd = str.str();
    int res = system(cmd.c_str());
    if (res) {
        LOG(WARN, "curl command failed");
    }
    return NULL;
}

void CurlRun() {
    pthread_attr_t attr;
    int ret;

    pthread_attr_init(&attr);
    if ((ret = pthread_create(&curl_thread, &attr, SendMultiGetRequest, NULL)) != 0) {
        LOG(ERROR, "pthread_create error : " <<  strerror(ret) );
        assert(0);
    }
}

void CurlStop() {
    pthread_join(curl_thread, NULL);
}

class TestInterfaceTable : public InterfaceTable {
public:
    TestInterfaceTable() :
        InterfaceTable(Agent::GetInstance()->db(), "test") {}
    ~TestInterfaceTable() {}

    bool FindVmUuidFromMetadataIp(const Ip4Address &ip,
                                  std::string *vm_ip,
                                  std::string *vm_uuid,
                                  std::string *vm_project_uuid) {
        *vm_ip = vm1_ip;
        *vm_uuid = "1234";
        *vm_project_uuid = "5678";
        return true;
    }
};

class MetadataTest : public ::testing::Test {
public:
    enum HttpMethod {
        GET_METHOD,
        HEAD_METHOD,
        POST_METHOD,
        PUT_METHOD,
        DELETE_METHOD
    };

    void SetUp() {
        agent_ = Agent::GetInstance();
        AddIPAM("vn1", ipam_info, 1);
        client->WaitForIdle();
    }

    void TearDown() {
        DelIPAM("vn1");
        client->WaitForIdle();
    }

    MetadataTest() : nova_api_proxy_(NULL), vm_http_client_(NULL),
                     done_(0), itf_count_(0), data_size_(0) {
        rid_ = Agent::GetInstance()->interface_table()->Register(
                boost::bind(&MetadataTest::ItfUpdate, this, _2));
        Agent::GetInstance()->set_router_id(Ip4Address::from_string("127.0.0.1"));
        Agent::GetInstance()->set_compute_node_ip(Ip4Address::from_string("127.0.0.1"));
    }

    ~MetadataTest() {
        Agent::GetInstance()->interface_table()->Unregister(rid_);
    }

    void ItfUpdate(DBEntryBase *entry) {
        Interface *itf = static_cast<Interface *>(entry);
        tbb::mutex::scoped_lock lock(mutex_);
        unsigned int i;
        for (i = 0; i < itf_id_.size(); ++i)
            if (itf_id_[i] == itf->id())
                break;
        if (entry->IsDeleted()) {
            if (itf_count_ && i < itf_id_.size()) {
                itf_count_--;
                LOG(DEBUG, "Metadata test : interface deleted " << itf_id_[0]);
                itf_id_.erase(itf_id_.begin()); // we delete in create order
            }
        } else {
            if (i == itf_id_.size()) {
                itf_count_++;
                itf_id_.push_back(itf->id());
                LOG(DEBUG, "Metadata test : interface added " << itf->id());
            }
        }
    }

    uint32_t GetItfCount() {
        tbb::mutex::scoped_lock lock(mutex_);
        return itf_count_;
    }

    void WaitForItfUpdate(unsigned int expect_count) {
        int count = 0;
        while (GetItfCount() != expect_count) {
            if (++count == MAX_WAIT_COUNT)
                assert(0);
            usleep(1000);
        }
    }

    void CheckSandeshResponse(Sandesh *sandesh) {
    }

    void SetupLinkLocalConfig() {
        std::stringstream global_config;
        global_config << "<linklocal-services>\n"
                      << "<linklocal-service-entry>\n"
	              << "<linklocal-service-name>metadata</linklocal-service-name>\n"
	              << "<linklocal-service-ip>169.254.169.254</linklocal-service-ip>\n"
	              << "<linklocal-service-port>80</linklocal-service-port>\n"
	              << "<ip-fabric-DNS-service-name></ip-fabric-DNS-service-name>\n"
	              << "<ip-fabric-service-port>"
                      << nova_api_proxy_->GetPort()
                      << "</ip-fabric-service-port>\n"
	              << "<ip-fabric-service-ip>127.0.0.1</ip-fabric-service-ip>\n"
	              << "</linklocal-service-entry>\n"
	              << "</linklocal-services>";

        char buf[BUF_SIZE];
        int len = 0;
        memset(buf, 0, BUF_SIZE);
        AddXmlHdr(buf, len);
        AddNodeString(buf, len, "global-vrouter-config",
                      "default-global-system-config:default-global-vrouter-config",
                      1024, global_config.str().c_str());
        AddXmlTail(buf, len);
        ApplyXmlString(buf);
    }

    void ClearLinkLocalConfig() {
        char buf[BUF_SIZE];
        int len = 0;
        memset(buf, 0, BUF_SIZE);
        AddXmlHdr(buf, len);
        AddNodeString(buf, len, "global-vrouter-config",
                      "default-global-system-config:default-global-vrouter-config",
                      1024, "");
        AddXmlTail(buf, len);
        ApplyXmlString(buf);
    }

    void StartHttpClient() {
        vm_http_client_ = new HttpClient(Agent::GetInstance()->event_manager());
        vm_http_client_->Init();
    }

    void StopHttpClient() {
        vm_http_client_->Shutdown();
        vm_http_client_ = NULL;
    }

    HttpConnection *SendHttpClientRequest(HttpMethod method) {
        std::string body("test body");
        return SendHttpClientRequest(method, body);
    }

    HttpConnection *SendHttpClientRequest(HttpMethod method,
                                          std::string &body) {
        boost::system::error_code ec;
        Ip4Address server(Ip4Address::from_string("127.0.0.1", ec));

        boost::asio::ip::tcp::endpoint http_ep;
        http_ep.address(server);
        http_ep.port(Agent::GetInstance()->metadata_server_port());

        data_size_ = 0;
        std::string uri("openstack");
        std::vector<std::string> header_options;
        header_options.push_back(std::string("Connection: close"));
        HttpConnection *conn = vm_http_client_->CreateConnection(http_ep);
        conn->RegisterEventCb(
              boost::bind(&MetadataTest::OnClientSessionEvent, this, _1, _2));
        switch (method) {
            case GET_METHOD:
                conn->HttpGet(uri, false, false, true, header_options,
                              boost::bind(&MetadataTest::HandleHttpResponse,
                                          this, conn, _1, _2));
                break;
            case HEAD_METHOD:
                conn->HttpHead(uri, false, false, true, header_options,
                               boost::bind(&MetadataTest::HandleHttpResponse,
                                           this, conn, _1, _2));
                break;
            case PUT_METHOD: {
                data_size_ = body.size();
                std::stringstream str;
                str << "Content-Length: ";
                str <HttpPut(body, uri, false, false, true, header_options,
                              boost::bind(&MetadataTest::HandleHttpResponse,
                                          this, conn, _1, _2));
                break;
            }
            case POST_METHOD:
                data_size_ = body.size();
                conn->HttpPost(body, uri, false, false, true, header_options,
                              boost::bind(&MetadataTest::HandleHttpResponse,
                                          this, conn, _1, _2));
                break;
            case DELETE_METHOD:
                conn->HttpDelete(uri, false, false, true, header_options,
                                 boost::bind(&MetadataTest::HandleHttpResponse,
                                             this, conn, _1, _2));
                break;
        }
        return conn;
    }

    void HandleHttpResponse(HttpConnection *conn, std::string &msg,
                            boost::system::error_code &ec) {
        if (ec)
            assert(0);
        CloseClientSession(conn);
        done_++;
    }

    void CloseClientSession(HttpConnection *conn) {
        HttpClient *client = conn->client();
        client->RemoveConnection(conn);
    }

    void OnClientSessionEvent(HttpClientSession *session, TcpSession::Event event) {
        switch (event) {
            case TcpSession::CLOSE: {
                session->Connection()->client()->RemoveConnection(
                                                 session->Connection());
                break;
            }

            default:
                break;
        }
    }

    void StartNovaApiProxy() {
        nova_api_proxy_ =
             new HttpServer(Agent::GetInstance()->event_manager());
        nova_api_proxy_->RegisterHandler(HTTP_WILDCARD_ENTRY,
             boost::bind(&MetadataTest::HandleNovaApiRequest, this, _1, _2));
        nova_api_proxy_->Initialize(0);
    }

    void StopNovaApiProxy() {
        nova_api_proxy_->ClearSessions();
        nova_api_proxy_->Shutdown();
        nova_api_proxy_ = NULL;
    }

    void HandleNovaApiRequest(HttpSession *session, const HttpRequest *request) {
        int found = 0;
        uint32_t data_received = 0;
        const HttpRequest::HeaderMap &req_header = request->Headers();
        for (HttpRequest::HeaderMap::const_iterator it = req_header.begin();
             it != req_header.end(); ++it) {
            std::string option = it->first;
            if (option == "X-Forwarded-For" || option == "X-Instance-ID" ||
                option == "X-Instance-ID-Signature" || option == "X-Tenant-ID")
                found++;
            if (option == "Content-Length") {
                std::stringstream str(it->second);
                str >> data_received;
            }
        }

        if (found != 4)
            assert(0);

        if ((data_received != data_size_ || request->Body().size() != data_size_))
            assert(0);

        const char body[] = "\n"
                            "\n"
                            " <title>Server Status Success</title>\n"
                            "\n"
                            "\n";
        char response[512];
        snprintf(response, sizeof(response),
                 "HTTP/1.1 200 OK\r\n"
                 "Content-Type: text/html; charset=UTF-8\r\n"
                 "Content-Length: %u\r\n"
                 "\r\n%s", (unsigned int)strlen(body), body);

        session->Send(reinterpret_cast<const u_int8_t *>(response),
                      strlen(response), NULL);
        delete request;
    }

    Agent *agent_;
private:
    HttpServer *nova_api_proxy_;
    HttpClient *vm_http_client_;
    uint32_t done_;
    uint32_t itf_count_;
    uint32_t data_size_;
    DBTableBase::ListenerId rid_;
    std::vector<std::size_t> itf_id_;
    tbb::mutex mutex_;
};

TEST_F(MetadataTest, MetadataReqTest) {
    int count = 0;
    MetadataProxy::MetadataStats stats;
    struct PortInfo input[] = {
        {"vnet1", 1, vm1_ip, "00:00:00:01:01:01", 1, 1},
    };

    StartNovaApiProxy();
    SetupLinkLocalConfig();

    CreateVmportEnv(input, 1, 0);
    client->WaitForIdle();
    client->Reset();

    VmInterface *intf = static_cast<VmInterface *>(VmPortGet(1));
    TxTcpPacket(intf->id(), vm1_ip, "169.254.169.254", 1000, 80, false);
    client->WaitForIdle();
    FlowEntry *flow = FlowGet(0, vm1_ip, "169.254.169.254", 6, 1000, 80,
                              intf->flow_key_nh()->id());
    EXPECT_TRUE(flow != NULL);
    FlowEntry *rflow = flow->reverse_flow_entry();
    EXPECT_TRUE(rflow != NULL);

    EXPECT_TRUE(flow->key().src_addr.to_v4() ==
                Ip4Address::from_string(vm1_ip));
    EXPECT_TRUE(flow->key().dst_addr.to_v4() ==
                Ip4Address::from_string("169.254.169.254"));
    EXPECT_TRUE(rflow->key().src_addr.to_v4() == agent_->router_id());
    EXPECT_TRUE(rflow->key().dst_addr.to_v4() == intf->mdata_ip_addr());

    StartHttpClient();

    // If the local address is not same as VM address in this request,
    // agent will have an internal error as it cannot find the VM
    HttpConnection *conn = SendHttpClientRequest(GET_METHOD);
    METADATA_CHECK (stats.internal_errors < 1);
    EXPECT_EQ(1U, stats.requests);
    EXPECT_EQ(0U, stats.responses);
    EXPECT_EQ(0U, stats.proxy_sessions);
    CloseClientSession(conn);
    client->WaitForIdle();

    // Introspect request
    MetadataInfo *sand = new MetadataInfo();
    Sandesh::set_response_callback(
             boost::bind(&MetadataTest::CheckSandeshResponse, this, _1));
    sand->HandleRequest();
    client->WaitForIdle();
    sand->Release();

    // for agent to identify the vm, the remote end should have vm's ip;
    // overload the FindVmUuidFromMetadataIp to return true
    InterfaceTable *intf_table = Agent::GetInstance()->interface_table();
    std::auto_ptr<InterfaceTable> interface_table(new TestInterfaceTable());
    Agent::GetInstance()->set_interface_table(interface_table.get());
    SendHttpClientRequest(GET_METHOD);
    METADATA_CHECK (stats.responses < 1);
    Agent::GetInstance()->set_interface_table(intf_table);
    EXPECT_EQ(2U, stats.requests);
    EXPECT_EQ(1U, stats.proxy_sessions);
    EXPECT_EQ(1U, stats.internal_errors);

    client->Reset();
    StopHttpClient();
    DeleteVmportEnv(input, 1, 1, 0);
    client->WaitForIdle();

    ClearLinkLocalConfig();
    StopNovaApiProxy();
    client->WaitForIdle();

    Agent::GetInstance()->services()->metadataproxy()->ClearStats();
}

// test multiple client requests in the same session
TEST_F(MetadataTest, MetadataMultiReqTest) {
    int count = 0;
    MetadataProxy::MetadataStats stats;
    struct PortInfo input[] = {
        {"vnet1", 1, vm1_ip, "00:00:00:01:01:01", 1, 1},
    };

    StartNovaApiProxy();
    SetupLinkLocalConfig();

    CreateVmportEnv(input, 1, 0);
    client->WaitForIdle();
    client->Reset();

    // for agent to identify the vm, the remote end should have vm's ip;
    // overload the FindVmUuidFromMetadataIp to return true
    InterfaceTable *intf_table = Agent::GetInstance()->interface_table();
    std::auto_ptr<InterfaceTable> interface_table(new TestInterfaceTable());
    Agent::GetInstance()->set_interface_table(interface_table.get());

    CurlRun();

    METADATA_CHECK (stats.responses < 12);
    Agent::GetInstance()->set_interface_table(intf_table);
    EXPECT_EQ(2U, stats.requests);
    EXPECT_EQ(12U, stats.responses);
    EXPECT_EQ(1U, stats.proxy_sessions);
    EXPECT_EQ(0U, stats.internal_errors);

    CurlStop();
    client->WaitForIdle();

    client->Reset();
    DeleteVmportEnv(input, 1, 1, 0);
    client->WaitForIdle();

    ClearLinkLocalConfig();
    StopNovaApiProxy();
    client->WaitForIdle();

    Agent::GetInstance()->services()->metadataproxy()->ClearStats();
}

TEST_F(MetadataTest, MetadataReqTest_services_ip) {
    int count = 0;
    MetadataProxy::MetadataStats stats;
    struct PortInfo input[] = {
        {"vnet1", 1, vm1_ip, "00:00:00:01:01:01", 1, 1},
    };

    // Change compute_node_ip
    Ip4Address old_ip = agent_->compute_node_ip();
    AgentParam param(agent_);
    param.BuildAddressList("5.5.5.5");
    TestAgentInit *init = static_cast<TestAgentInit *>(client->agent_init());
    init->ProcessComputeAddress(&param);

    StartNovaApiProxy();
    SetupLinkLocalConfig();

    CreateVmportEnv(input, 1, 0);
    client->WaitForIdle();
    client->Reset();

    StartHttpClient();

    // If the local address is not same as VM address in this request,
    // agent will have an internal error as it cannot find the VM
    HttpConnection *conn = SendHttpClientRequest(GET_METHOD);
    METADATA_CHECK (stats.internal_errors < 1);
    EXPECT_EQ(1U, stats.requests);
    EXPECT_EQ(0U, stats.responses);
    EXPECT_EQ(0U, stats.proxy_sessions);
    CloseClientSession(conn);
    client->WaitForIdle();

    VmInterface *intf = static_cast<VmInterface *>(VmPortGet(1));
    TxTcpPacket(intf->id(), vm1_ip, "169.254.169.254", 1000, 80, false);
    client->WaitForIdle();
    FlowEntry *flow = FlowGet(0, vm1_ip, "169.254.169.254", 6, 1000, 80,
                              intf->flow_key_nh()->id());
    EXPECT_TRUE(flow != NULL);
    FlowEntry *rflow = flow->reverse_flow_entry();
    EXPECT_TRUE(rflow != NULL);

    EXPECT_TRUE(flow->key().src_addr.to_v4() ==
                Ip4Address::from_string(vm1_ip));
    EXPECT_TRUE(flow->key().dst_addr.to_v4() ==
                Ip4Address::from_string("169.254.169.254"));

    EXPECT_TRUE(rflow->key().src_addr.to_v4() ==
                Ip4Address::from_string("5.5.5.5"));
    EXPECT_TRUE(rflow->key().dst_addr.to_v4() == intf->mdata_ip_addr());
    client->Reset();
    StopHttpClient();
    DeleteVmportEnv(input, 1, 1, 0);
    client->WaitForIdle();

    ClearLinkLocalConfig();
    StopNovaApiProxy();
    client->WaitForIdle();
    agent_->set_compute_node_ip(old_ip);

    Agent::GetInstance()->services()->metadataproxy()->ClearStats();
}

// Send PUT / POST / HEAD / DELETE requests
TEST_F(MetadataTest, MetadataOtherMethodsTest) {
    int count = 0;
    MetadataProxy::MetadataStats stats;
    struct PortInfo input[] = {
        {"vnet1", 1, vm1_ip, "00:00:00:01:01:01", 1, 1},
    };

    StartNovaApiProxy();
    SetupLinkLocalConfig();

    CreateVmportEnv(input, 1, 0);
    client->WaitForIdle();
    client->Reset();

    StartHttpClient();

    // for agent to identify the vm, the remote end should have vm's ip;
    // overload the FindVmUuidFromMetadataIp to return true
    InterfaceTable *intf_table = Agent::GetInstance()->interface_table();
    std::auto_ptr<InterfaceTable> interface_table(new TestInterfaceTable());
    Agent::GetInstance()->set_interface_table(interface_table.get());
    SendHttpClientRequest(POST_METHOD);
    METADATA_CHECK (stats.responses < 5);
    EXPECT_EQ(1U, stats.requests);
    EXPECT_EQ(5U, stats.responses);
    EXPECT_EQ(1U, stats.proxy_sessions);
    EXPECT_EQ(0U, stats.internal_errors);

    // test with large data, that goes in multiple packets
    std::string large_data;
    for (int i = 0; i < 100; i++) {
        large_data.append("add more data to be sent");
    }
    SendHttpClientRequest(POST_METHOD, large_data);
    METADATA_CHECK (stats.responses < 10);
    EXPECT_EQ(2U, stats.requests);
    EXPECT_EQ(10U, stats.responses);
    EXPECT_EQ(2U, stats.proxy_sessions);
    EXPECT_EQ(0U, stats.internal_errors);

    SendHttpClientRequest(PUT_METHOD);
    METADATA_CHECK (stats.responses < 15);
    EXPECT_EQ(3U, stats.requests);
    EXPECT_EQ(15U, stats.responses);
    EXPECT_EQ(3U, stats.proxy_sessions);
    EXPECT_EQ(0U, stats.internal_errors);

    SendHttpClientRequest(HEAD_METHOD);
    METADATA_CHECK (stats.responses < 20);
    EXPECT_EQ(4U, stats.requests);
    EXPECT_EQ(20U, stats.responses);
    EXPECT_EQ(4U, stats.proxy_sessions);
    EXPECT_EQ(0U, stats.internal_errors);

    SendHttpClientRequest(DELETE_METHOD);
    METADATA_CHECK (stats.responses < 25);
    Agent::GetInstance()->set_interface_table(intf_table);
    EXPECT_EQ(5U, stats.requests);
    EXPECT_EQ(25U, stats.responses);
    EXPECT_EQ(5U, stats.proxy_sessions);
    EXPECT_EQ(0U, stats.internal_errors);

    client->Reset();
    StopHttpClient();
    DeleteVmportEnv(input, 1, 1, 0);
    client->WaitForIdle();

    ClearLinkLocalConfig();
    StopNovaApiProxy();
    client->WaitForIdle();

    Agent::GetInstance()->services()->metadataproxy()->ClearStats();
}

// Send request without linklocal metadata
TEST_F(MetadataTest, MetadataNoLinkLocalTest) {
    int count = 0;
    MetadataProxy::MetadataStats stats;
    struct PortInfo input[] = {
        {"vnet1", 1, vm1_ip, "00:00:00:01:01:01", 1, 1},
    };

    StartNovaApiProxy();

    CreateVmportEnv(input, 1, 0);
    client->WaitForIdle();
    client->Reset();

    StartHttpClient();

    // for agent to identify the vm, the remote end should have vm's ip;
    // overload the FindVmUuidFromMetadataIp to return true
    InterfaceTable *intf_table = Agent::GetInstance()->interface_table();
    std::auto_ptr<InterfaceTable> interface_table(new TestInterfaceTable());
    Agent::GetInstance()->set_interface_table(interface_table.get());
    SendHttpClientRequest(GET_METHOD);
    METADATA_CHECK (stats.internal_errors < 1);
    Agent::GetInstance()->set_interface_table(intf_table);
    EXPECT_EQ(1U, stats.requests);
    EXPECT_EQ(0U, stats.responses);
    EXPECT_EQ(0U, stats.proxy_sessions);
    EXPECT_EQ(1U, stats.internal_errors);

    client->Reset();
    StopHttpClient();
    DeleteVmportEnv(input, 1, 1, 0);
    client->WaitForIdle();

    StopNovaApiProxy();
    client->WaitForIdle();

    Agent::GetInstance()->services()->metadataproxy()->ClearStats();
}

// Send message and close server connection while message is going
TEST_F(MetadataTest, MetadataCloseServerTest) {
    MetadataProxy::MetadataStats stats;
    struct PortInfo input[] = {
        {"vnet1", 1, vm1_ip, "00:00:00:01:01:01", 1, 1},
    };

    StartNovaApiProxy();
    SetupLinkLocalConfig();

    CreateVmportEnv(input, 1, 0);
    client->WaitForIdle();
    client->Reset();

    StartHttpClient();

    InterfaceTable *intf_table = Agent::GetInstance()->interface_table();
    std::auto_ptr<InterfaceTable> interface_table(new TestInterfaceTable());
    Agent::GetInstance()->set_interface_table(interface_table.get());
    std::string large_data;
    for (int i = 0; i < 200; i++) {
        large_data.append("add more data to be sent");
    }
    // stop server
    StopNovaApiProxy();
    client->WaitForIdle();

    Agent::GetInstance()->set_interface_table(intf_table);
    client->Reset();
    StopHttpClient();
    DeleteVmportEnv(input, 1, 1, 0);
    client->WaitForIdle();

    ClearLinkLocalConfig();
    client->WaitForIdle();

    Agent::GetInstance()->services()->metadataproxy()->ClearStats();
}

void RouterIdDepInit(Agent *agent) {
}

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

    client = TestInit(init_file, ksync_init, true, true, false);
    usleep(100000);
    client->WaitForIdle();

    int ret = RUN_ALL_TESTS();
    TestShutdown();
    delete client;
    return ret;
}

ext: cc

max_stars_repo_name: codilime/contrail-controller-arch

max_issues_repo_name: chnyda/contrail-controller

max_forks_repo_name: chnyda/contrail-controller


text:
// Warning! This file is autogenerated.
#include <boost/text/collation_table.hpp>
#include <boost/text/collate.hpp>
#include <boost/text/data/all.hpp>

#ifndef LIMIT_TESTING_FOR_CI
#include <boost/text/save_load_table.hpp>

#include <boost/filesystem.hpp>
#endif

#include <gtest/gtest.h>

using namespace boost::text;

auto const error = [](string const & s) { std::cout << s; };
auto const warning = [](string const & s) {};

collation_table make_save_load_table()
{
#ifdef LIMIT_TESTING_FOR_CI
    string const table_str(data::ha::standard_collation_tailoring());
    return tailored_collation_table(
        table_str,
        "ha::standard_collation_tailoring()", error, warning);
#else
    if (!exists(boost::filesystem::path("ha_standard.table"))) {
        string const table_str(data::ha::standard_collation_tailoring());
        collation_table table = tailored_collation_table(
            table_str,
            "ha::standard_collation_tailoring()", error, warning);
        save_table(table, "ha_standard.table.0");
        boost::filesystem::rename("ha_standard.table.0", "ha_standard.table");
    }
    return load_table("ha_standard.table");
#endif
}
collation_table const & table()
{
    static collation_table retval = make_save_load_table();
    return retval;
}
TEST(tailoring, ha_standard_000_001)
{
    {
    // greater than (or equal to, for =) preceeding cps
    auto const res = std::vector<uint32_t>(1, 0x0042);
    auto const rel = std::vector<uint32_t>(1, 0x0253);
    string const res_str = to_string(res);
    string const rel_str = to_string(rel);
    auto const res_view = as_utf32(res);
    auto const rel_view = as_utf32(rel);
    EXPECT_EQ(collate(
        res.begin(), res.end(),
        rel.begin(), rel.end(),
        table(), collation_strength::primary),
        -1);
    EXPECT_EQ(collate(
        res_view.begin(), res_view.end(),
        rel_view.begin(), rel_view.end(),
        table(), collation_strength::primary),
        -1);
    }
    {
    // greater than (or equal to, for =) preceeding cps
    auto const res = std::vector<uint32_t>(1, 0x0253);
    auto const rel = std::vector<uint32_t>(1, 0x0181);
    string const res_str = to_string(res);
    string const rel_str = to_string(rel);
    auto const res_view = as_utf32(res);
    auto const rel_view = as_utf32(rel);
    EXPECT_EQ(collate(
        res.begin(), res.end(),
        rel.begin(), rel.end(),
        table(), collation_strength::tertiary),
        -1);
    EXPECT_EQ(collate(
        res_view.begin(), res_view.end(),
        rel_view.begin(), rel_view.end(),
        table(), collation_strength::tertiary),
        -1);
    // equal to preceeding cps at next-lower strength
    EXPECT_EQ(collate(
        res.begin(), res.end(),
        rel.begin(), rel.end(),
        table(), collation_strength::secondary),
        0);
    // equal to preceeding cps at next-lower strength
    EXPECT_EQ(collate(
        res_view.begin(), res_view.end(),
        rel_view.begin(), rel_view.end(),
        table(), collation_strength::secondary),
        0);
    }
    {
    // greater than (or equal to, for =) preceeding cps
    auto const res = std::vector<uint32_t>(1, 0x0044);
    auto const rel = std::vector<uint32_t>(1, 0x0257);
    string const res_str = to_string(res);
    string const rel_str = to_string(rel);
    auto const res_view = as_utf32(res);
    auto const rel_view = as_utf32(rel);
    EXPECT_EQ(collate(
        res.begin(), res.end(),
        rel.begin(), rel.end(),
        table(), collation_strength::primary),
        -1);
    EXPECT_EQ(collate(
        res_view.begin(), res_view.end(),
        rel_view.begin(), rel_view.end(),
        table(), collation_strength::primary),
        -1);
    }
    {
    // greater than (or equal to, for =) preceeding cps
    auto const res = std::vector<uint32_t>(1, 0x0257);
    auto const rel = std::vector<uint32_t>(1, 0x018a);
    string const res_str = to_string(res);
    string const rel_str = to_string(rel);
    auto const res_view = as_utf32(res);
    auto const rel_view = as_utf32(rel);
    EXPECT_EQ(collate(
        res.begin(), res.end(),
        rel.begin(), rel.end(),
        table(), collation_strength::tertiary),
        -1);
    EXPECT_EQ(collate(
        res_view.begin(), res_view.end(),
        rel_view.begin(), rel_view.end(),
        table(), collation_strength::tertiary),
        -1);
    // equal to preceeding cps at next-lower strength
    EXPECT_EQ(collate(
        res.begin(), res.end(),
        rel.begin(), rel.end(),
        table(), collation_strength::secondary),
        0);
    // equal to preceeding cps at next-lower strength
    EXPECT_EQ(collate(
        res_view.begin(), res_view.end(),
        rel_view.begin(), rel_view.end(),
        table(), collation_strength::secondary),
        0);
    }
    {
    // greater than (or equal to, for =) preceeding cps
    auto const res = std::vector<uint32_t>(1, 0x004b);
    auto const rel = std::vector<uint32_t>(1, 0x0199);
    string const res_str = to_string(res);
    string const rel_str = to_string(rel);
    auto const res_view = as_utf32(res);
    auto const rel_view = as_utf32(rel);
    EXPECT_EQ(collate(
        res.begin(), res.end(),
        rel.begin(), rel.end(),
        table(), collation_strength::primary),
        -1);
    EXPECT_EQ(collate(
        res_view.begin(), res_view.end(),
        rel_view.begin(), rel_view.end(),
        table(), collation_strength::primary),
        -1);
    }
    {
    // greater than (or equal to, for =) preceeding cps
    auto const res = std::vector<uint32_t>(1, 0x0199);
    auto const rel = std::vector<uint32_t>(1, 0x0198);
    string const res_str = to_string(res);
    string const rel_str = to_string(rel);
    auto const res_view = as_utf32(res);
    auto const rel_view = as_utf32(rel);
    EXPECT_EQ(collate(
        res.begin(), res.end(),
        rel.begin(), rel.end(),
        table(), collation_strength::tertiary),
        -1);
    EXPECT_EQ(collate(
        res_view.begin(), res_view.end(),
        rel_view.begin(), rel_view.end(),
        table(), collation_strength::tertiary),
        -1);
    // equal to preceeding cps at next-lower strength
    EXPECT_EQ(collate(
        res.begin(), res.end(),
        rel.begin(), rel.end(),
        table(), collation_strength::secondary),
        0);
    // equal to preceeding cps at next-lower strength
    EXPECT_EQ(collate(
        res_view.begin(), res_view.end(),
        rel_view.begin(), rel_view.end(),
        table(), collation_strength::secondary),
        0);
    }
    {
    // greater than (or equal to, for =) preceeding cps
    auto const res = std::vector<uint32_t>(1, 0x0053);
    auto const rel = std::vector<uint32_t>{0x0073, 0x0068};
    string const res_str = to_string(res);
    string const rel_str = to_string(rel);
    auto const res_view = as_utf32(res);
    auto const rel_view = as_utf32(rel);
    EXPECT_EQ(collate(
        res.begin(), res.end(),
        rel.begin(), rel.end(),
        table(), collation_strength::primary),
        -1);
    EXPECT_EQ(collate(
        res_view.begin(), res_view.end(),
        rel_view.begin(), rel_view.end(),
        table(), collation_strength::primary),
        -1);
    }
    {
    // greater than (or equal to, for =) preceeding cps
    auto const res = std::vector<uint32_t>{0x0073, 0x0068};
    auto const rel = std::vector<uint32_t>{0x0053, 0x0068};
    string const res_str = to_string(res);
    string const rel_str = to_string(rel);
    auto const res_view = as_utf32(res);
    auto const rel_view = as_utf32(rel);
    EXPECT_EQ(collate(
        res.begin(), res.end(),
        rel.begin(), rel.end(),
        table(), collation_strength::tertiary),
        -1);
    EXPECT_EQ(collate(
        res_view.begin(), res_view.end(),
        rel_view.begin(), rel_view.end(),
        table(), collation_strength::tertiary),
        -1);
    // equal to preceeding cps at next-lower strength
    EXPECT_EQ(collate(
        res.begin(), res.end(),
        rel.begin(), rel.end(),
        table(), collation_strength::secondary),
        0);
    // equal to preceeding cps at next-lower strength
    EXPECT_EQ(collate(
        res_view.begin(), res_view.end(),
        rel_view.begin(), rel_view.end(),
        table(), collation_strength::secondary),
        0);
    }
    {
    // greater than (or equal to, for =) preceeding cps
    auto const res = std::vector<uint32_t>{0x0053, 0x0068};
    auto const rel = std::vector<uint32_t>{0x0053, 0x0048};
    string const res_str = to_string(res);
    string const rel_str = to_string(rel);
    auto const res_view = as_utf32(res);
    auto const rel_view = as_utf32(rel);
    EXPECT_EQ(collate(
        res.begin(), res.end(),
        rel.begin(), rel.end(),
        table(), collation_strength::tertiary),
        -1);
    EXPECT_EQ(collate(
        res_view.begin(), res_view.end(),
        rel_view.begin(), rel_view.end(),
        table(), collation_strength::tertiary),
        -1);
    // equal to preceeding cps at next-lower strength
    EXPECT_EQ(collate(
        res.begin(), res.end(),
        rel.begin(), rel.end(),
        table(), collation_strength::secondary),
        0);
    // equal to preceeding cps at next-lower strength
    EXPECT_EQ(collate(
        res_view.begin(), res_view.end(),
        rel_view.begin(), rel_view.end(),
        table(), collation_strength::secondary),
        0);
    }
    {
    // greater than (or equal to, for =) preceeding cps
    auto const res = std::vector<uint32_t>(1, 0x0054);
    auto const rel = std::vector<uint32_t>{0x0074, 0x0073};
    string const res_str = to_string(res);
    string const rel_str = to_string(rel);
    auto const res_view = as_utf32(res);
    auto const rel_view = as_utf32(rel);
    EXPECT_EQ(collate(
        res.begin(), res.end(),
        rel.begin(), rel.end(),
        table(), collation_strength::primary),
        -1);
    EXPECT_EQ(collate(
        res_view.begin(), res_view.end(),
        rel_view.begin(), rel_view.end(),
        table(), collation_strength::primary),
        -1);
    }
    {
    // greater than (or equal to, for =) preceeding cps
    auto const res = std::vector<uint32_t>{0x0074, 0x0073};
    auto const rel = std::vector<uint32_t>{0x0054, 0x0073};
    string const res_str = to_string(res);
    string const rel_str = to_string(rel);
    auto const res_view = as_utf32(res);
    auto const rel_view = as_utf32(rel);
    EXPECT_EQ(collate(
        res.begin(), res.end(),
        rel.begin(), rel.end(),
        table(), collation_strength::tertiary),
        -1);
    EXPECT_EQ(collate(
        res_view.begin(), res_view.end(),
        rel_view.begin(), rel_view.end(),
        table(), collation_strength::tertiary),
        -1);
    // equal to preceeding cps at next-lower strength
    EXPECT_EQ(collate(
        res.begin(), res.end(),
        rel.begin(), rel.end(),
        table(), collation_strength::secondary),
        0);
    // equal to preceeding cps at next-lower strength
    EXPECT_EQ(collate(
        res_view.begin(), res_view.end(),
        rel_view.begin(), rel_view.end(),
        table(), collation_strength::secondary),
        0);
    }
    {
    // greater than (or equal to, for =) preceeding cps
    auto const res = std::vector<uint32_t>{0x0054, 0x0073};
    auto const rel = std::vector<uint32_t>{0x0054, 0x0053};
    string const res_str = to_string(res);
    string const rel_str = to_string(rel);
    auto const res_view = as_utf32(res);
    auto const rel_view = as_utf32(rel);
    EXPECT_EQ(collate(
        res.begin(), res.end(),
        rel.begin(), rel.end(),
        table(), collation_strength::tertiary),
        -1);
    EXPECT_EQ(collate(
        res_view.begin(), res_view.end(),
        rel_view.begin(), rel_view.end(),
        table(), collation_strength::tertiary),
        -1);
    // equal to preceeding cps at next-lower strength
    EXPECT_EQ(collate(
        res.begin(), res.end(),
        rel.begin(), rel.end(),
        table(), collation_strength::secondary),
        0);
    // equal to preceeding cps at next-lower strength
    EXPECT_EQ(collate(
        res_view.begin(), res_view.end(),
        rel_view.begin(), rel_view.end(),
        table(), collation_strength::secondary),
        0);
    }
    {
    // greater than (or equal to, for =) preceeding cps
    auto const res = std::vector<uint32_t>(1, 0x0059);
    auto const rel = std::vector<uint32_t>(1, 0x01b4);
    string const res_str = to_string(res);
    string const rel_str = to_string(rel);
    auto const res_view = as_utf32(res);
    auto const rel_view = as_utf32(rel);
    EXPECT_EQ(collate(
        res.begin(), res.end(),
        rel.begin(), rel.end(),
        table(), collation_strength::primary),
        -1);
    EXPECT_EQ(collate(
        res_view.begin(), res_view.end(),
        rel_view.begin(), rel_view.end(),
        table(), collation_strength::primary),
        -1);
    }
    {
    // greater than (or equal to, for =) preceeding cps
    auto const res = std::vector<uint32_t>(1, 0x01b4);
    auto const rel = std::vector<uint32_t>{0x02bc, 0x0079};
    string const res_str = to_string(res);
    string const rel_str = to_string(rel);
    auto const res_view = as_utf32(res);
    auto const rel_view = as_utf32(rel);
    EXPECT_EQ(collate(
        res.begin(), res.end(),
        rel.begin(), rel.end(),
        table(), collation_strength::tertiary),
        -1);
    EXPECT_EQ(collate(
        res_view.begin(), res_view.end(),
        rel_view.begin(), rel_view.end(),
        table(), collation_strength::tertiary),
        -1);
    // equal to preceeding cps at next-lower strength
    EXPECT_EQ(collate(
        res.begin(), res.end(),
        rel.begin(), rel.end(),
        table(), collation_strength::secondary),
        0);
    // equal to preceeding cps at next-lower strength
    EXPECT_EQ(collate(
        res_view.begin(), res_view.end(),
        rel_view.begin(), rel_view.end(),
        table(), collation_strength::secondary),
        0);
    }
    {
    // greater than (or equal to, for =) preceeding cps
    auto const res = std::vector<uint32_t>{0x02bc, 0x0079};
    auto const rel = std::vector<uint32_t>{0x0027, 0x0079};
    string const res_str = to_string(res);
    string const rel_str = to_string(rel);
    auto const res_view = as_utf32(res);
    auto const rel_view = as_utf32(rel);
    EXPECT_EQ(collate(
        res.begin(), res.end(),
        rel.begin(), rel.end(),
        table(), collation_strength::tertiary),
        -1);
    EXPECT_EQ(collate(
        res_view.begin(), res_view.end(),
        rel_view.begin(), rel_view.end(),
        table(), collation_strength::tertiary),
        -1);
    // equal to preceeding cps at next-lower strength
    EXPECT_EQ(collate(
        res.begin(), res.end(),
        rel.begin(), rel.end(),
        table(), collation_strength::secondary),
        0);
    // equal to preceeding cps at next-lower strength
    EXPECT_EQ(collate(
        res_view.begin(), res_view.end(),
        rel_view.begin(), rel_view.end(),
        table(), collation_strength::secondary),
        0);
    }
    {
    // greater than (or equal to, for =) preceeding cps
    auto const res = std::vector<uint32_t>{0x0027, 0x0079};
    auto const rel = std::vector<uint32_t>(1, 0x01b3);
    string const res_str = to_string(res);
    string const rel_str = to_string(rel);
    auto const res_view = as_utf32(res);
    auto const rel_view = as_utf32(rel);
    EXPECT_EQ(collate(
        res.begin(), res.end(),
        rel.begin(), rel.end(),
        table(), collation_strength::tertiary),
        -1);
    EXPECT_EQ(collate(
        res_view.begin(), res_view.end(),
        rel_view.begin(), rel_view.end(),
        table(), collation_strength::tertiary),
        -1);
    // equal to preceeding cps at next-lower strength
    EXPECT_EQ(collate(
        res.begin(), res.end(),
        rel.begin(), rel.end(),
        table(), collation_strength::secondary),
        0);
    // equal to preceeding cps at next-lower strength
    EXPECT_EQ(collate(
        res_view.begin(), res_view.end(),
        rel_view.begin(), rel_view.end(),
        table(), collation_strength::secondary),
        0);
    }
    {
    // greater than (or equal to, for =) preceeding cps
    auto const res = std::vector<uint32_t>(1, 0x01b3);
    auto const rel = std::vector<uint32_t>{0x02bc, 0x0059};
    string const res_str = to_string(res);
    string const rel_str = to_string(rel);
    auto const res_view = as_utf32(res);
    auto const rel_view = as_utf32(rel);
    EXPECT_EQ(collate(
        res.begin(), res.end(),
        rel.begin(), rel.end(),
        table(), collation_strength::tertiary),
        -1);
    EXPECT_EQ(collate(
        res_view.begin(), res_view.end(),
        rel_view.begin(), rel_view.end(),
        table(), collation_strength::tertiary),
        -1);
    // equal to preceeding cps at next-lower strength
    EXPECT_EQ(collate(
        res.begin(), res.end(),
        rel.begin(), rel.end(),
        table(), collation_strength::secondary),
        0);
    // equal to preceeding cps at next-lower strength
    EXPECT_EQ(collate(
        res_view.begin(), res_view.end(),
        rel_view.begin(), rel_view.end(),
        table(), collation_strength::secondary),
        0);
    }
    {
    // greater than (or equal to, for =) preceeding cps
    auto const res = std::vector<uint32_t>{0x02bc, 0x0059};
    auto const rel = std::vector<uint32_t>{0x0027, 0x0059};
    string const res_str = to_string(res);
    string const rel_str = to_string(rel);
    auto const res_view = as_utf32(res);
    auto const rel_view = as_utf32(rel);
    EXPECT_EQ(collate(
        res.begin(), res.end(),
        rel.begin(), rel.end(),
        table(), collation_strength::tertiary),
        -1);
    EXPECT_EQ(collate(
        res_view.begin(), res_view.end(),
        rel_view.begin(), rel_view.end(),
        table(), collation_strength::tertiary),
        -1);
    // equal to preceeding cps at next-lower strength
    EXPECT_EQ(collate(
        res.begin(), res.end(),
        rel.begin(), rel.end(),
        table(), collation_strength::secondary),
        0);
    // equal to preceeding cps at next-lower strength
    EXPECT_EQ(collate(
        res_view.begin(), res_view.end(),
        rel_view.begin(), rel_view.end(),
        table(), collation_strength::secondary),
        0);
    }
}

ext: cpp

max_stars_repo_name: eightysquirrels/text

max_issues_repo_name: eightysquirrels/text

max_forks_repo_name: eightysquirrels/text


text:
/*
 Copyright (C) 2016 Quaternion Risk Management Ltd
 All rights reserved.

 This file is part of ORE, a free-software/open-source library
 for transparent pricing and risk analysis - http://opensourcerisk.org

 ORE is free software: you can redistribute it and/or modify it
 under the terms of the Modified BSD License.  You should have received a
 copy of the license along with this program.
 The license is also available online at <http://opensourcerisk.org>

 This program is distributed on the basis that it will form a useful
 contribution to risk analytics and model standardisation, but WITHOUT
 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 FITNESS FOR A PARTICULAR PURPOSE. See the license for more details.
*/

/*! \file test/scenariogenerator.hpp
    \brief scenario generation test
    \ingroup tests
*/

#pragma once

#include <boost/test/unit_test.hpp>

namespace testsuite {

//! Test scenario generation
/*!
  \ingroup tests
*/
class ScenarioGeneratorTest {
public:
    //! Test consistency of CPI Swap Exposure evolutions with Cpi Floor prices
    static void testCpiSwapExposure();
    //! Test Martingale properties for single currency LGM scenarios with Pseudo Random numbers
    static void testLgmMersenneTwister();
    //! Test Martingale properties for single currency LGM scenarios with Pseudo Random numbers and antithetic sampling
    static void testLgmMersenneTwisterAntithetic();
    //! Test Martingale properties for single currency LGM scenarios with Sobol sequences
    static void testLgmLowDiscrepancy();
    //! Test Martingale properties for single currency LGM scenarios with Sobol sequences and Brownian bridge
    static void testLgmLowDiscrepancyBrownianBridge();
    //! Test Martingale properties for cross asset IR/FX model with Pseudo Random numbers
    static void testCrossAssetMersenneTwister();
    //! Test Martingale properties for cross asset IR/FX model with Pseudo Random numbers and antithetic sampling
    static void testCrossAssetMersenneTwisterAntithetic();
    //! Test Martingale properties for cross asset IR/FX model with Sobol sequences
    static void testCrossAssetLowDiscrepancy();
    //! Test Martingale properties for cross asset IR/FX model with Sobol sequences and Brownian bridge
    static void testCrossAssetLowDiscrepancyBrownianBridge();
    //! Test Martingale properties for cross asset IR/FX model scenarios after applying them to the simulation market
    static void testCrossAssetSimMarket();
    //! Test cross asset IR/FX model scenarios after applying to simulation market pathwise to the direct model output
    static void testCrossAssetSimMarket2();
    //! Test consistency of EUR and USD swap exposure evolutions with European Swaption prices
    static void testVanillaSwapExposure();
    //! Test consistency of FX Forward Exposure evolutions with FX Option prices
    static void testFxForwardExposure();
    //! Test consistency of FX Forward Exposure evolutions with FX Forward Option prices when IR vols vanish
    static void testFxForwardExposureZeroIrVol();

    static boost::unit_test_framework::test_suite* suite();
};
} // namespace testsuite

ext: hpp

max_stars_repo_name: mrslezak/Engine

max_issues_repo_name: mrslezak/Engine

max_forks_repo_name: mrslezak/Engine


text:
#ifndef BOOST_METAPARSE_GETTING_STARTED_11_HPP
#define BOOST_METAPARSE_GETTING_STARTED_11_HPP

// Automatically generated header file

// Definitions before section 10
#include "10.hpp"

// Definitions of section 10
using lparen_token = token<lit_c<'('>>;

using rparen_token = token<lit_c<')'>>;

using plus_exp1 =
 foldl_start_with_parser<
   sequence<one_of<plus_token, minus_token>, mult_exp4>,
   mult_exp4,
   boost::mpl::quote2<binary_op>
 >;

using paren_exp1 = sequence<lparen_token, plus_exp1, rparen_token>;

#include <boost/metaparse/middle_of.hpp>

using paren_exp2 = middle_of<lparen_token, plus_exp1, rparen_token>;

using primary_exp1 = one_of<int_token, paren_exp2>;

struct plus_exp2;

using paren_exp3 = middle_of<lparen_token, plus_exp2, rparen_token>;

using primary_exp2 = one_of<int_token, paren_exp2>;

using unary_exp2 =
 foldr_start_with_parser<
   minus_token,
   primary_exp2,
   boost::mpl::lambda<boost::mpl::negate<boost::mpl::_1>>::type
 >;

using mult_exp5 =
 foldl_start_with_parser<
   sequence<one_of<times_token, divides_token>, unary_exp2>,
   unary_exp2,
   boost::mpl::quote2<binary_op>
 >;

struct plus_exp2 :
 foldl_start_with_parser<
   sequence<one_of<plus_token, minus_token>, mult_exp5>,
   mult_exp5,
   boost::mpl::quote2<binary_op>
 > {};

using exp_parser19 = build_parser<plus_exp2>;

// query:
//    exp_parser19::apply<BOOST_METAPARSE_STRING("(1 + 2) * 3")>::type

#endif


text:
/*
 * Copyright (c) 2020-2021, Marco Sánchez Beeckman
 * All rights reserved.
 *
 * This source code is licensed under the BSD-style license found in the
 * LICENSE file in the root directory of this source tree.
 */

#include <algorithm>
#include <cmath>
#include <chrono>
#include <iostream>

#include <Eigen/Dense>
#include <Eigen/Sparse>
#include <opencv2/core.hpp>

#include "RandomWalker.h"

namespace tfg {

/**
 * Construct an instance of a RandomWalker object. The random walk is performed over a graph
 * where the nodes are tracks, so a trackTable object is required to initialize the RandomWalker.
 * @param trackTable A trackTable containing information of the tracks where the random walk is going to be performed.
 */
RandomWalker::RandomWalker(std::shared_ptr<tfg::TrackTable> &trackTable) {
    this->trackTable = trackTable;
}

RandomWalker::~RandomWalker() {}

/**
 * Label nodes (tracks) of the graph as foreground (there can be various foreground labels) or background from a series of seeded images.
 * An initial labeling is needed so that the labels can be propagated to unlabeled nodes.
 * @param seedImages An unordered map of OpenCV Mat objects, whose key is an integer that indicates the frame the seed is for.
 */
void RandomWalker::seed(const std::unordered_map<int, cv::Mat> &seedImages) {
    std::cout << "Seeding tracks" << '\n';

    // For each matrix of seeds, label each track passing through a seed with its corresponding tag.
    // If a track passes through two or more different seeds, tag it as unreliable and don't use it for propagation.
    for(auto& [ frame, image ] : seedImages) {
        for(unsigned int t = 0; t < trackTable->numberOfTracks(); t++) {
            const int trackLabel = trackTable->labelOfTrack(t);
            if(trackLabel == -2) continue; // If unreliable, do nothing to the track
            if(frame < trackTable->firstFrameOfTrack(t) || frame > trackTable->firstFrameOfTrack(t) + trackTable->durationOfTrack(t) - 1) continue;

            const cv::Vec2f position = trackTable->pointsInTrack(t)[frame - trackTable->firstFrameOfTrack(t)];
            const int posX = std::round(position(0));
            const int posY = std::round(position(1));
            const cv::Vec3b& color = image.at<cv::Vec3b>(posY, posX);

            // The labeling in the images must be done using colors the following way:
            // a BGR value of 0-0-255 (full red) corresponds to the background
            // a BGR value of 0-0-0 (no color) means there is no label
            // another BGR value using 255 for any combination of colors (e. g. 0-255-255) is a different foreground label
            // non-255 values are considered 0.
            int label = color[2]/255 + 2*(color[1]/255) + 4*(color[0]/255) - 1;
            int labelWithValidityCheck = trackLabel < 0 || trackLabel == label ? label : -2; // If a track is already labeled and it passes through a seed with a different tag, mark it as unreliable
            trackTable->setLabelToTrack(labelWithValidityCheck, t);
        }
    }

    // Sort tracks by label. The only purpose of this is to separate labeled and unlabeled tracks, so it is easier to mount the
    // system of equations for the propagation
    trackTable->sortTracksByLabel();
    this->numberOfLabels = trackTable->labelOfTrack(trackTable->numberOfTracks() - 1) + 1;
    this->unlabeledTracks = trackTable->indexOfFirstLabel();
    this->labeledTracks = trackTable->numberOfTracks() - unlabeledTracks;

    std::cout << "Sorted tracks according to labels" << '\n';
    std::cout << "There are " << unlabeledTracks << " unlabeled tracks and " << labeledTracks << " labeled tracks" << '\n';
}

void RandomWalker::seedDavis(const std::unordered_map<int, cv::Mat> &seedImages) {
    std::cout << "Seeding tracks" << '\n';

    // For each matrix of seeds, label each track passing through a seed with its corresponding tag.
    // If a track passes through two or more different seeds, tag it as unreliable and don't use it for propagation.
    for(auto& [ frame, image ] : seedImages) {
        for(unsigned int t = 0; t < trackTable->numberOfTracks(); t++) {
            const int trackLabel = trackTable->labelOfTrack(t);
            if(trackLabel == -2) continue; // If unreliable, do nothing to the track
            if(frame < trackTable->firstFrameOfTrack(t) || frame > trackTable->firstFrameOfTrack(t) + trackTable->durationOfTrack(t) - 1) continue;

            const cv::Vec2f position = trackTable->pointsInTrack(t)[frame - trackTable->firstFrameOfTrack(t)];
            const int posX = std::round(position(0));
            const int posY = std::round(position(1));
            const cv::Vec3b& color = image.at<cv::Vec3b>(posY, posX);

            // The labeling in the images must be done using colors the following way:
            // a BGR value of 0-0-0 (black/no color) corresponds to the background
            // a BGR value of 0-0-128 (red) corresponds to the first object
            // any other combination of values using 128 or 255
            int label = color[2]/128 + 2*(color[1]/128) + 4*(color[0]/128);
            int labelWithValidityCheck = trackLabel < 0 || trackLabel == label ? label : -2; // If a track is already labeled and it passes through a seed with a different tag, mark it as unreliable
            trackTable->setLabelToTrack(labelWithValidityCheck, t);
        }
    }

    trackTable->sortTracksByLabel();
    this->numberOfLabels = trackTable->labelOfTrack(trackTable->numberOfTracks() - 1) + 1;
    this->unlabeledTracks = trackTable->indexOfFirstLabel();
    this->labeledTracks = trackTable->numberOfTracks() - unlabeledTracks;

    std::cout << "Sorted tracks according to labels" << '\n';
    std::cout << "There are " << unlabeledTracks << " unlabeled tracks and " << labeledTracks << " labeled tracks" << '\n';
}

/**
 * Propagate the labels of the graph. This is done by solving a system of equations, specifically
 * the one in eq. 10 in the following paper:
 *
 * L. Grady.
 * Random Walks for Image Segmentation
 * IEEE Transactions on Pattern Analysis and Machine Intelligence, 28(11): 1768-1783, 2006
 * DOI: 10.1109/TPAMI.2006.233
 */
void RandomWalker::propagateSeeds(float lambda) {
    std::cout << "Beginning propagation of seeds" << '\n';

    Eigen::SparseMatrix<float> laplacianUnlabeled(unlabeledTracks, unlabeledTracks);
    Eigen::SparseMatrix<float> minusBt(unlabeledTracks, labeledTracks);

    std::cout << "Filling matrices for the system of equations" << '\n';

    std::chrono::steady_clock::time_point beginFill = std::chrono::steady_clock::now();
    std::vector<Eigen::Triplet<float>> laplacianEntries;
    std::vector<Eigen::Triplet<float>> minusBtEntries;
    std::vector<float> unlabeledTracksDegrees(unlabeledTracks, 0.0f);
    laplacianEntries.reserve(unlabeledTracks * unlabeledTracks / 2);
    minusBtEntries.reserve(unlabeledTracks * labeledTracks / 2);

    // Fill the sparse matrices corresponding to L_U and -Bt (minus B transposed) in eq. 7 of Grady's paper
    // L_U is the graph laplacian of the subgraph formed by only the unlabeled tracks
    // B is the matrix of size |V_M| x |V_U| whose entry (i, j) corresponds to minus the weight of the edge from the ith labeled track to the jth unlabeled track
    for(unsigned int tA = 0; tA < unlabeledTracks; tA++) {
        for(unsigned int tB = tA + 1; tB < unlabeledTracks + labeledTracks; tB++) {

            // The weights between affinities are defined from a distance, defined in eq. 2
            // of the following paper:
            // P. Ochs, J. Malik, and T. Brox.
            // Segmentation of moving objects by long term video analysis.
            // IEEE Transactions on Pattern Analysis and Machine Intelligence, 36(6): 1187-1200, Jun 2014.
            // DOI: 10.1109/TPAMI.2013.242
            const float trackDistance2 = trackTable->distance2BetweenTracks(tA, tB);
            const float weightAB = std::exp(-lambda * trackDistance2);

            // If the weight between tracks is almost 0, consider it 0 and do not store it explicitly in the sparse matrix
            // if(weightAB < 1e-8) continue;
            if(weightAB == 0) continue;
            unlabeledTracksDegrees[tA] += weightAB;

            if(tB >= unlabeledTracks) {
                minusBtEntries.push_back(Eigen::Triplet<float>(tA, tB - unlabeledTracks, weightAB));
            } else {
                unlabeledTracksDegrees[tB] += weightAB;
                laplacianEntries.push_back(Eigen::Triplet<float>(tA, tB, -weightAB));
            }
        }
        laplacianEntries.push_back(Eigen::Triplet<float>(tA, tA, unlabeledTracksDegrees[tA]));
    }
    laplacianUnlabeled.setFromTriplets(laplacianEntries.begin(), laplacianEntries.end());
    minusBt.setFromTriplets(minusBtEntries.begin(), minusBtEntries.end());
    std::chrono::steady_clock::time_point endFill = std::chrono::steady_clock::now();
    std::cout << "Filled in " << (std::chrono::duration_cast<std::chrono::microseconds>(endFill-beginFill).count())/1000000.0 << " seconds" << '\n';

    std::cout << "Starting walking process" << '\n';
    std::vector<Eigen::VectorXf> labelProbabilities;
    labelProbabilities.reserve(numberOfLabels - 1);

    // Solve the system of equations L_U*x = -Bt*m_s given in eq. 10 in Grady's paper for each label except the last one
    // m_s (s for each label) is a vector of length |V_M| that serves as an indicator for each label: a 1 in the ith entry means that
    // the ith labeled node is marked with label s, and a 0 means it is not
    for(int k = 0; k < numberOfLabels - 1; k++) {

        Eigen::SparseMatrix<float> M(labeledTracks, 1);
        std::vector<Eigen::Triplet<float>> Mentries;
        Mentries.reserve(labeledTracks);
        // Create indicator for the label
        for(int i = 0; i < labeledTracks; i++) {
            float labelIsPresent = trackTable->labelOfTrack(i + unlabeledTracks) == k ? 1 : 0;
            if(labelIsPresent != 0) Mentries.push_back(Eigen::Triplet<float>(i, 0, labelIsPresent));
        }
        M.setFromTriplets(Mentries.begin(), Mentries.end());

        // Get product -Bt*m_s, right hand of eq. 10 in Grady's paper
        Eigen::SparseMatrix<float> minusBtByM = minusBt * M;
        std::cout << "Solving system of equations for label " << k << '\n';
        std::chrono::steady_clock::time_point beginSolve = std::chrono::steady_clock::now();

        // Solve the system of equations using an iterative conjugate gradient algorithm for self-adjoint sparse matrices
        Eigen::ConjugateGradient<Eigen::SparseMatrix<float>, Eigen::Upper> conjGrad;
        conjGrad.compute(laplacianUnlabeled);
        Eigen::VectorXf probabilityLabel = conjGrad.solve(minusBtByM);
        labelProbabilities.push_back(probabilityLabel);
        std::chrono::steady_clock::time_point endSolve = std::chrono::steady_clock::now();
        std::cout << "Solved in " << (std::chrono::duration_cast<std::chrono::microseconds>(endSolve-beginSolve).count())/1000000.0 << " seconds" << '\n';
    }

    probabilities.clear();
    // probabilities.reserve(unlabeledTracks + labeledTracks);
    probabilities.resize(unlabeledTracks + labeledTracks);

    // For every unlabeled track, save the probabilities for each label
    // and order them the original way
    for(unsigned int t = 0; t < unlabeledTracks; t++) {
        std::vector<float> trackLabelProbabilities;
        trackLabelProbabilities.reserve(numberOfLabels);
        float maxProbability = 0;
        float sumOfProbabilities = 0;
        unsigned int mostLikelyLabel = 0;
        for(unsigned int k = 0; k < numberOfLabels - 1; k++) {
            trackLabelProbabilities.push_back(labelProbabilities[k](t));
            sumOfProbabilities += labelProbabilities[k](t);
            if(labelProbabilities[k](t) > maxProbability) {
                maxProbability = labelProbabilities[k](t);
                mostLikelyLabel = k;
            }
        }
        float probabilityOfLastLabel = 1 - sumOfProbabilities;
        trackLabelProbabilities.push_back(probabilityOfLastLabel);
        if(probabilityOfLastLabel > maxProbability) {
            maxProbability = probabilityOfLastLabel;
            mostLikelyLabel = numberOfLabels - 1;
        }

        probabilities[trackTable->numberOfTrack(t)] = trackLabelProbabilities;
        trackTable->setLabelToTrack(mostLikelyLabel, t);
    }
    // Do the same for the labeled tracks
    for(unsigned int t = unlabeledTracks; t < unlabeledTracks + labeledTracks; t++) {
        std::vector<float> trackLabelProbabilities(numberOfLabels, 0.0f);
        trackLabelProbabilities[trackTable->labelOfTrack(t)] = 1.0f;
        probabilities[trackTable->numberOfTrack(t)] = trackLabelProbabilities;
    }
}

/**
 * Write in a file the probabilities that each track has for every label.
 * @param file The file where the probabilities are going to be stored.
 */
void RandomWalker::writeProbabilities(std::ofstream &file) {
    file << unlabeledTracks + labeledTracks << " " << numberOfLabels << '\n';
    for(unsigned int t = 0; t < unlabeledTracks + labeledTracks; t++) {
        // std::vector<float>& trackLabelProbabilities = probabilities.at(t);
        for(unsigned int k = 0; k < numberOfLabels - 1; k++) {
            // file << trackLabelProbabilities[k] << " ";
            file << probabilities[t][k] << " ";
        }
        // file << trackLabelProbabilities[numberOfLabels - 1] << '\n';
        file << probabilities[t][numberOfLabels - 1] << '\n';
    }
}

} // namespace tfg


text:
#pragma warning (disable: 4244 4018) //coversion from `long` to `float`, possible loss of data

#include <boost/spirit/include/qi.hpp>
#include <boost/spirit/include/karma.hpp>
#include <boost/thread/thread_pool.hpp>

#include <type_traits>

using namespace boost::spirit;

#include <System/Utils.h>

boost::basic_thread_pool g_ThreadPool;

template<typename T>
bool ConvertStrToData(const char *src, T &dest)
{
	const char
		*first_it = src,
		*last_it = first_it + strlen(src);

	return qi::parse(first_it, last_it,
		typename std::conditional<
			std::is_floating_point<T>::value,
				qi::real_parser<T>,
				qi::int_parser<T>
		>::type(),
		dest);
}

template bool ConvertStrToData(const char *src, int &dest);
template bool ConvertStrToData(const char *src, unsigned int &dest);
template bool ConvertStrToData(const char *src, short &dest);
template bool ConvertStrToData(const char *src, unsigned short &dest);
template bool ConvertStrToData(const char *src, char &dest);
template bool ConvertStrToData(const char *src, unsigned char &dest);
template bool ConvertStrToData(const char *src, long long &dest);
template bool ConvertStrToData(const char *src, unsigned long long &dest);
template bool ConvertStrToData(const char *src, float &dest);
template bool ConvertStrToData(const char *src, double &dest);

template<typename T>
bool ConvertDataToStr(T src, string &dest)
{
	return karma::generate(std::back_inserter(dest),
		typename std::conditional<
			std::is_floating_point<T>::value,
				karma::real_generator<T>,
				typename std::conditional<
					std::is_signed<T>::value,
						karma::int_generator<T>,
						karma::uint_generator<T>
				>::type
		>::type(),
		src);
}
template bool ConvertDataToStr(int src, string &dest);
template bool ConvertDataToStr(unsigned int src, string &dest);
template bool ConvertDataToStr(short src, string &dest);
template bool ConvertDataToStr(unsigned short src, string &dest);
template bool ConvertDataToStr(char src, string &dest);
template bool ConvertDataToStr(unsigned char src, string &dest);
template bool ConvertDataToStr(long long src, string &dest);
template bool ConvertDataToStr(unsigned long long src, string &dest);
template bool ConvertDataToStr(float src, string &dest);
template bool ConvertDataToStr(double src, string &dest);

void ExecuteThreaded(std::function<void()> &&func)
{
	g_ThreadPool.submit(std::move(func));
}


text:
//Author: Dr. Shantanu Shahane
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <float.h>
#include <string.h>
#include <iostream>
#include <vector>
#include <algorithm>
#include "class.hpp"
#include <unistd.h>
#include <limits.h>
#include <Eigen/Dense>
#include <unsupported/Eigen/SparseExtra>
#include <Eigen/SparseLU>
#include <Eigen/OrderingMethods>
#include <Eigen/Eigenvalues>
#include <Eigen/Core>
#include <Eigen/SparseCore>
#include <Spectra/GenEigsSolver.h>
#include <Spectra/MatOp/SparseGenMatProd.h>
#include <Spectra/GenEigsRealShiftSolver.h>
#include <Spectra/MatOp/SparseGenRealShiftSolve.h>
using namespace std;

PARAMETERS::PARAMETERS(string parameter_file, string gmsh_file)
{
    cout << "\n";
    check_mpi();
    meshfile = gmsh_file;
    cout << "PARAMETERS::PARAMETERS gmsh_file: " << meshfile << endl;
    does_file_exist(meshfile.c_str(), "Called from PARAMETERS::read_calc_parameters");
    read_calc_parameters(parameter_file);
    verify_parameters();
    get_problem_dimension_msh();
    calc_cloud_num_points();
    calc_polynomial_term_exponents();
    cout << "\n";
}

void PARAMETERS::read_calc_parameters(string parameter_file)
{
    does_file_exist(parameter_file.c_str(), "Called from PARAMETERS::read_calc_parameters");
    char ctemp[5000], ctemp_2[100];
    string stemp;
    int itemp;
    double dtemp;
    FILE *file;
    file = fopen(parameter_file.c_str(), "r");
    cout << endl;
    fscanf(file, "%[^,],%i\n", ctemp, &poly_deg);
    printf("PARAMETERS::read_calc_parameters Read %s = %i\n", ctemp, poly_deg);
    fscanf(file, "%[^,],%i\n", ctemp, &phs_deg);
    printf("PARAMETERS::read_calc_parameters Read %s = %i\n", ctemp, phs_deg);
    fscanf(file, "%[^,],%lf\n", ctemp, &cloud_size_multiplier);
    printf("PARAMETERS::read_calc_parameters Read %s = %g\n", ctemp, cloud_size_multiplier);

    fscanf(file, "%[^,],%i\n", ctemp, &nt);
    printf("PARAMETERS::read_calc_parameters Read %s = %i\n", ctemp, nt);
    fscanf(file, "%[^,],%lf\n", ctemp, &Courant);
    printf("PARAMETERS::read_calc_parameters Read %s = %g\n", ctemp, Courant);
    fscanf(file, "%[^,],%s\n", ctemp, ctemp_2);
    solver_type = ctemp_2;
    printf("PARAMETERS::read_calc_parameters Read %s = %s\n", ctemp, solver_type.c_str());
    fscanf(file, "%[^,],%lf\n", ctemp, &steady_tolerance);
    printf("PARAMETERS::read_calc_parameters Read %s = %g\n", ctemp, steady_tolerance);
    fscanf(file, "%[^,],%lf\n", ctemp, &solver_tolerance);
    printf("PARAMETERS::read_calc_parameters Read %s = %g\n", ctemp, solver_tolerance);
    fscanf(file, "%[^,],%i\n", ctemp, &euclid_precond_level_hypre);
    printf("PARAMETERS::read_calc_parameters Read %s = %i\n", ctemp, euclid_precond_level_hypre);
    fscanf(file, "%[^,],%i\n", ctemp, &gmres_kdim);
    printf("PARAMETERS::read_calc_parameters Read %s = %i\n", ctemp, gmres_kdim);
    fscanf(file, "%[^,],%lf\n", ctemp, &precond_droptol);
    printf("PARAMETERS::read_calc_parameters Read %s = %g\n", ctemp, precond_droptol);
    fscanf(file, "%[^,],%i\n", ctemp, &n_iter);
    printf("PARAMETERS::read_calc_parameters Read %s = %i\n", ctemp, n_iter);

    cout << endl;
}

void PARAMETERS::verify_parameters()
{
    if (poly_deg < 2 || poly_deg > 15)
    {
        printf("\n\nERROR from PARAMETERS::verify_parameters poly_deg should be in range [2, 15]; current value: %i\n\n", poly_deg);
        throw bad_exception();
    }
    if (phs_deg != 3 && phs_deg != 5 && phs_deg != 7 && phs_deg != 9 && phs_deg != 11)
    {
        printf("\n\nERROR from PARAMETERS::verify_parameters phs_deg should be 3, 5, 7, 9, or 11; current value: %i\n\n", phs_deg);
        throw bad_exception();
    }
    if (strcmp(solver_type.c_str(), "hypre_ilu_gmres") != 0 && strcmp(solver_type.c_str(), "eigen_direct") != 0 && strcmp(solver_type.c_str(), "eigen_ilu_bicgstab") != 0)
    {
        cout << "\n\nERROR from PARAMETERS::verify_parameters solver_type should be either hypre_ilu_gmres, eigen_ilu_bicgstab or eigen_direct; current value: " << solver_type << "\n\n";
        throw bad_exception();
    }

    if (meshfile.find_last_of(".") + 1 == meshfile.size())
    {
        cout << "\n\nERROR from PARAMETERS::verify_parameters extension of meshfile should be msh; \nUnable to find extension in the meshfile: " << meshfile << "\n\n";
        throw bad_exception();
    }
    string extension = meshfile.substr(meshfile.find_last_of(".") + 1, meshfile.size());
    if (extension != "msh")
    {
        cout << "\n\nERROR from PARAMETERS::verify_parameters extension of meshfile should be msh; \nCurrent extension: " << extension << " of the meshfile: " << meshfile << "\n\n";
        throw bad_exception();
    }

    output_file_prefix = meshfile;
    if (output_file_prefix.find_last_of("/") < output_file_prefix.size()) //remove location details from mesh file name if required
        output_file_prefix = output_file_prefix.substr(output_file_prefix.find_last_of("/") + 1, output_file_prefix.size());
    if (output_file_prefix.find_last_of(".") < output_file_prefix.size()) //remove extension from mesh file name if required
        output_file_prefix.replace(output_file_prefix.begin() + output_file_prefix.find_last_of("."), output_file_prefix.end(), "");
    output_file_prefix += "_polydeg_" + to_string(poly_deg);
    cout << "PARAMETERS::verify_parameters output_file_prefix: " << output_file_prefix << endl;
}

void PARAMETERS::calc_cloud_num_points()
{
    if (dimension == 2)
        num_poly_terms = ((int)(0.5 * (poly_deg + 1) * (poly_deg + 2)));
    else
        num_poly_terms = ((int)((poly_deg + 1) * (poly_deg + 2) * (poly_deg + 3) / 6)); //sum(0.5*(1:poly_deg+1).*(2:poly_deg+2)) = (poly_deg + 1) * (poly_deg + 2) * (poly_deg + 3) / 6
    cloud_size = (int)(ceil(cloud_size_multiplier * num_poly_terms));
    cout << "PARAMETERS::calc_cloud_num_points num_poly_terms: " << num_poly_terms << ", cloud_size: " << cloud_size << endl;
}

void PARAMETERS::get_problem_dimension_msh()
{ //identify whether its a 2D or 3D gmsh grid

    clock_t start = clock();
    int dim = 2;
    FILE *file;
    int itemp, ncv, cv_type, eof_flag = 0;
    double dtemp;
    char temp[50];
    file = fopen(meshfile.c_str(), "r");
    while (true)
    {
        eof_flag = fscanf(file, "%s ", temp);
        if (strcmp(temp, "$MeshFormat") == 0)
            break;
        if (eof_flag < 0)
        {
            printf("\n\nERROR from PARAMETERS::get_problem_dimension_msh $MeshFormat not found in %s file\n\n", meshfile.c_str());
            throw bad_exception();
        }
    }
    // fscanf(file, "%s", temp);
    fgets(temp, 50, file);
    if (strcmp(temp, "2.2 0 8\n") != 0)
    {
        printf("\n\nERROR from PARAMETERS::get_problem_dimension_msh MeshFormat should be 2.2 0 8; but %s file has %s format instead\n\n", meshfile.c_str(), temp);
        throw bad_exception();
    }
    eof_flag = 0;
    while (true)
    {
        eof_flag = fscanf(file, "%s ", temp);
        if (strcmp(temp, "$Nodes") == 0)
            break;
        if (eof_flag < 0)
        {
            printf("\n\nERROR from PARAMETERS::get_problem_dimension_msh $Nodes not found in %s file\n\n", meshfile.c_str());
            throw bad_exception();
        }
    }
    eof_flag = 0;
    while (true)
    {
        eof_flag = fscanf(file, "%s ", temp);
        if (strcmp(temp, "$Elements") == 0)
            break;
        if (eof_flag < 0)
        {
            printf("\n\nERROR from PARAMETERS::get_problem_dimension_msh $Elements not found in %s file\n\n", meshfile.c_str());
            throw bad_exception();
        }
    }
    fscanf(file, "%i ", &ncv);
    // cout << "get_problem_dimension_msh ncv = " << ncv << endl;
    //reference: http://www.manpagez.com/info/gmsh/gmsh-2.2.6/gmsh_63.php
    for (int icv = 0; icv < ncv; icv++)
    {
        fscanf(file, "%i ", &itemp);   //cv number
        fscanf(file, "%i ", &cv_type); //cv type
        if (cv_type == 4 || cv_type == 5 || cv_type == 6 || cv_type == 7)
        { //4-node tetrahedron || 8-node hexahedron || 6-node prism || 5-node pyramid
            dim = 3;
            break; //3D element found
        }
        else if (cv_type == 11 || cv_type == 12 || cv_type == 13 || cv_type == 14)
        { //10-node tetrahedron || 27-node hexahedron || 18-node prism || 14-node pyramid
            dim = 3;
            break; //3D element found
        }
        else if (cv_type == 17 || cv_type == 18 || cv_type == 19)
        { //20-node hexahedron || 15-node prism || 13-node pyramid
            dim = 3;
            break; //3D element found
        }
        else if (cv_type == 29 || cv_type == 30 || cv_type == 31)
        { //20-node tetrahedron || 35-node tetrahedron || 56-node tetrahedron
            dim = 3;
            break; //3D element found
        }
        fscanf(file, "%*[^\n]\n"); //skip reading remaining row
    }
    dimension = dim;
    cout << "PARAMETERS::get_problem_dimension_msh problem dimension: " << dimension << endl;
}

void PARAMETERS::calc_polynomial_term_exponents()
{
    polynomial_term_exponents.resize(num_poly_terms, dimension);
    Eigen::MatrixXi deg_temp;
    if (dimension == 2)
    {
        polynomial_term_exponents.row(0) = Eigen::MatrixXi::Zero(1, dimension);                                        //constant term
        polynomial_term_exponents.block(1, 0, dimension, dimension) = Eigen::MatrixXi::Identity(dimension, dimension); //two linear terms
        int previous_index = 1;
        for (int deg = 2; deg <= poly_deg; deg++)
        { //quadratic onwards terms
            deg_temp.resize(deg + 1, dimension);
            for (int i = previous_index; i < previous_index + deg; i++)
            {
                deg_temp(i - previous_index, 0) = 1 + polynomial_term_exponents(i, 0); //exponents of 'x'
                deg_temp(i - previous_index, 1) = polynomial_term_exponents(i, 1);     //exponents of 'y'
            }
            deg_temp(deg, 0) = 0;
            deg_temp(deg, 1) = deg;
            previous_index = previous_index + deg;
            polynomial_term_exponents.block(previous_index, 0, deg + 1, dimension) = deg_temp;
        }
    }
    else
    {
        polynomial_term_exponents.row(0) = Eigen::MatrixXi::Zero(1, dimension);                                        //constant term
        polynomial_term_exponents.block(1, 0, dimension, dimension) = Eigen::MatrixXi::Identity(dimension, dimension); //three linear terms
        int previous_index = 1;
        for (int deg = 2; deg <= poly_deg; deg++)
        { //quadratic onwards terms
            deg_temp.resize(0.5 * (deg + 1) * (deg + 2), dimension);
            for (int i = previous_index; i < previous_index + (deg * (deg + 1) / 2); i++)
            {
                deg_temp(i - previous_index, 0) = 1 + polynomial_term_exponents(i, 0); //exponents of 'x'
                deg_temp(i - previous_index, 1) = polynomial_term_exponents(i, 1);     //exponents of 'y'
                deg_temp(i - previous_index, 2) = polynomial_term_exponents(i, 2);     //exponents of 'z'
            }
            for (int i = previous_index + (0.5 * deg * (deg - 1)); i < previous_index + (0.5 * deg * (deg + 1)); i++)
            {
                deg_temp(i - previous_index + deg, 0) = polynomial_term_exponents(i, 0);     //exponents of 'x'
                deg_temp(i - previous_index + deg, 1) = 1 + polynomial_term_exponents(i, 1); //exponents of 'y'
                deg_temp(i - previous_index + deg, 2) = polynomial_term_exponents(i, 2);     //exponents of 'z'
            }
            deg_temp((0.5 * (deg + 1) * (deg + 2)) - 1, 0) = 0;
            deg_temp((0.5 * (deg + 1) * (deg + 2)) - 1, 1) = 0;
            deg_temp((0.5 * (deg + 1) * (deg + 2)) - 1, 2) = deg;
            previous_index = previous_index + (0.5 * deg * (deg + 1));
            polynomial_term_exponents.block(previous_index, 0, 0.5 * (deg + 1) * (deg + 2), dimension) = deg_temp;
        }
    }
}

void PARAMETERS::calc_dt(Eigen::SparseMatrix<double, Eigen::RowMajor> &grad_x, Eigen::SparseMatrix<double, Eigen::RowMajor> &grad_y, Eigen::SparseMatrix<double, Eigen::RowMajor> &grad_z, Eigen::SparseMatrix<double, Eigen::RowMajor> &laplacian, double u0, double v0, double w0, double alpha)
{
    Eigen::VectorXcd eigval_grad_x, eigval_grad_y, eigval_grad_z, eigval_laplacian;
    eigval_grad_x = calc_largest_magnitude_eigenvalue(grad_x);
    eigval_grad_y = calc_largest_magnitude_eigenvalue(grad_y);
    if (dimension == 3)
        eigval_grad_z = calc_largest_magnitude_eigenvalue(grad_z);
    eigval_laplacian = calc_largest_magnitude_eigenvalue(laplacian);
    grad_x_eigval_real = eigval_grad_x[0].real(), grad_x_eigval_imag = eigval_grad_x[0].imag();
    grad_y_eigval_real = eigval_grad_y[0].real(), grad_y_eigval_imag = eigval_grad_y[0].imag();
    if (dimension == 3)
        grad_z_eigval_real = eigval_grad_z[0].real(), grad_z_eigval_imag = eigval_grad_z[0].imag();
    laplace_eigval_real = eigval_laplacian[0].real(), laplace_eigval_imag = eigval_laplacian[0].imag();
    printf("\nPARAMETERS::calc_dt Eigenvalue with largest magnitude: grad_x: (%g, %g)\n", eigval_grad_x[0].real(), eigval_grad_x[0].imag());
    printf("PARAMETERS::calc_dt Eigenvalue with largest magnitude: grad_y: (%g, %g)\n", eigval_grad_y[0].real(), eigval_grad_y[0].imag());
    if (dimension == 3)
        printf("PARAMETERS::calc_dt Eigenvalue with largest magnitude: grad_z: (%g, %g)\n", eigval_grad_z[0].real(), eigval_grad_z[0].imag());
    printf("PARAMETERS::calc_dt Eigenvalue with largest magnitude: laplacian: (%g, %g)\n", eigval_laplacian[0].real(), eigval_laplacian[0].imag());

    double evalues = u0 * sqrt((eigval_grad_x[0].real() * eigval_grad_x[0].real()) + (eigval_grad_x[0].imag() * eigval_grad_x[0].imag()));
    evalues += v0 * sqrt((eigval_grad_y[0].real() * eigval_grad_y[0].real()) + (eigval_grad_y[0].imag() * eigval_grad_y[0].imag()));
    if (dimension == 3)
        evalues += w0 * sqrt((eigval_grad_z[0].real() * eigval_grad_z[0].real()) + (eigval_grad_z[0].imag() * eigval_grad_z[0].imag()));
    evalues += alpha * sqrt((eigval_laplacian[0].real() * eigval_laplacian[0].real()) + (eigval_laplacian[0].imag() * eigval_laplacian[0].imag()));
    dt = 2.0 / evalues; //forward Euler
    dt = dt * Courant;
    printf("PARAMETERS::calc_dt Courant: %g, dt: %g seconds\n\n", Courant, dt);
}

max_stars_repo_name: shahaneshantanu/memphys

max_issues_repo_name: shahaneshantanu/memphys

max_forks_repo_name: shahaneshantanu/memphys