text:
/*
 * PowerStateGraph.cpp
 *
 *  Created on: 3 Aug 2011
 *      Author: jack
 */

#include "PowerStateGraph.h"
#include "AggregateData.h"
#include <iostream>
#include <list>
#include <boost/graph/graphviz.hpp>
#include <cstdio> // sprintf

using namespace std;

PowerStateGraph::PowerStateGraph()
: totalCount(0), aggData(0)
{
    using namespace boost;

    // add a vertex to represent "off"
    // uses Statistic's default constructor to make a Statistic with all-zeros
    offVertex = add_vertex(powerStateGraph);

}

void PowerStateGraph::addItemToEdgeHistory(
        const PSGraph::edge_descriptor& edge
        )
{
    // Check if edge histories are disabled
    if (EDGE_HISTORY_SIZE == 0)
        return;

    if (edgeHistory.size() >= EDGE_HISTORY_SIZE)
        edgeHistory.erase( edgeHistory.begin() );

    edgeHistory.push_back(edge);
}

const Statistic< double >& PowerStateGraph::getEnergyConsumption() const
{
    return energyConsumption;
}

/**
 * @brief Update or initialise Power State Graph.
 */
void PowerStateGraph::update(
        const Signature& sig,
        const bool verbose
        )
{
    const double energyConsumptionFromSig = sig.getEnergyConsumption();
    energyConsumption.update( energyConsumptionFromSig );
    cout << "Energy consumption from sig" << sig.getID() << " = "
         << energyConsumptionFromSig / J_PER_KWH << " kWh" << endl;

    edgeHistory.clear();

    // get the gradient spikes for the signature
    list<Signature::Spike> spikes = sig.getDeltaSpikes();
    list<Signature::Spike>::iterator spike;
    PSGraph::vertex_descriptor targetVertex, sourceVertex=offVertex;
    size_t indexOfLastAcceptedSpike = 0;
    size_t start=0, end=0;

    // for each spike, locate the samples immediately before and immediately after the spike
    const size_t WINDOW = 8; // how far either side of the spike will we look?

    // take just the top TOP_SLICE_SIZE (whilst ordered by absolute value)
    const size_t TOP_SLICE_SIZE = 10;
    if (spikes.size() > TOP_SLICE_SIZE) {
        list<Signature::Spike>::iterator it = spikes.begin();
        advance( it, TOP_SLICE_SIZE );
        spikes.erase( it, spikes.end() );
    }

    // re-order by index (i.e. by time)
    spikes.sort( Signature::Spike::compareIndexAsc );

    for (spike = spikes.begin(); spike!=spikes.end(); spike++) {

        // calculate the start index for the 8 pre-spike samples
        start = ((spike->index > WINDOW) ? (spike->index - WINDOW) : 0 );

        // calculate the end index for the 8 post-spike samples
        end = spike->index + WINDOW + spike->n + 1;
        if (end > sig.getSize())
            end = sig.getSize();

        // Create statistics for pre- and post-spike
        Statistic<Sample_t> preSpikePowerState(sig, start, spike->index);
        Statistic<Sample_t> postSpikePowerState(sig, (spike->index + spike->n + 1), end);

        // check to see if spikes needs to be rejected and that the spikes aren't too close
        if (rejectSpike(preSpikePowerState, postSpikePowerState) || start < indexOfLastAcceptedSpike) {
            if (verbose) cout << " REJECT";
        } else {

            // Create inter-spike stats
            Statistic<Sample_t> betweenSpikesPowerState(
                    sig,
                    (spike->index + spike->n + 1),
                    indexOfNextSpike( spikes, spike, sig )
                    );

            targetVertex =
                    updateOrInsertVertex(
                            sig,
                            postSpikePowerState,
                            betweenSpikesPowerState
                            );

            if (sourceVertex != targetVertex) {
                updateOrInsertEdge( sourceVertex, targetVertex,
                        (spike->index - indexOfLastAcceptedSpike), spike->delta );
            }

            if (verbose)
                printSpikeInfo( spike, start, end, preSpikePowerState, postSpikePowerState, sig);

            sourceVertex = targetVertex;
            indexOfLastAcceptedSpike = spike->index;

        }
    }
}

/**
 * @brief Useful for diagnostics.
 */
void PowerStateGraph::printSpikeInfo(
        const list<Signature::Spike>::iterator spike,
        const size_t start,
        const size_t end,
        const Statistic<Sample_t>& preSpikePowerState,
        const Statistic<Sample_t>& postSpikePowerState,
        const Signature& sig
        ) const
{

    cout << endl
            << "SPIKE: index=" << spike->index
            << ", delta=" << spike->delta
            << ", duration=" << spike->n
            << endl
            << " KEEP" << endl
            << "Before=" << preSpikePowerState << endl
            << "After =" << postSpikePowerState  << endl;

    for (size_t i=start; i<end; i++ ) {
        cout << i << "\t" << sig[i];
        if (i < spike->index)
            cout << " before";
        else if (i == spike->index)
            cout << " <--SPIKE";
        else if (i < (spike->index + spike->n))
            cout << " <--Spike continues";
        else
            cout << " after";
        cout << endl;
    }

}

/**
 * @brief Determine the index of the next spike after @c spike
 */
const size_t PowerStateGraph::indexOfNextSpike(
        const list<Signature::Spike>& spikes,
        list<Signature::Spike>::iterator spike,
        const Signature& sig
        ) const
{
    list<Signature::Spike>::iterator spikePlusOne = spike;
    spikePlusOne++;
    if (spikePlusOne != spikes.end())
        return spikePlusOne->index;
    else
        return sig.getSize() - 1;
}

/**
 * @brief Decides whether or not to reject this spike based on the
 * @c before and @c after statistics.
 */
const bool PowerStateGraph::rejectSpike(
        const Statistic<Sample_t>& before,
        const Statistic<Sample_t>& after,
        const bool verbose
        ) const
{
    // check the stdev isn't too large
    if ( (before.stdev > fabs(before.mean)) ||
         ( after.stdev > fabs(after.mean))   ) {
        if (verbose) cout << "stdev too big";
        return true;
    }

    // check the means aren't too close
    if (Utils::within(before.mean, after.mean, Utils::highest(before.stdev, after.stdev)*2) ||
        Utils::within(before.mean, after.mean, Utils::highest(before.mean, after.mean)*0.2)     ) {
        if (verbose) cout << "means too close";
        return true;
    }

    // if we get to here then don't reject
    return false;
}

/**
 * @brief Attempts to find an existing vertex which is
 * statistically similar to @c stat.  If an existing similar
 * vertex is found then the vertex's stats are updated
 * with the new data points. If a similar vertex is not found,
 * and new vertex is inserted.  Either way, a vertex_descriptor
 * is returned to the new or existing similar vertex.
 */
PowerStateGraph::PSGraph::vertex_descriptor PowerStateGraph::updateOrInsertVertex(
        const Signature& sig, /**< source of the raw data */
        const Statistic<Sample_t>& postSpikePowerState,
        const Statistic<Sample_t>& betweenSpikesPowerState,
        const bool verbose  /**< cout debugging messages? */
    )
{
    bool foundSimilar; // return param for mostSimilarVertex()
    PSGraph::vertex_descriptor vertex = mostSimilarVertex( &foundSimilar, postSpikePowerState );
    // mostSimilarVertex() uses a T-Test to find the powerStateGraph
    // vertex with a mean most similar to postSpikePowerState.

    if (verbose)
        cout << " mostSimlar = " << powerStateGraph[vertex].postSpike << endl;

    if ( foundSimilar ) {
        if (vertex != offVertex) {
            // Update an existing vertex

            if (verbose) {
                cout << endl
                        << "Updating existing vertex: " << endl
                        << "    existing vertex = " << powerStateGraph[vertex].postSpike << endl
                        << "    new stat        = " << postSpikePowerState << endl;
            }

            powerStateGraph[vertex].postSpike.update( postSpikePowerState );
            powerStateGraph[vertex].betweenSpikes.update( betweenSpikesPowerState );

            if (verbose)
                cout << "    merged          = " << powerStateGraph[vertex].postSpike << endl;
        }
    } else {
        // Add a new vertex
        vertex = add_vertex(powerStateGraph);
        powerStateGraph[vertex].postSpike = postSpikePowerState;
        powerStateGraph[vertex].betweenSpikes = betweenSpikesPowerState;

        if (verbose) {
            cout << endl
                 << "Adding new vertex: " << powerStateGraph[vertex].postSpike << endl;
        }

    }
    return vertex;

}

/**
 * @brief Find the vertex statistically most similarto @c stat.
 *
 * @return vertex descriptor of best fit.
 *         @c success is also used as a return parameter.
 *
 * @todo remove commented-out code if it proves to be of no use
 */
PowerStateGraph::PSGraph::vertex_descriptor PowerStateGraph::mostSimilarVertex(
        bool * success, /**< return parameter.  Did we find a satisfactory match? */
        const Statistic<Sample_t>& stat, /**< stat to find in graph vertices */
        const double ALPHA /**< significance level (what constitutes as a "satisfactory" match?) */
    ) const
{
    PSGraph::vertex_descriptor vertex=0;
    std::pair<PSG_vertex_iter, PSG_vertex_iter> vp;
    double tTest, highestTTest=0;

    // Find the best fit
    for (vp = boost::vertices(powerStateGraph); vp.first != vp.second; ++vp.first) {
        // t test
        tTest = stat.tTest( powerStateGraph[*vp.first].postSpike );
        if (tTest > highestTTest) {
            highestTTest = tTest;
            vertex = *vp.first;
        }

    }

    // Check whether the best fit is satisfactory
    if ( (highestTTest > (ALPHA/2)) // T-Test
       ||  ( Utils::within
                 (
                 powerStateGraph[vertex].postSpike.mean,
                 stat.mean,
                 (Utils::highest(powerStateGraph[vertex].postSpike.mean, stat.mean) * 0.2)
                 ) // check if the means are within 20% of each other
         )
      )
        *success = true;
    else
        *success = false;

    return vertex;
}

/**
 * @brief Update or Insert a new edge into powerStateGraph.
 *
 * - If an edge already exists between beforeVertex and afterVertex then
 *     - check all out-edges from beforeVertex
 *     - if any out-edge has the same history as the current history
 *           and the same sign delta then
 *           and and is within 3 stdevs of the existing delta then:
 *         - update that edge's statistics
 *     - else
 *         - create a new edge.
 * - else
 *     - create a new edge
 */
void PowerStateGraph::updateOrInsertEdge(
        const PSGraph::vertex_descriptor& beforeVertex,
        const PSGraph::vertex_descriptor& afterVertex,
        const size_t samplesSinceLastSpike,
        const double spikeDelta,
        const bool verbose
        )
{
    totalCount++;

    if (verbose) cout << "-------------------" << endl;

    PSGraph::edge_descriptor existingEdge, newEdge;
    bool edgeExistsAlready;
    tie(existingEdge, edgeExistsAlready) = boost::edge(beforeVertex, afterVertex, powerStateGraph);

    if (verbose) {
        if ( edgeExistsAlready ) {
            cout << "edgeExistsAlready" << endl;
        } else {
            cout << "edge does not Exist Already" << endl;
        }
    }

    if ( edgeExistsAlready ) {

        if (verbose) {
            for (list<PSGraph::edge_descriptor>::iterator edge=edgeHistory.begin(); edge!=edgeHistory.end(); edge++) {
                cout << "                              edgeHistory = " << *edge << endl;
            }
            for (list<PSGraph::edge_descriptor>::iterator edge=powerStateGraph[existingEdge].edgeHistory.begin();
                    edge!=powerStateGraph[existingEdge].edgeHistory.end(); edge++) {
                cout << "powerStateGraph[existingEdge].edgeHistory = " << *edge << endl;
            }
            cout << endl;
        }

        // check if any of the out edges from beforeVertex have the same
        // history as our current history... if so, update that edge.
        PSGraph::out_edge_iterator out_e_i, out_e_end;
        tie(out_e_i, out_e_end) = out_edges(beforeVertex, powerStateGraph);
        for (; out_e_i != out_e_end; out_e_i++) {
            // check if the edge has the same history and same sign delta and is within 3 stdevs of the existing delta
            if ( edgeListsAreEqual(powerStateGraph[*out_e_i].edgeHistory, edgeHistory) &&
                    Utils::sameSign(powerStateGraph[*out_e_i].delta.mean, spikeDelta ) &&
                    Utils::within(powerStateGraph[*out_e_i].delta.mean,
                            spikeDelta, powerStateGraph[*out_e_i].delta.nonZeroStdev()*3 )) {

                if (verbose) {
                    cout << "edge histories the same. merging with" << *out_e_i
                            << powerStateGraph[*out_e_i].delta << endl;
                }

                // update existing edge's stats
                powerStateGraph[*out_e_i].delta.update( spikeDelta );
                powerStateGraph[*out_e_i].duration.update( samplesSinceLastSpike );
                powerStateGraph[*out_e_i].count++;

                addItemToEdgeHistory( *out_e_i );

                return;
            }
        }
    }

    // if we get to here then we know that either:
    //    an edge doesn't already exist between these 2 vertices
    // or
    //    an edge does exist but it doesn't share our history.
    // either way, we need to add a new edge:

    tie(newEdge, edgeExistsAlready) = boost::add_edge(beforeVertex, afterVertex, powerStateGraph);
    powerStateGraph[newEdge].delta    = Statistic<double>( spikeDelta );
    powerStateGraph[newEdge].duration = Statistic<size_t>( samplesSinceLastSpike );
    powerStateGraph[newEdge].count    = 1;
    powerStateGraph[newEdge].edgeHistory = edgeHistory;

    if (verbose) {
        cout << "adding new edge" << newEdge << endl
            << "   delta stats    = " << powerStateGraph[newEdge].delta << endl
            << "   duration stats = " << powerStateGraph[newEdge].duration << endl;
    }

    addItemToEdgeHistory( newEdge );
}

/**
 * @brief Check if edge list @c a and edge list @c b are equal.
 */
const bool PowerStateGraph::edgeListsAreEqual(
        const list< PSGraph::edge_descriptor >& a,
        const list< PSGraph::edge_descriptor >& b,
        const bool verbose
        ) const
{
    // Check if edge histories are disabled
    if (EDGE_HISTORY_SIZE == 0)
        return true;

    //******************* DIAGNOSTICS *******************
    if (verbose ) {
        cout << "************" << endl
                << "comparing lists:" << endl;

        for (list< PSGraph::edge_descriptor >::const_iterator i=a.begin(); i!=a.end(); i++) {
            cout << *i << endl;
        }
        cout << "----" << endl;
        for (list< PSGraph::edge_descriptor >::const_iterator i=b.begin(); i!=b.end(); i++) {
            cout << *i << endl;
        }
    }//___________________________________________________

    if (a.size() != b.size()) {
        if (verbose) cout << "sizes not equal" << endl;
        return false;
    }

    if (verbose) cout << "************" << endl;

    list<PSGraph::edge_descriptor>::const_iterator a_i, b_i;

    a_i = a.begin();
    b_i = b.begin();

    while (a_i != a.end() && b_i != b.end()) {

        if (verbose) cout << "*a_i=" << *a_i << " *b_i=" << *b_i << endl;

        if ( source(*a_i, powerStateGraph) != source(*b_i, powerStateGraph) &&
             target(*a_i, powerStateGraph) != target(*b_i, powerStateGraph)   ) {

            if (verbose) cout << "edge list are not equal" << endl << endl;

            return false;
        }

        a_i++;
        b_i++;
    }

    return true;
}

/**
 * @brief Update (or create) directional edges between power states (vertices)
 *        so the edges are consistent with the observed transitions between
 *        power states in @c sig.
 *
 * Basic strategy is to:
 *

 *
retrieve gradient spikes from @c sig, extract the most salient
 *      and put the spikes into temporal order.

 *
Create stats for the data points between each spike. Use mostSimilarVertex()
 *      to determine which of the existing power state vertices this belongs to.

 *
Check to see if there's an existing edge representing the power state transition,
 *      if so then check the @c delta and @c duration stats for the edge and update if necessary.
 *      If no edge exists then create one with the necessary @c delta and @c duration,
 *      with a hard-coded stdev.

 *

 *
 *  @deprecated just leaving this here to illustrate one of the strategies I tried.
 *  This has been superseded by update() (which used to be updateVertices) and updateOrInsertEdge()
 *
 */
void PowerStateGraph::updateEdges( const Signature& sig )
{
    Statistic<Sample_t> nextPowerState, prevPowerState;
    PSGraph::vertex_descriptor mostSimVertex, previousVertex;
    pair<PSGraph::edge_descriptor, bool> addedEdge;
    size_t previousIndex;

    // get the gradient spikes for the first signature
    list<Signature::Spike> spikes = sig.getDeltaSpikes();

    // take just the top ten (whilst ordered by absolute value)
    if (spikes.size() > 10) {
        list<Signature::Spike>::iterator it = spikes.begin();
        advance( it, 10 );
        spikes.erase( it, spikes.end() );
    }

    // re-order by index (i.e. by time)
    spikes.sort( Signature::Spike::compareIndexAsc );

    // calculate stats for the signature's values between start and first spike
    // as long as the first spike->index > 1
    previousVertex = offVertex;
    previousIndex  = 0;

    list<Signature::Spike>::iterator spike, prevSpike;
    prevSpike = spike = spikes.begin();
    advance(spike, 1);
    for (; spike!=spikes.end(); spike++ ) {

        nextPowerState = Statistic<Sample_t>( sig, prevSpike->index+1, prevSpike->index+10); // spike->index);

        bool foundSimilar;
        mostSimVertex = mostSimilarVertex( &foundSimilar, nextPowerState, 0.1 );

        if (foundSimilar && previousVertex!=mostSimVertex) {
            /* Add an edge from previousVertex to mostSimVertex.
             * If an edge already exists then boost::add_edge will
             * return an edge_descriptor to that edge. */
            addedEdge = boost::add_edge(previousVertex, mostSimVertex, powerStateGraph);

            if (addedEdge.second) { // then there was not an edge already in the graph
                powerStateGraph[addedEdge.first].delta =
                        Statistic<double>( prevSpike->delta );
                powerStateGraph[addedEdge.first].duration =
                        Statistic<size_t>( prevSpike->index - previousIndex );
                powerStateGraph[addedEdge.first].count = 1;
                totalCount++;
            } else { // there was an edge already in the graph so update that edge
                powerStateGraph[addedEdge.first].delta.update( prevSpike->delta );
                powerStateGraph[addedEdge.first].duration.update( prevSpike->index - previousIndex );
                powerStateGraph[addedEdge.first].count++;
                totalCount++;
            }

            previousIndex  = prevSpike->index;
            previousVertex = mostSimVertex;
        }
        prevSpike = spike;
    }

    // add the last edge back to "off"
    addedEdge = boost::add_edge(previousVertex, offVertex, powerStateGraph);
    if (addedEdge.second) { // then there was not an edge already in the graph
        powerStateGraph[addedEdge.first].delta =
                Statistic<double>( prevSpike->delta );
        powerStateGraph[addedEdge.first].duration =
                Statistic<size_t>( prevSpike->index - previousIndex );
        powerStateGraph[addedEdge.first].count = 1;
        totalCount++;
    } else { // there was an edge already in the graph so update that edge
        powerStateGraph[addedEdge.first].delta.update( prevSpike->delta );
        powerStateGraph[addedEdge.first].duration.update( prevSpike->index - previousIndex );
        powerStateGraph[addedEdge.first].count++;
        totalCount++;
    }

}

/**
 * @brief Produce a graphviz output
 */
void PowerStateGraph::writeGraphViz(ostream& out)
{
    write_graphviz(
            out,
            powerStateGraph,
            PSG_vertex_writer(powerStateGraph),
            PSG_edge_writer(powerStateGraph)
            );
}

/**
 * @brief This is the main public interface to the disaggregation algorithm
 *
 * Each candidate solution is represented as a
 * tree.  Each vertex
 * on this tree is a spike found in the aggregate data.
 * The edge weights are the mean of the probability density
 * functions for the spike size and the timing.  If an
 * edge described in PowerStateGraph (from training) cannot be found
 * in the aggregate data then this candidate is discarded.  If
 * we successfully get from vertex0 back to vertex0 then the
 * shortest path through the tree is calculated and saved.
 *
 *

 *
Retrieve edge @c e which connects @c vertex0 (offVertex)
 *     to vertex1.  Look through the AggregateData
 *     searching for any spike within a certain number of
 *     standard deviations of @c e.delta.

 *
When a spike is found, start a tree structure and
 *     start looking for the subsequent edges learnt during training.

 *
Look for the delta corresponding to each out edge from current vertex.
 *     Store the UNIX timestamp of each candidate.

 *

 *
 * @return a list of UNIX times when the device starts
 */
const list<PowerStateGraph::Fingerprint> PowerStateGraph::disaggregate(
        const AggregateData& aggregateData, /**< A populated array of AggregateData */
        const bool keep_overlapping, /**< Should we keep or remove overlapping candidates? */
        const bool verbose
        )
{
    cout << endl << "***** TRAINING FINISHED. DISAGGREGATION STARTING. *****" << endl << endl;
    cout << "Finding start deltas...";
    cout.flush();

    if (num_vertices( powerStateGraph ) < 2) {
        Utils::fatalError( "powerStateGraph is empty. Cannot continue with disaggregation." );
    }

    // Store a pointer to aggregateData for use later.
    aggData = &aggregateData;

    /* Fingerprint is a struct for bundling start
     * timestamp, duration, energy and avLikelihood */
    list<Fingerprint> fingerprintList; // what we return
    Fingerprint candidateFingerprint;

    /* Load the stats from the first PowerStateGraph edge.
     * The Boost Graph Lib provides an out_edges() function which returns
     * a pair of edge iterators: the first iterator points to the first edge;
     * the second iterator points 1 past the last edge. Dereferencing an
     * edge iterator returns an edge descriptor. tie() allows easy
     * access to each element of the pair of edge iterators. */
    PSG_out_edge_iter out_i, out_end;
    tie(out_i, out_end) = out_edges(offVertex, powerStateGraph);
    PowerStateEdge firstEdgeStats = powerStateGraph[*out_i];

    // Search through aggregateData for possible start spikes
    list<AggregateData::FoundSpike> posStartSpikes
        = aggregateData.findSpike(firstEdgeStats.delta);

    cout << " found " << posStartSpikes.size() << " possible start deltas. Following through... " << endl;

    // For each possible start spike in the delta aggregate, attempt
    // to find all the subsequent spikes specified by powerStateGraph edges.
    for (list<AggregateData::FoundSpike>::const_iterator posStartSpike=posStartSpikes.begin();
            posStartSpike!=posStartSpikes.end();
            posStartSpike++) {

        // Attempt to recursively "follow through" from this posStartSpike.
        // initTraceToEnd() calls the recursive function traceToEnd().
        candidateFingerprint = initTraceToEnd(
                *posStartSpike,
                ( posStartSpike->timestamp - firstEdgeStats.duration.mean ) // time the device probably started
                );

        // If this start spike was successfully traced all the way to
        // an off power state then add this item to fingerprintList.
        if ( candidateFingerprint.avLikelihood != -1 ) {
            fingerprintList.push_back( candidateFingerprint );
            if (verbose)
                cout << endl << "candidate found at " << endl << candidateFingerprint << endl;
            else {
                cout << ".";
                cout.flush();
            }
        }
    }

    cout << endl;

    if ( fingerprintList.empty() ) {
        cout << "No signatures found." << endl;
    } else {
        if (!keep_overlapping) {
            removeOverlapping( &fingerprintList );
        }

        displayAndPlotFingerprintList( fingerprintList, aggregateData.getFilename() );
    }

    return fingerprintList;
}

void PowerStateGraph::displayAndPlotFingerprintList(
        const list< Fingerprint >& fingerprintList,
        const string& aggDataFilename
        ) const
{
    cout << "Displaying disaggregation data and dumping to file for plotting. " << endl;

    fstream fs;
    Utils::openFile(fs, DATA_OUTPUT_PATH + "disagg.dat", fstream::out);
    list<Fingerprint>::const_iterator disagItem;
    size_t count = 0;
    for (disagItem=fingerprintList.begin(); disagItem!=fingerprintList.end(); disagItem++) {
        cout << endl << "Candidate fingerprint found: " << endl << *disagItem << endl;
        count++;
        for (list<TimeAndPower>::const_iterator tap_i=disagItem->timeAndPower.begin();
                tap_i!=disagItem->timeAndPower.end(); tap_i++) {
            fs << tap_i->timestamp << "\t" << tap_i->meanPower << endl;
        }
    }
    cout << endl
            << "Found " << count << " candidate" << (count==1 ? ". " : "s." ) << endl;

    // Calculate size of x-axis border
    const size_t begOfFirstFingerprint = fingerprintList.front().timestamp;
    const size_t endOfLastFingerprint  = fingerprintList.back().timestamp + fingerprintList.back().duration;
    const size_t BORDER = (endOfLastFingerprint - begOfFirstFingerprint) / 10;

    // Set Plot variables
    GNUplot::PlotVars pv;
    pv.inFilename  = "disagg";
    pv.outFilename = "disagg";
    pv.title       = "Automatic disaggregation for " + deviceName;
    pv.xlabel      = "time";
    pv.ylabel      = "power (kW)";

    // set sensible xrange, but only if we're not dealing with synthetic data with a low timecode
    if (fingerprintList.front().timestamp > BORDER) {
        const size_t dstOffset = 3600; // to correct for BST
        pv.plotArgs    = "[\"" + Utils::size_t_to_s( fingerprintList.front().timestamp - BORDER + dstOffset) + "\":\""
                + Utils::size_t_to_s( fingerprintList.back().timestamp + fingerprintList.back().duration + BORDER + dstOffset) + "\"]";
    }

    pv.data.push_back(
            GNUplot::PlotData(
                    "disagg", "Automatically determined device fingerprint", "DISAGG"));
    pv.data.push_back(
            GNUplot::PlotData(
                      aggDataFilename, "Aggregate data", "AGGDATA", false));

    GNUplot::plot( pv );
}

/**
 * Remove any overlapping list entries and leave the one with the highest likelihood
 */
void PowerStateGraph::removeOverlapping(
        list<Fingerprint> * fingerprintList, /**< Input and output parameter */
        const bool verbose
        )
{
    list<Fingerprint>::iterator currentDisagItem, prevDisagItem;
    size_t count = 0;

    currentDisagItem = prevDisagItem = fingerprintList->begin();
    advance( currentDisagItem, 1 );
    while ( currentDisagItem != fingerprintList->end() ) {

        // check to see whether the two disag items overlap
        if ( (prevDisagItem->timestamp + prevDisagItem->duration) >= currentDisagItem->timestamp ) {

            count++;

            if (verbose) {
                cout << "Overlap detected between items with start timestamps "
                    << prevDisagItem->timestamp << " and " << currentDisagItem->timestamp << " ...";
            }

            // they overlap.  so find which needs to be replaced
            if (prevDisagItem->avLikelihood < currentDisagItem->avLikelihood) {
                if (verbose) cout << "erasing item with timestamp " << prevDisagItem->timestamp << endl;
                fingerprintList->erase( prevDisagItem );
                prevDisagItem = currentDisagItem++;
            } else {
                if (verbose) cout << "erasing item with timestamp " << currentDisagItem->timestamp << endl;
                fingerprintList->erase( currentDisagItem++ );
                // don't change prevDisagItem
            }

        } else {
            prevDisagItem = currentDisagItem++;
        }
    }

    if (count == 0) {
        cout << "No candidate fingerprints overlap." << endl;
    } else {
        cout << "Removed " << count << " overlapping candidate fingerprints." << endl;
    }
}

/**
 *
 * @return DisaggregatedStruct.likelihood will be set to -1 if
 * this looks like it's not a good candidate match.
 */
const PowerStateGraph::Fingerprint PowerStateGraph::initTraceToEnd(
        const AggregateData::FoundSpike& spike,
        const size_t deviceStart, /**< The possible time the device started. */
        const bool verbose
        )
{
    DisagTree disagTree;

    // make the first vertex (which represents "off")
    DisagTree::vertex_descriptor disagOffVertex = add_vertex(disagTree);
    disagTree[disagOffVertex].timestamp = deviceStart;
    disagTree[disagOffVertex].meanPower = 0; // this is "off"
    disagTree[disagOffVertex].psgVertex = offVertex;

    // add a vertex to represent the first true power state
    DisagTree::vertex_descriptor firstVertex = add_vertex(disagTree);

    // retrieve info for firstVertex and for edge between disagOffVertex and firstVertex
    PSG_out_edge_iter out_e_i, out_e_end;
    PSG_vertex_iter v_i, v_end;
    tie(out_e_i, out_e_end) = out_edges(offVertex, powerStateGraph);
    PowerStateEdge firstEdgeStats = powerStateGraph[*out_e_i];
    tie(v_i, v_end) = vertices(powerStateGraph);
    v_i++;

    disagTree[firstVertex].timestamp = spike.timestamp;
    disagTree[firstVertex].meanPower = powerStateGraph[*v_i].betweenSpikes.mean;
    disagTree[firstVertex].psgVertex = *v_i;
    disagTree[firstVertex].psgEdge   = *out_e_i;
    if (EDGE_HISTORY_SIZE)
        disagTree[firstVertex].edgeHistory.push_back(*out_e_i);

    // add an edge between disagOffVertex and firstVertex
    DisagTree::edge_descriptor edge;
    bool existingEdge;
    tie(edge, existingEdge) = add_edge(
            disagOffVertex,   // source vertex
            firstVertex,      // target vertex
            spike.likelihood, // edge value
            disagTree);

    // now recursively trace from this edge to the end
    traceToEnd( &disagTree, firstVertex, deviceStart );

    if (verbose) {
        write_graphviz(cout, disagTree,
            Disag_vertex_writer(disagTree), Disag_edge_writer(disagTree));
    }

    // find route through the tree with highest average edge likelihoods
    listOfPaths.clear();

    LikelihoodAndVertex nextLAV;
    nextLAV.vertex = firstVertex;
    nextLAV.likelihood = disagTree[edge];

    findListOfPathsThroughDisagTree(
            disagTree,
            disagOffVertex,
            nextLAV);

    // Return the most confident path through the disagTree
    return findBestPath( disagTree, deviceStart );

}

/**
 * @brief Trace the tree downwards from @c vertex recursively finding
 * every path which successfully completes (i.e. reaches an off state).
 */
void PowerStateGraph::findListOfPathsThroughDisagTree(
        const DisagTree& disagTree,
        const DisagTree::vertex_descriptor vertex,
        const LikelihoodAndVertex lav,
        list<PowerStateGraph::LikelihoodAndVertex> path /**< Deliberately called-by-value because we want a copy. */
    )
{
    path.push_back( lav );

    // base case = we're at the end
    if ( vertex != 0 && // check we're not at the first vertex
            disagTree[ vertex ].meanPower == 0 ) {
        listOfPaths.push_back( path );
        return;
    }

    // iterate through each out-edge
    Disag_out_edge_iter out_e_i, out_e_end;
    tie(out_e_i, out_e_end) = out_edges(vertex, disagTree);

    for (; out_e_i!=out_e_end; out_e_i++) {

        DisagTree::vertex_descriptor downstreamVertex = target(*out_e_i, disagTree);

        LikelihoodAndVertex nextLav;
        nextLav.vertex = downstreamVertex;
        nextLav.likelihood = disagTree[*out_e_i];

        // we haven't hit the end yet so recursively follow tree downwards.
        findListOfPathsThroughDisagTree(
                disagTree,
                downstreamVertex,
                nextLav,
                path );

    }
    return;
}

/**
 * Trace from startVertex to the off state in PSGraph
 */
void PowerStateGraph::traceToEnd(
        DisagTree * disagTree_p, /**< input and output parameter */
        const DisagTree::vertex_descriptor& disagVertex,
        const size_t prevTimestamp, /**< timestamp of previous vertex */
        const bool verbose
        ) const
{

    if (verbose) cout << "***traceToEnd... prevTimestamp=" << prevTimestamp << " DisagTree startVertex=" << disagVertex << endl;

    list<AggregateData::FoundSpike> foundSpikes;

    // A handy reference to make the code more readable
    DisagTree& disagTree = *disagTree_p;

    // base case
    if ( disagTree[disagVertex].psgVertex == offVertex ) {
        return;
    }

    // For each out-edge from disagVertex.psgVertex, retrieve a list of
    // spikes which match and create a new DisagTree vertex for each match.
    PSG_out_edge_iter psg_out_i, psg_out_end;
    tie(psg_out_i, psg_out_end) =
            out_edges(disagTree[disagVertex].psgVertex, powerStateGraph);

    const size_t WINDOW_FRAME = 8; // number of seconds to widen window by

    for (; psg_out_i!=psg_out_end; psg_out_i++ ) {

        // Commented out alternative strategy for getting edge history. Almost certainly slower than
        // using disagTree[disagVertex].edgeHistory
//        std::list< PSGraph::edge_descriptor > eHistory = getEdgeHistoryForVertex(disagTree, disagVertex);

        if ( ! edgeListsAreEqual(
//                eHistory,
                disagTree[disagVertex].edgeHistory,
                powerStateGraph[*psg_out_i].edgeHistory ) ) {
            if (verbose) cout << "edge histories not equal" << endl;
            continue;
        }

        size_t begOfSearchWindow, endOfSearchWindow;

        size_t e = powerStateGraph[*psg_out_i].duration.nonZeroStdev();

        if (verbose)  cout << "disagTree[startVertex].timestamp=" << disagTree[disagVertex].timestamp << endl;

        begOfSearchWindow = (disagTree[disagVertex].timestamp + powerStateGraph[*psg_out_i].duration.min)
                - WINDOW_FRAME - e;

        endOfSearchWindow = (disagTree[disagVertex].timestamp + powerStateGraph[*psg_out_i].duration.max)
                + WINDOW_FRAME + e;

        //********************** DIAGNOSTICS *************************
        if (verbose) {
            cout << "endOfSearchWindow=" << endOfSearchWindow
                    << " disagTree[startVertex].timestamp=" << disagTree[disagVertex].timestamp
                    << " powerStateGraph[*psg_out_i].duration.max=" << powerStateGraph[*psg_out_i].duration.max
                    << " powerStateGraph[*psg_out_i].duration.mean=" << powerStateGraph[*psg_out_i].duration.mean
                    << " powerStateGraph[*psg_out_i].duration.stdev=" << powerStateGraph[*psg_out_i].duration.stdev
                    << endl;
        }//__________________________________________________________

        // ensure we're not looking backwards in time
        if (begOfSearchWindow <= disagTree[disagVertex].timestamp) {
            begOfSearchWindow =
                    disagTree[disagVertex].timestamp
                    + powerStateGraph[*psg_out_i].duration.min
                    - (powerStateGraph[*psg_out_i].duration.min / 10);
        }

        // check that we're not looking past the end of aggData
        if (endOfSearchWindow > (*aggData)[ (*aggData).getSize() - 1 ].timestamp ) {
            continue; // we can't process this if we're trying to look past the end of the aggData
        }

        //************************************************************//
        // get a list of candidate spikes matching this PSG-out-edge  //

        foundSpikes.clear();
        foundSpikes = aggData->findSpike(
                powerStateGraph[*psg_out_i].delta,  // spike stats
                begOfSearchWindow,
                endOfSearchWindow
                );

        //***************************************************************//
        // for each candidate spike, create a new vertex in disagTree   //
        // and recursively trace this to the end                         //

        for (list<AggregateData::FoundSpike>::const_iterator spike=foundSpikes.begin();
                spike!=foundSpikes.end();
                spike++) {

            // ensure that the absolute value of the aggregate data signal
            // does not drop below the minimum for our current power state
            if (aggData->readingGoesBelowPowerState(
                    disagTree[disagVertex].timestamp,
                    spike->timestamp,
                    powerStateGraph[ disagTree[disagVertex].psgVertex ].betweenSpikes ) ) {

                continue;
            }

            double normalisedLikelihoodForTime = powerStateGraph[*psg_out_i].duration.normalisedLikelihood(
                    spike->timestamp - disagTree[disagVertex].timestamp);

            // merge probability for time and for spike delta
            double avLikelihood = ( normalisedLikelihoodForTime + spike->likelihood ) / 2;

            // create new vertex
            DisagTree::vertex_descriptor newVertex=add_vertex( disagTree );

            // create new edge
            DisagTree::edge_descriptor newEdge;
            bool existingEdge;
            tie(newEdge, existingEdge) =
                    add_edge( disagVertex, newVertex, avLikelihood, disagTree );

            // add details to newVertex
            disagTree[newVertex].timestamp = spike->timestamp;
            // get vertex that *psg_out_i points to
            disagTree[newVertex].psgVertex = target(*psg_out_i, powerStateGraph);
            disagTree[newVertex].psgEdge   = *psg_out_i;
            disagTree[newVertex].meanPower =
                    powerStateGraph[disagTree[newVertex].psgVertex].betweenSpikes.mean;

            if (EDGE_HISTORY_SIZE) {
                disagTree[newVertex].edgeHistory = disagTree[disagVertex].edgeHistory;
                disagTree[newVertex].edgeHistory.push_back(*psg_out_i);
                if (disagTree[newVertex].edgeHistory.size() > EDGE_HISTORY_SIZE) {
                    disagTree[newVertex].edgeHistory.erase( disagTree[newVertex].edgeHistory.begin() );
                }
            }

            // recursively trace to end.
            traceToEnd(disagTree_p, newVertex, disagTree[disagVertex].timestamp);
        }
    }
}

/**
 * @brief Iterates through each path in @c listOfPaths
 *        to find the one with the highest likelihood.
 *
 * To be called after @c listOfPaths has been populated by findListOfPathsThroughDisagTree().
 *
 * @return details of the best path.  Confidence is set to -1 if no paths are available.
 */
const PowerStateGraph::Fingerprint PowerStateGraph::findBestPath(
        const DisagTree& disagTree,
        const size_t deviceStart,
        const bool verbose
        )
{
    Fingerprint fingerprint;
    fingerprint.timestamp = deviceStart;

    // first check that listOfPaths is populated
    if ( listOfPaths.empty() ) {
        fingerprint.avLikelihood = -1; // return error code
        return fingerprint;
    }

    fingerprint.avLikelihood = 0;

    list< list<LikelihoodAndVertex> >::const_iterator path_i, bestPath_i;
    list< LikelihoodAndVertex >::const_iterator lav_i;

    if (verbose) cout << "path dump:" << endl;

    double likelihoodAccumulator, avLikelihood;
    bool foundGoodPath = false;

    // iterate through each path to find the one with the highest likelihood
    for (path_i=listOfPaths.begin(); path_i!=listOfPaths.end(); path_i++) {

        likelihoodAccumulator = 0;
        for ( lav_i=path_i->begin(); lav_i!=path_i->end(); lav_i++ ) {
            if (verbose) cout << "vertex=" << lav_i->vertex << " conf=" << lav_i->likelihood << ", ";
            likelihoodAccumulator += lav_i->likelihood;
        }

        avLikelihood = likelihoodAccumulator / path_i->size();

        // if this is the most confident path we've seen yet then record its details.
        if ( avLikelihood >= fingerprint.avLikelihood ) {
            fingerprint.avLikelihood = avLikelihood;
            bestPath_i = path_i;
            foundGoodPath = true;
        }

        if (verbose) cout << endl << endl;
    }

    // now get energy usage and duration from bestPath_i
    size_t prevTimestamp, duration;
    double prevMeanPower = 0;

    fingerprint.energy = 0;
    fingerprint.duration = 0;
    duration = 0;
    prevTimestamp = deviceStart;
    size_t count = 0;
    if (foundGoodPath) {
        for (lav_i=bestPath_i->begin(); lav_i != bestPath_i->end(); lav_i++) {

            // Add to timeAndPower list for read-out later when we plot the power states
            if ( lav_i != bestPath_i->begin() ) {
                fingerprint.timeAndPower.push_back( TimeAndPower(
                        disagTree[ lav_i->vertex ].timestamp-1 ,
                        (count==1 ? 0 : prevMeanPower)
                ) );
                fingerprint.timeAndPower.push_back( TimeAndPower(
                        disagTree[ lav_i->vertex ].timestamp,
                        disagTree[ lav_i->vertex ].meanPower) );
            }

            // Now calculate duration and energy consumption
            duration = disagTree[ lav_i->vertex ].timestamp - prevTimestamp;
            fingerprint.energy += prevMeanPower * duration;

            fingerprint.duration += duration;

            prevTimestamp = disagTree[ lav_i->vertex ].timestamp;
            prevMeanPower = disagTree[ lav_i->vertex ].meanPower;
            count++;
        }

        fingerprint.avLikelihood = ( fingerprint.avLikelihood +
                energyConsumption.normalisedLikelihood(fingerprint.energy) ) / 2;
    } else {
        fingerprint.avLikelihood = -1;
    }

    return fingerprint;
}

/**
 * @brief Trace the disagTree backwards.
 */
list< PowerStateGraph::PSGraph::edge_descriptor > PowerStateGraph::getEdgeHistoryForVertex(
        const DisagTree& disagTree,
        const DisagTree::vertex_descriptor& startVertex
        ) const
{
    list< PSGraph::edge_descriptor > eHistory;
    DisagTree::vertex_descriptor disagVertex = startVertex;

    DisagTree::in_edge_iterator in_e_i, in_e_end;

    while (disagTree[disagVertex].psgVertex != offVertex && eHistory.size() < EDGE_HISTORY_SIZE) {
        eHistory.push_front( disagTree[disagVertex].psgEdge );
        tie(in_e_i, in_e_end) = in_edges(disagVertex, disagTree);
        disagVertex = source( *in_e_i, disagTree ); // the vertex upstream from in_e_i
    }

    return eHistory;
}

void PowerStateGraph::setDeviceName(const string& _deviceName)
{
    deviceName = _deviceName;
}

std::ostream& operator<<( std::ostream& o, const PowerStateGraph& psg )
{
    PowerStateGraph::PSG_vertex_index_map index = boost::get(boost::vertex_index, psg.powerStateGraph);

    o << "vertices(graph) = " << std::endl;
    std::pair<PowerStateGraph::PSG_vertex_iter, PowerStateGraph::PSG_vertex_iter> vp;
    for (vp = boost::vertices(psg.powerStateGraph); vp.first != vp.second; ++vp.first) {
        o << "vertex" << index[*vp.first] << " = {" << psg.powerStateGraph[*vp.first].postSpike <<  "}";
        if ( *vp.first == psg.offVertex ) {
            o << "\n          (offVertex)" << std::endl;
        }
        o << std::endl;
    }

    return o;
}


text:
#include "WindingNumbers.h"
#include "SPlisHSPlasH/Common.h"
#include <Eigen/Dense>

#define _USE_MATH_DEFINES
#include "math.h"

using namespace Utilities;
using namespace SPH;
using namespace Eigen;

Real WindingNumbers::computeGeneralizedWindingNumber(const Vector3r& p_, const Vector3r& a_, const Vector3r& b_, const Vector3r& c_)
{
 	const Vector3r& a = a_ - p_;
 	const Vector3r& b = b_ - p_;
 	const Vector3r& c = c_ - p_;

 	const Real normA = a.norm();
 	const Real normB = b.norm();
 	const Real normC = c.norm();

	Matrix3r A;
	A.row(0) = a;
	A.row(1) = b;
	A.row(2) = c;
 	const Real det = A.determinant();
 	const Real divisor = normA*normB*normC + (a.dot(b))*normC + (b.dot(c))*normA + (c.dot(a))*normB;

 	static const Real tau = 2.0 *  M_PI;
 	return std::atan2(det,divisor) / tau; // Only divide by 2*pi instead of 4*pi because there was a 2 out front
}

Real WindingNumbers::computeGeneralizedWindingNumber(const Vector3r& p, const TriangleMesh& mesh)
{
 	const unsigned int *faces = mesh.getFaces().data();
 	const unsigned int nFaces = mesh.numFaces();
	const Vector3r *v = mesh.getVertices().data();

 	// compute generalized winding number
 	Real w_p = 0;
 	#pragma omp parallel for reduction (+: w_p)
 	for (int idx = 0; idx < static_cast<int>(nFaces); idx ++)
 	{
 		const Vector3r& a = v[faces[idx*3+0]];
 		const Vector3r& b = v[faces[idx*3+1]];
 		const Vector3r& c = v[faces[idx*3+2]];

 		w_p += computeGeneralizedWindingNumber(p, a, b, c);
 	}

 	return w_p;
}

ext: cpp


text:
//  (C) Copyright Edward Diener 2011-2015
//  Use, modification and distribution are 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).

#include <boost/vmd/elem.hpp>
#include <boost/vmd/equal.hpp>
#include <boost/vmd/is_empty.hpp>
#include <boost/detail/lightweight_test.hpp>
#include <boost/preprocessor/tuple/elem.hpp>

int main()
  {

#if BOOST_PP_VARIADICS

  #define BOOST_VMD_REGISTER_ggh (ggh)

  #define ANIDENTIFIER ggh
  #define ANUMBER 249
  #define ASEQ (25)(26)(27)
  #define ATYPE BOOST_VMD_TYPE_SEQ
  #define ATYPE2 BOOST_VMD_TYPE_ARRAY
  #define ATYPE3 BOOST_VMD_TYPE_TYPE
  #define ATYPE5 BOOST_VMD_TYPE_EMPTY
  #define ALIST (0,(1,(2,(3,BOOST_PP_NIL))))
  #define ANARRAY (3,(a,b,38))
  #define ANARRAY2 (5,(c,d,133,22,15))
  #define ASEQUENCE ANUMBER ALIST ATYPE ANIDENTIFIER ANARRAY ASEQ
  #define ASEQUENCE2 ANIDENTIFIER ANARRAY2 ALIST ASEQ ATYPE2
  #define ASEQUENCE3 ASEQ ANIDENTIFIER ATYPE3 ALIST
  #define ASEQUENCE4
  #define ASEQUENCE5 ALIST ASEQ ATYPE5 ANIDENTIFIER

  BOOST_TEST(BOOST_VMD_IS_EMPTY(BOOST_PP_TUPLE_ELEM(0,BOOST_VMD_ELEM(0,ASEQUENCE,BOOST_VMD_RETURN_AFTER,BOOST_VMD_TYPE_TYPE))));
  BOOST_TEST(BOOST_VMD_IS_EMPTY(BOOST_PP_TUPLE_ELEM(1,BOOST_VMD_ELEM(0,ASEQUENCE2,BOOST_VMD_RETURN_AFTER,BOOST_VMD_TYPE_TYPE))));
  BOOST_TEST(BOOST_VMD_IS_EMPTY(BOOST_PP_TUPLE_ELEM(0,BOOST_VMD_ELEM(0,ASEQUENCE3,BOOST_VMD_RETURN_AFTER,BOOST_VMD_TYPE_TYPE))));

  BOOST_TEST(!BOOST_VMD_IS_EMPTY(BOOST_PP_TUPLE_ELEM(1,BOOST_VMD_ELEM(2,ASEQUENCE,BOOST_VMD_RETURN_AFTER,BOOST_VMD_TYPE_TYPE))));
  BOOST_TEST(BOOST_VMD_EQUAL(BOOST_PP_TUPLE_ELEM(0,BOOST_VMD_ELEM(2,ASEQUENCE,BOOST_VMD_RETURN_AFTER,BOOST_VMD_TYPE_TYPE)),BOOST_VMD_TYPE_SEQ,BOOST_VMD_TYPE_TYPE));
  BOOST_TEST(BOOST_VMD_EQUAL(BOOST_VMD_ELEM(2,ASEQUENCE,BOOST_VMD_TYPE_TYPE),BOOST_VMD_TYPE_SEQ,BOOST_VMD_TYPE_TYPE));

  BOOST_TEST(BOOST_VMD_IS_EMPTY(BOOST_PP_TUPLE_ELEM(1,BOOST_VMD_ELEM(4,ASEQUENCE2,BOOST_VMD_RETURN_AFTER,BOOST_VMD_TYPE_TYPE))));
  BOOST_TEST(BOOST_VMD_EQUAL(BOOST_VMD_TYPE_ARRAY,BOOST_VMD_ELEM(4,ASEQUENCE2,BOOST_VMD_TYPE_TYPE),BOOST_VMD_TYPE_TYPE));

  BOOST_TEST(BOOST_VMD_EQUAL(BOOST_VMD_TYPE_TYPE,BOOST_PP_TUPLE_ELEM(0,BOOST_VMD_ELEM(2,ASEQUENCE3,BOOST_VMD_RETURN_AFTER,BOOST_VMD_TYPE_TYPE)),BOOST_VMD_TYPE_TYPE));
  BOOST_TEST(!BOOST_VMD_IS_EMPTY(BOOST_PP_TUPLE_ELEM(1,BOOST_VMD_ELEM(2,ASEQUENCE3,BOOST_VMD_RETURN_AFTER,BOOST_VMD_TYPE_TYPE))));
  BOOST_TEST(BOOST_VMD_EQUAL(BOOST_VMD_TYPE_TYPE,BOOST_VMD_ELEM(2,ASEQUENCE3,BOOST_VMD_TYPE_TYPE),BOOST_VMD_TYPE_TYPE));

  BOOST_TEST(BOOST_VMD_IS_EMPTY(BOOST_PP_TUPLE_ELEM(0,BOOST_VMD_ELEM(0,ASEQUENCE4,BOOST_VMD_RETURN_AFTER,BOOST_VMD_TYPE_TYPE))));
  BOOST_TEST(BOOST_VMD_IS_EMPTY(BOOST_PP_TUPLE_ELEM(1,BOOST_VMD_ELEM(0,ASEQUENCE4,BOOST_VMD_RETURN_AFTER,BOOST_VMD_TYPE_TYPE))));
  BOOST_TEST(BOOST_VMD_IS_EMPTY(BOOST_VMD_ELEM(0,ASEQUENCE4,BOOST_VMD_TYPE_TYPE)));

  BOOST_TEST(!BOOST_VMD_IS_EMPTY(BOOST_PP_TUPLE_ELEM(1,BOOST_VMD_ELEM(2,ASEQUENCE5,BOOST_VMD_RETURN_AFTER,BOOST_VMD_TYPE_TYPE))));
  BOOST_TEST(BOOST_VMD_EQUAL(BOOST_VMD_TYPE_EMPTY,BOOST_VMD_ELEM(2,ASEQUENCE5,BOOST_VMD_TYPE_TYPE),BOOST_VMD_TYPE_TYPE));

#else

BOOST_ERROR("No variadic macro support");

#endif

  return boost::report_errors();

  }

ext: cxx


text:
// Copyright (C) 2001 Vladimir Prus <ghost@cs.msu.su>
// Copyright (C) 2001 Jeremy Siek <jsiek@cs.indiana.edu>
// Permission to copy, use, modify, sell and distribute this software is
// granted, provided this copyright notice appears in all copies and
// modified version are clearly marked as such. This software is provided
// "as is" without express or implied warranty, and with no claim as to its
// suitability for any purpose.

#include <boost/graph/vector_as_graph.hpp>
#include <boost/graph/transitive_closure.hpp>

#include <iostream>

#include <boost/graph/vector_as_graph.hpp>
#include <boost/graph/adjacency_list.hpp>
#include <boost/graph/graph_utility.hpp>

#include <cstdlib>
#include <ctime>
#include <boost/progress.hpp>
using namespace std;
using namespace boost;

void generate_graph(int n, double p, vector< vector<int> >& r1)
{
  static class {
  public:
    double operator()() {
      return double(rand())/RAND_MAX;
    }
  } gen;
  r1.clear();
  r1.resize(n);
  for (int i = 0; i < n; ++i)
    for (int j = 0; j < n; ++j)
      if (gen() < p)
        r1[i].push_back(j);
}

// (i,j) is in E' if j is reachable from i
// Hmm, is_reachable does not detect when there is a non-trivial path
// from i to i. It always returns true for is_reachable(i,i).
// This needs to be fixed/worked around.
template <typename Graph, typename GraphTC>
bool check_transitive_closure(Graph& g, GraphTC& tc)
{
  std::vector<default_color_type> color_map_vec(num_vertices(g));
  typename graph_traits<GraphTC>::vertex_iterator i, i_end;
  for (tie(i, i_end) = vertices(tc); i != i_end; ++i) {
    typename graph_traits<GraphTC>::out_edge_iterator j, j_end;
    for (tie(j, j_end) = out_edges(*i, tc); j != j_end; ++j)
      if (!is_reachable(source(*j, g), target(*j, g), g, &color_map_vec[0]))
        return false;
  }
  return true;
}

bool test(int n, double p)
{
  vector< vector<int> > g1, g1_tc;
  generate_graph(n, p, g1);
  cout << "Created graph with " << n << " vertices.\n";

  vector< vector<int> > g1_c(g1);

  {
    progress_timer t;
    cout << "transitive_closure" << endl;
    transitive_closure(g1, g1_tc, vertex_index_map(identity_property_map()));
  }

  if(check_transitive_closure(g1, g1_tc))
    return true;
  else {
    //cout << "Original graph was " << multiline << g1_c << endl;
    //cout << "Result is " << multiline << g1 << endl;
    return false;
  }
}

int main()
{
  srand(time(0));
  static class {
  public:
    double operator()() {
      return double(rand())/RAND_MAX;
    }
  } gen;

  for (size_t i = 0; i < 100; ++i) {
    int n = 0 + int(20*gen());
    double p = gen();
    if (!test(n, p)) {
      cout << "Failed." << endl;
      return 1;
    }
  }
  cout << "Passed." << endl;
}


text:
#include "Util.h"
//#include <Eigen/Dense>
#include "cryptoTools/Common/Log.h"
#include "cryptoTools/Common/Matrix.h"
#include "cryptoTools/Crypto/PRNG.h"
#include "cryptoTools/Common/BitIterator.h"
#include <numeric>
#include <iomanip>
#include <future>
#include <cmath>

#include "Mtx.h"
#include <deque>
#include "LdpcSampler.h"

#ifdef ENABLE_ALGO994
extern "C" {
#include "libOTe/Tools/LDPC/Algo994/utils.h"
#include "libOTe/Tools/LDPC/Algo994/generations.h"
#include "libOTe/Tools/LDPC/Algo994/data_defs.h"
}
#endif

#include <fstream>

namespace osuCrypto
{

    //using Mtx = Eigen::Matrix<int, Eigen::Dynamic, Eigen::Dynamic>;

#ifdef ENABLE_ALGO994
    int alg994 = ALG_SAVED_UNROLLED;
    int num_saved_generators =  5;
    int num_cores = 4;
    int num_permutations = 5;
    int print_matrices = 0;
#endif

    int minDist(std::string path, u64 numTHreads, bool verbose)
    {
#ifdef ENABLE_ALGO994

        char* inputMatrix;
        int   info, k, n, dist;
        num_cores = numTHreads;

        // Read input matrix.
        info = read_char_matrix((char*)path.c_str(), &inputMatrix, &k, &n);
        if (info != 0) {
            fprintf(stderr, "\n");
            fprintf(stderr, "ERROR in read_char_matrix: ");
            fprintf(stderr, "Error while reading input file with matrix.\n\n\n");
        }

        assert(
            alg994 == ALG_BASIC || // 1
            alg994 == ALG_OPTIMIZED || // 2
            alg994 == ALG_STACK || // 3
            alg994 == ALG_SAVED || // 4
            alg994 == ALG_SAVED_UNROLLED || // 5
            alg994 == ALG_GRAY  // 10
        );

        // Compute distance of input matrix.
        dist = compute_distance_of_matrix_impl(inputMatrix, k, n,
            alg994,
            num_saved_generators,
            num_cores,
            num_permutations,
            print_matrices);

        // Remove matrices.
        free(inputMatrix);

        return dist;
#else
        throw std::runtime_error("also 994 not enabled. " LOCATION);
#endif
    }

    int minDist2(const DenseMtx& mtx, u64 nt, bool verbose)
    {
        std::string outPath("./deleteMe");
        std::fstream out(outPath, std::fstream::trunc | std::fstream::out);

        if (out.is_open() == false)
        {
            std::cout << "failed to open: " << outPath << std::endl;
            return 0;
        }

        std::vector<std::pair<u64, u64>> swaps;
        auto G = computeGen(mtx, swaps);

        out << G.rows() << " " << G.cols() << " matrix dimensions\n"
            << G << std::endl;

        out.close();

        u64 d;

        d = minDist(outPath, nt, false);

        std::remove(outPath.c_str());
        return (int)d;
    }

    u64 numNonzeroRows(const DenseMtx& mtx)
    {
        u64 r = 0;
        for (u64 i = 0; i < mtx.rows(); ++i)
        {
            if (mtx.row(i).isZero())
                ++r;
        }

        return mtx.rows() - r;
    }

    //u64 rank(const DenseMtx& mtx)
    //{
    //    DenseMtx m2 = gaussianElim(mtx);
    //    return numNonzeroRows(m2);
    //}

    void ithCombination(u64 index, u64 n, std::vector<u64>& set)
    {
        //'''Yields the items of the single combination that would be at the provided
        //(0-based) index in a lexicographically sorted list of combinations of choices
        //of k items from n items [0,n), given the combinations were sorted in
        //descending order. Yields in descending order.
        //'''
        u64 nCk = 1;
        u64 nMinusI = n;
        u64 iPlus1 = 1;

        auto k = set.size();

        // nMinusI, iPlus1 in zip(range(n, n - k, -1), range(1, k + 1)):
        for (; nMinusI != n - k; --nMinusI, ++iPlus1)
        {
            nCk *= nMinusI;
            nCk /= iPlus1;
        }

        //std::cout << "nCk " << nCk << std::endl;

        auto curIndex = nCk;
        for (auto kk = k; kk != 0ull; --kk)//in range(k, 0, -1):
        {
            //std::cout << "kk " << kk << " " <<  nCk << std::endl;
            nCk *= kk;
            nCk /= n;
            while (curIndex - nCk > index) {
                curIndex -= nCk;
                nCk *= (n - kk);
                nCk -= nCk % kk;
                n -= 1;
                nCk /= n;
            }
            n -= 1;

            set[kk - 1] = n;
        }
    }

    std::vector<u64> ithCombination(u64 index, u64 n, u64 k)
    {
        std::vector<u64> set(k);
        ithCombination(index, n, set);
        return set;
    }

    u64 choose(u64 n, u64 k)
    {
        if (k == 0) return 1;
        return (n * choose(n - 1, k - 1)) / k;
    }

    struct NChooseK
    {
        u64 mN;
        u64 mK;
        u64 mI, mEnd;
        std::vector<u64> mSet;

        NChooseK(u64 n, u64 k, u64 begin = 0, u64 end = -1)
            : mN(n)
            , mK(k)
            , mI(begin)
            , mEnd(std::min<u64>(choose(n, k), end))
        {
            assert(k <= n);
            mSet = ithCombination(begin, n, k);
        }

        const std::vector<u64>& operator*() const
        {
            return mSet;
        }

        void operator++()
        {
            ++mI;
            assert(mI <= mEnd);

            u64 i = 0;
            while (i < mK - 1 && mSet[i] + 1 == mSet[i + 1])
                ++i;

            //if (i == mK - 1 && mSet.back() == mN - 1)
            //{
            //    mSet.clear();
            //    return;
            //    //assert(mSet.back() != mN - 1);
            //}

            ++mSet[i];
            for (u64 j = 0; j < i; ++j)
                mSet[j] = j;
        }

        explicit operator bool() const
        {
            return mI < mEnd;
        }

    };

    bool isZero(const span<block>& sum)
    {
        for (auto& b : sum)
            if (b != oc::ZeroBlock)
            {
                return false;
                break;
            }
        return true;
    }
    bool isEq(const span<const block>& u, const span<const block>& v)
    {
        assert(u.size() == v.size());
        return memcmp(u.data(), v.data(), v.size_bytes()) == 0;
    }

    template <typename T>
    class queue
    {
    private:
        std::mutex              d_mutex;
        std::condition_variable d_condition;
        std::deque<T>           d_queue;
    public:
        void push(T const& value) {
            {
                std::unique_lock<std::mutex> lock(this->d_mutex);
                d_queue.push_front(value);
            }
            this->d_condition.notify_one();
        }
        T pop() {
            std::unique_lock<std::mutex> lock(this->d_mutex);
            this->d_condition.wait(lock, [=] { return !this->d_queue.empty(); });
            T rc(std::move(this->d_queue.back()));
            this->d_queue.pop_back();
            return rc;
        }
    };

    std::pair<double, std::vector<u64>> minDist(const DenseMtx& mtx, bool verbose, u64 numThreads)
    {
        assert(mtx.rows() < mtx.cols());

        u64 percision = 2;
        u64 p = (u64)std::pow(10, percision);

        std::mutex mut;
        queue<std::function<void()>> queue;
        std::vector<std::thread> thrds(numThreads);
        for (u64 i = 0; i < thrds.size(); ++i)
        {
            thrds[i] = std::thread([&queue]() {

                while (true)
                {
                    std::function<void()> fn =queue.pop();

                    if (!fn)
                        return;
                    fn();
                }
                });
        }

        u64 dd = mtx.mData.cols();
#define ASSUME_DD_1
#ifdef ASSUME_DD_1
        if (dd != 1)
            throw RTE_LOC;
#endif
        for (u64 weight = 2; weight < mtx.rows(); ++weight)
        {
            auto total = choose(mtx.cols(), weight);
            u64 next = 0;
            std::atomic<u64> ii(0);
            //std::atomic<u64> rem = numThreads;

            bool done = false;
            using Ret = std::pair<double, std::vector<u64>>;
            std::vector<std::promise<Ret>> prom(numThreads);
            //auto fu = prom.get_future();

            for (u64 i = 0; i < numThreads; ++i)
            {
                queue.push([&, i, weight]() {

                    auto begin = i * total / numThreads;
                    auto end = (i + 1) * total / numThreads;
                    auto iter = NChooseK(mtx.cols(), weight, begin, end);
                    //u64& mI = iter.mI;
                    std::vector<u64> set;
#ifdef ASSUME_DD_1
                    block sum;
#else
                    std::vector<block> sum(dd);
#endif
                    while (begin++ != end && !done)
                    {
                        set = *iter;

                        if (verbose && ii >= next)
                        {
                            std::lock_guard<std::mutex> lock(mut);
                            if (verbose && ii >= next)
                            {
                                auto cur = ii * p / total;
                                next = (cur + 1) * total / p;;
                                std::cout << "\r" << weight << "." << std::setw(percision) << std::setfill('0') << cur << "       " << std::flush;
                            }
                        }
#ifdef ASSUME_DD_1
                        auto ptr = mtx.mData.data();
                        sum = ptr[set[0]];
                        for (u64 i = 1; i < weight; ++i)
                        {
                            auto col = ptr[set[i]];
                            sum = sum ^ col;
                        }

                        auto linDep = sum == oc::ZeroBlock;
#else

                        auto v = mtx.col(set[0]);
                        std::copy(v.data(), v.data() + dd, sum.data());

                        for (u64 i = 1; i < weight; ++i)
                        {
                            auto col = mtx.col(set[i]);
                            for (u64 j = 0; j < dd; ++j)
                                sum[j] = sum[j] ^ col[j];
                        }

                        auto linDep = isZero(sum);
#endif
                        if (linDep)
                        {
                            std::lock_guard<std::mutex> lock(mut);
                            if (verbose)
                            {
                                std::cout << std::endl;
                            }
                            done = true;

                            prom[i].set_value(std::make_pair(weight + ii / double(total), set));
                            return;
                        }

                        ++ii;

                        //++iter;
                        {
                            //++(mI);
                            //assert(mI <= iter.mEnd);

                            u64 i = 0;
                            while (i < iter.mK - 1 && iter.mSet[i] + 1 == iter.mSet[i + 1])
                                ++i;

                            ++iter.mSet[i];
                            for (u64 j = 0; j < i; ++j)
                                iter.mSet[j] = j;
                        }
                    }

                    prom[i].set_value({});
                    return;
                    });
            }

            Ret ret;
            for (u64 i = 0; i < numThreads; ++i)
            {
                auto cc = prom[i].get_future().get();
                if (cc.second.size())
                    ret = cc;
            }

            if (done)
            {
                for (u64 i = 0; i < numThreads; ++i)
                    queue.push({});
                for (u64 i = 0; i < numThreads; ++i)
                    thrds[i].join();
                return ret;
            }

        }
        assert(0);
        return {};
    }

    std::pair<double, std::vector<u64>> minDist(const DenseMtx& mtx, bool verbose)
    {
        assert(mtx.rows() < mtx.cols());

        u64 percision = 2;
        u64 p = (u64) pow(10, percision);

        for (u64 weight = 2; weight < mtx.rows(); ++weight)
        {
            auto iter = NChooseK(mtx.cols(), weight);
            auto total = choose(mtx.cols(), weight);
            u64 prev = -1;
            u64 ii = 0;
            std::vector<block> sum(mtx.mData.cols());

            while (iter)
            {
                auto& set = *iter;

                auto cur = ii * p / total;
                if (verbose && prev != cur)
                {
                    prev = cur;
                    std::cout << "\r" << weight << "." << std::setw(percision) << std::setfill('0') << cur << "       " << std::flush;
                }

                auto v = mtx.col(set[0]);
                std::copy(v.begin(), v.end(), sum.begin());

                for (u64 i = 1; i < set.size(); ++i)
                {
                    auto col = mtx.col(set[i]);
                    for (u64 j = 0; j < sum.size(); ++j)
                        sum[j] = sum[j] ^ col[j];
                }

                if (isZero(sum))
                {
                    if (verbose)
                        std::cout << std::endl;
                    return std::make_pair(weight + ii / double(total), set);
                }

                ++ii;
                ++iter;
            }
        }
        assert(0);
        return {};
    }

    DenseMtx computeGen(DenseMtx& H)
    {
        assert(H.rows() < H.cols());

        auto n = H.cols();
        auto m = H.rows();
        auto k = n - m;

        auto mtx = H.subMatrix(0, k, m, m);

        auto P = H.subMatrix(0, 0, m, k);

        for (u64 i = 0; i < m; ++i)
        {
            if (mtx(i, i) == 0)
            {
                for (u64 j = i + 1; j < m; ++j)
                {
                    if (mtx(j, i) == 1)
                    {
                        mtx.row(i).swap(mtx.row(j));
                        P.row(i).swap(P.row(j));
                        break;
                    }
                }

                if (mtx(i, i) == 0)
                {
                    //std::cout << mtx << std::endl;
                    return {};
                }
            }

            for (u64 j = 0; j < m; ++j)
            {
                if (j != i && mtx(j, i))
                {
                    for (u64 l = 0; l < m; ++l)
                    {
                        mtx(j, l) ^= mtx(i, l);
                    }

                    for (u64 l = 0; l < k; ++l)
                        P(j, l) ^= P(i, l);
                }
            }

        }

        DenseMtx G(k, n);
        for (u64 i = 0; i < k; ++i)
            G(i, i) = 1;

        for (u64 i = 0; i < m; ++i)
        {
            for (u64 j = 0; j < k; ++j)
            {
                G(j, i + k) = P(i, j);
            }
        }

        return G;
    }

    DenseMtx computeGen(DenseMtx H, std::vector<std::pair<u64, u64>>& colSwaps)
    {
        assert(H.rows() < H.cols());

        auto n = H.cols();
        auto m = H.rows();
        auto k = n - m;
        colSwaps.clear();

        for (u64 row = 0, col = k; row < m; ++row, ++col)
        {
            //std::cout << row << std::endl << H << std::endl;;

            if (H(row, col) == 0)
            {
                bool found = false;
                // look fow a row swap
                for (u64 row2 = row + 1; row2 < m && found == false; ++row2)
                {
                    if (H(row2, col) == 1)
                    {
                        H.row(row).swap(H.row(row2));
                        found = true;
                    }
                }

                if (found == false)
                {
                    // look for a col swap

                    for (u64 col2 = 0; col2 < k && found == false; ++col2)
                    {
                        for (u64 row2 = row; row2 < m && found == false; ++row2)
                        {
                            if (H(row2, col2) == 1)
                            {
                                H.row(row).swap(H.row(row2));

                                // swap columns.
                                colSwaps.push_back({ col,col2 });
                                auto c0 = H.col(col);
                                auto c1 = H.col(col2);
                                std::swap_ranges(c0.begin(), c0.end(), c1.begin());
                                found = true;
                            }
                        }
                    }
                }

                if (found == false)
                {
                    // can not be put in systematic form.
                    //std::cout  << H << std::endl;

                    return {};
                }
            }

            // clear all other ones from the current column.
            for (u64 row2 = 0; row2 < m; ++row2)
            {
                if (row2 != row && H(row2, col))
                {
                    // row2 = row ^ row2
                    for (u64 col2 = 0; col2 < n; ++col2)
                    {
                        H(row2, col2) ^= H(row, col2);
                    }
                }
            }

        }

        auto P = H.subMatrix(0, 0, m, k);

        DenseMtx G(k, n);
        for (u64 i = 0; i < k; ++i)
            G(i, i) = 1;

        for (u64 i = 0; i < m; ++i)
        {
            for (u64 j = 0; j < k; ++j)
            {
                G(j, i + k) = P(i, j);
            }
        }

        return G;
    }

    namespace tests
    {
        void computeGen_test(const oc::CLP& cmd)
        {
            u64 n = 30;
            u64 k = n / 4;
            u64 m = n - k;
            u64 t = 100;

            oc::PRNG prng(block(cmd.getOr("s", 0), 0));

            DenseMtx W(t, k);

            for (u64 i = 0; i < W.rows(); ++i)
                for (u64 j = 0; j < W.cols(); ++j)
                    W(i, j) = prng.getBit();

            DenseMtx H(m, n), G;

            while (G.rows() == 0)
            {
                for (u64 i = 0; i < H.rows(); ++i)
                {
                    //H(i, k + i) = 1;

                    for (u64 j = 0; j < H.cols(); ++j)
                        H(i, j) = prng.getBit();
                }

                G = computeGen(H);
            }

            //std::cout << H << std::endl;
            //std::cout << G << std::endl;

            auto C = W * G;

            auto S = H * C.transpose();

            for (u64 i = 0; i < S.rows(); ++i)
                for (u64 j = 0; j < S.cols(); ++j)
                    assert(S(i, j) == 0);
        }

        void computeGen_test2(const oc::CLP& cmd)
        {
            u64 n = 100;
            u64 k = n / 4;
            u64 m = n - k;
            u64 t = 100;

            oc::PRNG prng(block(cmd.getOr("s", 0), 0));

            DenseMtx W(t, k);

            for (u64 i = 0; i < W.rows(); ++i)
                for (u64 j = 0; j < W.cols(); ++j)
                    W(i, j) = prng.getBit();

            DenseMtx H(m, n), G;
            std::vector<std::pair<u64, u64>> swaps;

            while (G.rows() == 0)
            {
                for (u64 i = 0; i < H.rows(); ++i)
                {
                    //H(i, k + i) = 1;

                    for (u64 j = 0; j < H.cols(); ++j)
                        H(i, j) = prng.getBit();
                }

                G = computeGen(H, swaps);
            }

            //std::cout << H << std::endl;
            //std::cout << G << std::endl;

            auto C = W * G;

            for(auto s : swaps)
            {
                auto c0 = C.col(s.first);
                auto c1 = C.col(s.second);
                std::swap_ranges(c0.begin(), c0.end(), c1.begin());
            }

            auto S = H * C.transpose();

            for (u64 i = 0; i < S.rows(); ++i)
                for (u64 j = 0; j < S.cols(); ++j)
                    assert(S(i, j) == 0);
        }
    }

    struct selectPrt
    {
        const DenseMtx& mMtx;
        const std::vector<u64>& mCols;

        selectPrt(const DenseMtx& m, const std::vector<u64>& c)
            : mMtx(m)
            , mCols(c)
        {}
    };

    std::ostream& operator<<(std::ostream& o, const selectPrt& p)
    {
        for (u64 i = 0; i < p.mMtx.rows(); ++i)
        {
            auto iter = p.mCols.begin();
            for (u64 j = 0; j < p.mMtx.cols(); ++j)
            {
                if (iter != p.mCols.end() && *iter == j)
                    o << oc::Color::Green;

                o << p.mMtx(i, j) << " ";

                if (iter != p.mCols.end() && *iter == j)
                {
                    o << oc::Color::Default;
                    ++iter;
                }
            }

            o << std::endl;
        }
        return o;
    }

    //
    //DenseMtx uniformFixedColWeight(u64 rows, u64 cols, u64 w, PRNG& prng)
    //{
    //    DenseMtx mtx(rows, cols);
    //    mtx.setZero();
    //    std::vector<u64> rem; rem.reserve(cols * w);
    //    for (u64 i = 0; i < cols; ++i)
    //    {
    //        for (u64 j = 0; j < w; ++j)
    //            rem.push_back(i);
    //    }
    //
    //    std::shuffle(rem.begin(), rem.end(), prng);
    //
    //    while (rem.size())
    //    {
    //        auto i = prng.get<u64>() % rows;
    //        mtx(i, rem.back()) = 1;
    //        rem.pop_back();
    //    }
    //    return mtx;
    //}

    void rank(oc::CLP& cmd)
    {

        //u64 n = 6, k = 4;

        //NChooseK nCk(n, k);
        //u64 i = 0;
        //while (nCk)
        //{
        //    auto set0 = *nCk;
        //    auto set1 = ithCombination(i, n, k);

        //    std::cout << i << ":\n";
        //    for (u64 j = 0; j < k; ++j)
        //        std::cout << set0[j] << " ";
        //    std::cout << "  \n";
        //    for (u64 j = 0; j < k; ++j)
        //        std::cout << set1[j] << " ";
        //    std::cout << std::endl;

        //    ++i;
        //    ++nCk;
        //}
        //return;

        u64 rows = cmd.getOr("m", 20);
        u64 cols = static_cast<u64>(rows * cmd.getOr("e", 2.0));
        u64 weight = cmd.getOr("w", 4);
        auto gaps = cmd.getManyOr<u64>("g", { 0 });

        u64 trials = cmd.getOr("t", 1);
        bool verbose = cmd.isSet("v");
        u64 thrds = cmd.getOr("thrds", std::thread::hardware_concurrency());

        assert(cols > rows);
        assert(rows > weight);
        oc::PRNG prng(block(0, cmd.getOr("s", 0)));

        DenseMtx mtx(rows, cols);

        for (auto gap : gaps)
        {

            u64 dWeight = cmd.getOr("wd", (1 + gap) / 2);
            double avg = 0;

            for (u64 t = 0; t < trials; ++t)
            {
                mtx.setZero();

                if (gap)
                {
                    mtx = sampleTriangularBand(
                        rows, cols, weight, gap, dWeight, false, prng).dense();
                }
                else
                {
                    mtx = sampleFixedColWeight(rows, cols, weight, prng, false).dense();

                }

                auto d = minDist(mtx, verbose, thrds);

                if (verbose)
                {
                    std::cout << " " << d.first;
                    std::cout << "\n" << selectPrt(mtx, d.second) << std::endl;
                }

                avg += d.first;
            }

            std::cout << " ~~ " << gap << " " << avg / trials << std::endl;

        }
        //std::cout << "minDist = " << d.size() << std::endl
        //    << "[";
        //for (u64 i = 0; i < d.size(); ++i)
        //    std::cout << d[i] << " ";
        //std::cout << "]" << std::endl;

        //std::cout << selectPrt(mtx, d) << std::endl;

        //auto m2 = mtx.block(0, 0, rows, rows).transpose();
        //std::cout << m2 << std::endl << std::endl;

        //auto m3 = gaussianElim(m2);

        //std::cout << m3 << std::endl << std::endl;

        //auto r = rank(m3);
        //std::cout << "rank " << r << " / " << m3.rows() << std::endl;
        return;
        //Matrix5x3 m = Matrix5x3::Random();
        //std::cout << "Here is the matrix m:" << endl << m << endl;
        //Eigen::FullPivLU<Mtx> lu(mtx);
        //std::cout << "Here is, up to permutations, its LU decomposition matrix:"
        //    << std::endl << lu.matrixLU() << std::endl;
        //std::cout << "Here is the L part:" << std::endl;
        //Mtx l = Mtx::Identity();
        //l.block<5, 3>(0, 0).triangularView<StrictlyLower>() = lu.matrixLU();
        //cout << l << endl;
        //cout << "Here is the U part:" << endl;
        //Matrix5x3 u = lu.matrixLU().triangularView<Upper>();
        //cout << u << endl;
        //cout << "Let us now reconstruct the original matrix m:" << endl;
        //cout << lu.permutationP().inverse() * l * u * lu.permutationQ().inverse() << endl;
    }
}

max_stars_repo_name: WeDPR-Team/libOTe

max_issues_repo_name: WeDPR-Team/libOTe

max_forks_repo_name: WeDPR-Team/libOTe


text:
#include <boost/thread.hpp>
#include <boost/bind.hpp>
#include <boost/asio.hpp>
#include <iostream>

using namespace boost::asio;

io_service service;

void func(int i) {
    std::cout << "func called, i= " << i << std::endl;
}

void worker_thread() {
    service.run();
}

int main(int argc, char* argv[]) {
  for ( int i = 0; i < 10; ++i) {
      service.post(boost::bind(func, i));
  }

  boost::thread_group threads;
  for ( int i = 0; i < 3; ++i) {
      threads.create_thread(worker_thread);
  }
  // wait for all threads to be created
  //boost::this_thread::sleep( boost::posix_time::millisec(500));
  threads.join_all();

  return 0;
}

max_stars_repo_name: chenmingqiang/codesnippet

max_issues_repo_name: chenmingqiang/codesnippet

max_forks_repo_name: chenmingqiang/codesnippet


text:
#ifndef _HOMOMORPHINE_SEAL_BACKED_H_
#define _HOMOMORPHINE_SEAL_BACKED_H_

#include <iostream>
#include <utility>
#include <sstream>
#include <string>
#include <math.h>
#include <seal/seal.h>
#include <boost/log/trivial.hpp>
#include <boost/log/utility/setup.hpp>

#include "util.hpp"
#include "backend.hpp"
#include "arithmetic_backend.hpp"
#include "constants.hpp"

using namespace seal;
using namespace std;
using namespace boost::archive::iterators;

namespace homomorphine
{

  //! SEAL algorithms
  enum class SealAlgorithm {
    BFV,     /*!< SEAL BFV backend */
    CKKS,    /*!< SEAL CKKS backend */
    UNKNOWN  /*!< Unknown algorithm type - usually an error */
  };

  /*! /brief SEAL backend
   *
   * This class is an implementation of the SEAL backend.
   */
  class SealBackend : public ArithmeticBackend
  {
    private:
      std::shared_ptr<SEALContext> context;               /*!< SEAL context object */
      KeyGenerator *keygen = nullptr;                     /*!< SEAL key generator */
      SealAlgorithm algorithm;                            /*!< SEAL algorithm type (BFV, or CKKS) */
      EncryptionParameters *encryption_params = nullptr;  /*!< SEAL specific encryption parameters */
      PublicKey public_key;                               /*!< SEAL public key */
      SecretKey secret_key;                               /*!< SEAL secret key */
      RelinKeys relin_keys;                               /*!< SEAL relinearization keys */
      Ciphertext cipher;                                  /*!< cipher */
      double scale;                                       /*!< bit-precision of encoding, used by CKKS encryption implementation */

      /*!
       * Get SEAL algorithm type from algorithm name
       *
       * \param name algorithm name
       * \return algorithm type
       */
      SealAlgorithm getAlgorithmType(string name);

      /*!
       * Initialize BFV specific stuff of SEAL backend
       */
      void initBFV();

      /*!
       * Initialize CKKS specific stuff of SEAL backend
       */
      void initCKKS();

      /*!
       * Encode vector of values using BFV algorithm
       *
       * \param values vector of values
       * \return plaintext representation of values
       */
      Plaintext encodeWithBFV(vector<long> values);

      /*!
       * Encode single value using BFV algorithm
       *
       * \param value value
       * \return plaintext representation of value
       */
      Plaintext encodeWithBFV(long value);

      /*!
       * Encode vector of values using CKKS algorithm
       *
       * \param values vector of values
       * \return plaintext representation of values
       */
      Plaintext encodeWithCKKS(vector<long> values);

      /*!
       * Encode single value using CKKS algorithm
       *
       * \param value value
       * \return plaintext representation of value
       */
      Plaintext encodeWithCKKS(long value);

      /*!
       * Decode value from plaintext using BFV algorithm
       *
       * \param plain_result plaintext representation of value
       * \return value
       */
      long decodeWithBFV(Plaintext plain_result);

      /*!
       * Decode vector of values from plaintext using BFV algorithm
       *
       * \param plain_result plaintext representation of value
       * \return value
       */
      vector<long> decodeValuesWithBFV(Plaintext plain_result);

      /*!
       * Decode value from plaintext using CKKS algorithm
       *
       * \param plain_result plaintext representation of value
       * \return vector of values
       */
      long decodeWithCKKS(Plaintext plain_result);

      /*!
       * Decode vector of value from plaintext using CKKS algorithm
       *
       * \param plain_result plaintext representation of value
       * \return vector of values
       */
      vector<long> decodeValuesWithCKKS(Plaintext plain_result);

    public:

      /*!
       * SEAL backend cleanup
       */
      ~SealBackend();

      /*!
       * Initializes the SEAL backend
       */
      void init();

      /*!
       * Sets the specific SEAL algorithm implementation
       * that backend provides (currently, BFV and CKKS)
       *
       * \param algorithm homomorphic encryption algorithm
       */
      void setAlgorithm(string algorithm);

      /*!
       * Sets the specific SEAL algorithm implementation
       * that backend provides using the name of algorithm
       *
       * \param algorithm homomorphic encryption algorithm
       */
      void setAlgorithm(SealAlgorithm algorithm);

      /*!
       * Generates the public/secret key pair
       */
      void generateKeys();

      /*!
       * Returns the UUEncoded public key
       *
       * \return public key
       */
      string getPublicKey();

      /*!
       * Writes the public key to a stream
       *
       * \param stream public key stream
       */
      void writePublicKeyToStream(ostream& stream);

      /*!
       * Returns the UUEncoded secret key
       *
       * \return secret key
       */
      string getSecretKey();

      /*!
       * Writes the secret key to a stream
       *
       * \param stream secret key stream
       */
      void writeSecretKeyToStream(ostream& stream);

      /*!
       * Returns the pair of UUEncoded public and secret keys
       *
       * \return pair of public and secret keys
       */
      pair<string, string> getKeys();

      /*!
       * Sets the public key
       *
       * \param public_key UUEncoded public key
       */
      void setPublicKey(string public_key);

      /*!
       * Sets the public key from stream
       *
       * \param stream public key binary stream
       */
      void readPublicKeyFromStream(istream &stream);

      /*!
       * Sets the secret key
       *
       * \param secret_key UUEncoded secret key
       */
      void setSecretKey(string secret_key);

      /*!
       * Sets the secret key from stream
       *
       * \param stream secret key binary stream
       */
      void readSecretKeyFromStream(istream &stream);

      /*!
       * Sets the both public and secret keys
       *
       * \param public_key UUEncoded public key
       * \param secret_key UUEncoded secret key
       */
      void setKeys(string public_key, string secret_key);

       /*!
       * Returns the UUEncoded cipher containing ecrypted value, or vector of values
       *
       * \return UUEncoded cipher
       */
      string getCipher();

      /*!
       * Writes the cipher to output stream
       *
       * \param stream output stream
       */
      void writeCipherToStream(ostream& stream);

      /*!
       * Sets the UUEncoded cipher containing ecrypted value, or vector of values
       *
       * \param cipher UUEncoded cipher
       */
      void setCipher(string cipher);

      /*!
       * Reads the cipher from input stream
       *
       * \param stream cipher stream
       */
      void readCipherFromStream(istream &stream);

      /*!
       * Encrypts the vector of values using the public key
       *
       * \param values vector of values
       * \return UUEncoded cipher
       */
      void encrypt(vector<long> values);

      /*!
       * Encrypts the single value using the public key
       *
       * \param value value
       * \return UUEncoded cipher
       */
      void encrypt(long value);

      /*!
       * Decrypts the vector of values using the secret key
       *
       * \return vector of decrypted values
       */
      vector<long> decryptValues();

      /*!
       * Decrypts the single value using the secret key
       *
       * \return decrypted value
       */
      long decrypt();

      /*!
       * Adds the vector of values to encrypted vector of values
       *
       * \param values vector of values
       */
      void add(vector<long> values);

      /*!
       * Adds the value to encrypted value
       *
       * \param value value
       */
      void add(long value);

      /*!
       * Negates a single encrypted value, or a vector of encrypted values
       */
      void negate();

      /*!
       * Multiplies the vector of values with the encrypted vector of values
       *
       * \param values vector of values
       */
      void multiply(vector<long> values);

      /*!
       * Multiplies the value with the encrypted value
       *
       * \param value value
       */
      void multiply(long value);
  };
}

#endif

ext: hpp

max_stars_repo_name: caboom/homomorphine

max_issues_repo_name: caboom/homomorphine

max_forks_repo_name: caboom/homomorphine