AhoCorasickAmbiguous.h

Go to the documentation of this file.

// --------------------------------------------------------------------------
//                   OpenMS -- Open-Source Mass Spectrometry
// --------------------------------------------------------------------------
// Copyright The OpenMS Team -- Eberhard Karls University Tuebingen,
// ETH Zurich, and Freie Universitaet Berlin 2002-2022.
//
// This software is released under a three-clause BSD license:
//  * 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 any author or any participating institution
//    may be used to endorse or promote products derived from this software
//    without specific prior written permission.
// For a full list of authors, refer to the file AUTHORS.
// --------------------------------------------------------------------------
// 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 ANY OF THE AUTHORS OR THE CONTRIBUTING
// INSTITUTIONS 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.
//
// --------------------------------------------------------------------------
// $Maintainer: Chris Bielow $
// $Authors: Chris Bielow $
// --------------------------------------------------------------------------

#pragma once

#include <OpenMS/CONCEPT/Macros.h>
#include <OpenMS/DATASTRUCTURES/String.h>

#include <cassert>
#include <functional> // for std::hash
#include <limits>
#include <queue>
#include <string>
#include <unordered_map>
#include <vector>

namespace OpenMS
{
  constexpr char const AAtoChar[28] = {
    'A', // 00 Ala Alanine
    'Y', // 01 Tyr Tyrosine
    'C', // 02 Cys Cysteine
    'D', // 03 Asp Aspartic Acid   // B
    'N', // 04 Asn Asparagine      // B
    'F', // 05 Phe Phenylalanine
    'G', // 06 Gly Glycine
    'H', // 07 His Histidine
    'I', // 08 Ile Isoleucine      // J
    'L', // 09 Leu Leucine         // J
    'K', // 10 Lys Lysine
    'W', // 11 Trp Tryptophan
    'M', // 12 Met Methionine
    'O', // 13 Pyl Pyrrolysine
    'P', // 14 Pro Proline
    'E', // 15 Glu Glutamic Acid   // Z
    'Q', // 16 Gln Glutamine       // Z
    'R', // 17 Arg Arginine
    'S', // 18 Ser Serine
    'T', // 19 Thr Threonine
    'U', // 20 Selenocysteine
    'V', // 21 Val Valine
    // ambiguous AAs start here (index: 22...25)
    'B', // 22 Aspartic Acid, Asparagine  $   // the ambAA's need to be consecutive (B,J,Z,X,$)
    'J', // 23 Leucine, Isoleucine        $
    'Z', // 24 Glutamic Acid, Glutamine   $
    'X', // 25 Unknown 
    // non-regular AA's, which are special
    '$', // 26 superAA, i.e. it models a mismatch, which can be anything, including AAAs
    '?', // 27 invalid AA (will usually be skipped) -- must be the last AA (AA::operator++ and others rely on it)
  };

  constexpr char const CharToAA[128] = {
    // ASCII char (7-bit Int with values from 0..127) --> amino acid 
    27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, // 0
    27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, // 1
  //                 $
    27, 27, 27, 27, 26, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, // 2
    27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, // 3

  //  ,  A,  B,  C,  D,  E,  F,  G,  H,  I,  J,  K,  L,  M,  N,  O,
    27, 00, 22, 02, 03, 15, 05, 06, 07,  8, 23, 10,  9, 12, 04, 13, // 4

  // P,  Q,  R,  S,  T,  U,  V,  W,  X,  Y,  Z,   ,   ,   ,   ,   ,
    14, 16, 17, 18, 19, 20, 21, 11, 25, 01, 24, 27, 27, 27, 27, 27, // 5

  //  ,  a,  b,  c,  d,  e,  f,  g,  h,  i,  j,  k,  l,  m,  n,  o,
    27, 00, 22, 02, 03, 15, 05, 06, 07,  8, 23, 10,  9, 12, 04, 13,   // 6

  // p,  q,  r,  s,  t,  u,  v,  w,  x,  y,  z,   ,   ,   ,   ,   ,
    14, 16, 17, 18, 19, 20, 21, 11, 25, 01, 24, 27, 27, 27, 27, 27, // 7
  };

  struct OPENMS_DLLAPI Hit {
    using T = uint32_t;
    Hit() = default;
    Hit(T needle_index, T needle_length, T query_pos) : needle_index(needle_index), needle_length(needle_length), query_pos(query_pos) {};
    T needle_index;
    T needle_length;
    T query_pos;
    bool operator==(const Hit& rhs) const
    {
      return needle_index == rhs.needle_index && needle_length == rhs.needle_length && query_pos == rhs.query_pos;
    }
  };

  struct OPENMS_DLLAPI AA
  {
    constexpr AA() : aa_(AA('?').aa_)
    {
    }
    
    constexpr explicit AA(const char c) : aa_(CharToAA[(unsigned char)c])
    {
    }

    constexpr uint8_t operator()() const
    {
      return aa_;
    }

    constexpr bool operator==(const AA other) const
    {
      return aa_ == other.aa_;
    }

    constexpr bool operator<=(const AA other) const
    {
      return aa_ <= other.aa_;
    }

    constexpr bool isAmbiguous() const
    {
      return aa_ >= AA('B').aa_;
    }
    
    constexpr bool isValid() const
    {
      return aa_ != AA('?').aa_;
    }

    constexpr bool isValidForPeptide() const
    {
      return aa_ <= AA('X').aa_; // '$' or '?'
    }

    constexpr AA& operator++()
    {
      ++aa_;
      assert(aa_ <= AA('?').aa_); // make sure we don't overflow
      return *this;
    }

    constexpr AA operator++(int)
    {
      AA r(*this);
      ++aa_;
      assert(aa_ <= AA('?').aa_); // make sure we don't overflow
      return r;
    }

    constexpr AA operator-(const AA rhs) const
    {
      AA r(*this);
      r.aa_ -= rhs.aa_;
      return r;
    }

  private:
    uint8_t aa_; 
  };

  class OPENMS_DLLAPI Index
  {
  public:
    using T = uint32_t;
    Index() = default;
    
    Index(T val) : i_(val) {};

    bool isInvalid() const;

    bool isValid() const;

    T operator()() const;

    bool operator==(const Index other) const;

    T& pos();

  private:
    T i_ = std::numeric_limits<T>::max(); 
  };
 } // namespace OpenMS

// this needs to go into namespace std
template<> struct std::hash<OpenMS::Index>
{
  std::size_t operator()(OpenMS::Index const& s) const noexcept
  {
    return std::hash<OpenMS::Index::T> {}(s());
  }
};

namespace OpenMS
{
  struct OPENMS_DLLAPI ACNode
  {
    ACNode() {};

    ACNode(const AA label, const uint8_t depth) : edge(label)
    {
      depth_and_hits.depth = depth;
    }

    struct DepthHits {
      DepthHits()
      {
        memset(this, 0, sizeof *this); // make sure bitfields are 0; C++20 allows {} initialization ...
      };
      uint8_t has_hit : 1; 
      // we could add another bit here to distinguish between a local hit and suffix hit, but on Windows, this slows it down
      uint8_t depth : 7; 
    };

    using ChildCountType = uint8_t; 

    Index suffix {0};        
    Index first_child {0};   
    // there is room for optimization here, by pulling edge labels into a separate vector (allows for SSE-enabled search of children)
    AA edge {0};             
    ChildCountType nr_children = 0; 
    DepthHits depth_and_hits; 
  };
  
  // forward declaration
  struct ACTrieState;

  struct OPENMS_DLLAPI ACSpawn
  {
    ACSpawn() = delete;

    ACSpawn(std::string::const_iterator query_pos, Index tree_pos, uint8_t max_aa, uint8_t max_mm, uint8_t max_prefix_loss);

    size_t textPos(const ACTrieState& state) const;

    AA nextValidAA(const ACTrieState& state);

    std::string::const_iterator it_query; 
    Index tree_pos;                       
    uint8_t max_aaa_leftover {0};         
    uint8_t max_mm_leftover {0};          
    uint8_t max_prefix_loss_leftover {0}; 
  };

  OPENMS_DLLAPI AA nextValidAA(const std::string::const_iterator end, std::string::const_iterator& it_q);

  struct OPENMS_DLLAPI ACTrieState
  {
    friend ACSpawn;
    void setQuery(const std::string& haystack);

    size_t textPos() const;

    std::string::const_iterator textPosIt() const;

    const std::string& getQuery() const;

    AA nextValidAA();

    std::vector<Hit> hits;             
    Index tree_pos;                    
    std::queue<ACSpawn> spawns;        

  private:
    std::string query_;                
    std::string::const_iterator it_q_; 
  };

  class OPENMS_DLLAPI ACTrie
  {
  public:
    ACTrie(uint32_t max_aaa = 0, uint32_t max_mm = 0);

    ~ACTrie();

    void addNeedle(const std::string& needle);

    void addNeedles(const std::vector<std::string>& needles);

    void addNeedlesAndCompress(const std::vector<std::string>& needles);

    void compressTrie();

    size_t getNeedleCount() const;

    void setMaxAAACount(const uint32_t max_aaa);

    uint32_t getMaxAAACount() const;
    
    void setMaxMMCount(const uint32_t max_mm);

    uint32_t getMaxMMCount() const;

    bool nextHits(ACTrieState& state) const;

    void getAllHits(ACTrieState& state) const;

  private:
    bool nextHitsNoClear_(ACTrieState& state) const;

    Index add_(const Index from, const AA edge);

    bool addHits_(Index i, const size_t text_pos, std::vector<Hit>& hits) const;

    bool addHitsSpawn_(Index i, const ACSpawn& spawn, const size_t text_pos, std::vector<Hit>& hits, const int current_spawn_depths) const;

    Index follow_(const Index i, const AA edge) const;

    bool followSpawn_(ACSpawn& spawn, const AA edge, ACTrieState& state) const;

    Index stepMaster_(const Index i, const AA edge, ACTrieState& state) const;

    bool stepSpawn_(ACSpawn& spawn, ACTrieState& state) const;

    void createSpawns_(const Index i, const AA fromAA, const AA toAA, ACTrieState& state, const uint32_t aaa_left, const uint32_t mm_left) const;

    void createSubSpawns_(const ACSpawn& prototype, const AA fromAA, const AA toAA, ACTrieState& state) const;

    void createMMSpawns_(const Index i, const AA except_fromAA, const AA except_toAA, const AA except_edge, ACTrieState& state, const uint32_t aaa_left, const uint32_t mm_left) const;

    void createMMSubSpawns_(const ACSpawn& prototype, const AA except_fromAA, const AA except_toAA, const AA except_edge, ACTrieState& state) const;

    Index findChildNaive_(Index parent, AA child_label);

    Index findChildBFS_(const Index parent, const AA child_label) const;

    std::vector<ACNode> trie_;  
    uint32_t needle_count_ {0}; 
    uint32_t max_aaa_ {0};      
    uint32_t max_mm_ {0};       

    std::unordered_map<Index, std::vector<uint32_t>> umap_index2needles_;
    std::unordered_map<Index, std::vector<Index>> umap_index2children_naive_;
  };

} // namespace OpenMS