<commit_before>#ifndef VG_CLUSTER_HPP_INCLUDED
#define VG_CLUSTER_HPP_INCLUDED

#include "gcsa/gcsa.h"
#include "gcsa/lcp.h"
#include "position.hpp"
#include "gssw_aligner.hpp"
#include "utility.hpp"
#include "mem.hpp"
#include "xg.hpp"

#include <functional>
#include <string>
#include <vector>
#include <map>


/**
 * \file cluster.hpp
 *
 * Chaining and clustering tools to work with maximal exact matches.
 */

namespace vg {
    
using namespace std;

// Prime numbers spaced at approximately logarithmic intervals
static constexpr size_t spaced_primes[62] = {2ull, 5ull, 13ull, 29ull, 53ull, 127ull, 227ull, 487ull, 967ull, 2039ull, 4093ull, 8191ull, 16381ull, 32749ull, 65521ull, 131071ull, 262139ull, 524287ull, 1048573ull, 2097143ull, 4194301ull, 8388593ull, 16777213ull, 33554393ull, 67108859ull, 134217689ull, 268435399ull, 536870909ull, 1073741789ull, 2147483647ull, 4294967291ull, 8589934583ull, 17179869143ull, 34359738337ull, 68719476731ull, 137438953447ull, 274877906899ull, 549755813881ull, 1099511627689ull, 2199023255531ull, 4398046511093ull, 8796093022151ull, 17592186044399ull, 35184372088777ull, 70368744177643ull, 140737488355213ull, 281474976710597ull, 562949953421231ull, 1125899906842597ull, 2251799813685119ull, 4503599627370449ull, 9007199254740881ull, 18014398509481951ull, 36028797018963913ull, 72057594037927931ull, 144115188075855859ull, 288230376151711717ull, 576460752303423433ull, 1152921504606846883ull, 2305843009213693951ull, 4611686018427387847ull, 9223372036854775783ull};

// Precomputed primitive roots of unity paired with these primes (chosen randomly from the 20 smallest roots)
static constexpr size_t primitive_roots_of_unity[62] = {1ull, 3ull, 2ull, 21ull, 27ull, 56ull, 17ull, 45ull, 40ull, 28ull, 69ull, 70ull, 40ull, 31ull, 119ull, 75ull, 42ull, 61ull, 60ull, 46ull, 21ull, 13ull, 39ull, 13ull, 29ull, 15ull, 29ull, 32ull, 37ull, 73ull, 56ull, 45ull, 13ull, 90ull, 51ull, 12ull, 32ull, 11ull, 39ull, 24ull, 8ull, 39ull, 7ull, 51ull, 38ull, 67ull, 2ull, 34ull, 62ull, 19ull, 13ull, 30ull, 12ull, 45ull, 31ull, 57ull, 6ull, 57ull, 3ull, 37ull, 68ull, 54ull};

/**
 * Iterate over pairsets of integers in a pseudorandom but deterministic order.
 * We use the same permutation every time for a given number of items to pair
 * up.
 */
class ShuffledPairs {
public:
    /**
     * Make a new iterable pairing up the given number of items.
     */
    ShuffledPairs(size_t num_items);
    
    /**
     * Actual iterator class.
     */
    class iterator {
    public:
        /**
         * Advance to the next pair.
         */
        iterator& operator++();
        
        /**
         * Get the pair pointed to.
         */
        pair<size_t, size_t> operator*() const;
        
        /**
         * see if two iterators are equal.
         */
        bool operator==(const iterator& other) const;
        
        /**
         * see if two iterators are not equal.
         */
        bool operator!=(const iterator& other) const;
        
        friend class ShuffledPairs;
        
        // Default copy constructor and assignment operator.
        iterator(const iterator& other) = default;
        iterator& operator=(const iterator& other) = default;
        
    private:
        // What is the ordinal value of this element in the permutation?
        size_t permutation_idx;
        // Which permutation does this iterator mean?
        size_t permuted;
        // What are we iterating over?
        const ShuffledPairs& iteratee;
        
        // Make an iterator. Only the friend parent can do it.
        iterator(const ShuffledPairs& iteratee, size_t start_at);
    };
    
    using const_iterator = iterator;
    
    /**
     * Get an iterator to the first pair.
     */
    iterator begin() const;
    /**
     * Get an iterator to the past-the-end pair.
     */
    iterator end() const;
    
private:

    // How many items are we working with to pair up?
    size_t num_items;
    // How many pairs do we make?
    size_t num_pairs;
    // A prime number that is larger than the number of pairs
    size_t larger_prime;
    // A primitive root of unity for the prime
    size_t primitive_root;
};

class MEMChainModelVertex {
public:
    MaximalExactMatch mem;
    vector<pair<MEMChainModelVertex*, double> > next_cost; // for forward
    vector<pair<MEMChainModelVertex*, double> > prev_cost; // for backward
    double weight;
    double score;
    MEMChainModelVertex* prev;
    MEMChainModelVertex(void) = default;                                      // Copy constructor
    MEMChainModelVertex(const MEMChainModelVertex&) = default;               // Copy constructor
    MEMChainModelVertex(MEMChainModelVertex&&) = default;                    // Move constructor
    MEMChainModelVertex& operator=(const MEMChainModelVertex&) & = default;  // MEMChainModelVertexopy assignment operator
    MEMChainModelVertex& operator=(MEMChainModelVertex&&) & = default;       // Move assignment operator
    virtual ~MEMChainModelVertex() { }                     // Destructor
};

class MEMChainModel {
public:
    vector<MEMChainModelVertex> model;
    map<string, map<int64_t, vector<vector<MEMChainModelVertex>::iterator> > > positions;
    set<vector<MEMChainModelVertex>::iterator> redundant_vertexes;
    MEMChainModel(
        const vector<size_t>& aln_lengths,
        const vector<vector<MaximalExactMatch> >& matches,
        const function<int64_t(pos_t)>& approx_position,
        const function<map<string, vector<size_t> >(pos_t)>& path_position,
        const function<double(const MaximalExactMatch&, const MaximalExactMatch&)>& transition_weight,
        int band_width = 10,
        int position_depth = 8,
        int max_connections = 20);
    void score(const set<MEMChainModelVertex*>& exclude);
    MEMChainModelVertex* max_vertex(void);
    vector<vector<MaximalExactMatch> > traceback(int alt_alns, bool paired, bool debug);
    void display(ostream& out);
    void clear_scores(void);
};<commit_msg>allowing more than one MEM hit at a given approximate position very slightly hurts our precision<commit_after>#ifndef VG_CLUSTER_HPP_INCLUDED
#define VG_CLUSTER_HPP_INCLUDED

#include "gcsa/gcsa.h"
#include "gcsa/lcp.h"
#include "position.hpp"
#include "gssw_aligner.hpp"
#include "utility.hpp"
#include "mem.hpp"
#include "xg.hpp"

#include <functional>
#include <string>
#include <vector>
#include <map>


/**
 * \file cluster.hpp
 *
 * Chaining and clustering tools to work with maximal exact matches.
 */

namespace vg {
    
using namespace std;

// Prime numbers spaced at approximately logarithmic intervals
static constexpr size_t spaced_primes[62] = {2ull, 5ull, 13ull, 29ull, 53ull, 127ull, 227ull, 487ull, 967ull, 2039ull, 4093ull, 8191ull, 16381ull, 32749ull, 65521ull, 131071ull, 262139ull, 524287ull, 1048573ull, 2097143ull, 4194301ull, 8388593ull, 16777213ull, 33554393ull, 67108859ull, 134217689ull, 268435399ull, 536870909ull, 1073741789ull, 2147483647ull, 4294967291ull, 8589934583ull, 17179869143ull, 34359738337ull, 68719476731ull, 137438953447ull, 274877906899ull, 549755813881ull, 1099511627689ull, 2199023255531ull, 4398046511093ull, 8796093022151ull, 17592186044399ull, 35184372088777ull, 70368744177643ull, 140737488355213ull, 281474976710597ull, 562949953421231ull, 1125899906842597ull, 2251799813685119ull, 4503599627370449ull, 9007199254740881ull, 18014398509481951ull, 36028797018963913ull, 72057594037927931ull, 144115188075855859ull, 288230376151711717ull, 576460752303423433ull, 1152921504606846883ull, 2305843009213693951ull, 4611686018427387847ull, 9223372036854775783ull};

// Precomputed primitive roots of unity paired with these primes (chosen randomly from the 20 smallest roots)
static constexpr size_t primitive_roots_of_unity[62] = {1ull, 3ull, 2ull, 21ull, 27ull, 56ull, 17ull, 45ull, 40ull, 28ull, 69ull, 70ull, 40ull, 31ull, 119ull, 75ull, 42ull, 61ull, 60ull, 46ull, 21ull, 13ull, 39ull, 13ull, 29ull, 15ull, 29ull, 32ull, 37ull, 73ull, 56ull, 45ull, 13ull, 90ull, 51ull, 12ull, 32ull, 11ull, 39ull, 24ull, 8ull, 39ull, 7ull, 51ull, 38ull, 67ull, 2ull, 34ull, 62ull, 19ull, 13ull, 30ull, 12ull, 45ull, 31ull, 57ull, 6ull, 57ull, 3ull, 37ull, 68ull, 54ull};

/**
 * Iterate over pairsets of integers in a pseudorandom but deterministic order.
 * We use the same permutation every time for a given number of items to pair
 * up.
 */
class ShuffledPairs {
public:
    /**
     * Make a new iterable pairing up the given number of items.
     */
    ShuffledPairs(size_t num_items);
    
    /**
     * Actual iterator class.
     */
    class iterator {
    public:
        /**
         * Advance to the next pair.
         */
        iterator& operator++();
        
        /**
         * Get the pair pointed to.
         */
        pair<size_t, size_t> operator*() const;
        
        /**
         * see if two iterators are equal.
         */
        bool operator==(const iterator& other) const;
        
        /**
         * see if two iterators are not equal.
         */
        bool operator!=(const iterator& other) const;
        
        friend class ShuffledPairs;
        
        // Default copy constructor and assignment operator.
        iterator(const iterator& other) = default;
        iterator& operator=(const iterator& other) = default;
        
    private:
        // What is the ordinal value of this element in the permutation?
        size_t permutation_idx;
        // Which permutation does this iterator mean?
        size_t permuted;
        // What are we iterating over?
        const ShuffledPairs& iteratee;
        
        // Make an iterator. Only the friend parent can do it.
        iterator(const ShuffledPairs& iteratee, size_t start_at);
    };
    
    using const_iterator = iterator;
    
    /**
     * Get an iterator to the first pair.
     */
    iterator begin() const;
    /**
     * Get an iterator to the past-the-end pair.
     */
    iterator end() const;
    
private:

    // How many items are we working with to pair up?
    size_t num_items;
    // How many pairs do we make?
    size_t num_pairs;
    // A prime number that is larger than the number of pairs
    size_t larger_prime;
    // A primitive root of unity for the prime
    size_t primitive_root;
};

class MEMChainModelVertex {
public:
    MaximalExactMatch mem;
    vector<pair<MEMChainModelVertex*, double> > next_cost; // for forward
    vector<pair<MEMChainModelVertex*, double> > prev_cost; // for backward
    double weight;
    double score;
    MEMChainModelVertex* prev;
    MEMChainModelVertex(void) = default;                                      // Copy constructor
    MEMChainModelVertex(const MEMChainModelVertex&) = default;               // Copy constructor
    MEMChainModelVertex(MEMChainModelVertex&&) = default;                    // Move constructor
    MEMChainModelVertex& operator=(const MEMChainModelVertex&) & = default;  // MEMChainModelVertexopy assignment operator
    MEMChainModelVertex& operator=(MEMChainModelVertex&&) & = default;       // Move assignment operator
    virtual ~MEMChainModelVertex() { }                     // Destructor
};

class MEMChainModel {
public:
    vector<MEMChainModelVertex> model;
    map<string, map<int64_t, vector<vector<MEMChainModelVertex>::iterator> > > positions;
    set<vector<MEMChainModelVertex>::iterator> redundant_vertexes;
    MEMChainModel(
        const vector<size_t>& aln_lengths,
        const vector<vector<MaximalExactMatch> >& matches,
        const function<int64_t(pos_t)>& approx_position,
        const function<map<string, vector<size_t> >(pos_t)>& path_position,
        const function<double(const MaximalExactMatch&, const MaximalExactMatch&)>& transition_weight,
        int band_width = 10,
        int position_depth = 1,
        int max_connections = 20);
    void score(const set<MEMChainModelVertex*>& exclude);
    MEMChainModelVertex* max_vertex(void);
    vector<vector<MaximalExactMatch> > traceback(int alt_alns, bool paired, bool debug);
    void display(ostream& out);
    void clear_scores(void);
};
{'fraction_non_alphanumeric': 0.08935219657483247, 'fraction_numerical': 0.1272441466037892, 'mean_word_length': 3.2987197724039827, 'pattern_counts': {'":': 0, '<': 57, '<?xml version=': 0, '>': 57, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4690602', 'n_tokens_mistral': 4973, 'n_tokens_neox': 3764, 'n_words': 1065}