text:
#include "Bubbles.hpp"
#include "Assembler.hpp"
#include "computeLayout.hpp"
#include "orderPairs.hpp"
#include "prefixLength.hpp"
#include "shastaLapack.hpp"
#include "writeGraph.hpp"
using namespace shasta;

#include <boost/graph/connected_components.hpp>
#include <boost/graph/filtered_graph.hpp>
#include <boost/graph/iteration_macros.hpp>

#include <sstream>
#include <queue>

Bubbles::Bubbles(
    const Assembler& assembler, bool debug) :
    assembler(assembler),
    debug(debug)
{

    // Bubbles with discordant ratio greater than this value
    // are considered bad and will be removed.
    const double discordantRatioThreshold = 0.1;

    // Only used for graphics.
    const double minRelativePhase = 0.;

    findBubbles();
    cout << "Found " << bubbles.size() << " bubbles." << endl;

    fillOrientedReadsTable();
    writeOrientedReadsTable();

    createBubbleGraph();
    flagBadBubbles();
    removeBadBubbles(discordantRatioThreshold);
    if(debug) {
        writeBubbles();
        writeBubbleGraphGraphviz();
    }
    cout << "The bubble graph has " << num_vertices(bubbleGraph) <<
        " vertices and " << num_edges(bubbleGraph) << " edges." << endl;

    phase(minRelativePhase);

}

// For now we only consider diploid bubbles, defined using the
// following strict criteria:
// - Source vertex v0 has out-degree 2.
// - Target vertex v1 has in-degree 2.
// - There are two parallel edges eA and eB, both v0->v1.

void Bubbles::findBubbles()
{
    const bool debug = false;

    // Loop over possible choices for v0 (source vertex of the bubble).
    const AssemblyGraph& assemblyGraph = *assembler.assemblyGraphPointer;
    const AssemblyGraph::VertexId vertexCount = assemblyGraph.vertices.size();
    uint64_t suppressedDueToRepeats = 0;
    for(AssemblyGraph::VertexId v0=0; v0<vertexCount; v0++) {

        // Check the out-degree.
        const auto outEdges0 = assemblyGraph.edgesBySource[v0];
        if(outEdges0.size() != 2) {
            continue;
        }

        const AssemblyGraph::EdgeId eA = outEdges0[0];
        const AssemblyGraph::EdgeId eB = outEdges0[1];
        const AssemblyGraph::Edge& edgeA = assemblyGraph.edges[eA];
        const AssemblyGraph::Edge& edgeB = assemblyGraph.edges[eB];

        // Check the target vertex.
        const AssemblyGraph::VertexId v1 = edgeA.target;
        if(edgeB.target != v1) {
            continue;
        }
        const auto inEdges1 = assemblyGraph.edgesByTarget[v1];
        if(inEdges1.size() != 2) {
            continue;
        }

        // Assemble the RLE sequence for the two edges.
        vector<Base> sequenceA;
        vector<Base> sequenceB;
        assembler.assembleAssemblyGraphEdgeRleStrict(eA, sequenceA);
        assembler.assembleAssemblyGraphEdgeRleStrict(eB, sequenceB);

        // If the RLE sequences are identical, skip it.
        // This can happen due to missing alignments.
        if(sequenceA == sequenceB) {
            continue;
        }

        if(debug) {
            const MarkerGraph::VertexId u0 = assemblyGraph.vertices[v0];
            const MarkerGraph::VertexId u1 = assemblyGraph.vertices[v1];
            cout << "Bubble at marker graph " << u0 << "->" << u1 << "\n";
            copy(sequenceA.begin(), sequenceA.end(), ostream_iterator<Base>(cout));
            cout << "\n";
            copy(sequenceB.begin(), sequenceB.end(), ostream_iterator<Base>(cout));
            cout << "\n";
        }

        // If the two sequences differ by a 2- of 3-base repeat,
        // this is probably an erroneous bubble.
        if(isShortRepeatCopyNumberDifference(sequenceA, sequenceB)) {
            ++suppressedDueToRepeats;
            if(debug) {
                cout << "Is short repeat difference - skipping.\n";
            }
            continue;
        }

        // We have a diploid bubble.
        bubbles.push_back(Bubble(v0, v1, eA, eB, assembler.markerGraph, assemblyGraph));

    }

    cout << suppressedDueToRepeats << " bubbles were suppressed due to 2- or 3-base repeats." << endl;
}

void Bubbles::writeBubbles()
{
    ofstream csv("Bubbles.csv");
    csv << "BubbleId,VertexId0,VertexId1,CoverageA,CoverageB,ConcordantSum,DiscordantSum,DiscordantRatio\n";

    for(uint64_t bubbleId=0; bubbleId<bubbles.size(); bubbleId++) {
        const Bubble& bubble = bubbles[bubbleId];
        csv << bubbleId << ",";
        csv << bubble.mv0 << ",";
        csv << bubble.mv1 << ",";
        csv << bubble.orientedReadIds[0].size() << ",";
        csv << bubble.orientedReadIds[1].size() << ",";
        csv << bubble.concordantSum << ",";
        csv << bubble.discordantSum << ",";
        csv << bubble.discordantRatio() << "\n";
    }
}

// Figure out if two sequences differ only by copy numbers in
// a 2- or 3-base repeat.
bool Bubbles::isShortRepeatCopyNumberDifference(
    const vector<Base>& x,
    const vector<Base>& y)
{
    const bool debug = false;

    // Get the lengths and their differences.
    const uint64_t nx = x.size();
    const uint64_t ny = y.size();

    // When they have the same length, we can certainly return false.
    if(nx == ny) {
        return false;
    }

    // Recursive call so x is shorter than y.
    if(ny < nx) {
        return isShortRepeatCopyNumberDifference(y, x);
    }
    SHASTA_ASSERT(nx < ny);

    // If the length difference is not a multiple of 2 or 3,
    // we can certainly return false.
    const uint64_t dn = ny - nx;
    if( ((dn%2) != 0) and
        ((dn%3) != 0) ) {
        if(debug) {
            cout << "Length difference is not a multiple of 2 or 3\n";
            cout << "nx " << nx << ", ny " << ny << ", dn " << dn << "\n";
        }
        return false;
    }

    const uint64_t commonPrefixLength = shasta::commonPrefixLength(x, y);
    const uint64_t commonSuffixLength = shasta::commonSuffixLength(x, y);

    // Find the portion of y that is not in x.
    uint64_t ix = commonPrefixLength;
    uint64_t iy = commonPrefixLength;
    uint64_t jx = nx - commonSuffixLength;
    uint64_t jy = ny - commonSuffixLength;
    while((jx<ix) or (jy<iy)) {
        ++jx;
        ++jy;
    }

    if(debug) {
        cout << "commonPrefixLength " << commonPrefixLength << "\n";
        cout << "commonSuffixLength " << commonSuffixLength << "\n";
        cout << "nx " << nx << "\n";
        cout << "ny " << ny << "\n";
        cout << "dn " << dn << "\n";
        cout << "ix " << ix << "\n";
        cout << "jx " << jx << "\n";
        cout << "iy " << iy << "\n";
        cout << "jy " << jy << "\n";
    }

    if(ix != jx) {
        // There is more than just an insertion.
        return false;
    }

    // x and y differ by an insertion in iy of range [iy, jy).
    SHASTA_ASSERT(ix == jx);
    SHASTA_ASSERT(jy - iy == dn);

    // Check for k base repeat.
    // We kept the entire common prefix, so we can check just to the left of the insertion.
    for(uint64_t k=2; k<=3; k++) {
        if((dn % k) != 0) {
            continue;
        }
        if(debug) {
            cout << "Trying k = " << k << "\n";
        }

        // Check that the inserted bases are a repeat with period k.
        const uint64_t m = dn / k;
        bool repeatViolationFound = false;
        for(uint64_t i=0; i<m; i++) {
            for(uint64_t j=0; j<k; j++) {
                if(y[iy + i*k + j] != y[iy + j]) {
                    repeatViolationFound = true;
                    break;
                }
            }
        }
        if(repeatViolationFound) {
            // The inserted portion is not a repeat of period k.
            continue;
        }

        // Check the previous k bases in both x and y.
        if(ix < k) {
            continue;
        }
        if(iy < k) {
            continue;
        }
        for(uint64_t j=0; j<k; j++) {
            if(y[iy - k + j] != y[ix + j]) {
                repeatViolationFound = true;
                break;
            }
            if(x[ix - k + j] != y[ix + j]) {
                repeatViolationFound = true;
                break;
            }
        }
        if(repeatViolationFound) {
            continue;
        }

        // It is an insertion in y of a repeat of length k.
        return true;
    }

    // None of the k values we tried worked.
    return false;
}

Bubbles::Bubble::Bubble(
    AssemblyGraph::VertexId av0,
    AssemblyGraph::VertexId av1,
    AssemblyGraph::EdgeId eA,
    AssemblyGraph::EdgeId eB,
    const MarkerGraph& markerGraph,
    const AssemblyGraph& assemblyGraph) :
    av0(av0),
    av1(av1),
    aEdgeIds({eA, eB}),
    mv0(assemblyGraph.vertices[av0]),
    mv1(assemblyGraph.vertices[av1])
{
    fillInOrientedReadIds(markerGraph, assemblyGraph);
}

void Bubbles::Bubble::fillInOrientedReadIds(
    const MarkerGraph& markerGraph,
    const AssemblyGraph& assemblyGraph)
{

    // Loop over both sides of the bubble.
    array<std::set<OrientedReadId>, 2> orientedReadIdSets;
    for(uint64_t i=0; i<2; i++) {
        const AssemblyGraph::EdgeId aEdgeId = aEdgeIds[i];

        // Loop over marker graph edges of this side.
        const span<const MarkerGraph::EdgeId> mEdgeIds = assemblyGraph.edgeLists[aEdgeId];
        for(const MarkerGraph::EdgeId mEdgeId : mEdgeIds) {

            // Loop over MarkerIntervals of this marker graph edge.
            const span<const MarkerInterval> markerIntervals =
                markerGraph.edgeMarkerIntervals[mEdgeId];
            for(const MarkerInterval& markerInterval: markerIntervals) {
                orientedReadIdSets[i].insert(markerInterval.orientedReadId);
            }
        }
    }

    // Now store them, discarding the ones that appear on both sides.
    for(uint64_t i=0; i<2; i++) {
        const std::set<OrientedReadId>& x = orientedReadIdSets[i];
        const std::set<OrientedReadId>& y = orientedReadIdSets[1-i];
        for(const OrientedReadId orientedReadId: x) {
            if(y.find(orientedReadId) == y.end()) {
                orientedReadIds[i].push_back(orientedReadId);
            }
        }
    }
}

void Bubbles::fillOrientedReadsTable()
{
    orientedReadsTable.resize(2* assembler.getReads().readCount());
    for(uint64_t bubbleId=0; bubbleId<bubbles.size(); bubbleId++) {
        const Bubble& bubble = bubbles[bubbleId];
        for(uint64_t side=0; side<2; side++) {
            for(const OrientedReadId orientedReadId: bubble.orientedReadIds[side]) {
                orientedReadsTable[orientedReadId.getValue()].push_back(make_pair(bubbleId, side));
            }
        }
    }
}

void Bubbles::writeOrientedReadsTable()
{
    ofstream csv("OrientedReadsTable.csv");
    csv << "OrientedReadId,Bubble,Side\n";
    for(ReadId readId=0; readId<assembler.getReads().readCount(); readId++) {
        for(Strand strand=0; strand<2; strand++) {
            const OrientedReadId orientedReadId(readId, strand);
            for(const auto& v: orientedReadsTable[orientedReadId.getValue()]) {
                csv << orientedReadId << ",";
                csv << v.first << ",";
                csv << v.second << "\n";
            }
        }
    }
}

// Given two OrientedReadIds, use the orientedReadsTable
// to count the number of times they appear on the same
// side or opposite sides of the same bubble.
void Bubbles::findOrientedReadsRelativePhase(
    OrientedReadId orientedReadId0,
    OrientedReadId orientedReadId1,
    uint64_t& sameSideCount,
    uint64_t& oppositeSideCount
) const
{
    const vector<pair <uint64_t, uint64_t> >& v0 = orientedReadsTable[orientedReadId0.getValue()];
    const vector<pair <uint64_t, uint64_t> >& v1 = orientedReadsTable[orientedReadId1.getValue()];

    sameSideCount = 0;
    oppositeSideCount = 0;

    auto begin0 = v0.begin();
    auto begin1 = v1.begin();
    auto end0 = v0.end();
    auto end1 = v1.end();

    auto it0 = begin0;
    auto it1 = begin1;

    // The loop makes use of the fact that an OrientedReadId cannot appear
    // on both sides of a bubble.
    while((it0!=end0) and (it1!=end1)) {
        const uint64_t bubbleId0 = it0->first;
        const uint64_t bubbleId1 = it1->first;

        if(bubbleId0 < bubbleId1) {
            ++it0;
            continue;
        }

        if(bubbleId1 < bubbleId0) {
            ++it1;
            continue;
        }

        SHASTA_ASSERT(bubbleId0 == bubbleId1);

        if(not bubbles[bubbleId0].isBad) {
            const uint64_t side0 = it0->second;
            const uint64_t side1 = it1->second;
            if(side0 == side1) {
                ++sameSideCount;
            } else {
                ++oppositeSideCount;
            }
        }

        ++it0;
        ++it1;
    }
}

// Find OrientedReadIds that appear in at least one bubble
// together with a given OrientedReadId.
// They are returned sorted.
void Bubbles::findNeighborOrientedReadIds(
    OrientedReadId orientedReadId0,
    vector<OrientedReadId>& orientedReadIds
) const
{
    const vector <pair <uint64_t, uint64_t> >& v0 = orientedReadsTable[orientedReadId0.getValue()];
    std::set<OrientedReadId> orientedReadIdsSet;
    for(const auto& p: v0) {
        const uint64_t bubbleId = p.first;
        const Bubble& bubble = bubbles[bubbleId];
        if(bubble.isBad) {
            continue;
        }
        for(uint64_t side=0; side<2; side++) {
            for(const OrientedReadId orientedReadId1: bubble.orientedReadIds[side]) {
                if(orientedReadId1 != orientedReadId0) {
                    orientedReadIdsSet.insert(orientedReadId1);
                }
            }
        }
    }

    orientedReadIds.clear();
    copy(orientedReadIdsSet.begin(), orientedReadIdsSet.end(), back_inserter(orientedReadIds));
}

void Bubbles::createBubbleGraph()
{
    // Create the vertices.
    for(uint64_t bubbleId=0; bubbleId<bubbles.size(); bubbleId++) {
        bubbleGraph.vertexTable.push_back(add_vertex(BubbleGraphVertex(bubbleId), bubbleGraph));
    }

    // Create the edges.
    for(uint64_t bubbleIdA=0; bubbleIdA<bubbles.size(); bubbleIdA++) {
        const Bubble& bubbleA = bubbles[bubbleIdA];
        const BubbleGraph::vertex_descriptor vA = bubbleGraph.vertexTable[bubbleIdA];

        for(uint64_t sideA=0; sideA<2; sideA++) {
            const vector<OrientedReadId>& orientedReadIds = bubbleA.orientedReadIds[sideA];

            for(const OrientedReadId orientedReadId: orientedReadIds) {
                const auto& v = orientedReadsTable[orientedReadId.getValue()];

                for(const auto& p: v) {
                    const uint64_t bubbleIdB = p.first;
                    if(bubbleIdB <= bubbleIdA) {
                        continue;
                    }
                    const uint64_t sideB = p.second;
                    const BubbleGraph::vertex_descriptor vB = bubbleGraph.vertexTable[bubbleIdB];

                    // Locate the edge between these two bubbles,
                    // creating it if necessary.
                    BubbleGraph::edge_descriptor e;
                    bool edgeWasFound = false;
                    tie(e, edgeWasFound) = edge(vA, vB, bubbleGraph);
                    if(not edgeWasFound) {
                        tie(e, edgeWasFound) = add_edge(vA, vB, bubbleGraph);
                    }
                    SHASTA_ASSERT(edgeWasFound);

                    // Update the matrix.
                    BubbleGraphEdge& edge = bubbleGraph[e];
                    ++edge.matrix[sideA][sideB];
                }
            }
        }
    }

    ofstream csv("BubbleGraphEdges.csv");
    csv << "BubbleIdA,BubbleIdB,m00,m11,m01,m10\n";
    BGL_FORALL_EDGES(e, bubbleGraph, BubbleGraph) {
        const BubbleGraphEdge& edge = bubbleGraph[e];
        csv << bubbleGraph[source(e, bubbleGraph)].bubbleId << ",";
        csv << bubbleGraph[target(e, bubbleGraph)].bubbleId << ",";
        csv << edge.matrix[0][0] << ",";
        csv << edge.matrix[1][1] << ",";
        csv << edge.matrix[0][1] << ",";
        csv << edge.matrix[1][0] << "\n";
    }
}

void Bubbles::writeBubbleGraphGraphviz() const
{

    ofstream out("BubbleGraph.dot");
    out << "graph G{\n"
        "node [shape=point];\n";

    BGL_FORALL_VERTICES(v, bubbleGraph, BubbleGraph) {
        const BubbleGraphVertex& vertex = bubbleGraph[v];
        const uint64_t bubbleId = vertex.bubbleId;
        const Bubble& bubble = bubbles[bubbleId];
        const double discordantRatio = bubble.discordantRatio();
        SHASTA_ASSERT(discordantRatio <= 0.5);
        const double hue = (1. - 2. * discordantRatio) * 0.333333;  //  Good = green, bad = red, via yellow
        out << bubbleId <<
            " [color=\"" << hue << ",1.,0.8\""
            " tooltip=\""  << bubbleId << " " << discordantRatio << "\""
            "];\n";
    }

    BGL_FORALL_EDGES(e, bubbleGraph, BubbleGraph) {
        const BubbleGraphEdge& edge = bubbleGraph[e];
        const uint64_t bubbleIdA = bubbleGraph[source(e, bubbleGraph)].bubbleId;
        const uint64_t bubbleIdB = bubbleGraph[target(e, bubbleGraph)].bubbleId;

        const double hue = (1. - edge.ambiguity()) * 0.333333; //  Good = green, bad = red, via yellow

        out << bubbleIdA << "--" << bubbleIdB <<
            " ["
            // "penwidth=" << 0.1*double(total) <<
            " color=\"" << hue << ",1.,0.8\""
            "];\n";
    }

    out << "}\n";
}

// Write a single component of the BubbleGraph in svg format.
// To compute sfdp layout, only consider edges
// for which relativePhase() >= minRelativePhase.
void Bubbles::writeBubbleGraphComponentHtml(
    uint64_t componentId,
    const vector<BubbleGraph::vertex_descriptor>& component,
    double minRelativePhase) const
{
    using vertex_descriptor = BubbleGraph::vertex_descriptor;
    using edge_descriptor = BubbleGraph::edge_descriptor;

    // Map the vertices of this component to integers.
    std::map<vertex_descriptor, uint64_t> componentVertexMap;
    for(uint64_t i=0; i<component.size(); i++) {
        componentVertexMap.insert(make_pair(component[i], i));
    }

    // Create a graph with only this connected component and only
    // edges with relativePhase >= minRelativePhase.
    // This will be used to compute the sfdp layout.
    ComponentGraph componentGraph(component.size());
    for(uint64_t i0=0; i0<component.size(); i0++) {
        const vertex_descriptor v0 = component[i0];
        BGL_FORALL_OUTEDGES(v0, e, bubbleGraph, BubbleGraph) {
            if(bubbleGraph[e].relativePhase() < minRelativePhase) {
                continue;
            }
            const vertex_descriptor v1 = target(e, bubbleGraph);
            auto it = componentVertexMap.find(v1);
            SHASTA_ASSERT(it != componentVertexMap.end());
            const uint64_t i1 = it->second;
            if(i0 < i1) {
                add_edge(i0, i1, componentGraph);
            }
        }
    }

    // Compute the layout of this graph.
    std::map<ComponentGraph::vertex_descriptor, array<double, 2> > positionMap;
    SHASTA_ASSERT(computeLayout(componentGraph, "sfdp", 600., positionMap) == ComputeLayoutReturnCode::Success);
    for(uint64_t i=0; i<component.size(); i++) {
        componentGraph[i].position = positionMap[i];
        // cout << componentGraph[i].position[0] << " " << componentGraph[i].position[1] << "\n";
    }

    // Add all the remaining edges.
    for(uint64_t i0=0; i0<component.size(); i0++) {
        const vertex_descriptor v0 = component[i0];
        BGL_FORALL_OUTEDGES(v0, e, bubbleGraph, BubbleGraph) {
            if(bubbleGraph[e].relativePhase() >= minRelativePhase) {
                // We already added this edge.
                continue;
            }
            const vertex_descriptor v1 = target(e, bubbleGraph);
            auto it = componentVertexMap.find(v1);
            SHASTA_ASSERT(it != componentVertexMap.end());
            const uint64_t i1 = it->second;
            if(i0 < i1) {
                add_edge(i0, i1, componentGraph);
            }
        }
    }

    // Graphics scaling.
    double xMin = std::numeric_limits<double>::max();
    double xMax = std::numeric_limits<double>::min();
    double yMin = std::numeric_limits<double>::max();
    double yMax = std::numeric_limits<double>::min();
    BGL_FORALL_VERTICES_T(v, componentGraph, ComponentGraph) {
        const auto& position = componentGraph[v].position;
        xMin = min(xMin, position[0]);
        xMax = max(xMax, position[0]);
        yMin = min(yMin, position[1]);
        yMax = max(yMax, position[1]);
    }
    const double xyRange = max(xMax-xMin, yMax-yMin);
    const int svgSize = 1024;
    const double vertexRadiusPixels = 3.;
    const double vertexRadius = vertexRadiusPixels * xyRange / double(svgSize);
    const double edgeThicknessPixels = 1.;
    const double edgeThickness = edgeThicknessPixels * xyRange / double(svgSize);

    // Write the componentGraph in svg format.

    // Vertex attributes. Color by phase.
    std::map<ComponentGraph::vertex_descriptor, WriteGraph::VertexAttributes> vertexAttributes;
    BGL_FORALL_VERTICES(v, componentGraph, ComponentGraph) {
        const BubbleGraph::vertex_descriptor vb = component[v];
        const BubbleGraphVertex& vertex = bubbleGraph[vb];
        const uint64_t bubbleId = vertex.bubbleId;
        const int64_t phase = vertex.phase;
        auto& attributes = vertexAttributes[v];
        attributes.radius = vertexRadius;
        if(phase == +1) {
            attributes.color = "hsl(240,50%,50%)";
        }
        if(phase == -1) {
            attributes.color = "hsl(300,50%,50%)";
        }
        attributes.tooltip = to_string(bubbleId);
    }

    // Edge attributes. Color by relative phase.
    std::map<ComponentGraph::edge_descriptor, WriteGraph::EdgeAttributes> edgeAttributes;
    BGL_FORALL_EDGES(e, componentGraph, ComponentGraph) {
        const vertex_descriptor v0 = component[source(e, componentGraph)];
        const vertex_descriptor v1 = component[target(e, componentGraph)];
        edge_descriptor eb;
        bool edgeWasFound = false;
        tie(eb, edgeWasFound) = edge(v0, v1, bubbleGraph);
        const BubbleGraphEdge& edge = bubbleGraph[eb];
        const double relativePhase = edge.relativePhase();
        const double hue = (1. + relativePhase) * 60.; /// Goes from 0 (red) to 120 (green).
        auto& attributes = edgeAttributes[e];
        attributes.thickness = edgeThickness;
        if(relativePhase > 0.) {
            attributes.color = "hsla(" + to_string(int(hue)) + ",50%,50%,100%)";
        } else {
            attributes.color = "hsla(" + to_string(int(hue)) + ",50%,50%,20%)";
        }
    }

    // Draw the svg.
    ofstream out("BubbleGraph-Component-" + to_string(componentId) + ".html");
    out << "<html><body>";
    WriteGraph::writeSvg(
        componentGraph,
        "component" + to_string(componentId),
        svgSize, svgSize,
        vertexAttributes,
        edgeAttributes,
        out);
    out << "</body></html>";
}

// Write in html/svg format.
// To compute sfdp layout, only consider edges
// for which relativePhase() >= minRelativePhase.
void Bubbles::PhasingGraph::writeHtml(
    const string& fileName,
    double minRelativePhase) const
{
    using G = PhasingGraph;
    const G& g = *this;

    // Create a filtered PhasingGraph, containing only the edges
    // with relativePhase() >= minRelativePhase.
    using FilteredGraph = boost::filtered_graph<G, PhasingGraphEdgePredicate>;
    FilteredGraph filteredGraph(g, PhasingGraphEdgePredicate(g, minRelativePhase));

    // Compute the layout of the filtered graph.
    std::map<FilteredGraph::vertex_descriptor, array<double, 2> > positionMap;
    SHASTA_ASSERT(computeLayout(filteredGraph, "sfdp", 600., positionMap) == ComputeLayoutReturnCode::Success);
    BGL_FORALL_VERTICES(v, filteredGraph, FilteredGraph) {
        filteredGraph[v].position = positionMap[v];
    }

    // Graphics scaling.
    double xMin = std::numeric_limits<double>::max();
    double xMax = std::numeric_limits<double>::min();
    double yMin = std::numeric_limits<double>::max();
    double yMax = std::numeric_limits<double>::min();
    BGL_FORALL_VERTICES_T(v, g, G) {
        const auto& position = g[v].position;
        xMin = min(xMin, position[0]);
        xMax = max(xMax, position[0]);
        yMin = min(yMin, position[1]);
        yMax = max(yMax, position[1]);
    }
    const double xyRange = max(xMax-xMin, yMax-yMin);
    const int svgSize = 1024;
    const double vertexRadiusPixels = 3.;
    const double vertexRadius = vertexRadiusPixels * xyRange / double(svgSize);
    const double edgeThicknessPixels = 1.;
    const double edgeThickness = edgeThicknessPixels * xyRange / double(svgSize);

    // Write the componentGraph in svg format.

    // Vertex attributes. Color by phase.
    std::map<G::vertex_descriptor, WriteGraph::VertexAttributes> vertexAttributes;
    BGL_FORALL_VERTICES(v, g, G) {
        const PhasingGraphVertex& vertex = g[v];
        const int64_t phase = vertex.phase;
        auto& attributes = vertexAttributes[v];
        attributes.radius = vertexRadius;
        if(phase == +1) {
            attributes.color = "hsl(240,50%,50%)";
        }
        if(phase == -1) {
            attributes.color = "hsl(300,50%,50%)";
        }
        std::ostringstream s;
        s << vertex.eigenvectorComponent;
        attributes.tooltip = vertex.orientedReadId.getString() + " " +
            s.str();
    }

    // Edge attributes. Color by relative phase.
    std::map<G::edge_descriptor, WriteGraph::EdgeAttributes> edgeAttributes;
    BGL_FORALL_EDGES(e, g, G) {
        const PhasingGraphEdge& edge = g[e];
        const double relativePhase = edge.relativePhase();
        const double hue = (1. + relativePhase) * 60.; /// Goes from 0 (red) to 120 (green).
        auto& attributes = edgeAttributes[e];
        attributes.thickness = edgeThickness;
        if(relativePhase > 0.) {
            attributes.color = "hsla(" + to_string(int(hue)) + ",50%,50%,100%)";
        } else {
            attributes.color = "hsla(" + to_string(int(hue)) + ",50%,50%,20%)";
        }
    }

    // Draw the svg.
    ofstream out(fileName);
    out << "<html><body>";
    WriteGraph::writeSvg(
        g,
        "PhasingGraph",
        svgSize, svgSize,
        vertexAttributes,
        edgeAttributes,
        out);
    out << "</body></html>";
}

// Use the BubbleGraph to flag bad bubbles.
void Bubbles::flagBadBubbles()
{
    for(uint64_t bubbleId=0; bubbleId<bubbles.size(); bubbleId++) {
        Bubble& bubble = bubbles[bubbleId];
        bubble.concordantSum = 0;
        bubble.discordantSum = 0;
        BGL_FORALL_OUTEDGES(bubbleGraph.vertexDescriptor(bubbleId), e, bubbleGraph, BubbleGraph) {
            const BubbleGraphEdge& edge = bubbleGraph[e];
            const uint64_t diagonal = edge.matrix[0][0] + edge.matrix[1][1];
            const uint64_t offDiagonal = edge.matrix[0][1] + edge.matrix[1][0];
            const uint64_t concordant = max(diagonal, offDiagonal);
            const uint64_t discordant = min(diagonal, offDiagonal);
            bubble.concordantSum += concordant;
            bubble.discordantSum += discordant;
        }
    }
}

void Bubbles::removeBadBubbles(double discordantRatioThreshold)
{
    uint64_t removedCount = 0;
    for(uint64_t bubbleId=0; bubbleId<bubbles.size(); bubbleId++) {
        if(bubbles[bubbleId].discordantRatio() > discordantRatioThreshold) {
            const BubbleGraph::vertex_descriptor v = bubbleGraph.vertexDescriptor(bubbleId);
            clear_vertex(v, bubbleGraph);
            remove_vertex(v, bubbleGraph);
            bubbleGraph.vertexTable[bubbleId] = BubbleGraph::null_vertex();
            bubbles[bubbleId].isBad = true;
            ++removedCount;
        }
    }

    cout << "Removed " << removedCount << " bad bubbles out of " <<
        bubbles.size() << " total." << endl;

}

// Phase bubbles and reads.
// This should be called after bad bubbles were already removed
// from the bubble graph.
void Bubbles::phase(
    double minRelativePhase)
{
    bubbleGraph.computeConnectedComponents();
    cout << "Found " << bubbleGraph.connectedComponents.size() <<
        " connected components of the bubble graph." << endl;

    orientedReadsPhase.resize(orientedReadsTable.size());
        make_pair(
        std::numeric_limits<uint32_t>::max(),
        std::numeric_limits<uint32_t>::max());

    uint64_t phasedReadCount = 0;
    for(uint64_t componentId=0; componentId<bubbleGraph.connectedComponents.size(); componentId++) {
        if(debug) {
            cout << "Working on connected component " << componentId << endl;
        }
        const vector<BubbleGraph::vertex_descriptor>& component = bubbleGraph.connectedComponents[componentId];

        // Phase the bubbles.
        if(debug) {
            phaseComponentBubbles(component);
            writeBubbleGraphComponentHtml(componentId, component, minRelativePhase);
        }

        // Phase the oriented reads.
        vector<OrientedReadId> componentOrientedReadIds;
        findComponentOrientedReads(component, componentOrientedReadIds);
        phasedReadCount += componentOrientedReadIds.size();
        PhasingGraph phasingGraph;
        createPhasingGraph(componentOrientedReadIds, phasingGraph);
        phasingGraph.phaseSpectral(debug);

        if(true) {
            phasingGraph.writeHtml(
                "PhasingGraph-Component-" + to_string(componentId) + ".html",
                minRelativePhase);
        }

        // Copy the phasing for this connected component to the global
        // phasing vector orientedReadsPhase.
        BGL_FORALL_VERTICES(v, phasingGraph, PhasingGraph) {
            PhasingGraphVertex& vertex = phasingGraph[v];
            const OrientedReadId orientedReadId = vertex.orientedReadId;
            auto& p = orientedReadsPhase[orientedReadId.getValue()];
            p.first = uint32_t(componentId);
            if(vertex.phase == 1) {
                p.second = 0;
            } else {
                p.second = 1;
            }
        }

    }
    cout << phasedReadCount << " oriented reads were phased, out of " <<
        orientedReadsTable.size() << " total." << endl;

    // Write out the phasing.
    ofstream csv("Phasing.csv");
    for(uint64_t i=0; i<orientedReadsPhase.size(); i++) {
        const OrientedReadId orientedReadId = OrientedReadId::fromValue(ReadId(i));
        const auto& p = orientedReadsPhase[i];
        const uint32_t componentId = p.first;
        const uint32_t phase = 0;
        csv << orientedReadId << ",";
        if(componentId != std::numeric_limits<uint32_t>::max()) {
            csv << componentId << ",";
            csv << phase << "\n";
        } else {
            csv << ",\n";
        }

    }
}

// Phase the bubbles of a connected component of the BubbleGraph.
// This sets the phase in the component vertices.
void Bubbles::phaseComponentBubbles(
    const vector<BubbleGraph::vertex_descriptor>& component)
{
    using vertex_descriptor = BubbleGraph::vertex_descriptor;

    if(debug) {
        cout << "Phasing bubbles of a connected component with " << component.size() <<
            " bubbles." << endl;
    }

    // Map the vertices of this component to integers.
    std::map<vertex_descriptor, uint64_t> componentVertexMap;
    const uint64_t n = component.size();
    for(uint64_t i=0; i<n; i++) {
        componentVertexMap.insert(make_pair(component[i], i));
    }

    // Create the similarity matrix.
    // For the similarity of two bubbles use
    // diagonalCount() - offDiagonalCount().
    // This can be negative, but spectral clustering is still possible.
    vector<double> A(n * n, 0.);
    for(uint64_t i0=0; i0<component.size(); i0++) {
        const vertex_descriptor v0 = component[i0];
        BGL_FORALL_OUTEDGES(v0, e, bubbleGraph, BubbleGraph) {
            const BubbleGraphEdge& edge = bubbleGraph[e];
            const vertex_descriptor v1 = target(e, bubbleGraph);
            auto it = componentVertexMap.find(v1);
            SHASTA_ASSERT(it != componentVertexMap.end());
            const uint64_t i1 = it->second;
            SHASTA_ASSERT(i0 != i1);
            A[i0*n + i1] = double(edge.diagonalCount()) - double(edge.offDiagonalCount());
        }
    }

    // Check that the similarity matrix is symmetric.
    for(uint64_t i0=0; i0<n; i0++) {
        for(uint64_t i1=0; i1<n; i1++) {
            SHASTA_ASSERT(A[i0*n + i1] == A[i0 + i1*n]);
        }
    }

    // Compute the sum of each row. This will become the diagonal of the Laplacian matrix.
    vector<double> D(n, 0.);
    for(uint64_t i0=0; i0<n; i0++) {
        for(uint64_t i1=0; i1<n; i1++) {
            D[i0] += A[i0 + i1*n];
        }
    }

    // Now turn A from the similarity matrix into the Laplacian matrix.
    for(uint64_t i0=0; i0<n; i0++) {
        for(uint64_t i1=0; i1<n; i1++) {
            if(i0 == i1) {
                A[i0 + i1*n] = D[i0];
            } else {
                A[i0 + i1*n] = - A[i0 + i1*n];
            }
        }
    }

    // Compute eigenvalues and eigenvectors.
    int N = int(n);
    vector<double> W(n);
    const int LWORK = 10 * N;
    vector<double> WORK(LWORK);
    int INFO = 0;
    dsyev_("V", "L", N, &A[0], N, &W[0], &WORK[0], LWORK, INFO);
    SHASTA_ASSERT(INFO == 0);

    // Get the phase from the sign of the components of the first eigenvector.
    for(uint64_t i=0; i<n; i++) {
        if(A[i] >= 0.) {
            bubbleGraph[component[i]].phase = 1;
        } else {
            bubbleGraph[component[i]].phase = -1;
        }

    }

}

// Phase the oriented reads of a connected component of the BubbleGraph.
// This stores the phases in the orientedReadsPhase vector.
void Bubbles::PhasingGraph::phaseSpectral(bool debug)
{
    using G = PhasingGraph;
    G& g = *this;

    const uint64_t n = vertexMap.size();
    if(debug) {
        cout << "Phasing oriented reads of a PhasingGraph with " << n <<
            " oriented reads." << endl;
    }

    // Map vertices to integers.
    std::map<vertex_descriptor, uint64_t> vertexToIntegerMap;
    uint64_t i = 0;
    BGL_FORALL_VERTICES(v, g, G) {
        vertexToIntegerMap.insert(make_pair(v, i++));
    }

    // Create the similarity matrix.
    // For the similarity of two oriented reads use
    // sameSideCount - oppositeSideCount.
    // This can be negative, but spectral clustering is still possible.
    vector<double> A(n * n, 0.);
    BGL_FORALL_EDGES(e, g, G) {
        const PhasingGraphEdge& edge = g[e];
        const vertex_descriptor v0 = source(e, g);
        const vertex_descriptor v1 = target(e, g);

        // Find the matrix indexes for these two vertices.
        const auto it0 = vertexToIntegerMap.find(v0);
        SHASTA_ASSERT(it0 != vertexToIntegerMap.end());
        const uint64_t i0 = it0->second;
        const auto it1 = vertexToIntegerMap.find(v1);
        SHASTA_ASSERT(it1 != vertexToIntegerMap.end());
        const uint64_t i1 = it1->second;

        const double similarity = double(edge.sameSideCount) - double(edge.oppositeSideCount);
        A[i0*n + i1] = similarity;
        A[i0 + i1*n] = similarity;
    }

    // Check that the similarity matrix is symmetric.
    for(uint64_t i0=0; i0<n; i0++) {
        for(uint64_t i1=0; i1<n; i1++) {
            SHASTA_ASSERT(A[i0*n + i1] == A[i0 + i1*n]);
        }
    }

    // Compute the sum of each row. This will become the diagonal of the Laplacian matrix.
    vector<double> D(n, 0.);
    for(uint64_t i0=0; i0<n; i0++) {
        for(uint64_t i1=0; i1<n; i1++) {
            D[i0] += A[i0 + i1*n];
        }
    }

    // Now turn A from the similarity matrix into the Laplacian matrix.
    for(uint64_t i0=0; i0<n; i0++) {
        for(uint64_t i1=0; i1<n; i1++) {
            if(i0 == i1) {
                A[i0 + i1*n] = D[i0];
            } else {
                A[i0 + i1*n] = - A[i0 + i1*n];
            }
        }
    }

    // Compute eigenvalues and eigenvectors.
    int N = int(n);
    vector<double> W(n);
    const int LWORK = 10 * N;
    vector<double> WORK(LWORK);
    int INFO = 0;
    dsyev_("V", "L", N, &A[0], N, &W[0], &WORK[0], LWORK, INFO);
    SHASTA_ASSERT(INFO == 0);

    // Get the phase from the sign of the components of the first eigenvector.
    BGL_FORALL_VERTICES(v, g, G) {
        PhasingGraphVertex& vertex = g[v];
        const double eigenvectorComponent = A[vertexToIntegerMap[v]];
        vertex.eigenvectorComponent = eigenvectorComponent;

        if(eigenvectorComponent >= 0.) {
            vertex.phase = 1;
        } else {
            vertex.phase = -1;
        }
    }
}

// Phase a connected component using the SVD.
void Bubbles::phaseSvd(
    const vector<BubbleGraph::vertex_descriptor>& bubbleGraphComponent,
    PhasingGraph& phasingGraph)
{

    // Map global bubbleIds to indices in bubbleGraphComponent.
    std::map<uint64_t, uint64_t> bubbleMap;
    uint64_t i = 0;
    for(const BubbleGraph::vertex_descriptor v: bubbleGraphComponent) {
        bubbleMap.insert(make_pair(bubbleGraph[v].bubbleId, i++));
    }
    const uint64_t bubbleCount = bubbleMap.size();

    // Create a matrix with a row for each oriented read and a column for each bubble.
    // The matrix is stored in Fortran storage scheme (that is, by columns).
    // An element is set to +1 if the oriented read appears in the first side of the bubble
    // and to -1 if the oriented read appears on the second side of the bubble.
    // All other matrix elements are set to 0.
    const uint64_t orientedReadCount = num_vertices(phasingGraph);
    vector<double> A(orientedReadCount * bubbleCount, 0.);

    // To fill the matrix, loop over oriented reads.
    i = 0;
    BGL_FORALL_VERTICES(v, phasingGraph, PhasingGraph) {
        const OrientedReadId orientedReadId = phasingGraph[v].orientedReadId;

        // Loop over the bubbles that this OrientedReadId participates in.
        const vector <pair <uint64_t, uint64_t> >& orientedReadBubbles =
            orientedReadsTable[orientedReadId.getValue()];
        for(const pair <uint64_t, uint64_t>& p: orientedReadBubbles) {
            const uint64_t bubbleId = p.first;
            if(bubbleGraph.vertexTable[bubbleId] == BubbleGraph::null_vertex()) {
                continue; // This is a bad bubble.
            }
            const uint64_t side = p.second;

            // Find the index of this bubble in this connected component.
            const auto it = bubbleMap.find(bubbleId);
            SHASTA_ASSERT(it != bubbleMap.end());
            const auto j = it->second;
            SHASTA_ASSERT(j < bubbleCount);

            // Fill in the matrix element.
            double& matrixElement = A[j * orientedReadCount + i];
            if(side == 0) {
                matrixElement = 1.;
            } else if (side == 1) {
                matrixElement = -1.;
            } else {
                SHASTA_ASSERT(0);
            }
        }
        i++;
    }

    // Compute the SVD.
    const int M = int(orientedReadCount);
    const int N = int(bubbleCount);
    const string JOBU = "A";
    const string JOBVT = "A";
    const int LDA = M;
    vector<double> S(min(M, N));
    vector<double> U(M*M);
    const int LDU = M;
    vector<double> VT(N*N);
    const int LDVT = N;
    const int LWORK = 10 * max(M, N);
    vector<double> WORK(LWORK);
    int INFO = 0;
    cout << timestamp << "Computing the svd. Matrix is " << M << " by " << N << endl;
    dgesvd_(
        JOBU.data(), JOBVT.data(),
        M, N,
        &A[0], LDA, &S[0], &U[0], LDU, &VT[0], LDVT, &WORK[0], LWORK, INFO);
    cout << timestamp << "Done computing the svd." << endl;
    if(INFO != 0) {
        throw runtime_error("Error " + to_string(INFO) +
            " computing SVD decomposition of local read graph.");
    }
    for(uint64_t i=0; i<min(10UL, S.size()); i++) {
        cout << S[i] << endl;
    }

    // Get the phase from the sign of the components of the first left eigenvector.
    i = 0;
    BGL_FORALL_VERTICES(v, phasingGraph, PhasingGraph) {
        PhasingGraphVertex& vertex = phasingGraph[v];
        const double eigenvectorComponent = U[i++];
        vertex.eigenvectorComponent = eigenvectorComponent;

        if(eigenvectorComponent >= 0.) {
            vertex.phase = 1;
        } else {
            vertex.phase = -1;
        }
    }
}

void Bubbles::createPhasingGraph(
    const vector<OrientedReadId>& orientedReadIds,
    PhasingGraph& phasingGraph) const
{
    // Create the vertices and the vertexMap.
    phasingGraph.createVertices(orientedReadIds);

    // Create the edges.
    BGL_FORALL_VERTICES(v0, phasingGraph, PhasingGraph) {
        const OrientedReadId orientedReadId0 = phasingGraph[v0].orientedReadId;

        // Find the neighbors of v0 - that is, the ones that
        // have at least one bubble in common.
        vector<OrientedReadId> orientedReadIds1;
        findNeighborOrientedReadIds(orientedReadId0, orientedReadIds1);

        for(const OrientedReadId orientedReadId1: orientedReadIds1) {

            // Don't add an edge twice.
            if(orientedReadId1 <= orientedReadId0) {
                continue;
            }

            // Locate the vertex corresponding to orientedReadId1.
            auto it1 = phasingGraph.vertexMap.find(orientedReadId1);
            SHASTA_ASSERT(it1 != phasingGraph.vertexMap.end());
            const PhasingGraph::vertex_descriptor v1 = it1->second;

            uint64_t sameSideCount = 0;
            uint64_t oppositeSideCount = 0;
            findOrientedReadsRelativePhase(
                orientedReadId0, orientedReadId1,
                sameSideCount, oppositeSideCount);

            PhasingGraph::edge_descriptor e;
            bool edgeWasAdded = false;
            tie(e, edgeWasAdded) = add_edge(v0, v1, phasingGraph);
            SHASTA_ASSERT(edgeWasAdded);

            PhasingGraphEdge& edge = phasingGraph[e];
            edge.sameSideCount = sameSideCount;
            edge.oppositeSideCount = oppositeSideCount;
        }
    }

}

void Bubbles::PhasingGraph::createVertices(
    const vector<OrientedReadId>& orientedReadIds)
{
    PhasingGraph& g = *this;

    for(const OrientedReadId orientedReadId: orientedReadIds) {
        const vertex_descriptor v = add_vertex(g);
        g[v].orientedReadId = orientedReadId;
        vertexMap.insert(make_pair(orientedReadId, v));
    }
}

// Given a connected component of the BubbleGraph,
// find the OrientedReadIds that appear in it.
// The OrientedReadIds are returned sorted.
void Bubbles::findComponentOrientedReads(
    const vector<BubbleGraph::vertex_descriptor>& component,
    vector<OrientedReadId>& orientedReadIds
    ) const
{
#if 0
    cout << "Bubbles in this component:";
    for(const BubbleGraph::vertex_descriptor v: component) {
        const uint64_t bubbleId = bubbleGraph[v].bubbleId;
        cout << " " << bubbleId;
    }
    cout << endl;
#endif

    std::set<OrientedReadId> orientedReadIdsSet;
    for(const BubbleGraph::vertex_descriptor v: component) {
        const uint64_t bubbleId = bubbleGraph[v].bubbleId;
        const Bubble& bubble = bubbles[bubbleId];
        if(bubble.isBad) {
            continue;
        }
        for(uint64_t side=0; side<2; side++) {
            for(const OrientedReadId orientedReadId: bubble.orientedReadIds[side]) {
                orientedReadIdsSet.insert(orientedReadId);
            }
        }
    }

    orientedReadIds.clear();
    copy(orientedReadIdsSet.begin(), orientedReadIdsSet.end(),
        back_inserter(orientedReadIds));

#if 0
    cout << "OrientedReadIds in this component:" << endl;
    for(const OrientedReadId orientedReadId: orientedReadIds) {
        cout << orientedReadId << " ";
    }
    cout << endl;
#endif
}

void Bubbles::BubbleGraph::computeConnectedComponents()
{
    using boost::connected_components;
    using boost::get;
    using boost::color_map;

    using G = BubbleGraph;
    G& g = *this;

    connected_components(g,
        get(&BubbleGraphVertex::componentId, g),
        color_map(get(&BubbleGraphVertex::color, g)));

    // Gather the vertices in each connected component.
    std::map<uint64_t, vector<vertex_descriptor> > componentMap;
    BGL_FORALL_VERTICES(v, g, G) {
        componentMap[g[v].componentId].push_back(v);
    }
    connectedComponents.clear();
    for(const auto& p: componentMap) {
        connectedComponents.push_back(p.second);
    }
}

// Functions used to decide if an alignment should be used.
bool Bubbles::allowAlignment(
    const OrientedReadPair& orientedReadPair,
    bool useClustering) const
{
    const OrientedReadId orientedReadId0(orientedReadPair.readIds[0], 0);
    const OrientedReadId orientedReadId1(orientedReadPair.readIds[1], orientedReadPair.isSameStrand ? 0 : 1);

    OrientedReadId orientedReadId0Rc = orientedReadId0;
    OrientedReadId orientedReadId1Rc = orientedReadId1;
    orientedReadId0Rc.flipStrand();
    orientedReadId1Rc.flipStrand();

    return
        allowAlignment(orientedReadId0,   orientedReadId1, useClustering) and
        allowAlignment(orientedReadId0Rc, orientedReadId1Rc, useClustering);
}

bool Bubbles::allowAlignment(
    const OrientedReadId& orientedReadId0,
    const OrientedReadId& orientedReadId1,
    bool useClustering) const
{
    if(useClustering) {
        return allowAlignmentUsingClustering(orientedReadId0, orientedReadId1);
    } else {
        return allowAlignmentUsingBubbles(orientedReadId0, orientedReadId1);
    }
}

bool Bubbles::allowAlignmentUsingClustering(
    const OrientedReadId& orientedReadId0,
    const OrientedReadId& orientedReadId1) const
{
    const pair<uint64_t, uint64_t>& p0 = orientedReadsPhase[orientedReadId0.getValue()];
    const pair<uint64_t, uint64_t>& p1 = orientedReadsPhase[orientedReadId1.getValue()];

    const uint64_t componentId0 = p0.first;
    const uint64_t componentId1 = p1.first;
    const uint64_t unphased = std::numeric_limits<uint64_t>::max();
    const bool isPhased0 = not(componentId0 == unphased);
    const bool isPhased1 = not(componentId1 == unphased);

    // If both unphased, allow the alignment.
    if((not isPhased0) and (not isPhased1)) {
        return true;
    }

    // If only one is phased, don't allow the alignment.
    // This will have to be reviewed to allow alignments with unphased reads near
    // the beginning/end of a phased region.
    if(isPhased0 and (not isPhased1)) {
        return false;
    }
    if(isPhased1 and (not isPhased0)) {
        return false;
    }

    // If getting here, they are both phased.
    SHASTA_ASSERT(isPhased0 and isPhased1);

    // Only allow the alignment if they are in the same component and phase.
    return
        (componentId0 == componentId1) and (p0.second == p1.second);
}

bool Bubbles::allowAlignmentUsingBubbles(
    const OrientedReadId& orientedReadId0,
    const OrientedReadId& orientedReadId1) const
{
    uint64_t sameSideCount = 0;
    uint64_t oppositeSideCount = 0;
    findOrientedReadsRelativePhase(
        orientedReadId0, orientedReadId1,
        sameSideCount, oppositeSideCount);

    // Allow alignments in homozygous regions.
    if((sameSideCount==0) and (oppositeSideCount==0)) {
        return true;
    }

    // Thresholds to be exposed as options when code stabilizes.
    const uint64_t minSameSideCount = 6;
    const uint64_t maxOppositeSideCount = 1;

    return
        (sameSideCount >= minSameSideCount) and
        (oppositeSideCount <= maxOppositeSideCount);
}

max_stars_repo_name: rlorigro/shasta

max_issues_repo_name: rlorigro/shasta

max_forks_repo_name: rlorigro/shasta


text:
#include <test/test-models/good/variational/univariate_no_constraint.hpp>
#include <stan/variational/advi.hpp>
#include <stan/callbacks/stream_logger.hpp>
#include <gtest/gtest.h>
#include <test/unit/util.hpp>
#include <vector>
#include <string>
#include <iostream>
#include <boost/random/additive_combine.hpp> // L'Ecuyer RNG

typedef boost::ecuyer1988 rng_t;

class advi_test : public ::testing::Test {
public:
  advi_test()
    : logger(log_stream_, log_stream_, log_stream_, log_stream_, log_stream_),
      parameter_writer(parameter_stream_),
      diagnostic_writer(diagnostic_stream_) { }

  void SetUp() {
    err_msg1 = "Informational Message: The maximum number of iterations is reached!";
    err_msg2 = "Informational Message: The ELBO at a previous iteration is larger than the ELBO upon convergence!";

    // Create mock data_var_context
    static const std::string DATA = "";
    std::stringstream data_stream(DATA);
    stan::io::dump data_var_context(data_stream);
    // data_stream.close();

    model_ = new stan_model(data_var_context, &model_stream_);
    cont_params_ = Eigen::VectorXd::Zero(model_->num_params_r());
    base_rng_.seed(3021828106u);
    model_stream_.str("");
    log_stream_.str("");
    parameter_stream_.str("");
    diagnostic_stream_.str("");

    advi_meanfield_ = new stan::variational::advi<stan_model, stan::variational::normal_meanfield, rng_t>
      (*model_, cont_params_, base_rng_,
       1, 100,
       1, 1);
    advi_fullrank_ = new stan::variational::advi<stan_model, stan::variational::normal_fullrank, rng_t>
      (*model_, cont_params_, base_rng_,
       1, 100,
       1, 1);
  }

  void TearDown() {
    delete advi_meanfield_;
    delete advi_fullrank_;
    delete model_;
  }

  std::string err_msg1;
  std::string err_msg2;

  stan::variational::advi<stan_model, stan::variational::normal_meanfield, rng_t> *advi_meanfield_;
  stan::variational::advi<stan_model, stan::variational::normal_fullrank, rng_t> *advi_fullrank_;
  std::stringstream model_stream_;
  std::stringstream log_stream_;
  std::stringstream parameter_stream_;
  std::stringstream diagnostic_stream_;
  stan::callbacks::stream_logger logger;
  stan::callbacks::stream_writer parameter_writer;
  stan::callbacks::stream_writer diagnostic_writer;

private:
  stan_model *model_;
  rng_t base_rng_;
  Eigen::VectorXd cont_params_;
};

TEST_F(advi_test, max_iteration_warn_meanfield) {
  EXPECT_EQ(0, advi_meanfield_->run(10, 0, 50, 0.01, 1,
                                    logger,
                                    parameter_writer,
                                    diagnostic_writer));
  EXPECT_TRUE(log_stream_.str().find(err_msg1) != std::string::npos)
    << "The message should have err_msg1 inside it.";
}

TEST_F(advi_test, max_iteration_warn_fullrank) {
  EXPECT_EQ(0, advi_fullrank_->run(10, 0, 50, 0.01, 1,
                                   logger,
                                   parameter_writer,
                                   diagnostic_writer));
  EXPECT_TRUE(log_stream_.str().find(err_msg1) != std::string::npos)
    << "The message should have err_msg1 inside it.";
}

TEST_F(advi_test, prev_elbo_larger_meanfield) {
  EXPECT_EQ(0, advi_meanfield_->run(10, 0, 50, 0.1, 100,
                                    logger,
                                    parameter_writer,
                                    diagnostic_writer));
  EXPECT_TRUE(log_stream_.str().find(err_msg2) != std::string::npos)
    << "The message should have err_msg2 inside it.";
}

TEST_F(advi_test, prev_elbo_larger_fullrank) {
  EXPECT_EQ(0, advi_fullrank_->run(10, 0, 50, 0.2, 100,
                                   logger,
                                   parameter_writer,
                                   diagnostic_writer));
  EXPECT_TRUE(log_stream_.str().find(err_msg2) != std::string::npos)
    << "The message should have err_msg2 inside it.";
}

ext: cpp

max_stars_repo_name: alashworth/stan-monorepo

max_issues_repo_name: alashworth/stan-monorepo

max_forks_repo_name: alashworth/stan-monorepo


text:
/*
 * Copyright 2019 GridGain Systems, Inc. and Contributors.
 *
 * Licensed under the GridGain Community Edition License (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     https://www.gridgain.com/products/software/community-edition/gridgain-community-edition-license
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

#include <boost/test/unit_test.hpp>

#include <ignite/ignition.h>

#include <ignite/thin/ignite_client_configuration.h>
#include <ignite/thin/ignite_client.h>

#include <test_utils.h>

using namespace ignite::thin;
using namespace boost::unit_test;

class IgniteClientTestSuiteFixture
{
public:
    IgniteClientTestSuiteFixture()
    {
        serverNode = ignite_test::StartCrossPlatformServerNode("cache.xml", "ServerNode");
    }

    ~IgniteClientTestSuiteFixture()
    {
        ignite::Ignition::StopAll(false);
    }

private:
    /** Server node. */
    ignite::Ignite serverNode;
};

BOOST_FIXTURE_TEST_SUITE(IgniteClientTestSuite, IgniteClientTestSuiteFixture)

BOOST_AUTO_TEST_CASE(IgniteClientConnection)
{
    IgniteClientConfiguration cfg;

    cfg.SetEndPoints("127.0.0.1:11110");

    IgniteClient::Start(cfg);
}

BOOST_AUTO_TEST_CASE(IgniteClientConnectionFailover)
{
    IgniteClientConfiguration cfg;

    cfg.SetEndPoints("127.0.0.1:11109..11111");

    IgniteClient::Start(cfg);
}

BOOST_AUTO_TEST_SUITE_END()

max_stars_repo_name: FedorUporov/gridgain

max_issues_repo_name: FedorUporov/gridgain

max_forks_repo_name: FedorUporov/gridgain


text:
/*********************************************************************
 * BSD 3-Clause License
 *
 * Copyright (c) 2020 Northwestern University
 * 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 the copyright holder 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 HOLDER 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.
 *********************************************************************/
/**
 * @file collision.hpp
 * @author Boston Cleek
 * @date 27 Oct 2020
 * @brief Collision checking in 2D
 */
#ifndef COLLISION_HPP
#define COLLISION_HPP

#include <utility>
#include <unordered_map>
#include <tuple>
#include <limits>
#include <armadillo>

#include <ergodic_exploration/grid.hpp>

namespace ergodic_exploration
{
using arma::vec;
using std::tuple;

/** @brief State of the collision detector */
enum class CollisionMsg
{
  crash,
  obstacle,
  none
};

/** @brief Collision detection parameters */
struct CollisionConfig
{
  /**
   * @brief Constructor
   * @param rb - collision boundary
   * @param rc - collision threshold
   * @param rm - max search distance
   * @param cx - x-coordinate (jth column) center of circle
   * @param cy - y-coordinate (ith row) center of circle
   */
  CollisionConfig(int rb, int rc, int rm, int cx, int cy)
    : r_bnd(rb), r_col(rc), r_max(rm), cx(cx), cy(cy), dx(0), dy(0), sqrd_obs(-1)
  {
  }

  const int r_bnd, r_col, r_max;  // bounding , collision, and max search radii
  const int cx, cy;               // circle center
  int dx, dy;                     // displacement from circle center to obstacle
  int sqrd_obs;                   // squared distance circle center to obstacle
};

/** @brief 2D collision detection */
class Collision
{
public:
  /**
   * @brief Constructor
   * @param boundary_radius - collision boundary around robot represented as a circle
   * @param search_radius - radius outside of collision boundary to search for obstacles
   * @param obstacle_threshold - distance from boundary radius to obstacle to be considered a collision
   * @param occupied_threshold - min probaility [0 1] for a cell to be considered an obstacle
   */
  Collision(double boundary_radius, double search_radius, double obstacle_threshold,
            double occupied_threshold);

  /**
   * @brief Compose distance to closest obstacle
   * @param grid - grid map
   * @param pose - robot state [x, y, theta]
   * @return state of collision detector and distance to closest obstacle
   * @details if collision distance to obstacle is zero and if there are no
   * obstacles distance is max numerical value for a double
   */
  tuple<CollisionMsg, double> minDistance(const GridMap& grid, const vec& pose) const;

  /**
   * @brief Compose displacement vector from closest obstacle to robot
   * @param grid - grid map
   * @param pose - robot state [x, y, theta]
   * @return state of collision detector and vector from closest obstacle
   * @details if collision distance to vector is all zeros and if there are no
   * obstacles the vector contains the max numerical value for a double
   */
  tuple<CollisionMsg, vec> minDirection(const GridMap& grid, const vec& pose) const;

  /**
   * @brief Check if the given state is a collision
   * @param grid - grid map
   * @param pose - robot state [x, y, theta]
   * @return true if collision
   */
  bool collisionCheck(const GridMap& grid, const vec& pose) const;

  /** @brief Return the distance from robot center to collision boundary */
  double totalPadding() const;

private:
  /**
   * @brief Find occupied cells between collision boundary and max search radius
   * @param cfg[out] - collision configuration
   * @param grid - grid map
   * @return true if there is a collision
   */
  bool search(CollisionConfig& cfg, const GridMap& grid) const;

  /**
   * @brief Bresenham's circle algorithm
   * @param cfg[out] - collision configuration
   * @param grid - grid map
   * @param r - radius of circle
   * @return true if there is a collision
   * @details Perfroms Bresenham's circle algorithm to detect obstacle cells
   */
  bool bresenhamCircle(CollisionConfig& cfg, const GridMap& grid, int r) const;

  /**
   * @brief Check if cell is an obstacle
   * @param cfg[out] - collision configuration
   * @param grid - grid map
   * @param cj - x-coordinate (jth column) of cell on perimeter of circle
   * @param ci - y-coordinate (ith row) of cell on perimeter of circle
   * @return true if there is a collision
   */
  bool checkCell(CollisionConfig& cfg, const GridMap& grid, unsigned int cj,
                 unsigned int ci) const;

private:
  double boundary_radius_;     // circular radius around robot
  double search_radius_;       // search radius for obstacles
  double obstacle_threshold_;  // collision distance threshold
  double occupied_threshold_;  // probaility cell is occupied
};
}  // namespace ergodic_exploration
#endif

ext: hpp


text:
#ifndef __GRAPH__TYPES__H__
#define __GRAPH__TYPES__H__

#include <boost/graph/adjacency_list.hpp>
#include <boost/graph/graph_traits.hpp>

typedef boost::adjacency_list < boost::vecS, boost::vecS, boost::undirectedS > Graph;
typedef boost::graph_traits<Graph>::edge_descriptor Edge;
typedef boost::graph_traits<Graph>::vertex_descriptor Vertex;

#endif

ext: hh

max_stars_repo_name: JoeyEremondi/treewidth-memoization

max_issues_repo_name: JoeyEremondi/treewidth-memoization

max_forks_repo_name: JoeyEremondi/treewidth-memoization


text:
/*  Meru
 *  GraphicsResourceTexture.cpp
 *
 *  Copyright (c) 2009, Stanford University
 *  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 "CDNArchive.hpp"
#include "GraphicsResourceManager.hpp"
#include "GraphicsResourceTexture.hpp"
#include "ResourceDependencyTask.hpp"
#include "ResourceLoadTask.hpp"
#include "ResourceLoadingQueue.hpp"
#include "ResourceUnloadTask.hpp"
#include "SequentialWorkQueue.hpp"
#include <boost/bind.hpp>
#include <OgreResourceBackgroundQueue.h>

namespace Meru {

class TextureDependencyTask : public ResourceDependencyTask
{
public:
  TextureDependencyTask(DependencyManager* mgr, WeakResourcePtr resource, const String& hash);
  virtual ~TextureDependencyTask();

  virtual void run();
};

class TextureLoadTask : public ResourceLoadTask
{
public:
  TextureLoadTask(DependencyManager *mgr, SharedResourcePtr resource, const String &hash, unsigned int epoch);

  virtual void doRun();

protected:
  unsigned int mArchiveName;
};

class TextureUnloadTask : public ResourceUnloadTask
{
public:
  TextureUnloadTask(DependencyManager *mgr, WeakResourcePtr resource, const String &hash, unsigned int epoch);

  virtual void doRun();

protected:
  //bool mainThreadUnload(String name);
};

GraphicsResourceTexture::GraphicsResourceTexture(const RemoteFileId &resourceID)
  : GraphicsResourceAsset(resourceID, GraphicsResource::TEXTURE)
{

}

GraphicsResourceTexture::~GraphicsResourceTexture()
{
  if (mLoadState == LOAD_LOADED)
    doUnload();
}

ResourceDownloadTask* GraphicsResourceTexture::createDownloadTask(DependencyManager *manager, ResourceRequestor *resourceRequestor)
{
  return new ResourceDownloadTask(manager, mResourceID, resourceRequestor);
}

ResourceDependencyTask* GraphicsResourceTexture::createDependencyTask(DependencyManager *manager)
{
  return new TextureDependencyTask(manager, getWeakPtr(), mResourceID.toString());
}

ResourceLoadTask* GraphicsResourceTexture::createLoadTask(DependencyManager *manager)
{
  return new TextureLoadTask(manager, getSharedPtr(), mResourceID.toString(), mLoadEpoch);
}

ResourceUnloadTask* GraphicsResourceTexture::createUnloadTask(DependencyManager *manager)
{
  return new TextureUnloadTask(manager, getWeakPtr(), mResourceID.toString(), mLoadEpoch);
}

/***************************** TEXTURE DEPENDENCY TASK *************************/

TextureDependencyTask::TextureDependencyTask(DependencyManager *mgr, WeakResourcePtr resource, const String& hash)
  : ResourceDependencyTask(mgr, resource, hash)
{

}

TextureDependencyTask::~TextureDependencyTask()
{

}

void TextureDependencyTask::run()
{
  SharedResourcePtr resourcePtr = mResource.lock();
  if (!resourcePtr) {
    signalCompletion(false);
    return;
  }

  resourcePtr->setCost(mBuffer.size());
  resourcePtr->parsed(true);

  signalCompletion();
}

/***************************** TEXTURE LOAD TASK *************************/

TextureLoadTask::TextureLoadTask(DependencyManager *mgr, SharedResourcePtr resourcePtr, const String &hash, unsigned int epoch)
: ResourceLoadTask(mgr, resourcePtr, hash, epoch)
{
}

void TextureLoadTask::doRun()
{
  mArchiveName = CDNArchive::addArchive(CDNArchive::canonicalMhashName(mHash), mBuffer);
  Ogre::TextureManager::getSingleton().load(CDNArchive::canonicalMhashName(mHash), Ogre::ResourceGroupManager::DEFAULT_RESOURCE_GROUP_NAME);
  CDNArchive::removeArchive(mArchiveName);
  mResource->loaded(true, mEpoch);
}

/***************************** TEXTURE UNLOAD TASK *************************/

TextureUnloadTask::TextureUnloadTask(DependencyManager *mgr, WeakResourcePtr resource, const String &hash, unsigned int epoch)
: ResourceUnloadTask(mgr, resource, hash, epoch)
{

}
/*
bool TextureUnloadTask::mainThreadUnload(String name)
{
  Ogre::TextureManager::getSingleton().unload(name);
  operationCompleted(Ogre::BackgroundProcessTicket(), Ogre::BackgroundProcessResult());
  return false;
}
*/
void TextureUnloadTask::doRun()
{
  /*I REALLY wish this were true*/
  //    SequentialWorkQueue::getSingleton().queueWork(std::tr1::bind(&TextureUnloadTask::mainThreadUnload, this, mHash));

  Ogre::TextureManager* textureManager = Ogre::TextureManager::getSingletonPtr();
  textureManager->remove(mHash);

  Ogre::ResourcePtr textureResource = textureManager->getByName(mHash);
  assert(textureResource.isNull());

  SharedResourcePtr resource = mResource.lock();
  if (resource)
    resource->unloaded(true, mEpoch);
}

}

max_stars_repo_name: jordanmisiura/sirikata

max_issues_repo_name: jordanmisiura/sirikata

max_forks_repo_name: jordanmisiura/sirikata


text:
//| This file is a part of the sferes2 framework.
//| Copyright 2009, ISIR / Universite Pierre et Marie Curie (UPMC)
//| Main contributor(s): Jean-Baptiste Mouret, mouret@isir.fr
//|
//| This software is a computer program whose purpose is to facilitate
//| experiments in evolutionary computation and evolutionary robotics.
//|
//| This software is governed by the CeCILL license under French law
//| and abiding by the rules of distribution of free software.  You
//| can use, modify and/ or redistribute the software under the terms
//| of the CeCILL license as circulated by CEA, CNRS and INRIA at the
//| following URL "http://www.cecill.info".
//|
//| As a counterpart to the access to the source code and rights to
//| copy, modify and redistribute granted by the license, users are
//| provided only with a limited warranty and the software's author,
//| the holder of the economic rights, and the successive licensors
//| have only limited liability.
//|
//| In this respect, the user's attention is drawn to the risks
//| associated with loading, using, modifying and/or developing or
//| reproducing the software by the user in light of its specific
//| status of free software, that may mean that it is complicated to
//| manipulate, and that also therefore means that it is reserved for
//| developers and experienced professionals having in-depth computer
//| knowledge. Users are therefore encouraged to load and test the
//| software's suitability as regards their requirements in conditions
//| enabling the security of their systems and/or data to be ensured
//| and, more generally, to use and operate it in the same conditions
//| as regards security.
//|
//| The fact that you are presently reading this means that you have
//| had knowledge of the CeCILL license and that you accept its terms.

#ifndef _STAT_HPP_
#define _STAT_HPP_

#include <fstream>
#include <string>
#include <boost/shared_ptr.hpp>
#include <sferes/stc.hpp>

namespace sferes {
  namespace stat {
    template<typename Phen, typename Params, typename Exact = stc::Itself>
    class Stat {
     public:
      template<typename E>
      void refresh(const E& ea) {
        assert(!ea.pop().empty());
        stc::exact(this)->refresh(ea);
      }
      void show(std::ostream& os, size_t k) {
      }
      template<class Archive>
      void serialize(Archive & ar, const unsigned int version) {
      }
     protected:
      boost::shared_ptr<std::ofstream> _log_file;
      template<typename E>
      void _create_log_file(const E& ea, const std::string& name) {
        if (!_log_file && ea.dump_enabled()) {
          std::string log = ea.res_dir() + "/" + name;
          _log_file = boost::shared_ptr<std::ofstream>(new std::ofstream(log.c_str()));
        }
      }
    };

  }
}

#define SFERES_STAT(Class, Parent)					\
  template <typename Phen, typename Params, typename Exact = stc::Itself> \
  class Class : public Parent<Phen, Params, typename stc::FindExact<Class<Phen, Params, Exact>, Exact>::ret>

#define SFERES_STAT_PARENT(Class, Parent)				\
  Parent<Phen, Params,  typename stc::FindExact<Class<Phen, Params, Exact>, Exact>::ret>

#endif

max_stars_repo_name: Evolving-AI-Lab/innovation-engine

max_issues_repo_name: Evolving-AI-Lab/innovation-engine

max_forks_repo_name: Evolving-AI-Lab/innovation-engine


text:
//
// $Id$
//
//
// Original author: Matt Chambers <matt.chambers .@. vanderbilt.edu>
//
// Copyright 2010 Vanderbilt University - Nashville, TN 37232
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
//

#ifndef _PROTEINLIST_DECOYGENERATOR_HPP_
#define _PROTEINLIST_DECOYGENERATOR_HPP_

#include "pwiz/utility/misc/Export.hpp"
#include "pwiz/data/proteome/ProteinListWrapper.hpp"
#include <boost/cstdint.hpp>

namespace pwiz {
namespace analysis {

using namespace pwiz::proteome;

/// ProteinList decoy generator for creating decoy proteins on the fly
class PWIZ_API_DECL ProteinList_DecoyGenerator : public proteome::ProteinListWrapper
{
    public:

    /// client-implemented generator predicate -- called during construction of
    /// ProteinList_DecoyGenerator to create a decoy protein from a target protein
    struct PWIZ_API_DECL Predicate
    {
        /// return a decoy protein based on an input target protein
        virtual ProteinPtr generate(const Protein& protein) const = 0;

        /// return the string prefixed to a protein id to indicate it is a decoy
        virtual const std::string& decoyPrefix() const {return decoyPrefix_;}

        virtual ~Predicate() {}

        protected:
        std::string decoyPrefix_;
    };

    typedef boost::shared_ptr<Predicate> PredicatePtr;

    ProteinList_DecoyGenerator(const ProteinListPtr& original, const PredicatePtr& predicate);

    /// \name ProteinList interface
    //@{
    virtual size_t size() const;
    virtual size_t find(const std::string& id) const;
    virtual ProteinPtr protein(size_t index, bool getSequence = true) const;
    //@}

    private:
    struct Impl;
    boost::shared_ptr<Impl> impl_;
    ProteinList_DecoyGenerator(const ProteinList_DecoyGenerator&);
    ProteinList_DecoyGenerator& operator=(const ProteinList_DecoyGenerator&);
};

/// creates a reversed copy of every target protein with the specified decoy string prefixed to the id
class PWIZ_API_DECL ProteinList_DecoyGeneratorPredicate_Reversed : public ProteinList_DecoyGenerator::Predicate
{
    public:
    ProteinList_DecoyGeneratorPredicate_Reversed(const std::string& decoyPrefix);

    virtual ProteinPtr generate(const Protein& protein) const;
};

/// creates a randomly shuffled copy of every target protein with the specified decoy string prefixed to the id
class PWIZ_API_DECL ProteinList_DecoyGeneratorPredicate_Shuffled : public ProteinList_DecoyGenerator::Predicate
{
    public:
    ProteinList_DecoyGeneratorPredicate_Shuffled(const std::string& decoyPrefix, boost::uint32_t randomSeed = 0u);

    virtual ProteinPtr generate(const Protein& protein) const;

    private:
    struct Impl;
    boost::shared_ptr<Impl> impl_;
};

} // namespace analysis
} // namespace pwiz

#endif // _PROTEINLIST_DECOYGENERATOR_HPP_

max_stars_repo_name: austinkeller/pwiz

max_issues_repo_name: austinkeller/pwiz

max_forks_repo_name: austinkeller/pwiz


text:
/** @addtogroup Network
 *  @{
 */

/**
 * Represent an outgoing file transfer.
 * @class TransferReceiver
 */

#ifndef TRANSFER_RECEIVER_HPP
# define TRANSFER_RECEIVER_HPP

#include <iostream>
#include <string>
#include <functional>
#include <boost/filesystem.hpp>
#include <boost/filesystem/fstream.hpp>

#include <network/message.hpp>

class TLSSocket;

template <typename SocketType>
class MessageHandler;

class TransferReceiver
{
public:
  TransferReceiver(MessageHandler<TLSSocket>*, const std::string&, const std::string&, int,
                   std::function<void(const TransferReceiver*)>&&);
  ~TransferReceiver();
  void stop();

private:
  TransferReceiver(const TransferReceiver&);
  TransferReceiver& operator=(const TransferReceiver&);

  void get_next_chunk(Message*);

  MessageHandler<TLSSocket>* client;
  std::string id;
  boost::filesystem::path filename;
  int length;
  int received_length;
  std::function<void(const TransferReceiver*)> end_callback;
  boost::filesystem::ofstream file;
};

#endif // TRANSFER_RECEIVER_HPP

/**@}*/

max_stars_repo_name: louiz/batajelo

max_issues_repo_name: louiz/batajelo

max_forks_repo_name: louiz/batajelo


text: #include <boost/type_index/runtime_cast/pointer_cast.hpp>

max_stars_repo_name: miathedev/BoostForArduino

max_issues_repo_name: miathedev/BoostForArduino

max_forks_repo_name: miathedev/BoostForArduino


text:
#include "nmos/logging_api.h"

#include <boost/range/adaptor/transformed.hpp>
#include "nmos/api_utils.h"
#include "nmos/query_utils.h"
#include "nmos/slog.h"
#include "rql/rql.h"

namespace nmos
{
    namespace experimental
    {
        web::http::experimental::listener::api_router make_unmounted_logging_api(nmos::experimental::log_model& model, slog::base_gate& gate);

        web::http::experimental::listener::api_router make_logging_api(nmos::experimental::log_model& model, slog::base_gate& gate)
        {
            using namespace web::http::experimental::listener::api_router_using_declarations;

            api_router logging_api;

            logging_api.support(U("/?"), methods::GET, [](http_request req, http_response res, const string_t&, const route_parameters&)
            {
                set_reply(res, status_codes::OK, nmos::make_sub_routes_body({ U("log/") }, req, res));
                return pplx::task_from_result(true);
            });

            logging_api.support(U("/log/?"), methods::GET, [](http_request req, http_response res, const string_t&, const route_parameters&)
            {
                set_reply(res, status_codes::OK, nmos::make_sub_routes_body({ U("v1.0/") }, req, res));
                return pplx::task_from_result(true);
            });

            logging_api.mount(U("/log/v1.0"), make_unmounted_logging_api(model, gate));

            return logging_api;
        }

        bool match_logging_rql(const web::json::value& value, const web::json::value& query)
        {
            return query.is_null() || rql::evaluator
            {
                [&value](web::json::value& results, const web::json::value& key)
                {
                    return web::json::extract(value.as_object(), results, key.as_string());
                },
                rql::default_any_operators()
            }(query) == rql::value_true;
        }

        // Predicate to match events against a query
        struct log_event_query
        {
            typedef const nmos::experimental::log_event& argument_type;
            typedef bool result_type;

            log_event_query(const web::json::value& flat_query_params);

            result_type operator()(argument_type event) const;

            // the query/exemplar object for a Basic Query
            web::json::value basic_query;

            // a representation of the RQL abstract syntax tree for an Advanced Query
            web::json::value rql_query;
        };

        log_event_query::log_event_query(const web::json::value& flat_query_params)
            : basic_query(web::json::unflatten(flat_query_params))
        {
            if (basic_query.has_field(U("paging")))
            {
                basic_query.erase(U("paging"));
            }
            if (basic_query.has_field(U("query")))
            {
                auto& advanced = basic_query.at(U("query")).as_object();
                for (auto& field : advanced)
                {
                    if (field.first == U("rql"))
                    {
                        rql_query = rql::parse_query(field.second.as_string());
                    }
                    // an error is reported for unimplemented parameters
                    else
                    {
                        throw std::runtime_error("unimplemented parameter - query." + utility::us2s(field.first));
                    }
                }
                basic_query.erase(U("query"));
            }
        }

        log_event_query::result_type log_event_query::operator()(argument_type event) const
        {
            return web::json::match_query(event.data, basic_query, web::json::match_icase | web::json::match_substr)
                && match_logging_rql(event.data, rql_query);
        }

        // Cursor-based paging parameters
        struct log_event_paging
        {
            explicit log_event_paging(const web::json::value& flat_query_params, const nmos::tai& max_until = nmos::tai_max(), size_t default_limit = (std::numeric_limits<size_t>::max)(), size_t max_limit = (std::numeric_limits<size_t>::max)());

            // determine if the range [until, since) and limit are valid
            bool valid() const;

            // return only the results up until the time specified (inclusive)
            nmos::tai until;

            // return only the results since the time specified (non-inclusive)
            nmos::tai since;

            // restrict the response to the specified number of results
            size_t limit;

            // where both 'since' and 'until' parameters are specified, the 'since' value takes precedence
            // where a resulting data set is constrained by the server's value of 'limit'
            bool since_specified;

            template <typename Predicate>
            boost::any_range<const nmos::experimental::log_event, boost::bidirectional_traversal_tag, const nmos::experimental::log_event&, std::ptrdiff_t> page(const nmos::experimental::log_events& events, Predicate match)
            {
                return paging::cursor_based_page(events, match, until, since, limit, !since_specified);
            }
        };

        log_event_paging::log_event_paging(const web::json::value& flat_query_params, const nmos::tai& max_until, size_t default_limit, size_t max_limit)
            : until(max_until)
            , since(nmos::tai_min())
            , limit(default_limit)
            , since_specified(false)
        {
            auto query_params = web::json::unflatten(flat_query_params);

            // extract the supported paging options
            if (query_params.has_field(U("paging")))
            {
                auto& paging = query_params.at(U("paging")).as_object();
                for (auto& field : paging)
                {
                    if (field.first == U("until"))
                    {
                        until = nmos::parse_version(field.second.as_string());
                        // "These parameters may be used to construct a URL which would return the same set of bounded data on consecutive requests"
                        // therefore the response to a request with 'until' in the future should be fixed up...
                        if (until > max_until) until = max_until;
                    }
                    else if (field.first == U("since"))
                    {
                        since = nmos::parse_version(field.second.as_string());
                        since_specified = true;
                        // how should a request with 'since' in the future be fixed up?
                        if (since > until && until == max_until) until = since;
                    }
                    else if (field.first == U("limit"))
                    {
                        // "If the client had requested a page size which the server is unable to honour, the actual page size would be returned"
                        limit = (size_t)field.second.as_integer();
                        if (limit > max_limit) limit = max_limit;
                    }
                    // an error is reported for unimplemented parameters
                    else
                    {
                        throw std::runtime_error("unimplemented parameter - paging." + utility::us2s(field.first));
                    }
                }
            }
        }

        bool log_event_paging::valid() const
        {
            return since <= until;
        }

        namespace details
        {
            // make user error information (to be used with status_codes::BadRequest)
            utility::string_t make_valid_paging_error(const nmos::experimental::log_event_paging& paging)
            {
                return U("the value of the 'paging.since' parameter must be less than or equal to the value of the 'paging.until' parameter");
            }
        }

        // Cursor-based paging customisation points

        inline nmos::tai extract_cursor(const nmos::experimental::log_events&, nmos::experimental::log_events::const_iterator it)
        {
            return it->cursor;
        }

        inline nmos::experimental::log_events::const_iterator lower_bound(const nmos::experimental::log_events& index, const nmos::tai& cursor)
        {
            return std::lower_bound(index.begin(), index.end(), cursor, [](const nmos::experimental::log_event& event, const nmos::tai& cursor)
            {
                return event.cursor > cursor;
            });
        }

        namespace details
        {
            utility::string_t make_query_parameters(web::json::value flat_query_params)
            {
                // any non-string query parameters need serializing before encoding
                // but web::json::query_from_value does this...

                // special case, RQL is kept as the URI-encoded string
                const auto encoded_rql = nmos::fields::query_rql(flat_query_params);

                // all other string values need encoding
                nmos::details::encode_elements(flat_query_params);

                if (flat_query_params.has_field(nmos::fields::query_rql))
                {
                    flat_query_params[nmos::fields::query_rql] = web::json::value::string(encoded_rql);
                }

                return web::json::query_from_value(flat_query_params);
            }

            web::json::value parse_query_parameters(const utility::string_t& query)
            {
                auto flat_query_params = web::json::value_from_query(query);

                // special case, RQL is kept as the URI-encoded string
                const auto encoded_rql = nmos::fields::query_rql(flat_query_params);

                // all other string values need decoding
                nmos::details::decode_elements(flat_query_params);

                if (flat_query_params.has_field(nmos::fields::query_rql))
                {
                    flat_query_params[nmos::fields::query_rql] = web::json::value::string(encoded_rql);
                }

                // any non-string query parameters need parsing after decoding...
                if (flat_query_params.has_field(nmos::fields::paging_limit))
                {
                    flat_query_params[nmos::fields::paging_limit] = web::json::value::parse(nmos::fields::paging_limit(flat_query_params));
                }

                return flat_query_params;
            }

            web::uri make_query_uri_with_no_paging(const web::http::http_request& req, const nmos::settings& settings)
            {
                // could rebuild the query parameters from the decoded and parsed query string, rather than individually deleting the paging parameters from the request?
                auto query_params = parse_query_parameters(req.request_uri().query());
                if (query_params.has_field(U("paging.since")))
                {
                    query_params.erase(U("paging.since"));
                }
                if (query_params.has_field(U("paging.until")))
                {
                    query_params.erase(U("paging.until"));
                }
                if (query_params.has_field(U("paging.limit")))
                {
                    query_params.erase(U("paging.limit"));
                }

                // RFC 5988 allows relative URLs, but NMOS specification examples are all absolute URLs
                // See https://tools.ietf.org/html/rfc5988#section-5
                // See https://github.com/AMWA-TV/nmos-discovery-registration/blob/v1.2/docs/2.5.%20APIs%20-%20Query%20Parameters.md

                // get the request host and port (or use the primary host address, and port, from settings)
                auto req_host_port = web::http::get_host_port(req);
                if (req_host_port.first.empty())
                {
                    req_host_port.first = nmos::get_host(settings);
                }
                if (0 == req_host_port.second)
                {
                    req_host_port.second = nmos::experimental::fields::logging_port(settings);
                }

                return web::uri_builder()
                    .set_scheme(nmos::http_scheme(settings)) // no means to detect the API protocol from the request unless reverse proxy added X-Forwarded-Proto?
                    .set_host(req_host_port.first)
                    .set_port(req_host_port.second)
                    .set_path(req.request_uri().path())
                    .set_query(make_query_parameters(query_params))
                    .to_uri();
            }

            void add_paging_headers(web::http::http_headers& headers, const nmos::experimental::log_event_paging& paging, const web::uri& base_link)
            {
                // X-Paging-Limit "identifies the current limit being used for paging. This may not match the requested value if the requested value was too high for the implementation"
                headers.add(U("X-Paging-Limit"), paging.limit);

                // X-Paging-Since "identifies the current value of the query parameter 'paging.since' in use, or if not specified identifies what value it would have had to return this data set.
                // This value may be re-used as the paging.until value of a query to return the previous page of results."
                headers.add(U("X-Paging-Since"), make_version(paging.since));

                // X-Paging-Until "identifies the current value of the query parameter 'paging.until' in use, or if not specified identifies what value it would have had to return this data set.
                // This value may be re-used as the paging.since value of a query to return the next page of results."
                headers.add(U("X-Paging-Until"), make_version(paging.until));

                // Link header "provides references to cursors for paging. The 'rel' attribute may be one of 'next', 'prev', 'first' or 'last'"

                auto link = web::uri_builder(base_link)
                    .append_query(U("paging.limit=") + utility::ostringstreamed(paging.limit))
                    .to_string();

                // "The Link header identifies the 'next' and 'prev' cursors which an application may use to make its next requests."
                // (though note, both the 'next' and 'prev' cursors may return empty responses)

                headers.add(U("Link"), U("<") + link + U("&paging.until=") + make_version(paging.since) + U(">; rel=\"prev\""));

                headers.add(U("Link"), U("<") + link + U("&paging.since=") + make_version(paging.until) + U(">; rel=\"next\""));

                // "An implementation may also provide 'first' and 'last' cursors to identify the beginning and end of a collection of this can be addressed via a consistent URL."

                headers.add(U("Link"), U("<") + link + U("&paging.since=") + make_version(nmos::tai_min()) + U(">; rel=\"first\""));

                headers.add(U("Link"), U("<") + link + U(">; rel=\"last\""));
            }
        }

        web::http::experimental::listener::api_router make_unmounted_logging_api(nmos::experimental::log_model& model, slog::base_gate& gate_)
        {
            using namespace web::http::experimental::listener::api_router_using_declarations;

            api_router logging_api;

            logging_api.support(U("/?"), methods::GET, [](http_request req, http_response res, const string_t&, const route_parameters&)
            {
                set_reply(res, status_codes::OK, nmos::make_sub_routes_body({ U("events/") }, req, res));
                return pplx::task_from_result(true);
            });

            logging_api.support(U("/events/?"), methods::GET, [&model, &gate_](http_request req, http_response res, const string_t&, const route_parameters& parameters)
            {
                nmos::api_gate gate(gate_, req, parameters);
                auto lock = model.read_lock();

                // Extract and decode the query string

                auto flat_query_params = details::parse_query_parameters(req.request_uri().query());

                // Configure the query predicate

                log_event_query match(flat_query_params);

                // Configure the paging parameters

                // Use the paging limit (default and max) from the setings
                log_event_paging paging(flat_query_params, tai_now(), (size_t)nmos::experimental::fields::logging_paging_default(model.settings), (size_t)nmos::experimental::fields::logging_paging_limit(model.settings));

                if (paging.valid())
                {
                    // Get the payload and update the paging parameters
                    struct default_constructible_event_query_wrapper { const log_event_query* impl; bool operator()(const log_event& e) const { return (*impl)(e); } };
                    auto page = paging.page(model.events, default_constructible_event_query_wrapper{ &match });

                    size_t count = 0;

                    set_reply(res, status_codes::OK,
                        web::json::serialize(page,
                            [&count](const log_event& event){ ++count; return event.data; }),
                        web::http::details::mime_types::application_json);

                    slog::log<slog::severities::too_much_info>(gate, SLOG_FLF) << "Returning " << count << " matching log events";

                    details::add_paging_headers(res.headers(), paging, details::make_query_uri_with_no_paging(req, model.settings));
                }
                else
                {
                    set_error_reply(res, status_codes::BadRequest, U("Bad Request; ") + details::make_valid_paging_error(paging));
                }

                return pplx::task_from_result(true);
            });

            logging_api.support(U("/events/?"), methods::DEL, [&model](http_request req, http_response res, const string_t&, const route_parameters&)
            {
                auto lock = model.write_lock();

                if (req.request_uri().query().empty())
                {
                    model.events.clear();
                    set_reply(res, status_codes::NoContent);
                }
                else
                {
                    set_reply(res, status_codes::NotImplemented);
                }

                return pplx::task_from_result(true);
            });

            logging_api.support(U("/events/") + nmos::patterns::resourceId.pattern + U("/?"), methods::GET, [&model](http_request, http_response res, const string_t&, const route_parameters& parameters)
            {
                auto lock = model.read_lock();

                const string_t eventId = parameters.at(nmos::patterns::resourceId.name);

                auto& by_id = model.events.get<tags::id>();
                auto event = by_id.find(eventId);
                if (by_id.end() != event)
                {
                    set_reply(res, status_codes::OK, event->data);
                }
                else
                {
                    set_reply(res, status_codes::NotFound);
                }

                return pplx::task_from_result(true);
            });

            return logging_api;
        }
    }
}

max_stars_repo_name: rufusu-sony/nmos-cpp

max_issues_repo_name: rufusu-sony/nmos-cpp

max_forks_repo_name: rufusu-sony/nmos-cpp