text:
/*  Sirikata
 *  LibproxProximity.cpp
 *
 *  Copyright (c) 2009, Ewen Cheslack-Postava
 *  All rights reserved.
 *
 *  Redistribution and use in source and binary forms, with or without
 *  modification, are permitted provided that the following conditions are
 *  met:
 *  * Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *  * Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in
 *    the documentation and/or other materials provided with the
 *    distribution.
 *  * Neither the name of Sirikata nor the names of its contributors may
 *    be used to endorse or promote products derived from this software
 *    without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
 * IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
 * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
 * PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
 * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
 * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

#include "LibproxProximity.hpp"
#include "Options.hpp"
#include <sirikata/core/options/CommonOptions.hpp>

#include <algorithm>

#include <sirikata/space/QueryHandlerFactory.hpp>

#include "Protocol_Prox.pbj.hpp"
#include "Protocol_ServerProx.pbj.hpp"

#include <sirikata/core/network/IOServiceFactory.hpp>

#include <sirikata/space/AggregateManager.hpp>

// Property tree for old API for queries
#include <boost/property_tree/ptree.hpp>
#include <boost/property_tree/json_parser.hpp>

#define PROXLOG(level,msg) SILOG(prox,level,"[PROX] " << msg)

namespace Sirikata {

static SolidAngle NoUpdateSolidAngle = SolidAngle(0.f);
// Note that this is different from sentinals indicating infinite
// results. Instead, this indicates that we shouldn't even request the update,
// leaving it as is, because we don't actually have a new value.
static uint32 NoUpdateMaxResults = ((uint32)INT_MAX)+1;

namespace {

bool parseQueryRequest(const String& query, SolidAngle* qangle_out, uint32* max_results_out) {
    if (query.empty())
        return false;

    using namespace boost::property_tree;
    ptree pt;
    try {
        std::stringstream data_json(query);
        read_json(data_json, pt);
    }
    catch(json_parser::json_parser_error exc) {
        return false;
    }

    if (pt.find("angle") != pt.not_found())
        *qangle_out = SolidAngle( pt.get<float32>("angle") );
    else
        *qangle_out = SolidAngle::Max;

    if (pt.find("max_results") != pt.not_found())
        *max_results_out = pt.get<uint32>("max_results");
    else
        *max_results_out = 0;

    return true;
}

}

LibproxProximity::LibproxProximity(SpaceContext* ctx, LocationService* locservice, CoordinateSegmentation* cseg, SpaceNetwork* net, AggregateManager* aggmgr)
 : LibproxProximityBase(ctx, locservice, cseg, net, aggmgr),
   mServerQuerier(NULL),
   mDistanceQueryDistance(0.f),
   mMaxObject(0.0f),
   mMinObjectQueryAngle(SolidAngle::Max),
   mMaxMaxCount(1),
   mServerQueries(),
   mServerDistance(false),
   mServerHandlerPoller(mProxStrand, std::tr1::bind(&LibproxProximity::tickQueryHandler, this, mServerQueryHandler), Duration::milliseconds((int64)100)),
   mObjectQueries(),
   mObjectDistance(false),
   mObjectHandlerPoller(mProxStrand, std::tr1::bind(&LibproxProximity::tickQueryHandler, this, mObjectQueryHandler), Duration::milliseconds((int64)100)),
   mStaticRebuilderPoller(mProxStrand, std::tr1::bind(&LibproxProximity::rebuildHandler, this, OBJECT_CLASS_STATIC), Duration::seconds(3600.f)),
   mDynamicRebuilderPoller(mProxStrand, std::tr1::bind(&LibproxProximity::rebuildHandler, this, OBJECT_CLASS_DYNAMIC), Duration::seconds(3600.f))
{
    using std::tr1::placeholders::_1;
    using std::tr1::placeholders::_2;
    using std::tr1::placeholders::_3;
    using std::tr1::placeholders::_4;
    using std::tr1::placeholders::_5;

    // Server Querier (discover other servers)
    String pinto_type = GetOptionValue<String>(OPT_PINTO);
    String pinto_options = GetOptionValue<String>(OPT_PINTO_OPTIONS);
    mServerQuerier = PintoServerQuerierFactory::getSingleton().getConstructor(pinto_type)(mContext, pinto_options);
    mServerQuerier->addListener(this);

    // Generic query parameters
    mDistanceQueryDistance = GetOptionValue<float32>(OPT_PROX_QUERY_RANGE);

    // Server Queries
    String server_handler_type = GetOptionValue<String>(OPT_PROX_SERVER_QUERY_HANDLER_TYPE);
    String server_handler_options = GetOptionValue<String>(OPT_PROX_SERVER_QUERY_HANDLER_OPTIONS);
    for(int i = 0; i < NUM_OBJECT_CLASSES; i++) {
        if (i >= mNumQueryHandlers) {
            mServerQueryHandler[i] = NULL;
            continue;
        }
        mServerQueryHandler[i] = QueryHandlerFactory<ObjectProxSimulationTraits>(server_handler_type, server_handler_options);
        mServerQueryHandler[i]->setAggregateListener(this); // *Must* be before handler->initialize
        bool server_static_objects = (mSeparateDynamicObjects && i == OBJECT_CLASS_STATIC);
        mServerQueryHandler[i]->initialize(
            mLocCache, mLocCache, server_static_objects,
            std::tr1::bind(&LibproxProximity::handlerShouldHandleObject, this, server_static_objects, false, _1, _2, _3, _4, _5)
        );
    }
    if (server_handler_type == "dist" || server_handler_type == "rtreedist") mServerDistance = true;

    // Object Queries
    String object_handler_type = GetOptionValue<String>(OPT_PROX_OBJECT_QUERY_HANDLER_TYPE);
    String object_handler_options = GetOptionValue<String>(OPT_PROX_OBJECT_QUERY_HANDLER_OPTIONS);
    for(int i = 0; i < NUM_OBJECT_CLASSES; i++) {
        if (i >= mNumQueryHandlers) {
            mObjectQueryHandler[i] = NULL;
            continue;
        }
        mObjectQueryHandler[i] = QueryHandlerFactory<ObjectProxSimulationTraits>(object_handler_type, object_handler_options);
        mObjectQueryHandler[i]->setAggregateListener(this); // *Must* be before handler->initialize
        bool object_static_objects = (mSeparateDynamicObjects && i == OBJECT_CLASS_STATIC);
        mObjectQueryHandler[i]->initialize(
            mLocCache, mLocCache, object_static_objects,
            std::tr1::bind(&LibproxProximity::handlerShouldHandleObject, this, object_static_objects, true, _1, _2, _3, _4, _5)
        );
    }
    if (object_handler_type == "dist" || object_handler_type == "rtreedist") mObjectDistance = true;
}

LibproxProximity::~LibproxProximity() {
    for(int i = 0; i < NUM_OBJECT_CLASSES; i++) {
        delete mObjectQueryHandler[i];
        delete mServerQueryHandler[i];
    }

    delete mServerQuerier;
}

// MAIN Thread Methods: The following should only be called from the main thread.

void LibproxProximity::start() {
    Proximity::start();

    // Always initialize with CSeg's current size
    BoundingBoxList bboxes = mCSeg->serverRegion(mContext->id());
    BoundingBox3f bbox = aggregateBBoxes(bboxes);
    mServerQuerier->updateRegion(bbox);

    mContext->add(&mServerHandlerPoller);
    mContext->add(&mObjectHandlerPoller);
    mContext->add(&mStaticRebuilderPoller);
    mContext->add(&mDynamicRebuilderPoller);
}

void LibproxProximity::newSession(ObjectSession* session) {
    using std::tr1::placeholders::_1;
    using std::tr1::placeholders::_2;

    SST::Stream<SpaceObjectReference>::Ptr strm = session->getStream();
    SST::Connection<SpaceObjectReference>::Ptr conn = strm->connection().lock();
    assert(conn);

    SpaceObjectReference sourceObject = conn->remoteEndPoint().endPoint;

    conn->registerReadDatagramCallback( OBJECT_PORT_PROXIMITY,
        std::tr1::bind(
            &LibproxProximity::handleObjectProximityMessage, this,
            sourceObject.object().getAsUUID(),
            _1, _2
        )
    );
}

void LibproxProximity::sessionClosed(ObjectSession* session) {
    // Prox strand may  have some state to clean up
    mProxStrand->post(
        std::tr1::bind(&LibproxProximity::handleDisconnectedObject, this, session->id().getAsUUID())
    );

    ObjectProxStreamMap::iterator prox_stream_it = mObjectProxStreams.find(session->id().getAsUUID());
    if (prox_stream_it != mObjectProxStreams.end()) {
        prox_stream_it->second->disable();
        mObjectProxStreams.erase(prox_stream_it);
    }
}

void LibproxProximity::handleObjectProximityMessage(const UUID& objid, void* buffer, uint32 length) {
    Sirikata::Protocol::Prox::QueryRequest prox_update;
    bool parse_success = prox_update.ParseFromString( String((char*) buffer, length) );
    if (!parse_success) {
        LOG_INVALID_MESSAGE_BUFFER(prox, error, ((char*)buffer), length);
        return;
    }

    // Need to support old style (angle + max count) and new style.
    if (prox_update.has_query_parameters()) {
        SolidAngle sa;
        uint32 max_results;
        if (parseQueryRequest(prox_update.query_parameters(), &sa, &max_results))
            updateQuery(objid, mLocService->location(objid), mLocService->bounds(objid), sa, max_results);
    }
    else {
        if (!prox_update.has_query_angle()) return;

        SolidAngle query_angle(prox_update.query_angle());
        uint32 query_max_results = (prox_update.has_query_max_count() ? prox_update.query_max_count() : NoUpdateMaxResults);
        updateQuery(objid, mLocService->location(objid), mLocService->bounds(objid), query_angle, query_max_results);
    }
}

// Setup all known servers for a server query update
void LibproxProximity::addAllServersForUpdate() {
    boost::lock_guard<boost::mutex> lck(mServerSetMutex);
    for(ServerSet::const_iterator it = mServersQueried.begin(); it != mServersQueried.end(); it++)
        mNeedServerQueryUpdate.insert(*it);
}

void LibproxProximity::sendQueryRequests() {
    TimedMotionVector3f loc;

    // FIXME avoid computing this so much
    BoundingBoxList bboxes = mCSeg->serverRegion(mContext->id());
    BoundingBox3f bbox = aggregateBBoxes(bboxes);
    BoundingSphere3f bounds = bbox.toBoundingSphere();

    ServerSet sub_servers;
    {
        boost::lock_guard<boost::mutex> lck(mServerSetMutex);
        sub_servers.swap(mNeedServerQueryUpdate);
    }
    for(ServerSet::const_iterator it = sub_servers.begin(); it != sub_servers.end(); it++) {
        ServerID sid = *it;
        Sirikata::Protocol::Prox::Container container;
        Sirikata::Protocol::Prox::IServerQuery msg = container.mutable_query();
        msg.set_action(Sirikata::Protocol::Prox::ServerQuery::AddOrUpdate);
        Sirikata::Protocol::ITimedMotionVector msg_loc = msg.mutable_location();
        msg_loc.set_t(loc.updateTime());
        msg_loc.set_position(loc.position());
        msg_loc.set_position(loc.velocity());
        msg.set_bounds(bounds);
        msg.set_min_angle(mMinObjectQueryAngle.asFloat());
        msg.set_max_count(mMaxMaxCount);

        Message* smsg = new Message(
            mContext->id(),
            SERVER_PORT_PROX,
            sid,
            SERVER_PORT_PROX,
            serializePBJMessage(container)
        );
        bool sent = mProxServerMessageService->route(smsg);
        if (!sent) {
            delete smsg;
            {
                boost::lock_guard<boost::mutex> lck(mServerSetMutex);
                mNeedServerQueryUpdate.insert(sid);
            }
        }
    }
}

void LibproxProximity::receiveMessage(Message* msg) {
    assert(msg->dest_port() == SERVER_PORT_PROX);

    Sirikata::Protocol::Prox::Container prox_container;
    bool parsed = parsePBJMessage(&prox_container, msg->payload());
    if (!parsed) {
        PROXLOG(warn,"Couldn't parse message, ID=" << msg->id());
        delete msg;
        return;
    }

    ServerID source_server = msg->source_server();

    if (prox_container.has_query()) {
        Sirikata::Protocol::Prox::ServerQuery prox_query_msg = prox_container.query();

        if (prox_query_msg.action() == Sirikata::Protocol::Prox::ServerQuery::AddOrUpdate) {
            assert(
                prox_query_msg.has_location() &&
                prox_query_msg.has_bounds() &&
                prox_query_msg.has_min_angle()
            );

            Sirikata::Protocol::TimedMotionVector msg_loc = prox_query_msg.location();
            TimedMotionVector3f qloc(msg_loc.t(), MotionVector3f(msg_loc.position(), msg_loc.velocity()));
            SolidAngle minangle(prox_query_msg.min_angle());
            uint32 query_max_results = (prox_query_msg.has_max_count() ? prox_query_msg.max_count() : NoUpdateMaxResults);

            updateQuery(source_server, qloc, prox_query_msg.bounds(), minangle, query_max_results);
        }
        else if (prox_query_msg.action() == Sirikata::Protocol::Prox::ServerQuery::Remove) {
            removeQuery(source_server);
        }
        else {
            assert(false);
        }
    }

    if (prox_container.has_result()) {
        if (!mServerQueryResults[source_server])
            mServerQueryResults[source_server] = ObjectSetPtr(new ObjectSet());

        Sirikata::Protocol::Prox::ProximityResults prox_result_msg = prox_container.result();

        assert( prox_result_msg.has_t() );
        Time t = prox_result_msg.t();

        for(int32 idx = 0; idx < prox_result_msg.update_size(); idx++) {
            Sirikata::Protocol::Prox::ProximityUpdate update = prox_result_msg.update(idx);

            for(int32 aidx = 0; aidx < update.addition_size(); aidx++) {
                Sirikata::Protocol::Prox::ObjectAddition addition = update.addition(aidx);
                mServerQueryResults[source_server]->insert(addition.object());
                mLocService->addReplicaObject(
                    t,
                    addition.object(),
                    TimedMotionVector3f( addition.location().t(), MotionVector3f(addition.location().position(), addition.location().velocity()) ),
                    TimedMotionQuaternion( addition.orientation().t(), MotionQuaternion(addition.orientation().position(), addition.orientation().velocity()) ),
                    addition.bounds(),
                    (addition.has_mesh() ? addition.mesh() : ""),
                    (addition.has_physics() ? addition.physics() : "")
                );
            }

            for(int32 ridx = 0; ridx < update.removal_size(); ridx++) {
                Sirikata::Protocol::Prox::ObjectRemoval removal = update.removal(ridx);
                mLocService->removeReplicaObject(t, removal.object());
                mServerQueryResults[source_server]->erase(removal.object());
            }
        }
    }

    delete msg;
}

// MigrationDataClient Interface

std::string LibproxProximity::migrationClientTag() {
    return "prox";
}

std::string LibproxProximity::generateMigrationData(const UUID& obj, ServerID source_server, ServerID dest_server) {
    ObjectQueryAngleMap::iterator it = mObjectQueryAngles.find(obj);
    if (it == mObjectQueryAngles.end()) // no query registered, return nothing
        return std::string();
    else {
        SolidAngle query_angle = it->second;
        removeQuery(obj);

        Sirikata::Protocol::Prox::ObjectMigrationData migr_data;
        migr_data.set_min_angle( query_angle.asFloat() );
        if (mObjectQueryMaxCounts.find(obj) != mObjectQueryMaxCounts.end())
            migr_data.set_max_count( mObjectQueryMaxCounts[obj] );
        return serializePBJMessage(migr_data);
    }
}

void LibproxProximity::receiveMigrationData(const UUID& obj, ServerID source_server, ServerID dest_server, const std::string& data) {
    Sirikata::Protocol::Prox::ObjectMigrationData migr_data;
    bool parse_success = migr_data.ParseFromString(data);
    if (!parse_success) {
        LOG_INVALID_MESSAGE(prox, error, data);
        return;
    }

    SolidAngle obj_query_angle(migr_data.min_angle());
    uint32 obj_query_max_results = (migr_data.has_max_count() ? migr_data.max_count() : ObjectProxSimulationTraits::InfiniteResults);
    addQuery(obj, obj_query_angle, obj_query_max_results);
}

// PintoServerQuerierListener Interface

void LibproxProximity::addRelevantServer(ServerID sid) {
    if (sid == mContext->id()) return;

    // Potentially invoked from PintoServerQuerier IO thread
    boost::lock_guard<boost::mutex> lck(mServerSetMutex);
    mServersQueried.insert(sid);
    mNeedServerQueryUpdate.insert(sid);
}

void LibproxProximity::removeRelevantServer(ServerID sid) {
    if (sid == mContext->id()) return;

    // Potentially invoked from PintoServerQuerier IO thread
    boost::lock_guard<boost::mutex> lck(mServerSetMutex);
    mServersQueried.erase(sid);
}

void LibproxProximity::aggregateCreated(ProxAggregator* handler, const UUID& objid) {
    LibproxProximityBase::aggregateCreated(objid);
}

void LibproxProximity::aggregateChildAdded(ProxAggregator* handler, const UUID& objid, const UUID& child, const BoundingSphere3f& bnds) {
    LibproxProximityBase::aggregateChildAdded(objid, child, bnds);
}

void LibproxProximity::aggregateChildRemoved(ProxAggregator* handler, const UUID& objid, const UUID& child, const BoundingSphere3f& bnds) {
    LibproxProximityBase::aggregateChildRemoved(objid, child, bnds);
}

void LibproxProximity::aggregateBoundsUpdated(ProxAggregator* handler, const UUID& objid, const BoundingSphere3f& bnds) {
    LibproxProximityBase::aggregateBoundsUpdated(objid, bnds);
}

void LibproxProximity::aggregateDestroyed(ProxAggregator* handler, const UUID& objid) {
    LibproxProximityBase::aggregateDestroyed(objid);
}

void LibproxProximity::aggregateObserved(ProxAggregator* handler, const UUID& objid, uint32 nobservers) {
    LibproxProximityBase::aggregateObserved(objid, nobservers);
}

void LibproxProximity::onSpaceNetworkConnected(ServerID sid) {
    mProxStrand->post(
        std::tr1::bind(&LibproxProximity::handleConnectedServer, this, sid)
    );
}

void LibproxProximity::onSpaceNetworkDisconnected(ServerID sid) {
    mProxStrand->post(
        std::tr1::bind(&LibproxProximity::handleDisconnectedServer, this, sid)
    );
}

void LibproxProximity::updateQuery(ServerID sid, const TimedMotionVector3f& loc, const BoundingSphere3f& bounds, const SolidAngle& sa, uint32 max_results) {
    mProxStrand->post(
        std::tr1::bind(&LibproxProximity::handleUpdateServerQuery, this, sid, loc, bounds, sa, max_results)
    );
}

void LibproxProximity::removeQuery(ServerID sid) {
    mProxStrand->post(
        std::tr1::bind(&LibproxProximity::handleRemoveServerQuery, this, sid)
    );
}

void LibproxProximity::addQuery(UUID obj, SolidAngle sa, uint32 max_results) {
    updateQuery(obj, mLocService->location(obj), mLocService->bounds(obj), sa, max_results);
}

void LibproxProximity::addQuery(UUID obj, const String& params) {
    SolidAngle sa;
    uint32 max_results;
    if (parseQueryRequest(params, &sa, &max_results))
        updateQuery(obj, mLocService->location(obj), mLocService->bounds(obj), sa, max_results);
}

void LibproxProximity::updateQuery(UUID obj, const TimedMotionVector3f& loc, const BoundingSphere3f& bounds, SolidAngle sa, uint32 max_results) {
    SeqNoPtr obj_seqno = mContext->objectSessionManager()->getSession(ObjectReference(obj))->getSeqNoPtr();

    // Update the prox thread
    mProxStrand->post(
        std::tr1::bind(&LibproxProximity::handleUpdateObjectQuery, this, obj, loc, bounds, sa, max_results, obj_seqno)
    );

    bool update_remote_queries = false;
    if (sa != NoUpdateSolidAngle) {
        // Update the main thread's record
        mObjectQueryAngles[obj] = sa;

        // Update min query angle, and update remote queries if necessary
        if (sa < mMinObjectQueryAngle) {
            mMinObjectQueryAngle = sa;
            update_remote_queries = true;
        }
    }

    if (max_results != NoUpdateMaxResults) {
        // Update the main thread's record
        mObjectQueryMaxCounts[obj] = max_results;

        // Update min query angle, and update remote queries if necessary
        if (max_results > mMaxMaxCount) {
            mMaxMaxCount = max_results;
            update_remote_queries = true;
        }
    }

    if (update_remote_queries) {
        PROXLOG(debug,"Query addition initiated server query request.");
        addAllServersForUpdate();
        mServerQuerier->updateQuery(mMinObjectQueryAngle, mMaxMaxCount);
    }
}

void LibproxProximity::removeQuery(UUID obj) {
    // Update the main thread's record
    SolidAngle sa = mObjectQueryAngles[obj];
    mObjectQueryAngles.erase(obj);
    uint32 max_count = mObjectQueryMaxCounts[obj];
    mObjectQueryMaxCounts.erase(obj);

    // Update the prox thread
    mProxStrand->post(
        std::tr1::bind(&LibproxProximity::handleRemoveObjectQuery, this, obj, true)
    );

    // Update min query angle, and update remote queries if necessary
    if (sa == mMinObjectQueryAngle || max_count == mMaxMaxCount) {
        PROXLOG(debug,"Query removal initiated server query request.");
        SolidAngle minangle(SolidAngle::Max);
        for(ObjectQueryAngleMap::iterator it = mObjectQueryAngles.begin(); it != mObjectQueryAngles.end(); it++)
            if (it->second < minangle) minangle = it->second;
        uint32 maxcount = 1;
        for(ObjectQueryMaxCountMap::iterator it = mObjectQueryMaxCounts.begin(); it != mObjectQueryMaxCounts.end(); it++)
            if (it->second == ObjectProxSimulationTraits::InfiniteResults || (maxcount > ObjectProxSimulationTraits::InfiniteResults && it->second > maxcount))
                maxcount = it->second;

        // NOTE: Even if this condition is satisfied, we could only be increasing
        // the minimum angle, so we don't *strictly* need to update the query.
        // Some buffer timing might be in order here to avoid excessive updates
        // while still getting the benefit from reducing the query angle.
        if (minangle != mMinObjectQueryAngle || maxcount != mMaxMaxCount) {
            mMinObjectQueryAngle = minangle;
            mMaxMaxCount = maxcount;
            addAllServersForUpdate();
            mServerQuerier->updateQuery(mMinObjectQueryAngle, mMaxMaxCount);
        }
    }
}

void LibproxProximity::updateObjectSize(const UUID& obj, float rad) {
    mObjectSizes[obj] = rad;

    if (rad > mMaxObject) {
        mMaxObject = rad;
        mServerQuerier->updateLargestObject(mMaxObject);
    }
}

void LibproxProximity::removeObjectSize(const UUID& obj) {
    ObjectSizeMap::iterator it = mObjectSizes.find(obj);
    assert(it != mObjectSizes.end());

    float32 old_val = it->second;
    mObjectSizes.erase(it);

    if (old_val == mMaxObject) {
        // We need to recompute mMaxObject since it may not be accurate anymore
        mMaxObject = 0.0f;
        for(ObjectSizeMap::iterator oit = mObjectSizes.begin(); oit != mObjectSizes.end(); oit++)
            mMaxObject = std::max(mMaxObject, oit->second);

        if (mMaxObject != old_val)
            mServerQuerier->updateLargestObject(mMaxObject);
    }
}

void LibproxProximity::checkObjectClass(bool is_local, const UUID& objid, const TimedMotionVector3f& newval) {
    mProxStrand->post(
        std::tr1::bind(&LibproxProximity::handleCheckObjectClass, this, is_local, objid, newval)
    );
}

int32 LibproxProximity::objectQueries() const {
    return mObjectQueries[OBJECT_CLASS_STATIC].size();
}

int32 LibproxProximity::serverQueries() const {
    return mServerQueries[OBJECT_CLASS_STATIC].size();
}

void LibproxProximity::poll() {
    // Update server-to-server angles if necessary
    sendQueryRequests();

    // Get and ship server results
    std::deque<Message*> server_results_copy;
    mServerResults.swap(server_results_copy);
    mServerResultsToSend.insert(mServerResultsToSend.end(), server_results_copy.begin(), server_results_copy.end());

    bool server_sent = true;
    while(server_sent && !mServerResultsToSend.empty()) {
        Message* msg_front = mServerResultsToSend.front();
        server_sent = mProxServerMessageService->route(msg_front);
        if (server_sent)
            mServerResultsToSend.pop_front();
    }

    // Get and ship object results
    std::deque<Sirikata::Protocol::Object::ObjectMessage*> object_results_copy;
    mObjectResults.swap(object_results_copy);
    mObjectResultsToSend.insert(mObjectResultsToSend.end(), object_results_copy.begin(), object_results_copy.end());

    bool object_sent = true;
    while(object_sent && !mObjectResultsToSend.empty()) {
        Sirikata::Protocol::Object::ObjectMessage* msg_front = mObjectResultsToSend.front();
        sendObjectResult(msg_front);
        delete msg_front;
        mObjectResultsToSend.pop_front();
    }
}

void LibproxProximity::queryHasEvents(Query* query) {
    if (
        query->handler() == mServerQueryHandler[OBJECT_CLASS_STATIC] ||
        query->handler() == mServerQueryHandler[OBJECT_CLASS_DYNAMIC]
    )
        generateServerQueryEvents(query);
    else
        generateObjectQueryEvents(query);
}

// Note: LocationServiceListener interface is only used in order to get updates on objects which have
// registered queries, allowing us to update those queries as appropriate.  All updating of objects
// in the prox data structure happens via the LocationServiceCache
void LibproxProximity::localObjectAdded(const UUID& uuid, bool agg, const TimedMotionVector3f& loc, const TimedMotionQuaternion& orient, const BoundingSphere3f& bounds, const String& mesh, const String& physics) {
    updateObjectSize(uuid, bounds.radius());
}
void LibproxProximity::localObjectRemoved(const UUID& uuid, bool agg) {
    removeObjectSize(uuid);
}
void LibproxProximity::localLocationUpdated(const UUID& uuid, bool agg, const TimedMotionVector3f& newval) {
    updateQuery(uuid, newval, mLocService->bounds(uuid), NoUpdateSolidAngle, NoUpdateMaxResults);
    if (mSeparateDynamicObjects)
        checkObjectClass(true, uuid, newval);
}
void LibproxProximity::localBoundsUpdated(const UUID& uuid, bool agg, const BoundingSphere3f& newval) {
    updateQuery(uuid, mLocService->location(uuid), newval, NoUpdateSolidAngle, NoUpdateMaxResults);
    updateObjectSize(uuid, newval.radius());
}
void LibproxProximity::replicaLocationUpdated(const UUID& uuid, const TimedMotionVector3f& newval) {
    if (mSeparateDynamicObjects)
        checkObjectClass(false, uuid, newval);
}

void LibproxProximity::updatedSegmentation(CoordinateSegmentation* cseg, const std::vector<SegmentationInfo>& new_seg) {
    BoundingBoxList bboxes = mCSeg->serverRegion(mContext->id());
    BoundingBox3f bbox = aggregateBBoxes(bboxes);
    mServerQuerier->updateRegion(bbox);
}

// PROX Thread: Everything after this should only be called from within the prox thread.

void LibproxProximity::tickQueryHandler(ProxQueryHandler* qh[NUM_OBJECT_CLASSES]) {
    Time simT = mContext->simTime();
    for(int i = 0; i < NUM_OBJECT_CLASSES; i++) {
        if (qh[i] != NULL)
            qh[i]->tick(simT);
    }
}

void LibproxProximity::rebuildHandler(ObjectClass objtype) {
    if (mServerQueryHandler[objtype] != NULL)
        mServerQueryHandler[objtype]->rebuild();
    if (mObjectQueryHandler[objtype] != NULL)
        mObjectQueryHandler[objtype]->rebuild();
}

void LibproxProximity::generateServerQueryEvents(Query* query) {
    typedef std::deque<QueryEvent> QueryEventList;

    Time t = mContext->simTime();
    uint32 max_count = GetOptionValue<uint32>(PROX_MAX_PER_RESULT);

    assert(mInvertedServerQueries.find(query) != mInvertedServerQueries.end());
    ServerID sid = mInvertedServerQueries[query];
    SeqNoPtr seqNoPtr = getOrCreateSeqNoInfo(sid);

    QueryEventList evts;
    query->popEvents(evts);

    while(!evts.empty()) {
        Sirikata::Protocol::Prox::Container container;
        Sirikata::Protocol::Prox::IProximityResults contents = container.mutable_result();
        contents.set_t(t);
        uint32 count = 0;
        while(count < max_count && !evts.empty()) {
            const QueryEvent& evt = evts.front();
            Sirikata::Protocol::Prox::IProximityUpdate event_results = contents.add_update();
            // Each QueryEvent is made up of additions and
            // removals
            for(uint32 aidx = 0; aidx < evt.additions().size(); aidx++) {
                UUID objid = evt.additions()[aidx].id();
                if (mLocCache->tracking(objid)) { // If the cache already lost it, we can't do anything
                    count++;

                    mContext->mainStrand->post(
                        std::tr1::bind(&LibproxProximity::handleAddServerLocSubscription, this, sid, objid, seqNoPtr)
                    );

                    Sirikata::Protocol::Prox::IObjectAddition addition = event_results.add_addition();
                    addition.set_object( objid );

                    uint64 seqNo = (*seqNoPtr)++;
                    addition.set_seqno (seqNo);

                    TimedMotionVector3f loc = mLocCache->location(objid);
                    Sirikata::Protocol::ITimedMotionVector msg_loc = addition.mutable_location();
                    msg_loc.set_t(loc.updateTime());
                    msg_loc.set_position(loc.position());
                    msg_loc.set_velocity(loc.velocity());

                    TimedMotionQuaternion orient = mLocCache->orientation(objid);
                    Sirikata::Protocol::ITimedMotionQuaternion msg_orient = addition.mutable_orientation();
                    msg_orient.set_t(orient.updateTime());
                    msg_orient.set_position(orient.position());
                    msg_orient.set_velocity(orient.velocity());

                    addition.set_bounds( mLocCache->bounds(objid) );
                    const String& mesh = mLocCache->mesh(objid);
                    if (mesh.size() > 0)
                        addition.set_mesh(mesh);
                    const String& phy = mLocCache->physics(objid);
                    if (phy.size() > 0)
                        addition.set_physics(phy);
                }
            }
            for(uint32 ridx = 0; ridx < evt.removals().size(); ridx++) {
                UUID objid = evt.removals()[ridx].id();
                count++;

                mContext->mainStrand->post(
                    std::tr1::bind(&LibproxProximity::handleRemoveServerLocSubscription, this, sid, objid)
                );

                Sirikata::Protocol::Prox::IObjectRemoval removal = event_results.add_removal();
                removal.set_object(objid);
                uint64 seqNo = (*seqNoPtr)++;
                removal.set_seqno (seqNo);
                removal.set_type(
                    (evt.removals()[ridx].permanent() == QueryEvent::Permanent)
                    ? Sirikata::Protocol::Prox::ObjectRemoval::Permanent
                    : Sirikata::Protocol::Prox::ObjectRemoval::Transient
                );
            }

            evts.pop_front();
        }

        //PROXLOG(insane,"Reporting " << contents.addition_size() << " additions, " << contents.removal_size() << " removals to server " << sid);

        Message* msg = new Message(
            mContext->id(),
            SERVER_PORT_PROX,
            sid,
            SERVER_PORT_PROX,
            serializePBJMessage(container)
        );
        mServerResults.push(msg);
    }
}

void LibproxProximity::generateObjectQueryEvents(Query* query) {
    typedef std::deque<QueryEvent> QueryEventList;

    uint32 max_count = GetOptionValue<uint32>(PROX_MAX_PER_RESULT);

    assert(mInvertedObjectQueries.find(query) != mInvertedObjectQueries.end());
    UUID query_id = mInvertedObjectQueries[query];
    SeqNoPtr seqNoPtr = getSeqNoInfo(query_id);

    QueryEventList evts;
    query->popEvents(evts);

    while(!evts.empty()) {
        Sirikata::Protocol::Prox::ProximityResults prox_results;
        prox_results.set_t(mContext->simTime());

        uint32 count = 0;
        while(count < max_count && !evts.empty()) {
            const QueryEvent& evt = evts.front();
            Sirikata::Protocol::Prox::IProximityUpdate event_results = prox_results.add_update();

            for(uint32 aidx = 0; aidx < evt.additions().size(); aidx++) {
                UUID objid = evt.additions()[aidx].id();
                if (mLocCache->tracking(objid)) { // If the cache already lost it, we can't do anything
                    count++;

                    mContext->mainStrand->post(
                        std::tr1::bind(&LibproxProximity::handleAddObjectLocSubscription, this, query_id, objid)
                    );

                    Sirikata::Protocol::Prox::IObjectAddition addition = event_results.add_addition();
                    addition.set_object( objid );

                    //query_id contains the uuid of the object that is receiving
                    //the proximity message that obj_id has been added.
                    uint64 seqNo = (*seqNoPtr);
                    addition.set_seqno (seqNo);

                    Sirikata::Protocol::ITimedMotionVector motion = addition.mutable_location();
                    TimedMotionVector3f loc = mLocCache->location(objid);
                    motion.set_t(loc.updateTime());
                    motion.set_position(loc.position());
                    motion.set_velocity(loc.velocity());

                    TimedMotionQuaternion orient = mLocCache->orientation(objid);
                    Sirikata::Protocol::ITimedMotionQuaternion msg_orient = addition.mutable_orientation();
                    msg_orient.set_t(orient.updateTime());
                    msg_orient.set_position(orient.position());
                    msg_orient.set_velocity(orient.velocity());

                    addition.set_bounds( mLocCache->bounds(objid) );
                    const String& mesh = mLocCache->mesh(objid);
                    if (mesh.size() > 0)
                        addition.set_mesh(mesh);
                    const String& phy = mLocCache->physics(objid);
                    if (phy.size() > 0)
                        addition.set_physics(phy);
                }
            }
            for(uint32 ridx = 0; ridx < evt.removals().size(); ridx++) {
                UUID objid = evt.removals()[ridx].id();
                count++;
                // Clear out seqno and let main strand remove loc
                // subcription
                mContext->mainStrand->post(
                    std::tr1::bind(&LibproxProximity::handleRemoveObjectLocSubscription, this, query_id, objid)
                );

                Sirikata::Protocol::Prox::IObjectRemoval removal = event_results.add_removal();
                removal.set_object( objid );
                uint64 seqNo = (*seqNoPtr)++;
                removal.set_seqno (seqNo);
                removal.set_type(
                    (evt.removals()[ridx].permanent() == QueryEvent::Permanent)
                    ? Sirikata::Protocol::Prox::ObjectRemoval::Permanent
                    : Sirikata::Protocol::Prox::ObjectRemoval::Transient
                );
            }
            evts.pop_front();
        }

        Sirikata::Protocol::Object::ObjectMessage* obj_msg = createObjectMessage(
            mContext->id(),
            UUID::null(), OBJECT_PORT_PROXIMITY,
            query_id, OBJECT_PORT_PROXIMITY,
            serializePBJMessage(prox_results)
        );
        mObjectResults.push(obj_msg);
    }
}

SeqNoPtr LibproxProximity::getOrCreateSeqNoInfo(const ServerID server_id)
{
    // server_id == querier
    ServerSeqNoInfoMap::iterator proxSeqNoIt = mServerSeqNos.find(server_id);
    if (proxSeqNoIt == mServerSeqNos.end())
        proxSeqNoIt = mServerSeqNos.insert( ServerSeqNoInfoMap::value_type(server_id, SeqNoPtr(new SeqNo())) ).first;
    return proxSeqNoIt->second;
}

void LibproxProximity::eraseSeqNoInfo(const ServerID server_id)
{
    // server_id == querier
    ServerSeqNoInfoMap::iterator proxSeqNoIt = mServerSeqNos.find(server_id);
    if (proxSeqNoIt == mServerSeqNos.end()) return;
    mServerSeqNos.erase(proxSeqNoIt);
}

SeqNoPtr LibproxProximity::getSeqNoInfo(const UUID& obj_id)
{
    // obj_id == querier
    ObjectSeqNoInfoMap::iterator proxSeqNoIt = mObjectSeqNos.find(obj_id);
    assert(proxSeqNoIt != mObjectSeqNos.end());
    return proxSeqNoIt->second;
}

void LibproxProximity::eraseSeqNoInfo(const UUID& obj_id)
{
    // obj_id == querier
    ObjectSeqNoInfoMap::iterator proxSeqNoIt = mObjectSeqNos.find(obj_id);
    if (proxSeqNoIt == mObjectSeqNos.end()) return;
    mObjectSeqNos.erase(proxSeqNoIt);
}

void LibproxProximity::handleUpdateServerQuery(const ServerID& server, const TimedMotionVector3f& loc, const BoundingSphere3f& bounds, const SolidAngle& angle, const uint32 max_results) {
    BoundingSphere3f region(bounds.center(), 0);
    float ms = bounds.radius();

    for(int i = 0; i < NUM_OBJECT_CLASSES; i++) {
        if (mServerQueryHandler[i] == NULL) continue;

        ServerQueryMap::iterator it = mServerQueries[i].find(server);
        if (it == mServerQueries[i].end()) {
            PROXLOG(debug,"Add server query from " << server << ", min angle " << angle.asFloat() << ", object class " << ObjectClassToString((ObjectClass)i));

            Query* q = mServerDistance ?
                mServerQueryHandler[i]->registerQuery(loc, region, ms, SolidAngle::Min, mDistanceQueryDistance) :
                mServerQueryHandler[i]->registerQuery(loc, region, ms, angle) ;
            if (max_results != NoUpdateMaxResults && max_results > 0)
                q->maxResults(max_results);
            mServerQueries[i][server] = q;
            mInvertedServerQueries[q] = server;
            q->setEventListener(this);
        }
        else {
            PROXLOG(debug,"Update server query from " << server << ", min angle " << angle.asFloat() << ", object class " << ObjectClassToString((ObjectClass)i));

            Query* q = it->second;
            q->position(loc);
            q->region( region );
            q->maxSize( ms );
            q->angle(angle);
            if (max_results != NoUpdateMaxResults && max_results > 0)
                q->maxResults(max_results);
        }
    }
}

void LibproxProximity::handleRemoveServerQuery(const ServerID& server) {
    PROXLOG(debug,"Remove server query from " << server);

    for(int i = 0; i < NUM_OBJECT_CLASSES; i++) {
        if (mServerQueryHandler[i] == NULL) continue;

        ServerQueryMap::iterator it = mServerQueries[i].find(server);
        if (it == mServerQueries[i].end()) continue;

        Query* q = it->second;
        mServerQueries[i].erase(it);
        mInvertedServerQueries.erase(q);
        delete q; // Note: Deleting query notifies QueryHandler and unsubscribes.
    }

    // Clear out sequence numbers
    eraseSeqNoInfo(server);

    mContext->mainStrand->post(
        std::tr1::bind(&LibproxProximity::handleRemoveAllServerLocSubscription, this, server)
    );
}

void LibproxProximity::handleConnectedServer(ServerID sid) {
    // Sometimes we may get forcefully disconnected from a server. To
    // reestablish our previous setup, if that server appears in our queried
    // set (we were told it was relevant to us by some higher level pinto
    // service), we mark it as needing another update. In the case that we
    // are just getting an initial connection, this shouldn't change anything
    // since it would already be markedas needing an update if we had been told
    // by pinto that it needed it.
    boost::lock_guard<boost::mutex> lck(mServerSetMutex);
    if (mServersQueried.find(sid) != mServersQueried.end())
        mNeedServerQueryUpdate.insert(sid);
}

void LibproxProximity::handleDisconnectedServer(ServerID sid) {
    // When we lose a connection, we need to clear out everything related to
    // that server.
    PROXLOG(debug, "Handling unexpected disconnection from " << sid << " by clearing out proximity state.");

    // Clear out remote server's query to us
    handleRemoveServerQuery(sid);

    // And our query (and associated state) to them
    // They clear out our query (as we do for them above), but we need to also
    // remove the state we maintain locally, i.e. replicated objects.
    ObjectSetPtr objects = mServerQueryResults[sid];
    if (!objects) return;
    mServerQueryResults.erase(sid);
    Time t = mContext->simTime();
    for(ObjectSet::const_iterator it = objects->begin(); it != objects->end(); it++)
        mLocService->removeReplicaObject(t, *it);
}

void LibproxProximity::handleUpdateObjectQuery(const UUID& object, const TimedMotionVector3f& loc, const BoundingSphere3f& bounds, const SolidAngle& angle, uint32 max_results, SeqNoPtr seqno) {
    BoundingSphere3f region(bounds.center(), 0);
    float ms = bounds.radius();

    PROXLOG(detailed,"Update object query from " << object.toString() << ", min angle " << angle.asFloat() << ", max results " << max_results);

    if (mObjectSeqNos.find(object) == mObjectSeqNos.end())
        mObjectSeqNos.insert( ObjectSeqNoInfoMap::value_type(object, seqno) );

    for(int i = 0; i < NUM_OBJECT_CLASSES; i++) {
        if (mObjectQueryHandler[i] == NULL) continue;

        ObjectQueryMap::iterator it = mObjectQueries[i].find(object);

        if (it == mObjectQueries[i].end()) {
            // We only add if we actually have all the necessary info, most importantly a real minimum angle.
            // This is necessary because we get this update for all location updates, even those for objects
            // which don't have subscriptions.
            if (angle != NoUpdateSolidAngle) {
                Query* q = mObjectDistance ?
                    mObjectQueryHandler[i]->registerQuery(loc, region, ms, SolidAngle::Min, mDistanceQueryDistance) :
                    mObjectQueryHandler[i]->registerQuery(loc, region, ms, angle);
                if (max_results != NoUpdateMaxResults && max_results > 0)
                    q->maxResults(max_results);
                mObjectQueries[i][object] = q;
                mInvertedObjectQueries[q] = object;
                q->setEventListener(this);
            }
        }
        else {
            Query* query = it->second;
            query->position(loc);
            query->region( region );
            query->maxSize( ms );
            if (angle != NoUpdateSolidAngle)
                query->angle(angle);
            if (max_results != NoUpdateMaxResults && max_results > 0)
                query->maxResults(max_results);
        }
    }
}

void LibproxProximity::handleRemoveObjectQuery(const UUID& object, bool notify_main_thread) {
    // Clear out queries
    for(int i = 0; i < NUM_OBJECT_CLASSES; i++) {
        if (mObjectQueryHandler[i] == NULL) continue;

        ObjectQueryMap::iterator it = mObjectQueries[i].find(object);
        if (it == mObjectQueries[i].end()) continue;

        Query* q = it->second;
        mObjectQueries[i].erase(it);
        mInvertedObjectQueries.erase(q);
        delete q; // Note: Deleting query notifies QueryHandler and unsubscribes.
    }

    // Clear out sequence numbers
    eraseSeqNoInfo(object);

    // Optionally let the main thread know to clear its communication state
    if (notify_main_thread) {
        // There's no corresponding removeAllSeqNoPtr because we
        // should have erased it above.
        mContext->mainStrand->post(
            std::tr1::bind(&LibproxProximity::handleRemoveAllObjectLocSubscription, this, object)
        );
    }
}

void LibproxProximity::handleDisconnectedObject(const UUID& object) {
    // Clear out query state if it exists
    handleRemoveObjectQuery(object, false);
}

bool LibproxProximity::handlerShouldHandleObject(bool is_static_handler, bool is_global_handler, const UUID& obj_id, bool is_local, const TimedMotionVector3f& pos, const BoundingSphere3f& region, float maxSize) {
    // We just need to decide whether the query handler should handle
    // the object. We need to consider local vs. replica and static
    // vs. dynamic.  All must 'vote' for handling the object for us to
    // say it should be handled, so as soon as we find a negative
    // response we can return false.

    // First classify by local vs. replica. Only say no on a local
    // handler looking at a replica.
    if (!is_local && !is_global_handler) return false;

    // If we're not doing the static/dynamic split, then this is a non-issue
    if (!mSeparateDynamicObjects) return true;

    // If we are splitting them, check velocity against is_static_handler. The
    // value here as arbitrary, just meant to indicate such small movement that
    // the object is effectively
    bool is_static = velocityIsStatic(pos.velocity());
    if ((is_static && is_static_handler) ||
        (!is_static && !is_static_handler))
        return true;
    else
        return false;
}

void LibproxProximity::handleCheckObjectClassForHandlers(const UUID& objid, bool is_static, ProxQueryHandler* handlers[NUM_OBJECT_CLASSES]) {
    if ( (is_static && handlers[OBJECT_CLASS_STATIC]->containsObject(objid)) ||
        (!is_static && handlers[OBJECT_CLASS_DYNAMIC]->containsObject(objid)) )
        return;

    // Validate that the other handler has the object.
    assert(
        (is_static && handlers[OBJECT_CLASS_DYNAMIC]->containsObject(objid)) ||
        (!is_static && handlers[OBJECT_CLASS_STATIC]->containsObject(objid))
    );

    // If it wasn't in the right place, switch it.
    int swap_out = is_static ? OBJECT_CLASS_DYNAMIC : OBJECT_CLASS_STATIC;
    int swap_in = is_static ? OBJECT_CLASS_STATIC : OBJECT_CLASS_DYNAMIC;
    PROXLOG(debug, "Swapping " << objid.toString() << " from " << ObjectClassToString((ObjectClass)swap_out) << " to " << ObjectClassToString((ObjectClass)swap_in));
    handlers[swap_out]->removeObject(objid);
    handlers[swap_in]->addObject(objid);
}

void LibproxProximity::handleCheckObjectClass(bool is_local, const UUID& objid, const TimedMotionVector3f& newval) {
    assert(mSeparateDynamicObjects == true);

    // Basic approach: we need to check if the object has switched between
    // static/dynamic. We need to do this for both the local (object query) and
    // global (server query) handlers.
    bool is_static = velocityIsStatic(newval.velocity());
    handleCheckObjectClassForHandlers(objid, is_static, mObjectQueryHandler);
    if (is_local)
        handleCheckObjectClassForHandlers(objid, is_static, mServerQueryHandler);
}

} // namespace Sirikata

max_stars_repo_name: pathorn/sirikata

max_issues_repo_name: pathorn/sirikata

max_forks_repo_name: pathorn/sirikata


text:
// Distributed under the MIT License.
// See LICENSE.txt for details.

#pragma once

#include <boost/preprocessor/arithmetic/dec.hpp>
#include <boost/preprocessor/arithmetic/inc.hpp>
#include <boost/preprocessor/control/expr_iif.hpp>
#include <boost/preprocessor/list/adt.hpp>
#include <boost/preprocessor/repetition/for.hpp>
#include <boost/preprocessor/repetition/repeat.hpp>
#include <boost/preprocessor/tuple/to_list.hpp>
#include <limits>
#include <pup.h>
#include <vector>

#include "DataStructures/Tensor/TypeAliases.hpp"
#include "Options/Options.hpp"
#include "Parallel/CharmPupable.hpp"
#include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp"  // IWYU pragma: keep
#include "Utilities/Math.hpp"
#include "Utilities/TMPL.hpp"

/// \cond
class DataVector;
/// \endcond

// IWYU pragma: no_forward_declare Tensor

namespace EquationsOfState {
/*!
 * \ingroup EquationsOfStateGroup
 * \brief A spectral equation of state
 *
 * This equation of state is determined as a function of \f$x =
 * \ln(\rho/\rho_0)\f$ where \f$\rho\f$ is the rest mass density and
 * \f$\rho_0\f$ is the provided reference density.  The adiabatic
 * index \f$\Gamma(x)\f$ is defined such that
 * \f{equation}{
 * \frac{d \ln p}{dx} = \Gamma(x) = \sum_{n=0}^N
 * \gamma_n x^n
 * \f}
 *
 * for the set of spectral coefficinets \f$\gamma_n\f$ when
 * \f$0 < x < x_u = \ln(\rho_u/\rho_0)\f$, where \f$\rho_u\f$ is the provided
 * upper density.
 *
 * For \f$ x < 0 \f$, \f$ \Gamma(x) = \gamma_0 \f$.
 *
 * For \f$ x > x_u \f$, \f$ \Gamma(x) = \Gamma(x_u) \f$
 *
 *
 */
class Spectral : public EquationOfState<true, 1> {
 public:
  static constexpr size_t thermodynamic_dim = 1;
  static constexpr bool is_relativistic = true;

  struct ReferenceDensity {
    using type = double;
    static constexpr Options::String help = {"Reference density rho_0"};
    static double lower_bound() { return 0.0; }
  };

  struct ReferencePressure {
    using type = double;
    static constexpr Options::String help = {"Reference pressure p_0"};
    static double lower_bound() { return 0.0; }
  };

  struct Coefficients {
    using type = std::vector<double>;
    static constexpr Options::String help = {"Spectral coefficients gamma_i"};
  };

  struct UpperDensity {
    using type = double;
    static constexpr Options::String help = {"Upper density rho_u"};
    static double lower_bound() { return 0.0; }
  };

  static constexpr Options::String help = {
      "A spectral equation of state.  Defining x = log(rho/rho_0), Gamma(x) = "
      "Sum_i gamma_i x^i, then the pressure is determined from d(log P)/dx = "
      "Gamma(x) for x > 0.  For x < 0 the EOS is a polytrope with "
      "Gamma(x)=Gamma(0).  For x > x_u = log(rho_u/rho_0), Gamma(x) = "
      "Gamma(x_u).\n"
      "To get smooth equations of state, it is recommended that the second "
      "and third supplied coefficient should be 0. It is up to the user to "
      "choose coefficients that are physically reasonable, e.g. that "
      "satisfy causality."};

  using options = tmpl::list<ReferenceDensity, ReferencePressure, Coefficients,
                             UpperDensity>;

  Spectral() = default;
  Spectral(const Spectral&) = default;
  Spectral& operator=(const Spectral&) = default;
  Spectral(Spectral&&) = default;
  Spectral& operator=(Spectral&&) = default;
  ~Spectral() override = default;

  Spectral(double reference_density, double reference_pressure,
           std::vector<double> coefficients, double upper_density);

  EQUATION_OF_STATE_FORWARD_DECLARE_MEMBERS(Spectral, 1)

  WRAPPED_PUPable_decl_base_template(  // NOLINT
      SINGLE_ARG(EquationOfState<true, 1>), Spectral);

  /// The lower bound of the rest mass density that is valid for this EOS
  double rest_mass_density_lower_bound() const override { return 0.0; }

  /// The upper bound of the rest mass density that is valid for this EOS
  double rest_mass_density_upper_bound() const override {
    return std::numeric_limits<double>::max();
  }

  /// The lower bound of the specific internal energy that is valid for this EOS
  /// at the given rest mass density \f$\rho\f$
  double specific_internal_energy_lower_bound(
      const double /* rest_mass_density */) const override {
    return 0.0;
  }

  /// The upper bound of the specific internal energy that is valid for this EOS
  /// at the given rest mass density \f$\rho\f$
  double specific_internal_energy_upper_bound(
      const double /* rest_mass_density */) const override {
    return std::numeric_limits<double>::max();
  }

  /// The lower bound of the specific enthalpy that is valid for this EOS
  double specific_enthalpy_lower_bound() const override { return 1.0; }

 private:
  EQUATION_OF_STATE_FORWARD_DECLARE_MEMBER_IMPLS(1)

  double gamma(const double x) const;
  double integral_of_gamma(const double x) const;
  double chi_from_density(const double density) const;
  double specific_internal_energy_from_density(const double density) const;
  double specific_enthalpy_from_density(const double density) const;
  double pressure_from_density(const double density) const;
  double pressure_from_log_density(const double x) const;
  double rest_mass_density_from_enthalpy(const double specific_enthalpy) const;

  double reference_density_ = std::numeric_limits<double>::signaling_NaN();
  double reference_pressure_ = std::numeric_limits<double>::signaling_NaN();
  std::vector<double> integral_coefficients_{};
  std::vector<double> gamma_coefficients_{};
  double x_max_ = std::numeric_limits<double>::signaling_NaN();
  double gamma_of_x_max_ = std::numeric_limits<double>::signaling_NaN();
  double integral_of_gamma_of_x_max_ =
      std::numeric_limits<double>::signaling_NaN();
  std::vector<double> table_of_specific_energies_{};
  // Information for Gaussian quadrature
  size_t number_of_quadrature_coefs_ =
      std::numeric_limits<size_t>::signaling_NaN();
  std::vector<double> quadrature_weights_{};
  std::vector<double> quadrature_points_{};
};

}  // namespace EquationsOfState

ext: hpp

max_stars_repo_name: Shabibti/spectre

max_issues_repo_name: Shabibti/spectre

max_forks_repo_name: Shabibti/spectre


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 portfolio/builders/swap.hpp
    \brief
    \ingroup builders
*/

#pragma once

#include <boost/make_shared.hpp>
#include <ored/portfolio/builders/cachingenginebuilder.hpp>
#include <ored/portfolio/enginefactory.hpp>
#include <ored/utilities/log.hpp>
#include <ql/pricingengines/swap/discountingswapengine.hpp>
#include <qle/pricingengines/discountingcurrencyswapengine.hpp>
#include <qle/pricingengines/discountingswapenginemulticurve.hpp>

namespace ore {
namespace data {

//! Engine Builder base class for Single Currency Swaps
/*! Pricing engines are cached by currency
    \ingroup builders
*/
class SwapEngineBuilderBase : public CachingPricingEngineBuilder<string, const Currency&> {
public:
    SwapEngineBuilderBase(const std::string& model, const std::string& engine)
        : CachingEngineBuilder(model, engine, {"Swap"}) {}

protected:
    virtual string keyImpl(const Currency& ccy) override { return ccy.code(); }
};

//! Engine Builder for Single Currency Swaps
/*! This builder uses QuantLib::DiscountingSwapEngine
    \ingroup builders
*/
class SwapEngineBuilder : public SwapEngineBuilderBase {
public:
    SwapEngineBuilder() : SwapEngineBuilderBase("DiscountedCashflows", "DiscountingSwapEngine") {}

protected:
    virtual boost::shared_ptr<PricingEngine> engineImpl(const Currency& ccy) override {

        Handle<YieldTermStructure> yts = market_->discountCurve(ccy.code(), configuration(MarketContext::pricing));
        return boost::make_shared<QuantLib::DiscountingSwapEngine>(yts);
    }
};

//! Engine Builder for Single Currency Swaps
/*! This builder uses QuantExt::DiscountingSwapEngineMultiCurve
    \ingroup builders
*/
class SwapEngineBuilderOptimised : public SwapEngineBuilderBase {
public:
    SwapEngineBuilderOptimised() : SwapEngineBuilderBase("DiscountedCashflows", "DiscountingSwapEngineOptimised") {}

protected:
    virtual boost::shared_ptr<PricingEngine> engineImpl(const Currency& ccy) override {

        Handle<YieldTermStructure> yts = market_->discountCurve(ccy.code(), configuration(MarketContext::pricing));
        return boost::make_shared<QuantExt::DiscountingSwapEngineMultiCurve>(yts);
    }
};

//! Engine Builder for Cross Currency Swaps
/*! Pricing engines are cached by currencies (represented as a string list)

    \ingroup builders
*/
class CrossCurrencySwapEngineBuilder
    : public CachingPricingEngineBuilder<string, const std::vector<Currency>&, const Currency&> {
public:
    CrossCurrencySwapEngineBuilder()
        : CachingEngineBuilder("DiscountedCashflows", "DiscountingCrossCurrencySwapEngine", {"CrossCurrencySwap"}) {}

protected:
    virtual string keyImpl(const std::vector<Currency>& ccys, const Currency& base) override {
        std::ostringstream ccyskey;
        ccyskey << base << "/";
        for (Size i = 0; i < ccys.size(); ++i)
            ccyskey << ccys[i] << ((i < ccys.size() - 1) ? "-" : "");
        return ccyskey.str();
    }

    virtual boost::shared_ptr<PricingEngine> engineImpl(const std::vector<Currency>& ccys,
                                                        const Currency& base) override {

        std::vector<Handle<YieldTermStructure>> discountCurves;
        std::vector<Handle<Quote>> fxQuotes;
        for (Size i = 0; i < ccys.size(); ++i) {
            discountCurves.push_back(market_->discountCurve(ccys[i].code(), configuration(MarketContext::pricing)));
            string pair = ccys[i].code() + base.code();
            fxQuotes.push_back(market_->fxSpot(pair, configuration(MarketContext::pricing)));
        }

        return boost::make_shared<QuantExt::DiscountingCurrencySwapEngine>(discountCurves, fxQuotes, ccys, base);
    }
};

} // namespace data
} // namespace ore

max_stars_repo_name: PiotrSiejda/Engine

max_issues_repo_name: PiotrSiejda/Engine

max_forks_repo_name: PiotrSiejda/Engine


text:
#include "dpmeans.hpp"
#include <iostream>
#include <limits>
#include <random>
#include <cmath>
#include <ctime>
#include <unordered_set>
#include <fstream>
#include <vector>

#include <Eigen/Dense>

using namespace std;
using namespace Eigen;

mt19937 rng;

dpmeans::dpmeans(const Ref<const MatrixXf>& X)
{
	K = 1;
	K_init = 4;
	n = X.rows();
	d = X.cols();
	z = VectorXi::Zero(n);
	mu = MatrixXf::Zero(K, d);

	nk = VectorXf::Zero(K);
	pik = VectorXf::Ones(K);

    //init mu
	mu.row(0) = X.colwise().sum() / (float) n;

	cout << "mu = " << endl << mu << endl;
    //init lambda
	lambda = kpp_init(X, K_init);

	max_iter = 100;
	vector<double> obj(max_iter, 0);
	vector<double> em_time(max_iter, 0);
}

float dpmeans::kpp_init(const Ref<const MatrixXf>& X, int k)
{
	// k++ init
	// lambda is max distance to k++ means
	float lambda = 0.0;

	int n = X.rows();
	int d = X.cols();
	MatrixXf mu = MatrixXf::Zero(k, d);

	VectorXf dist = VectorXf::Ones(n);
	VectorXf pdist = VectorXf::Zero(n);
	dist = dist * numeric_limits<float>::max();

	uniform_int_distribution<int> dis0n(0,n-1);
	uniform_real_distribution<float> dis01(0, 1);
	int idx = dis0n(rng);

	mu.row(0) = X.row(idx);
	MatrixXf D = MatrixXf::Zero(n, d);

	for (int i = 1; i < k; ++i)
	{
		D = X - mu.row(i - 1).replicate(n, 1);
		dist = dist.cwiseMin(D.cwiseProduct(D).rowwise().sum());
		//cout << "X = " << endl << X << endl;
		//cout << "mu = " << endl << mu << endl;
		//cout << "dist = " << endl << dist << endl;

		//sample discrete
		pdist = dist / dist.sum();

		for (int j = 1; j < pdist.size(); ++j)
		{
			pdist[j] = pdist[j] + pdist[j - 1];  //cumsum
		}

		int pidx = 0;
		float z01 = dis01(rng);
		for (int item = 0; item < pdist.size(); ++item)
		{
			if (z01 < pdist[item])
			{
				pidx = item;
				break;
			}
		}
		mu.row(i) = X.row(pidx);
		lambda = dist.maxCoeff();
	}

	return lambda;
}

VectorXi dpmeans::dpmeans_fit(const Ref<const MatrixXf>& X)
{

	int n = X.rows();
	int d = X.cols();
	int K = this->K;
	int max_iter = this->max_iter;

	double obj_tol = 1e-3;
	cout << "running dp-means..." << endl;

	//cout << "X = " << endl << X << endl;
	for (int iter = 0; iter < max_iter; ++iter)
	{
		clock_t tic = clock();
		MatrixXf dist = MatrixXf::Zero(n, K);

		//assignment step
		MatrixXf Xm = MatrixXf::Zero(n, d);
		for (int k = 0; k < K; ++k)
		{
			Xm = X - mu.row(k).replicate(n, 1);
			dist.col(k) = Xm.cwiseProduct(Xm).rowwise().sum();
		}
		//cout << "mu = " << endl << mu << endl;
		//cout << "dist = " << endl << dist << endl;

		//update labels
		VectorXf dmin = VectorXf::Zero(n);
		for (int ridx = 0; ridx < n; ++ridx)
		{
			MatrixXf::Index minIndex;
			dmin(ridx) = dist.row(ridx).minCoeff(&minIndex);
			z(ridx) = minIndex;
		}
		//cout << "dmin = " << endl << dmin << endl;
		//cout << "z = " << endl << z << endl;
		//cout << "lambda = " << endl << lambda << endl;

		VectorXi dmin_idx = VectorXi::Zero(n);
		for (int ridx = 0; ridx < n; ++ridx)
		{
			if (dmin(ridx) > lambda)
				dmin_idx(ridx) = 1;
		}
		int num_new = dmin_idx.sum();
		//cout << "num_new = " << endl << num_new << endl;

		if (num_new > 0)
		{
			//create a new cluster for points
			//with dmin > lambda
			K = K + 1;
			//cout << "K = " << endl << K << endl;

			VectorXf new_mean = VectorXf::Zero(d);
			for (int ridx = 0; ridx < n; ++ridx)
			{
				// if assigned to new cluster
				if (dmin_idx(ridx) > 0)
				{
					z(ridx) = K-1;  // cluster labels: [0,...,K-1]
					new_mean = new_mean + X.row(ridx).transpose();
				}
			}
			//cout << "z = " << endl << z << endl;
			mu.conservativeResize(mu.rows() + 1, NoChange);
			//cout << "new mean: " << endl << new_mean << endl;
			//cout << "num_new: " << endl << num_new << endl;
			mu.row(K-1) = new_mean / (float) num_new;  // add new mean

			MatrixXf Xm = MatrixXf::Zero(n, d);
			Xm = X - mu.row(K-1).replicate(n, 1);
			dist.conservativeResize(NoChange, dist.cols() + 1);
			dist.col(K-1) = Xm.cwiseProduct(Xm).rowwise().sum(); //add dist to new mean
			//cout << "dist = " << endl << dist << endl;

		}

		//update step
		nk = VectorXf::Zero(K);
		mu = MatrixXf::Zero(K, d);
		//for (int ridx = 0; ridx < z.size(); ++ridx) { nk(z(ridx))++; } //histogram
		for (int k = 0; k < K; ++k)
		{
			nk(k) = (float) (z.array() == k).count();

			for (int ridx = 0; ridx < z.size(); ++ridx)
			{
				if (z(ridx) == k)
				{
					mu.row(k) = mu.row(k) + X.row(ridx);
				}

			}
			mu.row(k) = mu.row(k) / nk(k);
		}
		pik = nk / nk.sum();
		cout << "mu = " << endl << mu << endl;
		cout << "nk = " << endl << nk << endl;
		cout << "pik = " << endl << pik << endl;

		//compute objective
		int kidx = 0;
		double suml2sq = 0.0;
		for (int ridx = 0; ridx < n; ++ridx)
		{
			kidx = z(ridx);
			suml2sq = suml2sq + dist(ridx, kidx);
		}
		suml2sq += lambda * K;
		obj.push_back(suml2sq);
		cout << "obj = " << endl << suml2sq << endl;

		//check convergence
		if (iter > 0 && abs(obj.at(iter) - obj.at(iter - 1)) < obj_tol *obj.at(iter))
		{
			cout << "converged in " << iter << " iterations." << endl;
			break;
		}

		clock_t toc = clock();
		double elapsed_sec = double(toc - tic) / (double) CLOCKS_PER_SEC;
		cout << "elapsed sec = " << endl << elapsed_sec << endl;
		em_time.push_back(elapsed_sec);
	}

	this->K = K;

	return z;
}

float dpmeans::compute_nmi(const Ref<const VectorXi>& z1, const Ref<const VectorXi>& z2)
{
	// compute normalized mutual information
	float nmi = 0.0;

	int n = z1.size();

	//compute number of unique
	//labels in z1 and z2
	unordered_set<int> s1;
	unordered_set<int> s2;

	int item;
	item = 0;
	for (int idx = 0; idx < n; ++idx)
	{
		item = z1(idx);
		s1.insert(item);
	}

	item = 0;
	for (int idx = 0; idx < n; ++idx)
	{
		item = z2(idx);
		s2.insert(item);
	}

	int k1 = s1.size();
	int k2 = s2.size();

	VectorXf nk1 = VectorXf::Zero(k1);
	VectorXf nk2 = VectorXf::Zero(k2);
	for (int idx = 0; idx < k1; ++idx)
	{
		nk1(idx) = (float) (z1.array() == idx).count();
	}
	for (int idx = 0; idx < k2; ++idx)
	{
		nk2(idx) = (float) (z2.array() == idx).count();
	}

	VectorXf pk1 = nk1 / (float)nk1.sum();
	VectorXf pk2 = nk2 / (float)nk2.sum();

	MatrixXi nk12 = MatrixXi::Zero(k1, k2);
	for (int idx1 = 0; idx1 < k1; ++idx1)
	{
		for (int idx2 = 0; idx2 < k2; ++idx2)
		{
			VectorXi b1 = VectorXi::Zero(n);
			b1 = (z1.array() == idx1).select(VectorXi::Ones(n), VectorXi::Zero(n));

			VectorXi b2 = VectorXi::Zero(n);
			b2 = (z2.array() == idx2).select(VectorXi::Ones(n), VectorXi::Zero(n));

			VectorXi b1b2 = b1.cwiseProduct(b2);

			nk12(idx1, idx2) = b1b2.sum();
		}
	}
	MatrixXf pk12 = nk12.cast<float>() / n;
	cout << "nk12 = " << endl << nk12 << endl;
	cout << "pk12 = " << endl << pk12 << endl;

	VectorXf logpk1 = (pk1 + numeric_limits<float>::epsilon()*VectorXf::Ones(k1)).array().log();
	VectorXf logpk2 = (pk2 + numeric_limits<float>::epsilon()*VectorXf::Ones(k2)).array().log();
	MatrixXf logpk12 = (pk12 + numeric_limits<float>::epsilon()*MatrixXf::Ones(k1, k2)).array().log();

	float Hx = -pk1.dot(logpk1);
	float Hy = -pk2.dot(logpk2);
	float Hxy = -pk12.cwiseProduct(logpk12).sum();

	float MI = Hx + Hy - Hxy;
	nmi = MI / (float) (0.5*(Hx + Hy));

	return nmi;
}

void dpmeans::display_params()
{
	cout << "K = " << this->K << endl;
	cout << "K_init = " << this->K_init << endl;
	//cout << "Labels: " << endl << this->z << endl;
	cout << "Means: " << endl << this->mu << endl;
	cout << "Counts: " << endl << this->nk << endl;
	cout << "Proportions: " << endl << this->pik << endl;
	cout << "Lambda: " << this->lambda << endl;
}

void load_iris(Ref<MatrixXf> X, Ref<VectorXi> y)
{

	//read-in labels
	int cnt = 0;
	string line;
	ifstream fin2("./data/iris_labels.txt");
	if (fin2.is_open())
	{
		while (getline(fin2, line))
		{
			y(cnt) = stoi(line);
			cnt++;
		}
		fin2.close();
	}
	else cout << "Unable to open fin2\n";
	//cout << "y = " << endl << y << endl;
	//cout << "y.size = " << endl << y.size() << endl;

	//read-in data
	cnt = 0;
	int nrows = 0;
	int ncols = 0;
	ifstream fin3("./data/iris_data.txt");
	if (fin3.is_open())
	{
		fin3 >> nrows >> ncols;
		for (int row = 0; row < nrows; row++)
			for (int col = 0; col < ncols; col++)
			{
				float num = 0.0;
				fin3 >> num;
				X(row, col) = num;
			}
		fin3.close();
	}
	else cout << "Unable to open fin3\n";
	//cout << "X = " << endl << X << endl;

}

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

	//generate data
	//int rows = 150, cols = 4;
	//MatrixXf X = MatrixXf::Random(rows, cols);

	//load iris
	int nrows = 0;
	int ncols = 0;
	string line;
	ifstream fin1("./data/iris_data.txt");
	if (fin1.is_open())
	{
		fin1 >> nrows >> ncols;
		fin1.close();
	}
	else cout << "Unable to open data file\n";
	cout << "nrows = " << nrows << endl;
	cout << "ncols = " << ncols << endl;

	VectorXi y = VectorXi::Zero(nrows);
	MatrixXf X = MatrixXf::Zero(nrows, ncols);

	load_iris(X,y);
	//cout << "y = " << endl << y << endl;
	//cout << "X = " << endl << X << endl;

	dpmeans dp(X);
	dp.display_params();

	VectorXi z;
	z = dp.dpmeans_fit(X);
	//cout << "z = " << endl << z << endl;

	float nmi;
	nmi = dp.compute_nmi(z, y);
	cout << "nmi = " << endl << nmi << endl;

	return 0;
}

ext: cpp

max_stars_repo_name: vishalbelsare/cpp