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

#include <sys/types.h>
#include <net/ethernet.h>
#include <boost/uuid/uuid_io.hpp>

#include "base/address_util.h"
#include "base/logging.h"
#include "db/db.h"
#include "db/db_entry.h"
#include "db/db_table.h"
#include "ifmap/ifmap_node.h"

#include <cfg/cfg_init.h>
#include <cmn/agent.h>
#include <init/agent_param.h>
#include <oper/operdb_init.h>
#include <oper/route_common.h>
#include <oper/vm.h>
#include <oper/vn.h>
#include <oper/vrf.h>
#include <oper/nexthop.h>
#include <oper/mpls.h>
#include <oper/mirror_table.h>
#include <oper/interface_common.h>
#include <oper/vrf_assign.h>
#include <oper/vxlan.h>

#include <vnc_cfg_types.h>
#include <oper/agent_sandesh.h>
#include <oper/sg.h>
#include "sandesh/sandesh_trace.h"
#include "sandesh/common/vns_types.h"
#include "sandesh/common/vns_constants.h"
#include <resource_manager/resource_manager.h>
#include <resource_manager/resource_table.h>
#include <resource_manager/mpls_index.h>

using namespace std;
using namespace boost::uuids;

/////////////////////////////////////////////////////////////////////////////
// Inet Interface Key and Data methods
/////////////////////////////////////////////////////////////////////////////
InetInterfaceKey::InetInterfaceKey(const std::string &name) :
    InterfaceKey(AgentKey::ADD_DEL_CHANGE, Interface::INET, nil_uuid(), name,
                 false) {
}

Interface *InetInterfaceKey::AllocEntry(const InterfaceTable *table) const {
    return new InetInterface(name_);
}

Interface *InetInterfaceKey::AllocEntry(const InterfaceTable *table,
                                        const InterfaceData *data)const {
    const InetInterfaceData *vhost_data  =
        static_cast<const InetInterfaceData *>(data);

    VrfKey key(data->vrf_name_);
    VrfEntry *vrf = static_cast<VrfEntry *>
        (table->agent()->vrf_table()->FindActiveEntry(&key));
    assert(vrf);

    Interface *xconnect = NULL;
    if (vhost_data->sub_type_ == InetInterface::VHOST) {
        PhysicalInterfaceKey key(vhost_data->xconnect_);
        xconnect = static_cast<Interface *>
            (table->agent()->interface_table()->FindActiveEntry(&key));
        assert(xconnect != NULL);
    }

    InetInterface *intf = new InetInterface(name_, vhost_data->sub_type_, vrf,
                                            vhost_data->ip_addr_,
                                            vhost_data->plen_, vhost_data->gw_,
                                            xconnect, vhost_data->vn_name_);
    return intf;
}

InterfaceKey *InetInterfaceKey::Clone() const {
    return new InetInterfaceKey(name_);
}

InetInterfaceData::InetInterfaceData(InetInterface::SubType sub_type,
                                     const std::string &vrf_name,
                                     const Ip4Address &addr, int plen,
                                     const Ip4Address &gw,
                                     const std::string &xconnect,
                                     const std::string vn_name,
                                     Interface::Transport transport) :
    InterfaceData(NULL, NULL, transport), sub_type_(sub_type), ip_addr_(addr),
    plen_(plen), gw_(gw), xconnect_(xconnect), vn_name_(vn_name) {
    InetInit(vrf_name);
}

/////////////////////////////////////////////////////////////////////////////
// SIMPLE GATEWAY Utility functions
/////////////////////////////////////////////////////////////////////////////
void InetInterface::ActivateSimpleGateway() {
    InterfaceTable *table = static_cast<InterfaceTable *>(get_table());
    Agent *agent = table->agent();

    // Create InterfaceNH before MPLS is created
    InterfaceNH::CreateInetInterfaceNextHop(name(), vrf()->GetName(),
                                            agent->pkt_interface_mac());

    //There is no policy enabled nexthop created for VGW interface,
    //hence use interface nexthop without policy as flow key index
    InterfaceNHKey key(new InetInterfaceKey(name()),
                       false, InterfaceNHFlags::INET4,
                       agent->pkt_interface_mac());
    flow_key_nh_ = static_cast<const NextHop *>(
            agent->nexthop_table()->FindActiveEntry(&key));
    assert(flow_key_nh_);
    label_ = flow_key_nh_->mpls_label()->label();
}

void InetInterface::DeActivateSimpleGateway() {
    InetUnicastAgentRouteTable *uc_rt_table =
        (VrfTable::GetInstance()->GetInet4UnicastRouteTable(vrf()->GetName()));

    InterfaceTable *table = static_cast<InterfaceTable *>(get_table());
    Agent *agent = table->agent();

    // Delete routes
    Ip4Address addr = Address::GetIp4SubnetAddress(ip_addr_, plen_);
    uc_rt_table->DeleteReq(agent->local_vm_peer(), agent->fabric_vrf_name(),
                           addr, plen_, NULL);

    uc_rt_table->DeleteReq(agent->local_vm_peer(),
                           vrf()->GetName(), Ip4Address(0), 0, NULL);

    // Delete NH
    InterfaceNH::DeleteInetInterfaceNextHop(name(), agent->pkt_interface_mac());

    // Reset Label
    label_ = MplsTable::kInvalidLabel;
    flow_key_nh_ = NULL;
}

/////////////////////////////////////////////////////////////////////////////
// VHOST / LL Utility functions
/////////////////////////////////////////////////////////////////////////////

// Add default route with given gateway
static void AddDefaultRoute(Agent *agent, InetUnicastAgentRouteTable *table,
                            const VrfEntry *vrf, const Interface *xconnect,
                            const Ip4Address &gw, const string &vn_name) {

    if (xconnect) {
        const PhysicalInterface *physical_intf =
            static_cast<const PhysicalInterface *>(xconnect);
        if (physical_intf->no_arp() ||
            (physical_intf->encap_type() == PhysicalInterface::RAW_IP)) {
            table->AddInterfaceRouteReq(agent, agent->local_peer(),
                                        vrf->GetName(), Ip4Address(0), 0,
                                        xconnect, vn_name);
            return;
        }
    }

    VnListType vn_list;
    vn_list.insert(vn_name);
    table->AddGatewayRoute(agent->local_peer(),
                           vrf->GetName(), Ip4Address(0), 0, gw, vn_list,
                           MplsTable::kInvalidLabel, SecurityGroupList(),
                           TagList(), CommunityList(), true);
    return;
}

static void DeleteDefaultRoute(Agent *agent, InetUnicastAgentRouteTable *table,
                               const VrfEntry *vrf, const Ip4Address &addr) {
    table->Delete(agent->local_peer(), vrf->GetName(), Ip4Address(0), 0);
}

// Following routes are added due o an inet interface
// - Receive route for IP address assigned
// - Receive route for the sub-net broadcast address
// - Resolve route for the subnet address
static void AddHostRoutes(Agent *agent, InetUnicastAgentRouteTable *table,
                          const VrfEntry *vrf, const string &interface,
                          const Interface *xconnect, const Ip4Address &addr,
                          int plen, const string &vn_name) {

    InetInterfaceKey intf_key(interface);
    table->AddVHostRecvRoute(agent->local_peer(), vrf->GetName(), intf_key,
                             addr, 32, vn_name, false, true);

    const PhysicalInterface *physical_intf =
        static_cast<const PhysicalInterface *>(xconnect);
    if (physical_intf) {
       if (physical_intf->no_arp())
           return;
       if (physical_intf->encap_type() != PhysicalInterface::ETHERNET)
           return;
    }

    table->AddVHostSubnetRecvRoute(agent->local_peer(), vrf->GetName(),
                                   intf_key,
                                   GetIp4SubnetBroadcastAddress(addr, plen),
                                   32, vn_name, false);
    if (agent->params()->subnet_hosts_resolvable() == true) {
        InetInterfaceKey intf_key(interface);
        table->AddResolveRoute(agent->local_peer(), vrf->GetName(),
                Address::GetIp4SubnetAddress(addr, plen), plen,
                intf_key, MplsTable::kInvalidLabel, false, vn_name,
                SecurityGroupList(), TagList());
    }

    if ((agent->tsn_enabled() == true) &&
        (agent->forwarding_enabled() == false)) {
        InetUnicastAgentRouteTable *inet_table =
            static_cast<InetUnicastAgentRouteTable *>
            (agent->fabric_policy_vrf()->GetInet4UnicastRouteTable());
        VmInterfaceKey vmi_key(AgentKey::ADD_DEL_CHANGE, nil_uuid(),
                               agent->params()->vhost_name());
        inet_table->AddVHostRecvRouteReq(agent->local_peer(),
                                         agent->fabric_policy_vrf_name(),
                                         vmi_key,
                                         agent->params()->vhost_addr(),
                                         32,
                                         agent->fabric_vn_name(),
                                         false, true);
    }
}

static void DeleteHostRoutes(Agent *agent, InetUnicastAgentRouteTable *table,
                             const VrfEntry *vrf, const Interface *xconnect,
                             const Ip4Address &addr, int plen) {
    table->Delete(agent->local_peer(), vrf->GetName(), addr, 32);
    const PhysicalInterface *physical_intf =
        static_cast<const PhysicalInterface *>(xconnect);
    if (physical_intf && physical_intf->no_arp()) {
        return;
    }

    table->Delete(agent->local_peer(), vrf->GetName(),
                  Address::GetIp4SubnetAddress(addr, plen), plen);
    table->Delete(agent->local_peer(), vrf->GetName(),
                  GetIp4SubnetBroadcastAddress(addr, plen), 32);
}

void InetInterface::AddHostMulticastRoutes() {
    VrfTable *vrf_table = static_cast<VrfTable *>(vrf()->get_table());
    InetUnicastAgentRouteTable *uc_rt_table =
        (vrf_table->GetInet4UnicastRouteTable(vrf()->GetName()));
    boost::system::error_code ec;
    InetInterfaceKey intf_key(name());
    // Add v4 route for covering multicast
    uc_rt_table->
        AddVHostRecvRoute(uc_rt_table->agent()->local_peer(),
                          vrf()->GetName(),
                          intf_key,
                          Ip4Address::from_string(IPV4_MULTICAST_BASE_ADDRESS,
                                                  ec),
                          MULTICAST_BASE_ADDRESS_PLEN,
                          vn_name_,
                          false, true);
    // Add v6 route for covering multicast
    uc_rt_table->
        AddVHostRecvRoute(uc_rt_table->agent()->local_peer(),
                          vrf()->GetName(),
                          intf_key,
                          Ip6Address::from_string(IPV6_MULTICAST_BASE_ADDRESS,
                                                  ec),
                          MULTICAST_BASE_ADDRESS_PLEN,
                          vn_name_,
                          false, true);
}

void InetInterface::DelHostMulticastRoutes() {
    VrfTable *vrf_table = static_cast<VrfTable *>(vrf()->get_table());
    InetUnicastAgentRouteTable *uc_rt_table =
        (vrf_table->GetInet4UnicastRouteTable(vrf()->GetName()));
    boost::system::error_code ec;
    // Del v4 route for covering multicast
    uc_rt_table->Delete(uc_rt_table->agent()->local_peer(),
                        vrf()->GetName(),
                        Ip4Address::from_string(IPV4_MULTICAST_BASE_ADDRESS,
                                                ec),
                        MULTICAST_BASE_ADDRESS_PLEN);
    // Del v6 route for covering multicast
    uc_rt_table->Delete(uc_rt_table->agent()->local_peer(),
                        vrf()->GetName(),
                        Ip6Address::from_string(IPV6_MULTICAST_BASE_ADDRESS,
                                                ec),
                        MULTICAST_BASE_ADDRESS_PLEN);
}

// Things to do to activate VHOST/LL interface
// 1. Create the receive next-hops for interface (with policy and witout policy)
// 2. Add routes needed to manage the IP address on interface
// 3. Add default route for the gateway
// 4. Add broadcast route in multicast route table
void InetInterface::ActivateHostInterface() {
    InterfaceTable *table = static_cast<InterfaceTable *>(get_table());
    Agent *agent = table->agent();

    // Create receive nexthops
    ReceiveNH::Create(agent->nexthop_table(), this, false);

    VrfTable *vrf_table = static_cast<VrfTable *>(vrf()->get_table());
    InetUnicastAgentRouteTable *uc_rt_table =
        (vrf_table->GetInet4UnicastRouteTable(vrf()->GetName()));
    if (ip_addr_.to_ulong()) {
        AddHostRoutes(agent, uc_rt_table, vrf(), name(), xconnect_.get(),
                      ip_addr_, plen_, vn_name_);
    }
    AddHostMulticastRoutes();

    if (gw_.to_ulong()) {
        AddDefaultRoute(agent, uc_rt_table, vrf(), xconnect_.get(), gw_,
                        vn_name_);
    }

    ReceiveNHKey nh_key(new InetInterfaceKey(name()), false);
    flow_key_nh_ = static_cast<const NextHop *>(
            agent->nexthop_table()->FindActiveEntry(&nh_key));
}

void InetInterface::DeActivateHostInterface() {
    ipv4_active_ = false;

    Agent *agent = static_cast<InterfaceTable *>(get_table())->agent();
    VrfTable *vrf_table = static_cast<VrfTable *>(vrf()->get_table());
    InetUnicastAgentRouteTable *uc_rt_table =
        (vrf_table->GetInet4UnicastRouteTable(vrf()->GetName()));
    if (ip_addr_.to_ulong()) {
        DeleteHostRoutes(agent, uc_rt_table, vrf(), xconnect_.get(), ip_addr_,
                         plen_);
    }
    DelHostMulticastRoutes();

    if (gw_.to_ulong()) {
        DeleteDefaultRoute(agent, uc_rt_table, vrf(), gw_);
    }

    // Delete receive nexthops
    ReceiveNH::Delete(agent->nexthop_table(), this, false);
    flow_key_nh_ = NULL;
}

/////////////////////////////////////////////////////////////////////////////
// Inet Interface methods
/////////////////////////////////////////////////////////////////////////////
InetInterface::InetInterface(const std::string &name) :
    Interface(Interface::INET, nil_uuid(), name, NULL, true, nil_uuid()),
    sub_type_(VHOST), ip_addr_(0), plen_(0), gw_(0), xconnect_(NULL),
    vn_name_("") {
    ipv4_active_ = false;
    l2_active_ = false;
}

InetInterface::InetInterface(const std::string &name, SubType sub_type,
                             VrfEntry *vrf, const Ip4Address &ip_addr, int plen,
                             const Ip4Address &gw, Interface *xconnect,
                             const std::string &vn_name) :
    Interface(Interface::INET, nil_uuid(), name, vrf, true, nil_uuid()),
    sub_type_(sub_type),
    ip_addr_(ip_addr),
    plen_(plen),
    gw_(gw),
    xconnect_(xconnect),
    vn_name_(vn_name) {
    ipv4_active_ = false;
    l2_active_ = false;
}

bool InetInterface::CmpInterface(const DBEntry &rhs) const {
    const InetInterface &intf = static_cast<const InetInterface &>(rhs);
    return name() < intf.name();
}

DBEntryBase::KeyPtr InetInterface::GetDBRequestKey() const {
    InterfaceKey *key = new InetInterfaceKey(name());
    return DBEntryBase::KeyPtr(key);
}

void InetInterface::Activate() {
    ipv4_active_ = true;
    if (sub_type_ == SIMPLE_GATEWAY) {
        InetInterface::ActivateSimpleGateway();
    } else {
        InetInterface::ActivateHostInterface();
    }

    return;
}

void InetInterface::DeActivate() {
    ipv4_active_ = false;

    if (sub_type_ == SIMPLE_GATEWAY) {
        InetInterface::DeActivateSimpleGateway();
    } else {
        InetInterface::DeActivateHostInterface();
    }
}

bool InetInterface::Delete(const DBRequest *req) {
    DeActivate();
    return true;
}

// Interface Activate cannot be done in AllocEntry. It must be done in PostAdd
// Activating an interface results in adding Interface Nexthops. Creating of
// Interface NH from AllocEntry will fail since the interface is not yet added
// in DB Table PartitionInterface. So, Activate an interface in PostAdd
void InetInterface::PostAdd() {
    Activate();
}

bool InetInterface::OnChange(InetInterfaceData *data) {
    bool ret = false;

    // A Delete followed by Add will result in OnChange callback. Interface is
    // Deactivated on delete. Activate it again to create Interface NH etc...
    if (ipv4_active_ != true) {
        Activate();
        ret = true;
    }

    if (sub_type_ == SIMPLE_GATEWAY) {
        return ret;
    }

    InterfaceTable *table = static_cast<InterfaceTable *>(get_table());
    Agent *agent = table->agent();

    VrfTable *vrf_table = static_cast<VrfTable *>(vrf()->get_table());
    InetUnicastAgentRouteTable *uc_rt_table =
        (vrf_table->GetInet4UnicastRouteTable(vrf()->GetName()));

    if (ip_addr_ != data->ip_addr_ || plen_ != data->plen_) {
        // Delete routes based on old ip-addr and prefix
        if (ip_addr_.to_ulong()) {
            DeleteHostRoutes(agent, uc_rt_table, vrf(), xconnect_.get(),
                             ip_addr_, plen_);
        }

        ip_addr_ = data->ip_addr_;
        plen_ = data->plen_;
        vn_name_ = data->vn_name_;
        // Add routes for new ip-address and prefix
        if (data->ip_addr_.to_ulong()) {
            AddHostRoutes(agent, uc_rt_table, vrf(), name(), xconnect_.get(),
                          ip_addr_, plen_, vn_name_);
        }
        ret = true;
    } else if (vn_name_ != data->vn_name_) {
        // Change in vn_name, update route with new VN Name
        vn_name_ = data->vn_name_;
        AddHostRoutes(agent, uc_rt_table, vrf(), name(), xconnect_.get(),
                      ip_addr_, plen_, vn_name_);
        ret = true;
    }

    if (gw_ != data->gw_) {
        // Delete routes based on old gateway
        if (gw_.to_ulong()) {
            DeleteDefaultRoute(agent, uc_rt_table, vrf(), gw_);
        }

        gw_ = data->gw_;
        // Add route for new gateway
        if (gw_.to_ulong()) {
            AddDefaultRoute(agent, uc_rt_table, vrf(), xconnect_.get(), gw_,
                            vn_name_);
        }

        ret = true;
    }

    if (transport_ != data->transport_) {
        transport_ = data->transport_;
    }
    return ret;
}

// Helper to Inet Interface
void InetInterface::CreateReq(InterfaceTable *table, const std::string &ifname,
                              SubType sub_type, const std::string &vrf_name,
                              const Ip4Address &addr, int plen,
                              const Ip4Address &gw,
                              const std::string &xconnect,
                              const std::string &vn_name,
                              Interface::Transport transport) {
    DBRequest req(DBRequest::DB_ENTRY_ADD_CHANGE);
    req.key.reset(new InetInterfaceKey(ifname));
    req.data.reset(new InetInterfaceData(sub_type, vrf_name, Ip4Address(addr),
                                         plen, Ip4Address(gw), xconnect,
                                         vn_name, transport));
    table->Enqueue(&req);
}

// Helper to Inet Interface
void InetInterface::Create(InterfaceTable *table, const std::string &ifname,
                           SubType sub_type, const std::string &vrf_name,
                           const Ip4Address &addr, int plen,
                           const Ip4Address &gw,
                           const std::string &xconnect,
                           const std::string &vn_name,
                           Interface::Transport transport) {
    DBRequest req(DBRequest::DB_ENTRY_ADD_CHANGE);
    req.key.reset(new InetInterfaceKey(ifname));
    req.data.reset(new InetInterfaceData(sub_type, vrf_name, Ip4Address(addr),
                                         plen, Ip4Address(gw), xconnect,
                                         vn_name, transport));
    table->Process(req);
}

// Helper to delete Inet Interface
void InetInterface::DeleteReq(InterfaceTable *table, const string &ifname) {
    DBRequest req(DBRequest::DB_ENTRY_DELETE);
    req.key.reset(new InetInterfaceKey(ifname));
    req.data.reset(NULL);
    table->Enqueue(&req);
}

void InetInterface::Delete(InterfaceTable *table, const string &ifname) {
    DBRequest req(DBRequest::DB_ENTRY_DELETE);
    req.key.reset(new InetInterfaceKey(ifname));
    req.data.reset(NULL);
    table->Process(req);
}

max_stars_repo_name: jnpr-pranav/contrail-controller

max_issues_repo_name: jnpr-pranav/contrail-controller

max_forks_repo_name: jnpr-pranav/contrail-controller


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

#include <boost/multiprecision/cpp_int.hpp>
#include <memory>

#if defined(__apple_build_version__) && (__clang_major__ < 9)
//
// Apples clang fails with:
// error: no matching function for call to '__implicit_conversion_to'
// Which is nothing to do with us really...
//
#define DISABLE_TEST
#endif

typedef boost::multiprecision::cpp_int mp_int;

#if !defined(BOOST_NO_CXX11_EXPLICIT_CONVERSION_OPERATORS) && !defined(BOOST_NO_CXX11_DELETED_FUNCTIONS) && !defined(DISABLE_TEST)

class Int1
{
 public:
   Int1(const mp_int& ) {}
   Int1(const Int1& ) {}
};

class Int2
{
 public:
   Int2(const mp_int& ) {}
   Int2(const Int2& ) = delete;
};

int main()
{
   using namespace boost::multiprecision;

   mp_int i(10);
   Int1   a(i + 10);
   Int2   b(i + 20);

#ifndef BOOST_NO_CXX11_SMART_PTR

   std::shared_ptr<Int1> p1 = std::make_shared<Int1>(i + 10);
   std::shared_ptr<Int2> p2 = std::make_shared<Int2>(i + 10);

#endif

   return 0;
}

#else

int main()
{
   return 0;
}

#endif

ext: cpp

max_issues_repo_name: qingkouwei/mediaones

max_forks_repo_name: qingkouwei/mediaones


text:
/*=============================================================================
    Copyright (c) 2015 Paul Fultz II
    constexpr_deduce.h
    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_HOF_GUARD_FUNCTION_CONSTEXPR_DEDUCE_H
#define BOOST_HOF_GUARD_FUNCTION_CONSTEXPR_DEDUCE_H

#include <boost/hof/config.hpp>

#define BOOST_HOF_CONST_FOLD(x) (__builtin_constant_p(x) ? (x) : (x))

#ifndef BOOST_HOF_HAS_CONST_FOLD
#if defined(__GNUC__) && !defined (__clang__) && __GNUC__ == 4 && __GNUC_MINOR__ < 7
#define BOOST_HOF_HAS_CONST_FOLD 0
#elif defined(__clang__) || defined (__GNUC__)
#define BOOST_HOF_HAS_CONST_FOLD 1
#else
#define BOOST_HOF_HAS_CONST_FOLD 0
#endif
#endif

namespace boost { namespace hof {

namespace detail {

struct constexpr_deduce
{
    constexpr constexpr_deduce()
    {}
    template<class F>
    constexpr operator F() const
    {
        return F();
    }
};

template<class T>
struct constexpr_deduce_unique
{
    constexpr constexpr_deduce_unique()
    {}
#if BOOST_HOF_HAS_CONST_FOLD
    template<class F>
    constexpr operator const F&() const
    {
        static_assert(BOOST_HOF_IS_EMPTY(F), "Function or lambda expression must be empty");
        return BOOST_HOF_CONST_FOLD(reinterpret_cast<const F&>(static_const_var<T>()));
    }
#else
    template<class F>
    constexpr operator F() const
    {
        // static_assert(std::is_default_constructible<F>::value, "Function not default constructible");
        return F();
    }
#endif
};

}}} // namespace boost::hof

#define BOOST_HOF_DETAIL_CONSTEXPR_DEDUCE true ? boost::hof::detail::constexpr_deduce() :
#define BOOST_HOF_DETAIL_CONSTEXPR_DEDUCE_UNIQUE(T) true ? boost::hof::detail::constexpr_deduce_unique<T>() :

#ifdef _MSC_VER
#define BOOST_HOF_DETAIL_MSVC_CONSTEXPR_DEDUCE BOOST_HOF_DETAIL_CONSTEXPR_DEDUCE
#else
#define BOOST_HOF_DETAIL_MSVC_CONSTEXPR_DEDUCE
#endif

#endif

ext: hpp


text:
// Copyright David Abrahams 2003.
// 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 COPY_CTOR_MUTATES_RHS_DWA2003219_HPP
# define COPY_CTOR_MUTATES_RHS_DWA2003219_HPP

#include <boost/python/detail/is_auto_ptr.hpp>
#include <boost/mpl/bool.hpp>

namespace boost { namespace python { namespace detail {

template <class T>
struct copy_ctor_mutates_rhs
    : is_auto_ptr<T>
{
};

}}} // namespace boost::python::detail

#endif // COPY_CTOR_MUTATES_RHS_DWA2003219_HPP


text:
/*
 * Copyright (c) Facebook, Inc. and its affiliates.
 *
 * This source code is licensed under the MIT license found in the
 * LICENSE file in the root directory of this source tree.
 */

#include "LocalDcePass.h"

#include <array>
#include <iostream>
#include <unordered_set>
#include <vector>

#include <boost/dynamic_bitset.hpp>

#include "ControlFlow.h"
#include "DexClass.h"
#include "DexUtil.h"
#include "GraphUtil.h"
#include "IRCode.h"
#include "IRInstruction.h"
#include "MethodOverrideGraph.h"
#include "Purity.h"
#include "Resolver.h"
#include "Transform.h"
#include "TypeSystem.h"
#include "Walkers.h"

namespace {

constexpr const char* METRIC_NPE_INSTRUCTIONS = "num_npe_instructions";
constexpr const char* METRIC_DEAD_INSTRUCTIONS = "num_dead_instructions";
constexpr const char* METRIC_UNREACHABLE_INSTRUCTIONS =
    "num_unreachable_instructions";
constexpr const char* METRIC_COMPUTED_NO_SIDE_EFFECTS_METHODS =
    "num_computed_no_side_effects_methods";
constexpr const char* METRIC_COMPUTED_NO_SIDE_EFFECTS_METHODS_ITERATIONS =
    "num_computed_no_side_effects_methods_iterations";

} // namespace

void LocalDcePass::run_pass(DexStoresVector& stores,
                            ConfigFiles& conf,
                            PassManager& mgr) {
  auto scope = build_class_scope(stores);
  auto pure_methods = get_pure_methods();
  auto configured_pure_methods = conf.get_pure_methods();
  pure_methods.insert(configured_pure_methods.begin(),
                      configured_pure_methods.end());
  auto override_graph = method_override_graph::build_graph(scope);
  std::unordered_set<const DexMethod*> computed_no_side_effects_methods;
  auto computed_no_side_effects_methods_iterations =
      compute_no_side_effects_methods(scope, override_graph.get(), pure_methods,
                                      &computed_no_side_effects_methods);
  for (auto m : computed_no_side_effects_methods) {
    pure_methods.insert(const_cast<DexMethod*>(m));
  }

  bool may_allocate_registers = !mgr.regalloc_has_run();
  if (!may_allocate_registers) {
    // compute_no_side_effects_methods might have found methods that have no
    // implementors. Let's not silently remove invocations to those in LocalDce,
    // as invoking them *will* unconditionally cause an exception.
    for (auto it = pure_methods.begin(); it != pure_methods.end();) {
      auto m = *it;
      if (m->is_def() && !has_implementor(override_graph.get(), m->as_def())) {
        it = pure_methods.erase(it);
      } else {
        it++;
      }
    }
  }

  auto stats =
      walk::parallel::methods<LocalDce::Stats>(scope, [&](DexMethod* m) {
        auto* code = m->get_code();
        if (code == nullptr || m->rstate.no_optimizations()) {
          return LocalDce::Stats();
        }

        LocalDce ldce(pure_methods, override_graph.get(),
                      may_allocate_registers);
        ldce.dce(code);
        return ldce.get_stats();
      });
  mgr.incr_metric(METRIC_NPE_INSTRUCTIONS, stats.npe_instruction_count);
  mgr.incr_metric(METRIC_DEAD_INSTRUCTIONS, stats.dead_instruction_count);
  mgr.incr_metric(METRIC_UNREACHABLE_INSTRUCTIONS,
                  stats.unreachable_instruction_count);
  mgr.incr_metric(METRIC_COMPUTED_NO_SIDE_EFFECTS_METHODS,
                  computed_no_side_effects_methods.size());
  mgr.incr_metric(METRIC_COMPUTED_NO_SIDE_EFFECTS_METHODS_ITERATIONS,
                  computed_no_side_effects_methods_iterations);

  TRACE(DCE, 1, "instructions removed -- npe: %d, dead: %d, unreachable: %d",
        stats.npe_instruction_count, stats.dead_instruction_count,
        stats.unreachable_instruction_count);
}

static LocalDcePass s_pass;

max_stars_repo_name: agampe/redex

max_issues_repo_name: agampe/redex

max_forks_repo_name: agampe/redex


text:
/*
 * TypeDefs.hpp
 *
 *  Created on: March 18, 2014
 *      Author: Péter Fankhauser
 *	 Institute: ETH Zurich, ANYbotics
 */

// Eigen
#pragma once

#include <Eigen/Core>

namespace grid_map {

  using Matrix = Eigen::MatrixXf;
  using DataType = Matrix::Scalar;
  using Position = Eigen::Vector2d;
  using Vector = Eigen::Vector2d;
  using Position3 = Eigen::Vector3d;
  using Vector3 = Eigen::Vector3d;
  using Index = Eigen::Array2i;
  using Size = Eigen::Array2i;
  using Length = Eigen::Array2d;
  using Time = uint64_t;

  /*
   * Interpolations are ordered in the order
   * of increasing accuracy and computational complexity.
   * INTER_NEAREST - fastest, but least accurate,
   * INTER_CUBIC - slowest, but the most accurate.
   * see:
   * https://en.wikipedia.org/wiki/Bicubic_interpolation
   * https://web.archive.org/web/20051024202307/http://www.geovista.psu.edu/sites/geocomp99/Gc99/082/gc_082.htm
   * for more info. Cubic convolution algorithm is also known as piecewise cubic
   * interpolation and in general does not guarantee continuous
   * first derivatives.
   */
  enum class InterpolationMethods{
      INTER_NEAREST, // nearest neighbor interpolation
      INTER_LINEAR,   // bilinear interpolation
      INTER_CUBIC_CONVOLUTION, //piecewise bicubic interpolation using convolution algorithm
      INTER_CUBIC // standard bicubic interpolation
  };

}  // namespace grid_map


text:
/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\
|  Phycas: Python software for phylogenetic analysis                          |
|  Copyright (C) 2006 Mark T. Holder, Paul O. Lewis and David L. Swofford     |
|                                                                             |
|  This program is free software; you can redistribute it and/or modify       |
|  it under the terms of the GNU General Public License as published by       |
|  the Free Software Foundation; either version 2 of the License, or          |
|  (at your option) any later version.                                        |
|                                                                             |
|  This program is distributed in the hope that it will be useful,            |
|  but WITHOUT ANY WARRANTY; without even the implied warranty of             |
|  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the              |
|  GNU General Public License for more details.                               |
|                                                                             |
|  You should have received a copy of the GNU General Public License along    |
|  with this program; if not, write to the Free Software Foundation, Inc.,    |
|  51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.                |
\~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/

#ifndef SLICE_SAMPLER_HPP
#define SLICE_SAMPLER_HPP

#include <cstdlib>
#include <climits>
#include <cfloat>
#include "basic_lot.hpp"
#include "probability_distribution.hpp"
#include "xprobdist.hpp"
#include <boost/shared_ptr.hpp>
#include <boost/weak_ptr.hpp>
#if defined(PYTHON_ONLY)
#	include <boost/python/call_method.hpp>
#endif
typedef std::pair<double, double> ParamAndLnProb;
typedef std::pair<double, double> SliceInterval;

namespace phycas
{

struct SliceStats
	{
	SliceStats() : nsamples(0), value(0.0), width(0.0), diff(0.0), failed(0.0), evals(0.0) {}

	unsigned nsamples;
	double value;
	double width;
	double diff;
	double failed;
	double evals;
	};

// Note: AdHocDensity is defined in probablity_distribution.hpp
typedef boost::shared_ptr<AdHocDensity> FuncToSampleShPtr;

#undef WEAK_FUNCTOSAMPLE
#if defined(WEAK_FUNCTOSAMPLE)
typedef boost::weak_ptr<AdHocDensity> FuncToSampleWkPtr;
#endif

/*----------------------------------------------------------------------------------------------------------------------
|	Implements the univariate slice sampler described in Neal, Radford M. 2003. Slice sampling. Annals of Statistics
|	31:705-741.
*/
class SliceSampler
	{
	public:
								SliceSampler();
#if defined(WEAK_FUNCTOSAMPLE)
								SliceSampler(LotShPtr rnd, FuncToSampleWkPtr f);
#else
								SliceSampler(LotShPtr rnd, FuncToSampleShPtr f);
#endif
		virtual					~SliceSampler();

		double					Sample();
		std::vector<double>					DebugSample();

		double					OverrelaxedSample();
		std::vector<double>					DebugOverrelaxedSample();

#if defined(WEAK_FUNCTOSAMPLE)
		void					AttachFunc(FuncToSampleWkPtr f);
#else
		void					AttachFunc(FuncToSampleShPtr f);
#endif
		void					AttachRandomNumberGenerator(LotShPtr rnd);

		void					SetXValue(double x);
		void					SetMaxUnits(unsigned umax);

		void					SetSliceUnitWidth(double uwidth);
		double					GetSliceUnitWidth() const;

		double					AdaptSimple(double multiplier);
		double					AdaptNeal(double multiplier);
		void					AdaptYConditional(double from_ends, double multiplier);
		SliceInterval			FindSliceInterval(ParamAndLnProb x0, const double ln_y0, double tol, unsigned max_steps = UINT_MAX) const;
		double					CalcW(double y0) const;

		double					GetMode() const;
		double					GetLnDensityAtMode() const;

		double					GetLastSampledXValue();
		double					GetLastSampledYValue();
		double					GetSliceYValue();

		double					GetLnZero() const {return -DBL_MAX;}
        void                    UseDoublingMethod(bool d);

		// For diagnosing inefficiency
		//
		double					GetMinX();
		double					GetMaxX();
		double					GetOrigLeftEdgeOfSlice();
		double					GetOrigRightEdgeOfSlice();
		double					GetLeftEdgeOfSlice();
		double					GetRightEdgeOfSlice();
		unsigned				GetNumFuncEvals();
		unsigned				GetNumFailedSamples();
		unsigned				GetNumUnitsRequired();
		unsigned				GetNumSamples();

		SliceStats				SummarizeDiagnostics();
		void					ResetDiagnostics();

	protected:

		void					Init();
		ParamAndLnProb			GetNextSample(const ParamAndLnProb);
		ParamAndLnProb			GetNextOverrelaxedSample(const ParamAndLnProb);
		SliceInterval			BisectionSqueeze(double left, double lnf_left, double right, double lnf_right, const double ln_y0, double tol, unsigned max_steps) const;

#if defined(WEAK_FUNCTOSAMPLE)
		FuncToSampleWkPtr		func;				/**< is a functor representing the probability distribution to be sampled */
#else
		FuncToSampleShPtr		func;				/**< is a functor representing the probability distribution to be sampled */
#endif
		LotShPtr				r;					/**< is the random number generator */
		ParamAndLnProb			lastSampled;		/**< most recent valid sample and its relative density */
		double					w;					/**< unit size for interval I */
		unsigned				maxUnits;			/**< maximum number of units of size w to use for interval I (set to UINT_MAX for unlimited) */

		// These quantities are used for diagnostic purposes in GetNextSample()
		//
		double					orig_left_edge;		/**< x coordinate of left edge of most recent slice (before cropping) */
		double					orig_right_edge;	/**< x coordinate of right edge of most recent slice (before cropping) */
		double					left_edge;			/**< x coordinate of left edge of most recent slice (after cropping) */
		double					right_edge;			/**< x coordinate of right edge of most recent slice (after cropping) */
		double					ln_y;				/**< log of f(x) representing most recent slice */

		// For diagnosing inefficiency, all reset with call to ResetDiagnostics()
		//
		double					min_x;				/**< minimum x value tried */
		double					max_x;				/**< maximum x value tried */
		double					sumValues;			/**< sum of last `num_samples' sampled values */
		double					sumWidths;			/**< sum of last `num_samples' cropped slice widths (where slice width = right_edge - left_edge) */
		double					sumDiffs;			/**< sum of last `num_samples' differences between successively sampled values */
		unsigned				func_evals;			/**< counts number of function evaluations */
		unsigned				failed_samples;		/**< counts number of samples that failed because they were not in the slice */
		unsigned				realized_m;			/**< counts number of units, each of size w, that were required to bracket the slice */
		unsigned				num_samples;		/**< counts number of times GetNextSample called */
		std::vector<ParamAndLnProb>	most_recent;	/**< vector of all (x, lnx) pairs evaluated in most recent sampling effort */

		// These are needed for y-conditional adaptation
		bool					ycond_on;			/**< if true, w chosen anew for each sample based on y-coordinate of slice */
		double					ycond_a;
		double					ycond_b;
		double					ycond_multiplier;
		ParamAndLnProb			mode;

		// These are for overrelaxed sampling
		unsigned				num_overrelaxed_samples;	/**< counts number of times GetNextOverrelaxedSample called */

        bool                    doubling;           /**< if true, doubling method will be used to increase slice interval */
	};

typedef boost::shared_ptr<SliceSampler> SliceSamplerShPtr;

} // namespace phycas

#endif

max_stars_repo_name: plewis/phycas

max_issues_repo_name: plewis/phycas

max_forks_repo_name: plewis/phycas


text:
#include <cslibs_kdl/external_forces.h>
#include <cslibs_kdl/kdl_conversion.h>
#include <Eigen/Jacobi>
using namespace cslibs_kdl;

ExternalForcesSerialChain::ExternalForcesSerialChain() :
    use_fingers_(false),
    set_model_(false),
    init_sensor_mat_(false),
    n_joints_(0),
    end_eff_frame_("")//,
{

}

bool ExternalForcesSerialChain::initialize()
{
    if(!set_model_){
        return false;
    }
    std::size_t nj_sq = n_joints_ * n_joints_;

    if(!init_sensor_mat_){
        sensor_mat_.setZero(n_joints_, nj_sq);

        std::size_t row = 0;

        for(auto link : link_names_){
            KDL::Vector rot_axis;
            model_->getRotationAxis(link, rot_axis);
//            Eigen::Matrix<double,1,6> si;
//            si << rot_axis(0),rot_axis(1),rot_axis(2),0,0,0;
//            sensor_mat_.block(row,6*row,1,6) = si;
//            sensor_mat_.block(row, 6*row,1,6) << rot_axis(0),rot_axis(1),rot_axis(2),0,0,0;
            sensor_mat_(row, 6*row)     = rot_axis(0);
            sensor_mat_(row, 6*row + 1) = rot_axis(1);
            sensor_mat_(row, 6*row + 2) = rot_axis(2);
            ++row;
        }
    }
    init_sensor_mat_ = true;

    return true;
}

void ExternalForcesSerialChain::setModel(const std::string &robot_model,
                              const std::string &chain_root,
                              const std::string &chain_tip,
                              const std::string &finger_1_tip,
                              const std::string &finger_2_tip,
                              const std::string &finger_3_tip)
{
    robot_model_ = robot_model;
    chain_root_ = chain_root;
    chain_tip_ = chain_tip;
    finger_1_tip_ = finger_1_tip;
    finger_2_tip_ = finger_2_tip;
    finger_3_tip_ = finger_3_tip;

    model_.reset(new cslibs_kdl::KinematicModel(robot_model, chain_root, chain_tip));
    if(finger_1_tip_ != ""){
      model_f1_.reset(new cslibs_kdl::KinematicModel(robot_model, chain_tip, finger_1_tip));
      use_fingers_ = true;
    }
    if(finger_2_tip_ != ""){
      model_f2_.reset(new cslibs_kdl::KinematicModel(robot_model, chain_tip, finger_2_tip));
      use_fingers_ = true;
    }
    if(finger_3_tip_ != ""){
      model_f3_.reset(new cslibs_kdl::KinematicModel(robot_model, chain_tip, finger_3_tip));
      use_fingers_ = true;
    }

    link_names_ = model_->getLinkNames();
    n_joints_ = static_cast<std::size_t>(model_->getNrOfJoints());
    set_model_ = true;
    init_sensor_mat_ = false;
}
void ExternalForcesSerialChain::getGeometricJacobianTransposed(const std::vector<double>& pos,
                                    std::size_t n_joints,
                                    Eigen::MatrixXd& result ) const
{
    if(!init_sensor_mat_){
        return;
    }
    std::size_t max = n_joints;
    if(n_joints > n_joints_){
        max = n_joints_;
    }

    result.setZero(max, 6);
    std::string frame_end = link_names_[max -1];
    for(std::size_t row = 0; row < max; ++row){
        std::string link = link_names_[row];
        KDL::Frame T0i;
        model_->getFKPose(pos, T0i, link);
//        KDL::Vector rot_axis = rot_axis_[row];
        KDL::Frame T0j;
        model_->getFKPose(pos, T0j, frame_end);
        KDL::Frame Tij = T0i.Inverse() * T0j;
        Eigen::Matrix<double, 6, 6> T = cslibs_kdl::convert2EigenWrenchTransform(Tij);
//        Eigen::Matrix<double,1,6> si = sensor_mat_.block(row,6*row,1,6).eval();
        Eigen::Matrix<double,1,6>  column = sensor_mat_.block(row,6*row,1,6).eval() * T;
        result.block(row,0,1,6) =  column;
    }
}

void ExternalForcesSerialChain::getGeometricJacobianTransposed(const std::vector<double> &pos, std::string segment, Eigen::MatrixXd &result) const
{

    std::string frame = segment;
    if(use_fingers_ && segment.find("finger") != std::string::npos ){
        frame = chain_tip_;
    }
    int index = model_->getKDLSegmentIndex(frame);
    if(index < 0){
        throw std::runtime_error("Cannot find link id for " + segment);
    }

    getGeometricJacobianTransposed(pos,index +1, result);
}

void ExternalForcesSerialChain::getGeometricJacobianTransposed(const cslibs_kdl_data::JointStateData &state, std::size_t n_joints, Eigen::MatrixXd &result) const
{
    getGeometricJacobianTransposed(state.position, n_joints, result);
}

void ExternalForcesSerialChain::getWrenchProjetion(const cslibs_kdl_data::JointStateData& state, Eigen::MatrixXd& result ) const
{
    std::size_t nj_sq = n_joints_ * n_joints_;
    result.setZero(nj_sq, nj_sq);

    std::size_t row = 0;

    for(const::std::string& link : link_names_){
        KDL::Frame T0i;
        model_->getFKPose(state.position, T0i, link);
        for(std::size_t col = row; col < n_joints_; ++col){
            KDL::Frame T0j;
            model_->getFKPose(state.position, T0j, link_names_[col]);
            KDL::Frame Tij = T0i.Inverse() * T0j;
            result.block<6,6>(6*row,6*col) = cslibs_kdl::convert2EigenWrenchTransform(Tij);
        }
        ++row;
    }
}

void ExternalForcesSerialChain::getWrenchProjetion(const cslibs_kdl_data::JointStateData& state,
                                      std::size_t up_to_joint,
                                      Eigen::MatrixXd& result ) const
{
    std::size_t nj_sq = n_joints_ * n_joints_;
    result.setZero(nj_sq, nj_sq);

    for(std::size_t row = 0; row <= up_to_joint; ++row){
        std::string link = link_names_[row];
        KDL::Frame T0i;
        model_->getFKPose(state.position, T0i, link);
        for(std::size_t col = row; col <= up_to_joint; ++col){
            KDL::Frame T0j;
            model_->getFKPose(state.position, T0j, link_names_[col]);
            KDL::Frame Tij = T0i.Inverse() * T0j;
            result.block<6,6>(6*row,6*col) = cslibs_kdl::convert2EigenWrenchTransform(Tij);
        }
        ++row;
    }
}

Eigen::MatrixXd ExternalForcesSerialChain::getJacobian(const std::vector<double> &pos, int id)
{

    return model_->getJacobian(pos, id).data;
}

Eigen::MatrixXd ExternalForcesSerialChain::getJacobian(const cslibs_kdl_data::JointStateData &state, int id)
{
    return model_->getJacobian(state.position, id).data;
}

bool ExternalForcesSerialChain::getJacobian(const cslibs_kdl_data::JointStateData &state, const std::string &frame, Eigen::MatrixXd& jacobian)
{
    KDL::Chain chain;
    bool succ = model_->getTree().getChain(chain_root_, frame, chain);
    if(succ){
        KDL::ChainJntToJacSolver jac_solver(chain);
        if(state.position.size() > n_joints_ ){
            throw std::logic_error("Dimension mismatch. More joint values expected");
        }
        KDL::Jacobian jac(chain.getNrOfJoints());

        KDL::JntArray qvals;
        cslibs_kdl::convert(state.position, qvals, state.position.size() - chain.getNrOfJoints());

        jac_solver.JntToJac(qvals, jac);
        jacobian = jac.data;
    }
    return succ;
}

Eigen::VectorXd ExternalForcesSerialChain::getExternalTorques(const std::vector<double>& pos,
                                   std::string frame_id,
                                   const KDL::Wrench &w) const
{
    Eigen::MatrixXd J;
    KDL::Wrench wt = w;

    std::string frame = frame_id;
    if(use_fingers_ && frame_id.find("finger") != std::string::npos ){
        frame = chain_tip_;
        KDL::Frame hTf = getFKPose(pos, chain_tip_);
        wt = hTf * wt;
    }

    Eigen::VectorXd F = cslibs_kdl::convert2Eigen(wt);

    int index = model_->getKDLSegmentIndex(frame);
    if(index < 0){
        throw std::runtime_error("Cannot find link id for " + frame_id);
    }

    getGeometricJacobianTransposed(pos,index +1, J);

    //        geometry_msgs::WrenchStamped msg;
    //        msg.wrench.force.x = n_i.force.x();
    //        msg.wrench.force.y = n_i.force.y();
    //        msg.wrench.force.z = n_i.force.z();
    //        msg.wrench.torque.x = n_i.torque.x();
    //        msg.wrench.torque.y = n_i.torque.y();
    //        msg.wrench.torque.z = n_i.torque.z();
    //        msg.header.frame_id = link_names_[id];
    //        msg.header.stamp = ros::Time::now();
    //        pubTorques.publish(msg);
    //        ros::spinOnce();

//    res = sensor_mat_ * J * F;
    Eigen::VectorXd res = J * F;
    return res;
}

Eigen::VectorXd ExternalForcesSerialChain::getExternalTorquesKDL(const std::vector<double>& pos, std::string frame_id, KDL::Wrench &w_local) const
{
    std::string frame = frame_id;
    KDL::Wrench wchain = w_local;
    // if local frame is a finger frame project wrench to finger frame
    if( use_fingers_ && frame_id.find("finger") != std::string::npos ){
        frame = chain_tip_;
        KDL::Frame hTf = getFKPose(pos, chain_tip_);
        wchain = hTf * wchain;
    }

    int index = model_->getKDLSegmentIndex(frame);
    if(index < 0){
        throw std::runtime_error("Cannot find link id for " + frame_id);
    }

    // rotate wrench in to base frame coordinates
    KDL::Frame baseTframe = getFKPose(pos, frame_id);
    KDL::Wrench w_in(baseTframe.M * wchain.force, baseTframe.M * wchain.torque);

    KDL::Jacobian jac = model_->getJacobian(pos, index +1);
    Eigen::MatrixXd jac_t = jac.data.transpose();
    Eigen::Matrix<double, 6, 1> ft;
    ft << w_in.force(0) ,  w_in.force(1) ,w_in.force(2),
            w_in.torque(0) , w_in.torque(1) , w_in.torque(2);
    Eigen::VectorXd tau = jac_t * ft;
    return tau;
}

Eigen::VectorXd ExternalForcesSerialChain::getExternalTorques(const cslibs_kdl_data::JointStateData &state,
                                   std::string frame_id,
                                   const KDL::Wrench& w) const
{

    return getExternalTorques(state.position, frame_id, w);
}

Eigen::VectorXd ExternalForcesSerialChain::getExternalTorquesKDL(const cslibs_kdl_data::JointStateData &state, std::string frame_id, KDL::Wrench& w_local) const
{
    return getExternalTorquesKDL(state.position, frame_id, w_local);
}

KDL::Frame ExternalForcesSerialChain::getFKPose(const std::vector<double>& position, const std::string& link) const
{
    KDL::Frame res;
    int ec = 0;
    auto getQvec = [](const std::vector<double>& pos, const std::size_t& n_joints, const std::size_t offset){
         std::vector<double> q(2,0);
        if(pos.size() > n_joints + offset){
            q[0] = pos[n_joints + offset];
        }
        return q;
    };
    if(link.find("finger") == std::string::npos){
        ec = model_->getFKPose(position, res, link);
    }
    else if(use_fingers_ && link.find("1") != std::string::npos ){
        KDL::Frame hand;
        KDL::Frame f1;
        ec = model_->getFKPose(position, hand, model_->getTipLink());
        std::vector<double> q = getQvec(position, n_joints_, 0);
        ec *= model_f1_->getFKPose(q, f1, link);
        res = hand * f1;
    }
    else if(use_fingers_ && link.find("2") != std::string::npos ){
        KDL::Frame hand;
        KDL::Frame f2;
        ec = model_->getFKPose(position, hand, model_->getTipLink());
        std::vector<double> q = getQvec(position, n_joints_, 1);
        ec *= model_f2_->getFKPose(q, f2, link);
        res = hand * f2;
    }
    else if(use_fingers_ && link.find("3") != std::string::npos ){
        KDL::Frame hand;
        KDL::Frame f3;
        ec = model_->getFKPose(position, hand, model_->getTipLink());
        std::vector<double> q = getQvec(position, n_joints_, 2);
        ec *= model_f3_->getFKPose(q, f3, link);
        res = hand * f3;
    }
    if(ec < 0){
        std::stringstream ss;
        ss << "can not compute forward kinematics. Frame id: " << link;
        throw std::runtime_error(ss.str());
    }
    return res;
}

KDL::Frame ExternalForcesSerialChain::getFKPose(const std::vector<double> &pos, const std::string& parent, const std::string& link) const
{
    // FK returns pose in base frame coordinates!
    KDL::Frame b_T_p = getFKPose(pos, parent);
    KDL::Frame b_T_l = getFKPose(pos, link);
    KDL::Frame p_T_l = b_T_p.Inverse() * b_T_l;
    return p_T_l;
}

KDL::Frame ExternalForcesSerialChain::getFKPose(const cslibs_kdl_data::JointStateData &state, const std::string& link) const
{
    return getFKPose(state.position, link);
}
std::vector<std::string> ExternalForcesSerialChain::getAllLinkNames() const
{
    std::vector<std::string> res = model_->getLinkNames();
    if(use_fingers_){
        std::vector<std::string> f1 = model_f1_->getLinkNames();
        std::vector<std::string> f2 = model_f2_->getLinkNames();
        std::vector<std::string> f3 = model_f3_->getLinkNames();

        for(auto s : f1){
            res.emplace_back(s);
        }
        for(auto s : f2){
            res.emplace_back(s);
        }
        for(auto s : f3){
            res.emplace_back(s);
        }
    }
    return res;
}

KDL::Wrench ExternalForcesSerialChain::createWrench(const KDL::Vector& position, const KDL::Vector& direction_vector, const KDL::Vector& force)
{
    KDL::Vector x(1,0,0);
    KDL::Vector axis = x * direction_vector;
    KDL::Wrench w(force, KDL::Vector::Zero());
    double alpha = std::acos(dot(x, direction_vector));
    KDL::Frame T(KDL::Rotation::Rot(axis, alpha), position);
    w = T * w;
    return w;
}

KDL::Wrench ExternalForcesSerialChain::createWrench(const KDL::Vector& position, const KDL::Vector& direction_vector, const double theta, const double phi)
{
    KDL::Vector z(0,0,1);
    KDL::Vector axis = z * direction_vector;
    KDL::Vector force(std::sin(theta)*std::cos(phi),
                      std::sin(theta)*std::sin(phi),
                      std::cos(theta));
    KDL::Wrench w(-force, KDL::Vector::Zero());
    double alpha = std::acos(dot(z, direction_vector));
    KDL::Frame T(KDL::Rotation::Rot(axis, alpha), position);
    w = T * w;
    return w;
}

void ExternalForcesSerialChain::forceRegression(const cslibs_kdl_data::JointStateData &state, std::vector<cslibs_kdl_data::ContactPoint> &contacts) const
{
    std::size_t n_c = contacts.size();
    std::size_t offset = state.torque.size() - n_joints_;
    Eigen::MatrixXd reg_mat = Eigen::MatrixXd::Zero(n_joints_, n_c);
    if(state.torque.size() < n_joints_){
        throw std::runtime_error("Dimension missmatch! Expected " + std::to_string(n_joints_) + " torques");
    }
    Eigen::MatrixXd r = state.getEigenVector(cslibs_kdl_data::JointStateData::DataType::JOINT_TORQUE, offset);

    for(std::size_t i = 0; i < n_c; ++i){
        const cslibs_kdl_data::ContactPoint& p = contacts[i];
        KDL::Vector pos(p.position(0), p.position(1), p.position(2));
        KDL::Vector dir(p.direction(0), p.direction(1), p.direction(2));
        KDL::Wrench wi = createWrench(pos, dir);
        Eigen::VectorXd tau_i = getExternalTorques(state.position, p.header.frame_id, wi);
        reg_mat.block(0, i, tau_i.rows(), 1) = tau_i.eval();
    }

     Eigen::JacobiSVD<Eigen::MatrixXd> svd(reg_mat, Eigen::ComputeThinU | Eigen::ComputeThinV | Eigen::FullPivHouseholderQRPreconditioner);
     Eigen::VectorXd forces = svd.solve(r);
     //createWrench supposes the force is antiparallel to the normal!
     for(std::size_t i = 0; i < n_c; ++i){
         cslibs_kdl_data::ContactPoint& p = contacts[i];
         p.force = -forces(i);
     }
}


text:
// Copyright Abel Sinkovics (abel@sinkovics.hu) 2010.
// 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)

#include <boost/metaparse/util/is_letter.hpp>

#include "common.hpp"

#include <boost/mpl/apply_wrap.hpp>
#include <boost/mpl/not.hpp>
#include <boost/mpl/assert.hpp>

#include "test_case.hpp"

BOOST_METAPARSE_TEST_CASE(util_is_letter)
{
  using boost::metaparse::util::is_letter;

  using boost::mpl::apply_wrap1;
  using boost::mpl::not_;

  // test_lcase_letter
  BOOST_MPL_ASSERT((apply_wrap1<is_letter<>, char_k>));

  // test_ucase_letter
  BOOST_MPL_ASSERT((apply_wrap1<is_letter<>, char_K>));

  // test_not_letter
  BOOST_MPL_ASSERT((not_<apply_wrap1<is_letter<>, char_7> >));
}


text:
#include "dynet/dynet.h"
#include "dynet/training.h"
#include "dynet/expr.h"
#include "dynet/lstm.h"
#include "dynet/mp.h"
#include "rnnlm.h"
#include <boost/algorithm/string.hpp>

#include <iostream>
#include <fstream>
#include <vector>
/*
TODO:
- The shadow params in the trainers need to be shared.
*/

using namespace std;
using namespace dynet;
using namespace dynet::mp;
using namespace boost::interprocess;

typedef vector<int> Datum;

vector<Datum> ReadData(string filename) {
  vector<Datum> data;
  ifstream fs(filename);
  if (!fs.is_open()) {
    cerr << "ERROR: Unable to open " << filename << endl;
    exit(1);
  }
  string line;
  while (getline(fs, line)) {
    data.push_back(read_sentence(line, d));
  }
  return data;
}

template<class T, class D>
class Learner : public ILearner<D, dynet::real> {
public:
  explicit Learner(RNNLanguageModel<T>& rnnlm, unsigned data_size) : rnnlm(rnnlm) {}
  ~Learner() {}

  dynet::real LearnFromDatum(const D& datum, bool learn) {
    ComputationGraph cg;
    Expression loss_expr = rnnlm.BuildLMGraph(datum, cg);
    dynet::real loss = as_scalar(cg.forward(loss_expr));
    if (learn) {
      cg.backward(loss_expr);
    }
    return loss;
  }

  void SaveModel() {}

private:
  RNNLanguageModel<T>& rnnlm;
};

int main(int argc, char** argv) {
  if (argc < 4) {
    cerr << "Usage: " << argv[0] << " cores corpus.txt dev.txt [iterations]" << endl;
    return 1;
  }
  srand(time(NULL));
  unsigned num_children = atoi(argv[1]);
  assert (num_children <= 64);
  vector<Datum> data = ReadData(argv[2]);
  vector<Datum> dev_data = ReadData(argv[3]);
  unsigned num_iterations = (argc >= 5) ? atoi(argv[4]) : UINT_MAX;
  unsigned dev_frequency = 5000;
  unsigned report_frequency = 10;

  dynet::initialize(argc, argv, true);

  Model model;
  SimpleSGDTrainer sgd(model, 0.2);
  //AdagradTrainer sgd(model, 0.0);
  //AdamTrainer sgd(model, 0.0);

  RNNLanguageModel<LSTMBuilder> rnnlm(model);

  Learner<LSTMBuilder, Datum> learner(rnnlm, data.size());
  run_multi_process<Datum>(num_children, &learner, &sgd, data, dev_data, num_iterations, dev_frequency, report_frequency);
}

ext: cc

max_stars_repo_name: MalcolmSun/dynet

max_issues_repo_name: MalcolmSun/dynet

max_forks_repo_name: MalcolmSun/dynet


text:
//---------------------------------------------------------------------------//
//!
//! \file   Utility_RotationCartesianCoordinateConversionPolicy.hpp
//! \author Alex Robinson
//! \brief  Rotation Cartesian coordinate conversion policy declaration
//!
//---------------------------------------------------------------------------//

#ifndef UTILITY_ROTATION_CARTESIAN_COORDINATE_CONVERSION_POLICY_HPP
#define UTILITY_ROTATION_CARTESIAN_COORDINATE_CONVERSION_POLICY_HPP

// Boost Includes
#include <boost/serialization/array_wrapper.hpp>

// FRENSIE Includes
#include "Utility_CartesianSpatialCoordinateConversionPolicy.hpp"
#include "Utility_CartesianDirectionalCoordinateConversionPolicy.hpp"

namespace Utility{

//! The rotation Cartesian coordinate conversion policy class
class RotationCartesianCoordinateConversionPolicy : public CartesianSpatialCoordinateConversionPolicy,
                                                    public CartesianDirectionalCoordinateConversionPolicy
{

public:

  //! Constructor
  RotationCartesianCoordinateConversionPolicy( const double axis[3] );

  //! Destructor
  ~RotationCartesianCoordinateConversionPolicy()
  { /* ... */ }

  //! Convert the spatial coordinates to cartesian coordinates
  void convertToCartesianSpatialCoordinates(
                                      const double primary_spatial_coord,
                                      const double secondary_spatial_coord,
                                      const double tertiary_spatial_coord,
                                      double& x_spatial_coord,
                                      double& y_spatial_coord,
                                      double& z_spatial_coord ) const override;

  //! Convert the cartesian coordinates to the spatial coordinate system
  void convertFromCartesianSpatialCoordinates(
                               const double x_spatial_coord,
                               const double y_spatial_coord,
                               const double z_spatial_coord,
                               double& primary_spatial_coord,
                               double& secondary_spatial_coord,
                               double& tertiary_spatial_coord ) const override;

  //! Convert the directional coordinates to cartesian coordinates
  void convertToCartesianDirectionalCoordinates(
                                  const double primary_directional_coord,
                                  const double secondary_directional_coord,
                                  const double tertiary_directional_coord,
                                  double& x_directional_coord,
                                  double& y_directional_coord,
                                  double& z_directional_coord ) const override;

  //! Convert the cartesian coordinates to the directional coordinate system
  void convertFromCartesianDirectionalCoordinates(
                           const double x_directional_coord,
                           const double y_directional_coord,
                           const double z_directional_coord,
                           double& primary_directional_coord,
                           double& secondary_directional_coord,
                           double& tertiary_directional_coord ) const override;

  //! Convert the spatial coordinates to cartesian coordinates
  using CartesianSpatialCoordinateConversionPolicy::convertToCartesianSpatialCoordinates;

  //! Convert the cartesian coordinates to the spatial coordinate system
  using CartesianSpatialCoordinateConversionPolicy::convertFromCartesianSpatialCoordinates;

  //! Convert the directional coordinates to cartesian coordinates
  using CartesianDirectionalCoordinateConversionPolicy::convertToCartesianDirectionalCoordinates;

  //! Convert the cartesian coordinates to the directional coordinate system
  using CartesianDirectionalCoordinateConversionPolicy::convertFromCartesianDirectionalCoordinates;

private:

  // The default constructor should not be used - if the z-axis of the local
  // coordinate system w.r.t. the global coordinate system aligns with the
  // z-axis of the global coordinate system use the basic conversion policy
  RotationCartesianCoordinateConversionPolicy()
  { /* ... */ }

  // We have C-arrays as members - hide the copy constructor and assignment
  // operator
  RotationCartesianCoordinateConversionPolicy( const RotationCartesianCoordinateConversionPolicy& that );
  RotationCartesianCoordinateConversionPolicy& operator=( const RotationCartesianCoordinateConversionPolicy& that );

  // Save the policy to an archive
  template<typename Archive>
  void save( Archive& ar, const unsigned version ) const;

  // Load the policy from an archive
  template<typename Archive>
  void load( Archive& ar, const unsigned version );

  BOOST_SERIALIZATION_SPLIT_MEMBER();

  // Declare the boost serialization access object as a friend
  friend class boost::serialization::access;

  // The z-axis (unit vector) of the local Cartesian coordinate system w.r.t.
  // the global Cartesian coordinate system
  double d_axis[3];
};

// Save the policy to an archive
template<typename Archive>
void RotationCartesianCoordinateConversionPolicy::save( Archive& ar, const unsigned version ) const
{
  // Save the base classes
  ar & BOOST_SERIALIZATION_BASE_OBJECT_NVP( CartesianSpatialCoordinateConversionPolicy );
  ar & BOOST_SERIALIZATION_BASE_OBJECT_NVP( CartesianDirectionalCoordinateConversionPolicy );

  // Save the local data
  ar & boost::serialization::make_nvp( "d_axis", boost::serialization::make_array( d_axis, 3 ) );
}

// Load the policy from an archive
template<typename Archive>
void RotationCartesianCoordinateConversionPolicy::load( Archive& ar, const unsigned version )
{
  // Load the base class
  ar & BOOST_SERIALIZATION_BASE_OBJECT_NVP( CartesianSpatialCoordinateConversionPolicy );
  ar & BOOST_SERIALIZATION_BASE_OBJECT_NVP( CartesianDirectionalCoordinateConversionPolicy );

  // Load the local data
  ar & boost::serialization::make_nvp( "d_axis", boost::serialization::make_array( d_axis, 3 ) );
}

} // end Utility namespace

BOOST_SERIALIZATION_CLASS_VERSION( RotationCartesianCoordinateConversionPolicy, Utility, 0 );
BOOST_SERIALIZATION_CLASS_EXPORT_STANDARD_KEY( RotationCartesianCoordinateConversionPolicy, Utility );
EXTERN_EXPLICIT_CLASS_SAVE_LOAD_INST( Utility, RotationCartesianCoordinateConversionPolicy );

#endif // end UTILITY_ROTATION_CARTESIAN_COORDINATE_CONVERSION_POLICY_HPP

//---------------------------------------------------------------------------//
// end Utility_RotationCartesianCoordinateConversionPolicy.hpp
//---------------------------------------------------------------------------//


text:
#ifndef BOOST_ARCHIVE_CODECVT_NULL_HPP
#define BOOST_ARCHIVE_CODECVT_NULL_HPP

// MS compatible compilers support #pragma once
#if defined(_MSC_VER) && (_MSC_VER >= 1020)
# pragma once
#endif

/////////1/////////2/////////3/////////4/////////5/////////6/////////7/////////8
// codecvt_null.hpp:

// (C) Copyright 2004 Robert Ramey - http://www.rrsd.com .
// Use, modification and distribution is subject to the Boost Software
// License, Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)

//  See http://www.boost.org for updates, documentation, and revision history.

#include <locale>
#include <cstddef>

#include <boost/config.hpp>
#include <boost/archive/detail/auto_link_archive.hpp>

namespace std{
    #if defined(__LIBCOMO__)
        using ::mbstate_t;
    #elif defined(__QNXNTO__)
        //using std::mbstate_t;
    #elif defined(BOOST_DINKUMWARE_STDLIB) && BOOST_DINKUMWARE_STDLIB == 1
        using ::mbstate_t;
    #elif defined(__SGI_STL_PORT)
    #elif defined(BOOST_NO_STDC_NAMESPACE)
        using ::codecvt;
        using ::mbstate_t;
    #elif defined(BOOST_RWSTD_VER)
        using ::mbstate_t;
   #endif
} // namespace std

namespace boost {
namespace archive {

template<class Ch>
class codecvt_null;

template<>
class codecvt_null<char> : public std::codecvt<char, char, std::mbstate_t>
{
    virtual bool do_always_noconv() const throw() {
        return true;
    }
public:
    explicit codecvt_null(std::size_t no_locale_manage = 0) :
        std::codecvt<char, char, std::mbstate_t>(no_locale_manage)
    {}
};

template<>
class codecvt_null<wchar_t> : public std::codecvt<wchar_t, char, std::mbstate_t>
{
    virtual BOOST_ARCHIVE_DECL(std::codecvt_base::result)
    do_out(
        std::mbstate_t & state,
        const wchar_t * first1,
        const wchar_t * last1,
        const wchar_t * & next1,
        char * first2,
        char * last2,
        char * & next2
    ) const;
    virtual BOOST_ARCHIVE_DECL(std::codecvt_base::result)
    do_in(
        std::mbstate_t & state,
        const char * first1,
        const char * last1,
        const char * & next1,
        wchar_t * first2,
        wchar_t * last2,
        wchar_t * & next2
    ) const;
    virtual int do_encoding( ) const throw( ){
        return sizeof(wchar_t) / sizeof(char);
    }
    virtual int do_max_length( ) const throw( ){
        return do_encoding();
    }
};

} // namespace archive
} // namespace boost

#endif //BOOST_ARCHIVE_CODECVT_NULL_HPP


text:
#include <boost/filesystem.hpp>
#include <boost/program_options.hpp>
#include <iostream>

#include <Eigen/Dense>
#include <algorithm>

#include "aux/eigen2hdf.hpp"
#include "quadrature/qhermite.hpp"
#include "spectral/basis/spectral_basis_factory_hermite.hpp"
#include "spectral/basis/spectral_basis_factory_ks.hpp"
#include "spectral/basis/toolbox/spectral_basis.hpp"
#include "spectral/h2n_1d.hpp"
#include "spectral/hermite_to_nodal.hpp"
#include "spectral/polar_to_hermite.hpp"

//#include <deal.II/base/exceptions.h>

using namespace std;
using namespace boltzmann;

namespace po = boost::program_options;

int main(int argc, char* argv[])
{
  typedef Eigen::VectorXd vec_t;

  int K = atoi(argv[1]);

  // create basis files
  typedef typename SpectralBasisFactoryHN::basis_type hermite_basis_t;
  hermite_basis_t hermite_basis;
  SpectralBasisFactoryHN::create(hermite_basis, K, 2);
  cout << "Hermite basis size: " << hermite_basis.n_dofs() << endl;

  // hermite to nodal
  typedef Eigen::MatrixXd mat_t;
  Hermite2Nodal<hermite_basis_t> h2n(
      hermite_basis, K, [K](mat_t& m1, mat_t& m2) { H2N_1d<>::create(m1, m2, K); });
  Hermite2Nodal<hermite_basis_t> h2ng(
      hermite_basis, K, [K](mat_t& m1, mat_t& m2) { H2NG_1d::create(m1, m2, K, 1.0); });

  auto& M = h2n.get_h2n();

  auto& M2 = h2ng.get_h2n();

  auto diff = M - M2;
  diff = diff.select(abs(diff) < 1e-15, Eigen::MatrixXd::Zero(diff.rows(), diff.cols()), diff);

  cout << "difference matrix\n";
  cout << diff << endl;
  // cout << "M := H2N-matrix\n";
  // cout << "M.T * M\n";
  // cout << M.transpose()* M  << endl;

  // fname = "coefficients.h5";
  // if(!boost::filesystem::exists(fname)) {
  //   cout << endl << fname << " does not exist. Abort!\n";
  //   return 1;
  // }
  // vec_t cp(N); // polar coefficients
  // hid_t h5_init = H5Fopen(fname.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
  // eigen2hdf::load(h5_init, "coeffs", cp);
  // H5Fclose(h5_init);

  // AssertDimension(polar_basis.n_dofs(), cp.size());

  hid_t h5f = H5Fcreate("test.h5", H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
  eigen2hdf::save(h5f, "M", M);
  eigen2hdf::save(h5f, "Mg", M2);
  // eigen2hdf::save(h5f, "coeffs_nodal", cn);
  // eigen2hdf::save(h5f, "coeffs_hermite", ch);
  // eigen2hdf::save(h5f, "coeffs", cp2);
  H5Fclose(h5f);

  return 0;
}


text:
#include <boost/make_shared.hpp>
#include <boost/algorithm/string.hpp>
#include <boost/algorithm/string/predicate.hpp>

#include "TextDisplay.h"
#include "TextRenderer.h"

#include "../MessageSystem.h"

namespace Kriti {
namespace GUI {

TextDisplay::TextDisplay(Math::Vector minSize, Math::Vector stretch,
    boost::shared_ptr<Font::Instance> font, Math::Vector colour,
    VerticalAlignment valign) : Widget(stretch), m_minSize(minSize),
    m_font(font), m_colour(colour), m_valign(valign) {

    m_regen = false;
    m_reflow = false;
}

Math::Vector TextDisplay::minSize() {
    if(m_reflow) {
        flowText(m_minSize.x());
        m_minSize.setY(m_font->lineSpacing() * m_lines.size());
    }
    return m_minSize;
}

void TextDisplay::fill(
    boost::shared_ptr<Render::RenderableContainer> container) {

    if(m_renderable) container->addTransparent(m_renderable);
}

void TextDisplay::flush(
    boost::shared_ptr<Render::RenderableContainer> container) {

    if(m_renderable) container->removeTransparent(m_renderable);
}

void TextDisplay::updated(
    boost::shared_ptr<OutlineRegistry> __attribute__((unused)) registry,
        Math::Vector clipStart, Math::Vector clipEnd) {

    Math::Vector outlineStart, outlineEnd;
    standardOutlineUpdate(nullptr, clipStart, clipEnd, outlineStart,
        outlineEnd);

    double spacing = m_font->lineSpacing();

    if(!m_renderable) m_renderable = boost::make_shared<Render::Renderable>();
    if(m_regen || m_lastScale != scale()) {
        m_renderable->clearRenderSequences();
        m_regen = false;
        // XXX: text should have been flowed already by call to minSize().
        if(m_reflow) {
            flowText(size().x());
            m_reflow = false;
        }

        Math::Vector sv(0,-spacing);
        Math::Vector offset(0,
            spacing * (m_lines.size()-1) - m_font->baselineOffset());

        for(auto &line : m_lines) {
            auto rline = TextRenderer().renderString(m_font, line, m_colour,
                scale());

            auto seq = rline->renderSequence(0);
            Math::AffineTransformation at;
            at.translate(offset);
            seq->sequenceTransform() *= at.matrix();

            m_renderable->addRenderSequence(rline->renderSequence(0));

            offset += sv;
        }
    }

    // update clip regions
    for(int i = 0; i < m_renderable->renderSequenceCount(); i ++) {
        m_renderable->renderSequence(i)->extraParams().setParam(
            "gui_clip_start", outlineStart);
        m_renderable->renderSequence(i)->extraParams().setParam(
            "gui_clip_end", outlineEnd);
    }

    m_lastScale = scale();

    double extra = size().y() - spacing*m_lines.size();
    switch(m_valign) {
    case Top:
        m_renderable->location() = pos();
        break;
    case VCentre:
        m_renderable->location() = pos() - Math::Vector(0,extra/2);
        break;
    case Bottom:
        m_renderable->location() = pos() - Math::Vector(0,extra);
        break;
    }
}

void TextDisplay::flowText(double overall_width) {
    m_lines.clear();
    std::vector<std::string> paragraphs;
    boost::split(paragraphs, m_text, boost::algorithm::is_any_of("\n\r"),
        boost::algorithm::token_compress_off);

    for(auto &paragraph : paragraphs) {
        std::vector<std::string> words;
        boost::split(words, paragraph,
            boost::algorithm::is_any_of("\t "),
            boost::algorithm::token_compress_on);

        std::string line;
        double running_width = 0;
        for(auto word : words) {
            if(running_width != 0) word = " " + word;
            Math::Vector ul, lr;
            TextRenderer().sizeString(m_font, word, ul, lr);
            double width = lr.x() - ul.x();
            // can we fit it on the current line?
            if(running_width + width <= overall_width) {
                line += word;
                running_width += width;
                continue;
            }
            // remove the space
            word = word.substr(1);
            TextRenderer().sizeString(m_font, word, ul, lr);
            running_width = width;
            m_lines.push_back(line);
            line = word;
        }
        m_lines.push_back(line);
    }
}

}  // namespace GUI
}  // namespace Kriti

max_stars_repo_name: etherealvisage/kriti

max_issues_repo_name: etherealvisage/kriti

max_forks_repo_name: etherealvisage/kriti


text:
/*
 * LaplaceM2L.cpp
 *
 *  Created on: Oct 12, 2016
 *      Author: wyan
 */

#include "SVD_pvfmm.hpp"

#include <Eigen/Dense>

#include <algorithm>
#include <chrono>
#include <iomanip>
#include <iostream>

#define DIRECTLAYER 2
#define PI314 (static_cast<double>(3.1415926535897932384626433))

namespace Laplace3D3D {

using EVec3 = Eigen::Vector3d;
using EVec4 = Eigen::Vector4d;
constexpr double eps = 1e-10;

// real and wave sum of 2D Laplace kernel Ewald

inline double freal(double xi, double r) { return std::erfc(xi * r) / r; }

inline double frealp(double xi, double r) {
    return -(2. * exp(-r * r * (xi * xi)) * xi) / (sqrt(M_PI) * r) - std::erfc(r * xi) / (r * r);
}

// xm: target, xn: source
inline double realSum(const double xi, const EVec3 &xn, const EVec3 &xm) {
    EVec3 rmn = xm - xn;
    double rnorm = rmn.norm();
    if (rnorm < eps) {
        return 0;
    }
    return freal(xi, rnorm);
}

inline double potEwald(const EVec3 &xm, const EVec3 &xn) {
    const double xi = 2; // recommend for box=1 to get machine precision
    EVec3 target = xm;
    EVec3 source = xn;
    target[0] = target[0] - floor(target[0]); // periodic BC
    target[1] = target[1] - floor(target[1]);
    target[2] = target[2] - floor(target[2]);
    source[0] = source[0] - floor(source[0]);
    source[1] = source[1] - floor(source[1]);
    source[2] = source[2] - floor(source[2]);

    // real sum
    int rLim = 4;
    double Kreal = 0;
    for (int i = -rLim; i <= rLim; i++) {
        for (int j = -rLim; j <= rLim; j++) {
            for (int k = -rLim; k <= rLim; k++) {
                Kreal += realSum(xi, target, source - EVec3(i, j, k));
            }
        }
    }

    // wave sum
    int wLim = 4;
    double Kwave = 0;
    EVec3 rmn = target - source;
    const double xi2 = xi * xi;
    const double rmnnorm = rmn.norm();
    for (int i = -wLim; i <= wLim; i++) {
        for (int j = -wLim; j <= wLim; j++) {
            for (int k = -wLim; k <= wLim; k++) {
                if (i == 0 && j == 0 && k == 0) {
                    continue;
                }
                EVec3 kvec = EVec3(i, j, k) * (2 * PI314);
                double k2 = kvec.dot(kvec);
                Kwave += 4 * PI314 * cos(kvec.dot(rmn)) * exp(-k2 / (4 * xi2)) / k2;
            }
        }
    }

    double Kself = rmnnorm < 1e-10 ? -2 * xi / sqrt(PI314) : 0;

    return Kreal + Kwave + Kself - PI314 / xi2;
}

inline void realGradSum(double xi, const EVec3 &target, const EVec3 &source, EVec3 &v) {
    EVec3 rvec = target - source;
    double rnorm = rvec.norm();
    if (rnorm < eps) {
        v.setZero();
    } else {
        v = (frealp(xi, rnorm) / rnorm) * rvec;
    }
}

inline void gradkernel(const EVec3 &target, const EVec3 &source, EVec3 &answer) {
    EVec3 rst = target - source;
    double rnorm = rst.norm();
    if (rnorm < eps) {
        answer.setZero();
        return;
    }
    double rnorm3 = rnorm * rnorm * rnorm;
    answer = -rst / rnorm3;
}

// grad of Laplace potential, without 1/4pi prefactor, periodic of -r_k/r^3
inline void gradEwald(const EVec3 &target_, const EVec3 &source_, EVec3 &answer) {
    EVec3 target = target_;
    EVec3 source = source_;
    target[0] = target[0] - floor(target[0]); // periodic BC
    target[1] = target[1] - floor(target[1]);
    target[2] = target[2] - floor(target[2]);
    source[0] = source[0] - floor(source[0]);
    source[1] = source[1] - floor(source[1]);
    source[2] = source[2] - floor(source[2]);

    double xi = 0.54;

    // real sum
    int rLim = 10;
    EVec3 Kreal = EVec3::Zero();
    for (int i = -rLim; i < rLim + 1; i++) {
        for (int j = -rLim; j < rLim + 1; j++) {
            for (int k = -rLim; k < rLim + 1; k++) {
                EVec3 v = EVec3::Zero();
                realGradSum(xi, target, source + EVec3(i, j, k), v);
                Kreal += v;
            }
        }
    }

    // wave sum
    int wLim = 10;
    EVec3 rmn = target - source;
    double xi2 = xi * xi;
    EVec3 Kwave(0., 0., 0.);
    for (int i = -wLim; i < wLim + 1; i++) {
        for (int j = -wLim; j < wLim + 1; j++) {
            for (int k = -wLim; k < wLim + 1; k++) {
                if (i == 0 && j == 0 && k == 0)
                    continue;
                EVec3 kvec = EVec3(i, j, k) * (2 * M_PI);
                double k2 = kvec.dot(kvec);
                double knorm = kvec.norm();
                Kwave += -kvec * (sin(kvec.dot(rmn)) * exp(-k2 / (4 * xi2)) / k2);
            }
        }
    }

    answer = Kreal + Kwave;
}

inline double pot(const EVec3 &target, const EVec3 &source) {
    EVec3 rst = target - source;
    double rnorm = rst.norm();
    return rnorm < eps ? 0 : 1 / rnorm;
}

inline EVec4 gKernel(const EVec3 &target, const EVec3 &source) {
    EVec3 rst = target - source;
    EVec4 pgrad = EVec4::Zero();
    double rnorm = rst.norm();
    if (rnorm < eps) {
        pgrad.setZero();
    } else {
        pgrad[0] = 1 / rnorm;
        double rnorm3 = rnorm * rnorm * rnorm;
        pgrad.block<3, 1>(1, 0) = -rst / rnorm3;
    }
    return pgrad;
}

inline EVec4 gKernelEwald(const EVec3 &target, const EVec3 &source) {
    EVec4 pgrad = EVec4::Zero();
    pgrad[0] = potEwald(target, source);
    EVec3 grad;
    gradEwald(target, source, grad);
    pgrad[1] = grad[0];
    pgrad[2] = grad[1];
    pgrad[3] = grad[2];
    return pgrad;
}

inline EVec4 gKernelNF(const EVec3 &target, const EVec3 &source, int N = DIRECTLAYER) {
    EVec4 gNF = EVec4::Zero();
    for (int i = -N; i < N + 1; i++) {
        for (int j = -N; j < N + 1; j++) {
            for (int k = -N; k < N + 1; k++) {
                EVec4 gFree = gKernel(target, source + EVec3(i, j, k));
                gNF += gFree;
            }
        }
    }
    return gNF;
}

// Out of Direct Sum Layer, far field part
inline EVec4 gKernelFF(const EVec3 &target, const EVec3 &source) {
    EVec4 fEwald = gKernelEwald(target, source);
    fEwald -= gKernelNF(target, source);
    return fEwald;
}

inline double potNF(const EVec3 &target, const EVec3 &source, int N = DIRECTLAYER) {
    double gNF = 0;
    for (int i = -N; i < N + 1; i++) {
        for (int j = -N; j < N + 1; j++) {
            for (int k = -N; k < N + 1; k++) {
                gNF += pot(target, source + EVec3(i, j, k));
            }
        }
    }
    return gNF;
}

// Out of Direct Sum Layer, far field part
inline double potFF(const EVec3 &target, const EVec3 &source) {
    double fEwald = potEwald(target, source);
    fEwald -= potNF(target, source);
    return fEwald;
}

/**
 * \brief Returns the coordinates of points on the surface of a cube.
 * \param[in] p Number of points on an edge of the cube is (n+1)
 * \param[in] c Coordinates to the centre of the cube (3D array).
 * \param[in] alpha Scaling factor for the size of the cube.
 * \param[in] depth Depth of the cube in the octree.
 * \return Vector with coordinates of points on the surface of the cube in the
 * format [x0 y0 z0 x1 y1 z1 .... ].
 */

template <class Real_t>
std::vector<Real_t> surface(int p, Real_t *c, Real_t alpha, int depth) {
    size_t n_ = (6 * (p - 1) * (p - 1) + 2); // Total number of points.

    std::vector<Real_t> coord(n_ * 3);
    coord[0] = coord[1] = coord[2] = -1.0;
    size_t cnt = 1;
    for (int i = 0; i < p - 1; i++)
        for (int j = 0; j < p - 1; j++) {
            coord[cnt * 3] = -1.0;
            coord[cnt * 3 + 1] = (2.0 * (i + 1) - p + 1) / (p - 1);
            coord[cnt * 3 + 2] = (2.0 * j - p + 1) / (p - 1);
            cnt++;
        }
    for (int i = 0; i < p - 1; i++)
        for (int j = 0; j < p - 1; j++) {
            coord[cnt * 3] = (2.0 * i - p + 1) / (p - 1);
            coord[cnt * 3 + 1] = -1.0;
            coord[cnt * 3 + 2] = (2.0 * (j + 1) - p + 1) / (p - 1);
            cnt++;
        }
    for (int i = 0; i < p - 1; i++)
        for (int j = 0; j < p - 1; j++) {
            coord[cnt * 3] = (2.0 * (i + 1) - p + 1) / (p - 1);
            coord[cnt * 3 + 1] = (2.0 * j - p + 1) / (p - 1);
            coord[cnt * 3 + 2] = -1.0;
            cnt++;
        }
    for (size_t i = 0; i < (n_ / 2) * 3; i++)
        coord[cnt * 3 + i] = -coord[i];

    Real_t r = 0.5 * pow(0.5, depth);
    Real_t b = alpha * r;
    for (size_t i = 0; i < n_; i++) {
        coord[i * 3 + 0] = (coord[i * 3 + 0] + 1.0) * b + c[0];
        coord[i * 3 + 1] = (coord[i * 3 + 1] + 1.0) * b + c[1];
        coord[i * 3 + 2] = (coord[i * 3 + 2] + 1.0) * b + c[2];
    }
    return coord;
}

int main(int argc, char **argv) {
    Eigen::initParallel();
    Eigen::setNbThreads(1);

    std::cout << std::scientific << std::setprecision(18);

    std::chrono::high_resolution_clock::time_point t1 = std::chrono::high_resolution_clock::now();
    const int pEquiv = atoi(argv[1]); // (8-1)^2*6 + 2 points
    const int pCheck = atoi(argv[1]);
    const double scaleEquiv = 1.05;
    const double scaleCheck = 2.95;
    const double pCenterEquiv[3] = {-(scaleEquiv - 1) / 2, -(scaleEquiv - 1) / 2, -(scaleEquiv - 1) / 2};
    const double pCenterCheck[3] = {-(scaleCheck - 1) / 2, -(scaleCheck - 1) / 2, -(scaleCheck - 1) / 2};

    const double scaleLEquiv = 1.05;
    const double scaleLCheck = 2.95;
    const double pCenterLEquiv[3] = {-(scaleLEquiv - 1) / 2, -(scaleLEquiv - 1) / 2, -(scaleLEquiv - 1) / 2};
    const double pCenterLCheck[3] = {-(scaleLCheck - 1) / 2, -(scaleLCheck - 1) / 2, -(scaleLCheck - 1) / 2};

    auto pointMEquiv = surface(pEquiv, (double *)&(pCenterEquiv[0]), scaleEquiv, 0);
    // center at 0.5,0.5,0.5, periodic box 1,1,1, scale 1.05, depth = 0
    auto pointMCheck = surface(pCheck, (double *)&(pCenterCheck[0]), scaleCheck, 0);
    // center at 0.5,0.5,0.5, periodic box 1,1,1, scale 1.05, depth = 0

    auto pointLEquiv = surface(pEquiv, (double *)&(pCenterLCheck[0]), scaleLCheck, 0);
    // center at 0.5,0.5,0.5, periodic box 1,1,1, scale 1.05, depth = 0
    auto pointLCheck = surface(pCheck, (double *)&(pCenterLEquiv[0]), scaleLEquiv, 0);
    // center at 0.5,0.5,0.5, periodic box 1,1,1, scale 1.05, depth = 0

    // calculate the operator M2L with least square
    const int equivN = pointMEquiv.size() / 3;
    const int checkN = pointLCheck.size() / 3;
    Eigen::MatrixXd M2L(equivN, equivN); // Laplace, 1->1

    // Aup for solving MEquiv
    Eigen::MatrixXd Aup(checkN, equivN);
    Eigen::MatrixXd AuppinvU(Aup.cols(), Aup.rows());
    Eigen::MatrixXd AuppinvVT(Aup.cols(), Aup.rows());
    for (int k = 0; k < checkN; k++) {
        Eigen::Vector3d Cpoint(pointMCheck[3 * k], pointMCheck[3 * k + 1], pointMCheck[3 * k + 2]);
        for (int l = 0; l < equivN; l++) {
            const Eigen::Vector3d Lpoint(pointMEquiv[3 * l], pointMEquiv[3 * l + 1], pointMEquiv[3 * l + 2]);
            Aup(k, l) = pot(Cpoint, Lpoint);
        }
    }
    pinv(Aup, AuppinvU, AuppinvVT);

    // condition number
    Eigen::JacobiSVD<Eigen::MatrixXd> svd(Aup);
    double cond = svd.singularValues()(0) / svd.singularValues()(svd.singularValues().size() - 1);
    std::cout << cond << std::endl;
    std::cout << "s:" << svd.singularValues() << std::endl;

    // Adown for solving LEquiv
    Eigen::MatrixXd Adown(checkN, equivN);
    Eigen::MatrixXd AdownpinvU(Adown.cols(), Adown.rows());
    Eigen::MatrixXd AdownpinvVT(Adown.cols(), Adown.rows());
    for (int k = 0; k < checkN; k++) {
        Eigen::Vector3d Cpoint(pointLCheck[3 * k], pointLCheck[3 * k + 1], pointLCheck[3 * k + 2]);
        for (int l = 0; l < equivN; l++) {
            const Eigen::Vector3d Lpoint(pointLEquiv[3 * l], pointLEquiv[3 * l + 1], pointLEquiv[3 * l + 2]);
            Adown(k, l) = pot(Cpoint, Lpoint);
        }
    }
    pinv(Adown, AdownpinvU, AdownpinvVT);

#pragma omp parallel for
    for (int i = 0; i < equivN; i++) {
        const Eigen::Vector3d Mpoint(pointMEquiv[3 * i], pointMEquiv[3 * i + 1], pointMEquiv[3 * i + 2]);
        const EVec3 Npoint(0.5, 0.5, 0.5);
        Eigen::VectorXd f(checkN);
        f.setZero();
        for (int k = 0; k < checkN; k++) {
            EVec3 Cpoint(pointLCheck[3 * k], pointLCheck[3 * k + 1], pointLCheck[3 * k + 2]);
            f[k] = potFF(Cpoint, Mpoint) - potFF(Cpoint, Npoint);
        }

        M2L.col(i) = (AdownpinvU.transpose() * (AdownpinvVT.transpose() * f));
    }
    std::chrono::high_resolution_clock::time_point t2 = std::chrono::high_resolution_clock::now();
    auto duration = std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count();

    // dump M2L
    for (int i = 0; i < equivN; i++) {
        for (int j = 0; j < equivN; j++) {
            std::cout << i << " " << j << " " << std::scientific << std::setprecision(18) << M2L(i, j) << std::endl;
        }
    }

    std::cout << "Precomputing time:" << duration / 1e6 << std::endl;

    // Test
    EVec3 center(0.6, 0.5, 0.5);
    std::vector<Eigen::Vector3d, Eigen::aligned_allocator<Eigen::Vector3d>> chargePoint(2);
    std::vector<double> chargeValue(2);
    chargePoint[0] = center + EVec3(0.1, 0, 0);
    chargeValue[0] = 1;
    chargePoint[1] = center + EVec3(-0.1, 0., 0.);
    chargeValue[1] = -1;
    // chargePoint[2] = center + EVec3(-0.2, 0., 0.);
    // chargeValue[2] = 1;

    // solve M
    Eigen::VectorXd f(checkN);
    for (int k = 0; k < checkN; k++) {
        double temp = 0;
        Eigen::Vector3d Cpoint(pointMCheck[3 * k], pointMCheck[3 * k + 1], pointMCheck[3 * k + 2]);
        for (size_t p = 0; p < chargePoint.size(); p++) {
            temp = temp + pot(Cpoint, chargePoint[p]) * (chargeValue[p]);
        }
        f[k] = temp;
    }
    Eigen::VectorXd Msource = (AuppinvU.transpose() * (AuppinvVT.transpose() * f));

    std::cout << "Msource " << Msource << std::endl;

    std::cout << "backward error: " << f - Aup * Msource << std::endl;

    // check dipole moment
    {
        EVec3 dipole = EVec3::Zero();
        for (int i = 0; i < chargePoint.size(); i++) {
            dipole += chargeValue[i] * chargePoint[i];
        }
        std::cout << "charge dipole " << dipole.transpose() << std::endl;
    }
    {
        EVec3 dipole = EVec3::Zero();
        for (int i = 0; i < equivN; i++) {
            Eigen::Vector3d Mpoint(pointMEquiv[3 * i], pointMEquiv[3 * i + 1], pointMEquiv[3 * i + 2]);
            dipole += Mpoint * Msource[i];
        }
        std::cout << "Mequiv dipole " << dipole.transpose() << std::endl;
    }

    Eigen::VectorXd M2Lsource = M2L * (Msource);

    for (int is = 0; is < 5; is++) {

        Eigen::Vector3d samplePoint = EVec3::Random() * 0.2 + EVec3(0.5, 0.5, 0.5);

        EVec4 UFFL2T = EVec4::Zero();
        EVec4 UFFS2T = EVec4::Zero();
        EVec4 UFFM2T = EVec4::Zero();

#pragma omp sections
        {
#pragma omp section
            for (int p = 0; p < chargePoint.size(); p++) {
                UFFS2T += gKernelFF(samplePoint, chargePoint[p]) * chargeValue[p];
            }

#pragma omp section
            for (int p = 0; p < equivN; p++) {
                Eigen::Vector3d Lpoint(pointLEquiv[3 * p], pointLEquiv[3 * p + 1], pointLEquiv[3 * p + 2]);
                UFFL2T += gKernel(samplePoint, Lpoint) * M2Lsource[p];
            }

#pragma omp section
            for (int p = 0; p < equivN; p++) {
                Eigen::Vector3d Mpoint(pointMEquiv[3 * p], pointMEquiv[3 * p + 1], pointMEquiv[3 * p + 2]);
                UFFM2T += gKernelFF(samplePoint, Mpoint) * Msource[p];
            }
        }
        std::cout << std::scientific << std::setprecision(10);
        std::cout << "-----------------------------------------------" << std::endl;
        std::cout << "samplePoint:" << samplePoint.transpose() << std::endl;
        std::cout << "UFF S2T: " << UFFS2T.transpose() << std::endl;
        std::cout << "UFF M2T: " << UFFM2T.transpose() << std::endl;
        std::cout << "UFF L2T: " << UFFL2T.transpose() << std::endl;
        std::cout << "Error M2T: " << (UFFM2T - UFFS2T).transpose() << std::endl;
        std::cout << "Error L2T: " << (UFFL2T - UFFS2T).transpose() << std::endl;
    }

    return 0;
}

} // namespace Laplace3D3D

#undef DIRECTLAYER
#undef PI314