reduce.hh

Go to the documentation of this file.

// This file is part of Awali.
// Copyright 2016-2022 Sylvain Lombardy, Victor Marsault, Jacques Sakarovitch
//
// Awali is a free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
 
 
#ifndef AWALI_ALGOS_REDUCE_HH
# define AWALI_ALGOS_REDUCE_HH
 
# include <map>
# include <unordered_map>
# include <vector>
# include <cmath>
 
#include <awali/sttc/algos/copy.hh>
#include <awali/sttc/algos/transpose.hh>
#include <awali/sttc/core/mutable_automaton.hh>
#include <awali/sttc/weightset/q.hh>
#include <awali/sttc/weightset/r.hh>
#include <awali/sttc/weightset/z.hh>
 
namespace awali { namespace sttc {
 
  namespace internal
  {
    /*
      The core algorithm computes a scaled basis from a list of
      vectors.  There is a permutation on entries in such a way that
      the first non-zero entry of the i-th vector of the basis has
      index i.  (in the field case, this entry is equal to 1)
 
      If v is a new vector, for every vector b_i of the bases:
 
      - either v[i] is zero and there is nothing to do
 
      - or v[i] is different from 0 in this case, v is modified as v:=
        v - v[i].b which makes v[i]=0 (in the Z case, b may also be
        modified) This is performed by the reduce_vector method
 
      After the iteration, if v is null, it means that the vector was in the
      vector space (or the Z-module) generated by the basis.
      Otherwise, a non zero entry of v is selected (find_pivot),
      v is normalized (normalization_vector) and the permutation is
      updated in such a way that the pivot becomes the first  non zero
      entry of v, which is inserted in the basis.
 
      Once the basis is stabilized, a bottom-up reduction is applied
      to the scaled basis in order to make it more "diagonal".
      Finally, the new automaton is built w.r.t. the computed basis.
      The construction of the automaton can not be done on-the-fly due
      to both the Z case where the basis is modified and the bottom-up
      reduction.
    */
 
 
    /*
 
     The implementation of elementary methods slightly differ
     w.r.t. the weightset.
 
     Generically,
 
       find_pivot returns the first non zero entry
 
       reduce_vector computes current := current - current[i].b_i
 
       normalisation_vector applied a ratio in such a way that
       the pivot is 1
 
       bottom_up_reduction applies, for each vector of the
       basis, reduce_vector to every preceding vector in the
       basis.
 
 
       In Q : find_pivot searchs for the entry where
       |numerator|+|denominator| is minimal
 
       In R : find_pivot searchs for the entry x where |x|+1/|x|
       is minimal
 
       In Z : find_pivot searchs for the (non zero) entry x where
       |x| is minimal
       reduce_vector both reduce the current and the basis vector
         (w.r.t. the gcd of current[i] and b_i[i])
       normalisation_vector does not apply
       bottom_up_reduction does not apply
 
    */
 
    template<typename Weightset>
    struct select
    {
      template<typename Reduc, typename Vector>
      static unsigned
      find_pivot(Reduc* that, const Vector& v,
                 unsigned begin, unsigned* permutation)
      {
        return that->find_pivot(v, begin, permutation);
      }
 
      template<typename Reduc, typename Vector>
      static void
      reduce_vector(Reduc* that, Vector& vbasis,
                    Vector& current, unsigned b, unsigned* permutation)
      {
        that->reduce_vector(vbasis, current, b, permutation);
      }
 
      template<typename Reduc, typename Vector>
      static void
      normalisation_vector(Reduc* that, Vector& v,
                           unsigned pivot, unsigned* permutation)
      {
        that->normalisation_vector(v, pivot, permutation);
      }
 
      template<typename Reduc, typename Basis>
      static void
      bottom_up_reduction(Reduc* that, Basis& basis, unsigned* permutation)
      {
        that->bottom_up_reduction(basis, permutation);
      }
 
      template<typename Reduc, typename Basis, typename Vector>
      static void
      vector_in_new_basis(Reduc* that, Basis& basis,
                          Vector& current, Vector& new_vector,
                          unsigned* permutation)
      {
        that->vector_in_new_basis(basis, current, new_vector, permutation);
      }
    };
 
    template <>
    struct select<q> : select<void>
    {
      template<typename Reduc, typename Vector>
      static unsigned
      find_pivot(Reduc* that, const Vector& v,
                 unsigned begin, unsigned* permutation)
      {
        return that->find_pivot_by_norm(v, begin, permutation);
      }
    };
 
    template <>
    struct select<r> : select<void>
    {
      template<typename Reduc, typename Vector>
      static unsigned
      find_pivot(Reduc* that, const Vector& v,
                 unsigned begin, unsigned* permutation)
      {
        return that->find_pivot_by_norm(v, begin, permutation);
      }
    };
 
    template <>
    struct select<z> : select<void>
    {
      template<typename Reduc, typename Vector>
      static unsigned
      find_pivot(Reduc* that, const Vector& v,
                 unsigned begin, unsigned* permutation)
      {
        return that->find_pivot_by_norm(v, begin, permutation);
      }
 
      template <typename Reduc, typename Vector>
      static void
      reduce_vector(Reduc* that, Vector& vbasis,
                    Vector& current, unsigned b, unsigned* permutation)
      {
        that->z_reduce_vector(vbasis, current, b, permutation);
      }
 
      template <typename Reduc, typename Vector>
      static void normalisation_vector(Reduc*, Vector&, unsigned, unsigned*)
      {}
 
      template<typename Reduc, typename Basis>
      static void bottom_up_reduction(Reduc*, Basis&, unsigned*)
      {}
 
      template<typename Reduc, typename Basis, typename Vector>
      static void
      vector_in_new_basis(Reduc* that, Basis& basis,
                          Vector& current, Vector& new_vector,
                          unsigned* permutation)
      {
        that->z_vector_in_new_basis(basis, current, new_vector, permutation);
      }
    };
 
    template <typename Aut, typename AutOutput>
    class reductioner
    {
      static_assert(context_t_of<Aut>::is_lal,
                    "requires free labelset");
 
      using automaton_t = Aut;
      using context_t = context_t_of<automaton_t>;
      using weightset_t = typename context_t::weightset_t;
      using output_automaton_t = AutOutput;
      using label_t = label_t_of<automaton_t>;
      using weight_t = typename context_t::weight_t;
      using vector_t = std::vector<weight_t>;
      using matrix_t = std::vector<std::map<std::size_t, weight_t> > ;
      using matrix_set_t = std::map<label_t, matrix_t>;
 
    public:
      reductioner(const automaton_t& input)
        : input_(input)
        , res_(make_shared_ptr<output_automaton_t>(input_->context()))
      {}
 
      void linear_representation()
      {
        std::unordered_map<state_t, unsigned> state_to_index;
        unsigned i = 0;
        for (auto s: input_->states())
          state_to_index[s] = i++;
        dimension = i;
        if (dimension == 0)
          return;
        init.resize(i);
        final.resize(i);
        // Computation of the initial vector.
        for (auto t : input_->initial_transitions())
          init[state_to_index[input_->dst_of(t)]] = input_->weight_of(t);
        // Computation of the final vector.
        for (auto t : input_->final_transitions())
          final[state_to_index[input_->src_of(t)]] = input_->weight_of(t);
        // For each letter, we define an adjency matrix.
        for (auto t : input_->transitions())
          {
            auto it = letter_matrix_set.find(input_->label_of(t));
            if (it == letter_matrix_set.end())
              it = letter_matrix_set.emplace(input_->label_of(t),
                                             matrix_t(dimension)).first;
            it->second[state_to_index[input_->src_of(t)]]
              [state_to_index[input_->dst_of(t)]] = input_->weight_of(t);
          }
      }
 
      //utility methods
 
      void product_vector_matrix(const vector_t& v,
                                 const matrix_t& m,
                                 vector_t& res)
      {
        for (unsigned i = 0; i < dimension; i++)
          for (auto it : m[i])
            {
              unsigned j = it.first;
              res[j] = ws_.add(res[j], ws_.mul(v[i], it.second));
            }
      }
 
      weight_t scalar_product(const vector_t& v,
                              const vector_t& w)
      {
        weight_t res = ws_.zero();
        for (unsigned i = 0; i < dimension; ++i)
          res = ws_.add(res, ws_.mul(v[i], w[i]));
        return res;
      }
 
      using z_weight_t = sttc::z::value_t; // int or long
      using q_weight_t = sttc::q::value_t;
      using r_weight_t = sttc::r::value_t; // = double
 
      /*
        The pivot is the entry x of the vector such that norm(x) is minimal.
        This "norm" highly depends on the semiring.
 
        For rational numbers, it is the sum of the numerator and the
        denominator (absolute value)
 
        For "real" numbers, it is |x|+|1/x| .
      */
 
      static unsigned norm(const q_weight_t& w)
      {
        return w.den+abs(w.num);
      }
 
      static double norm(const r_weight_t& w)
      {
        return std::fabs(w)+std::fabs(1/w);
      }
 
      static unsigned norm(const z_weight_t& w)
      {
        return abs(w);
      }
 
      // Works for both Q and R.
      unsigned
      find_pivot_by_norm(const vector_t& v, unsigned begin,
                         unsigned* permutation)
      {
        unsigned i=begin;
        for(; i < dimension && ws_.is_zero(v[permutation[i]]); ++i)
          ;
        if(i==dimension)
          return dimension;
        unsigned pivot=i;
        auto min = norm(v[permutation[i++]]);
        for (; i < dimension; ++i)
          if (!ws_.is_zero(v[permutation[i]])
              && (pivot == dimension
                  || norm(v[permutation[i]])<min))
            {
              pivot = i;
              min = norm(v[permutation[i]]);
            }
//      std::cerr << v[permutation[pivot]] << std::endl;
        return pivot;
      }
 
      // End of Q/R specializations.
 
 
      // Specialization for Z.
 
      // Gcd function that also computes the Bezout coefficients.
      // Used in the z reduction.
      static z_weight_t
      gcd(z_weight_t x, z_weight_t y, z_weight_t& a, z_weight_t& b)
      {
        //z_weight_t res = 0;
        //gcd = ax + by
        int sx=(x<0)?-1:1, sy=(y<0)?-1:1;
        bool inv=false;
        if(x<0) x=-x;
        if(y<0) y=-y;
        if(x>y){
          z_weight_t t=x; x=y; y=t;
          inv=true;
        }
        a=1; b=0;
        int a0=0, b0=1;
        while(y!=0) {
          z_weight_t q=x/y, r=x%y;
          z_weight_t a1=a0, b1=b0;
          a0=a-q*a0; b0=b-q*b0;
          a=a1;      b=b1;
          x=y;  y=r;
        }
        if(inv) {
          int c=a; a=b; b=c;
        }
        a*=sx; b*=sy;
        return x;
      }
 
      /*
        gcd= a.vbasis[pivot] + b.current[pivot]
        current[pivot] is made zero with a unimodular linear transformation
        This also modifies the basis vector
        [ vbasis  ] := [        a                      b       ] . [ vbasis   ]
        [ current ]    [-current[pivot]/gcd   vbasis[pivot]/gcd]   [ current ]
      */
      void z_reduce_vector(vector_t& vbasis, vector_t& current,
                           unsigned nb, unsigned* permutation)
      {
        unsigned pivot = permutation[nb]; //pivot of vector vbasis
        if (ws_.is_zero(current[pivot]))
          return;
        weight_t bp = vbasis[pivot];
        weight_t cp = current[pivot];
        weight_t a,b;
        weight_t g = gcd(bp, cp, a, b);
        bp /= g;
        cp /= g;
        for (unsigned i = nb; i < dimension; ++i)
          {
            weight_t tmp = current[permutation[i]];
            current[permutation[i]] = bp*tmp - cp*vbasis[permutation[i]];
            vbasis[permutation[i]] = a*vbasis[permutation[i]] + b*tmp;
          }
      }
 
      void z_vector_in_new_basis(std::vector<vector_t>& basis,
                                 vector_t& current, vector_t& new_vector,
                                 unsigned* permutation)
      {
        for (unsigned b = 0; b < basis.size(); ++b)
          {
            vector_t& vbasis = basis[b];
            unsigned pivot = permutation[b]; //pivot of vector vbasis
            new_vector[b] = current[pivot] / vbasis[pivot];
            if (!ws_.is_zero(new_vector[b]))
              {
                current[pivot] = ws_.zero();
                for (unsigned i = b+1; i < dimension; ++i)
                  current[permutation[i]]
                    -= new_vector[b] * vbasis[permutation[i]];
              }
          }
      }
      // End of Z specializations.
 
     /* Generic subroutines.
         These methods are written for any (skew) field.
         Some are specialized for Q and R for stability issues.
         Some are specialized for Z which is an Euclidean domain.
      */
 
      unsigned find_pivot(const vector_t& v, unsigned begin,
                          unsigned* permutation)
      {
        for (unsigned i = begin; i < dimension; ++i)
          if (!ws_.is_zero(v[permutation[i]]))
            return i;
        return dimension;
      }
 
      weight_t reduce_vector(vector_t& vbasis,
                             vector_t& current, unsigned b,
                             unsigned* permutation)
      {
        unsigned pivot = permutation[b]; //pivot of vector vbasis
        weight_t ratio = current[pivot];//  vbasis[pivot] is one
        if (ws_.is_zero(ratio))
          return ratio;
        // This is safer than current[p] = current[p]-ratio*vbasis[p];
        current[pivot] = ws_.zero();
        for (unsigned i = b+1; i < dimension; ++i)
          current[permutation[i]]
            = ws_.sub(current[permutation[i]],
                             ws_.mul(ratio, vbasis[permutation[i]]));
        return ratio;
      }
 
      void normalisation_vector(vector_t& v, unsigned pivot,
                                unsigned* permutation)
      {
        weight_t k = v[permutation[pivot]];
        if (!ws_.is_one(k))
          {
            v[permutation[pivot]] = ws_.one();
            for (unsigned r = pivot + 1; r < dimension; ++r)
              v[permutation[r]] = ws_.rdiv(v[permutation[r]], k);
          }
      }
 
      void bottom_up_reduction(std::vector<vector_t>& basis,
                               unsigned* permutation)
      {
        for (unsigned b = basis.size()-1; 0 < b; --b)
          for (unsigned c = 0; c < b; ++c)
            reduce_vector(basis[b], basis[c], b, permutation);
      }
 
      void vector_in_new_basis(std::vector<vector_t>& basis,
                               vector_t& current, vector_t& new_vector,
                               unsigned* permutation)
      {
        for (unsigned b = 0; b < basis.size(); ++b)
          new_vector[b] = reduce_vector(basis[b], current, b, permutation);
      }
 
    public:
      void left_reduce()
      {
        // Used to select the proper overload depending on weightset_t.
        using type_t = select<weightset_t>;
 
        linear_representation();
        // The basis is a list of vectors, each vector is associated with
        // a state of the output
        std::vector<vector_t> basis;
        // The permutation array corresponds to a permutation of the indices
        // (i.e. the columns of the matrices)
        // such that permtuation[k] is the pivot of the k-th vector of the basis.
        unsigned* permutation = new unsigned[dimension];
        for (unsigned i = 0; i < dimension; ++i)
          permutation[i] = i;
        // If the initial vector is null, the function immediatly returns
 
        // A non zero entry is chosen as pivot
        unsigned pivot = type_t::find_pivot(this, init, 0, permutation);
        if (pivot == dimension) //all components of init are 0
          return;
        // The pivot of the first basis vector is permutation[0];
        permutation[0] = pivot;
        permutation[pivot] = 0;
        // The initial vector is the first element of the new basis
        // (up to the normalisation w.r.t the pivot)
        vector_t first(init);
        type_t::normalisation_vector(this, first, 0, permutation);
        basis.push_back(first);
        // To each vector of the basis, all the successor vectors are
        // computed, reduced to respect to the basis, and finally, if
        // linearly independant, pushed at the end of the basis
        // itself.
        for (unsigned nb = 0; nb < basis.size(); ++nb)
          // All the vectors basis[nb].mu(a) are processed
          for (auto mu : letter_matrix_set) //mu is a pair (letter,matrix)
            {
              vector_t current(dimension);
              product_vector_matrix(basis[nb], mu.second, current);
              //reduction of current w.r.t each basis vector;
              for (unsigned b = 0; b < basis.size(); ++b)
                type_t::reduce_vector(this, basis[b],
                                      current, b, permutation);
              // After reduction, we put current in the basis if it is
              // not null and we search for the pivot of current.
              pivot = type_t::find_pivot(this, current, basis.size(),
                                         permutation);
              if (pivot != dimension) //otherwise, current is null
                {
                  if (pivot != basis.size())
                    std::swap(permutation[pivot], permutation[basis.size()]);
                  type_t::normalisation_vector(this, current,
                                               basis.size(), permutation);
                  basis.push_back(current);
                }
            }
 
        // now, we use each vector to reduce the preceding vectors in
        // the basis.  If weightset=Z we do not do it.
        type_t::bottom_up_reduction(this, basis, permutation);
 
        // Construction of the output automaton
        // 1. States
        std::vector<state_t> states(basis.size());
        for (unsigned b = 0; b < basis.size(); ++b)
          states[b] = res_->add_state();
        // 2. Initial vector
        vector_t vect_new_basis(basis.size());
        type_t::vector_in_new_basis(this, basis, init,
                                    vect_new_basis, permutation);
        for (unsigned b = 0; b < basis.size(); ++b)
          res_->set_initial(states[b], vect_new_basis[b]);
        // 3. Each vector of the basis is a state; computation of the
        // final function and the successor function
        for (unsigned v = 0; v < basis.size(); ++v)
          {
            weight_t k = scalar_product(basis[v],final);
            if(!ws_.is_zero(k))
              res_->set_final(states[v],k);
            for (auto mu : letter_matrix_set)
              {
                // mu is a pair (letter,matrix).
                vector_t current(dimension);
                product_vector_matrix(basis[v], mu.second, current);
                type_t::vector_in_new_basis(this, basis, current,
                                            vect_new_basis, permutation);
                for (unsigned b = 0; b < basis.size(); ++b)
                  res_->new_transition(states[v], states[b],
                                       mu.first, vect_new_basis[b]);
              }
        }
        delete[] permutation;
      }
 
 
      output_automaton_t get_output() const
      {
        return res_;
      }
 
    private:
      automaton_t input_;
      const weightset_t_of<automaton_t> ws_ = *input_->weightset();
 
      output_automaton_t res_;
 
      // Linear representation of the input.
      unsigned dimension;
      vector_t init;
      vector_t final;
      matrix_set_t letter_matrix_set;
 
    };
 
  }
 
  template<typename Aut>
  Aut reduce(const Aut& input)
  {
    auto tmp = transpose_view(input);
    internal::reductioner<decltype(tmp), Aut> algo(tmp);
    algo.left_reduce();
    auto tmp2=transpose_view(algo.get_output());
    internal::reductioner<decltype(tmp2), Aut> algo2(tmp2);
    algo2.left_reduce();
    auto ret=copy(algo2.get_output());
    if(ret->num_states() >= input->num_states())
      ret= copy(input);
    if(!input->get_name().empty()) {
      ret->set_desc("Reduction of "+input->get_name());
      ret->set_name("red-"+input->get_name());
    }
    else {
      ret->set_desc("Reduction");
      ret->set_name("red");
    }
    return ret;
  }
 
  template<typename Aut>
  Aut left_reduce(const Aut& input)
  {
    internal::reductioner<decltype(input), Aut> algo(input);
    algo.left_reduce();
    return copy(algo.get_output());
  }
 
}}//end of ns awali::stc
 
#endif