sparsity_interface.hpp

/*
 *    This file is part of CasADi.
 *
 *    CasADi -- A symbolic framework for dynamic optimization.
 *    Copyright (C) 2010-2023 Joel Andersson, Joris Gillis, Moritz Diehl,
 *                            KU Leuven. All rights reserved.
 *    Copyright (C) 2011-2014 Greg Horn
 *
 *    CasADi is free software; you can redistribute it and/or
 *    modify it under the terms of the GNU Lesser General Public
 *    License as published by the Free Software Foundation; either
 *    version 3 of the License, or (at your option) any later version.
 *
 *    CasADi 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
 *    Lesser General Public License for more details.
 *
 *    You should have received a copy of the GNU Lesser General Public
 *    License along with CasADi; if not, write to the Free Software
 *    Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
 *
 */
 
 
#ifndef CASADI_SPARSITY_INTERFACE_HPP
#define CASADI_SPARSITY_INTERFACE_HPP
 
#include "casadi_misc.hpp"
 
namespace casadi {
  struct CASADI_EXPORT SparsityInterfaceCommon {};
 
#ifndef SWIG
  template<typename MatType>
  class SparsityInterface : public SparsityInterfaceCommon {
  protected:
    // Helper functions
    inline const MatType& self() const { return static_cast<const MatType&>(*this); }
    inline MatType& self() { return static_cast<MatType&>(*this); }
  public:
 
    static std::vector< std::vector< MatType > >
      blocksplit(const MatType& x, const std::vector<casadi_int>& vert_offset,
                 const std::vector<casadi_int>& horz_offset);
    static std::vector< std::vector< MatType > >
      blocksplit(const MatType& x, casadi_int vert_incr, casadi_int horz_incr);
    static MatType veccat(const std::vector< MatType >& x);
    static MatType vec(const MatType& x);
    static MatType repmat(const MatType& x, casadi_int n, casadi_int m=1);
    static std::vector<casadi_int> offset(const std::vector< MatType > &v, bool vert=true);
    static std::vector< MatType > diagsplit(const MatType& x,
                                            const std::vector<casadi_int>& output_offset);
    static std::vector< MatType > diagsplit(const MatType& x, casadi_int incr);
    static std::vector< MatType > diagsplit(const MatType& x, casadi_int incr1, casadi_int incr2);
    static MatType mtimes(const std::vector<MatType> &args,
                          const std::string& blas = "reference");
    static std::vector<MatType > horzsplit(const MatType& x, casadi_int incr);
    static std::vector<MatType > vertsplit(const MatType& x, casadi_int incr);
    static std::vector<MatType > horzsplit_n(const MatType& x, casadi_int n);
    static std::vector<MatType > vertsplit_n(const MatType& x, casadi_int n);
    static MatType repmat(const MatType &A, const std::pair<casadi_int, casadi_int>& rc) {
      return MatType::repmat(A, rc.first, rc.second);
    }
 
    inline friend MatType horzcat(const std::vector<MatType> &v) {
      return MatType::horzcat(v);
    }
 
    inline friend MatType vertcat(const std::vector<MatType> &v) {
      return MatType::vertcat(v);
    }
 
    inline friend std::vector<MatType >
      horzsplit(const MatType &x, const std::vector<casadi_int>& offset) {
      return MatType::horzsplit(x, offset);
    }
 
    inline friend std::vector<MatType > horzsplit(const MatType& x, casadi_int incr=1) {
      return MatType::horzsplit(x, incr);
    }
 
    inline friend std::vector<MatType > horzsplit_n(const MatType& x, casadi_int n) {
      return MatType::horzsplit_n(x, n);
    }
 
    inline friend std::vector<MatType >
      vertsplit(const MatType& x, const std::vector<casadi_int>& offset) {
      return MatType::vertsplit(x, offset);
    }
 
    inline friend std::vector<casadi_int > offset(const std::vector<MatType> &v, bool vert=true) {
      return MatType::offset(v, vert);
    }
 
    inline friend std::vector<MatType > vertsplit(const MatType &x, casadi_int incr=1) {
      return MatType::vertsplit(x, incr);
    }
 
    inline friend std::vector<MatType > vertsplit_n(const MatType& x, casadi_int n) {
      return MatType::vertsplit_n(x, n);
    }
 
    inline friend MatType blockcat(const std::vector< std::vector<MatType > > &v) {
      return MatType::blockcat(v);
    }
 
    inline friend MatType
      blockcat(const MatType &A, const MatType &B, const MatType &C, const MatType &D) {
      return vertcat(horzcat(A, B), horzcat(C, D));
    }
 
    inline friend std::vector< std::vector< MatType > >
      blocksplit(const MatType& x,
                 const std::vector<casadi_int>& vert_offset,
                 const std::vector<casadi_int>& horz_offset) {
      return MatType::blocksplit(x, vert_offset, horz_offset);
    }
 
    inline friend std::vector< std::vector< MatType > >
      blocksplit(const MatType& x, casadi_int vert_incr=1, casadi_int horz_incr=1) {
      return MatType::blocksplit(x, vert_incr, horz_incr);
    }
 
    inline friend MatType diagcat(const std::vector<MatType> &A) {
      return MatType::diagcat(A);
    }
 
    friend std::vector< MatType >
      diagsplit(const MatType& x,
                const std::vector<casadi_int>& output_offset1,
                const std::vector<casadi_int>& output_offset2) {
      return MatType::diagsplit(x, output_offset1, output_offset2);
    }
 
    inline friend std::vector< MatType >
      diagsplit(const MatType& x, const std::vector<casadi_int>& output_offset) {
      return MatType::diagsplit(x, output_offset);
    }
 
    inline friend std::vector< MatType >
      diagsplit(const MatType& x, casadi_int incr=1) {
      return MatType::diagsplit(x, incr);
    }
 
    inline friend std::vector< MatType >
      diagsplit(const MatType& x, casadi_int incr1, casadi_int incr2) {
      return MatType::diagsplit(x, incr1, incr2);
    }
 
    inline friend MatType veccat(const std::vector< MatType >& x) {
      return MatType::veccat(x);
    }
 
    inline friend MatType mtimes(const MatType &x, const MatType &y,
                                 const std::string& blas = "reference") {
      return MatType::mtimes(x, y, blas);
    }
 
    inline friend MatType mtimes(const std::vector<MatType> &args,
                                 const std::string& blas = "reference") {
      return MatType::mtimes(args, blas);
    }
 
    inline friend MatType
      mac(const MatType &x, const MatType &y, const MatType &z,
          const std::string& blas = "reference") {
      return MatType::mac(x, y, z, blas);
    }
 
    inline friend MatType transpose(const MatType& X) {
      return X.T();
    }
 
    inline friend MatType vec(const MatType& x) {
      return MatType::vec(x);
    }
 
    inline friend MatType reshape(const MatType& x, casadi_int nrow, casadi_int ncol) {
      return MatType::reshape(x, nrow, ncol);
    }
 
    inline friend MatType reshape(const MatType& x, std::pair<casadi_int, casadi_int> rc) {
      return MatType::reshape(x, rc.first, rc.second);
    }
 
    inline friend MatType reshape(const MatType& x, const Sparsity& sp) {
      return MatType::reshape(x, sp);
    }
 
    inline friend MatType sparsity_cast(const MatType& x, const Sparsity& sp) {
      return MatType::sparsity_cast(x, sp);
    }
 
    inline friend casadi_int sprank(const MatType& x) {
      return MatType::sprank(x);
    }
 
    inline friend casadi_int norm_0_mul(const MatType &x, const MatType &y) {
      return MatType::norm_0_mul(x, y);
    }
 
    inline friend MatType triu(const MatType& x, bool includeDiagonal=true) {
      return MatType::triu(x, includeDiagonal);
    }
 
    inline friend MatType tril(const MatType& x, bool includeDiagonal=true) {
      return MatType::tril(x, includeDiagonal);
    }
 
    inline friend MatType kron(const MatType& a, const MatType& b) {
      return MatType::kron(a, b);
    }
 
    inline friend MatType repmat(const MatType &A, casadi_int n, casadi_int m=1) {
      return MatType::repmat(A, n, m);
    }
 
    inline friend MatType repmat(const MatType &A, const std::pair<casadi_int, casadi_int>& rc) {
      return MatType::repmat(A, rc);
    }
 
    inline friend MatType horzcat(const MatType &x, const MatType &y) {
      return horzcat(std::vector<MatType>{x, y});
    }
 
    inline friend MatType horzcat(const MatType &x, const MatType &y, const MatType &z) {
      return horzcat(std::vector<MatType>{x, y, z});
    }
 
    inline friend MatType horzcat(const MatType &x, const MatType &y, const MatType &z,
                                  const MatType &w) {
      return horzcat(std::vector<MatType>{x, y, z, w});
    }
 
    inline friend MatType horzcat(const MatType &x, const MatType &y, const MatType &z,
                                  const MatType &w, const MatType &v) {
      return horzcat(std::vector<MatType>{x, y, z, w, v});
    }
 
    inline friend MatType horzcat(const MatType &x, const MatType &y, const MatType &z,
                                  const MatType &w, const MatType &v, const MatType &u) {
      return horzcat(std::vector<MatType>{x, y, z, w, v, u});
    }
 
    inline friend MatType vertcat(const MatType &x, const MatType &y) {
      return vertcat(std::vector<MatType>{x, y});
    }
 
    inline friend MatType vertcat(const MatType &x, const MatType &y, const MatType &z) {
      return vertcat(std::vector<MatType>{x, y, z});
    }
 
    inline friend MatType vertcat(const MatType &x, const MatType &y, const MatType &z,
                                  const MatType &w) {
      return vertcat(std::vector<MatType>{x, y, z, w});
    }
 
    inline friend MatType vertcat(const MatType &x, const MatType &y, const MatType &z,
                                  const MatType &w, const MatType &v) {
      return vertcat(std::vector<MatType>{x, y, z, w, v});
    }
 
    inline friend MatType vertcat(const MatType &x, const MatType &y, const MatType &z,
                                  const MatType &w, const MatType &v, const MatType &u) {
      return vertcat(std::vector<MatType>{x, y, z, w, v, u});
    }
 
    inline friend MatType diagcat(const MatType &x, const MatType &y) {
      return diagcat(std::vector<MatType>{x, y});
    }
 
    inline friend MatType diagcat(const MatType &x, const MatType &y, const MatType &z) {
      return diagcat(std::vector<MatType>{x, y, z});
    }
 
    inline friend MatType diagcat(const MatType &x, const MatType &y, const MatType &z,
                                  const MatType &w) {
      return diagcat(std::vector<MatType>{x, y, z, w});
    }
 
    inline friend MatType diagcat(const MatType &x, const MatType &y, const MatType &z,
                                  const MatType &w, const MatType &v) {
      return diagcat(std::vector<MatType>{x, y, z, w, v});
    }
 
    inline friend MatType diagcat(const MatType &x, const MatType &y, const MatType &z,
                                  const MatType &w, const MatType &v, const MatType &u) {
      return diagcat(std::vector<MatType>{x, y, z, w, v, u});
    }
 
    inline friend MatType sum1(const MatType &x) { return MatType::sum1(x);}
 
    inline friend MatType sum2(const MatType &x) { return MatType::sum2(x);}
 
    inline friend MatType sum(const MatType &x) {
      if (x.is_vector()) {
        if (x.is_column()) {
          return MatType::sum1(x);
        } else {
          return MatType::sum2(x);
        }
      }
      if (x.size1()>x.size2()) {
        return MatType::sum2(MatType::sum1(x));
      } else {
        return MatType::sum1(MatType::sum2(x));
      }
    }
 
  };
#endif // SWIG
 
#ifndef SWIG
  template<typename MatType>
  MatType SparsityInterface<MatType>::vec(const MatType& x) {
    if (x.size2()==1) return x;
    return reshape(x, x.numel(), 1);
  }
 
  template<typename MatType>
  MatType SparsityInterface<MatType>::repmat(const MatType& x, casadi_int n, casadi_int m) {
    casadi_assert(n>=0 && m>=0,
      "repmat(x,n,m): n (" + str(n) + ") and m (" + str(m) + ") must be non-negative");
    if (n==1 && m==1) return x;
    MatType allrows = vertcat(std::vector<MatType>(n, x));
    if (n==0) allrows = MatType(0, x.size2());
    MatType ret = horzcat(std::vector<MatType>(m, allrows));
    if (m==0) ret = MatType(allrows.size1(), 0);
    return ret;
  }
 
  template<typename MatType>
  std::vector< std::vector< MatType > >
  SparsityInterface<MatType>::blocksplit(const MatType& x,
                                         const std::vector<casadi_int>& vert_offset,
                                         const std::vector<casadi_int>& horz_offset) {
    std::vector<MatType> rows = MatType::vertsplit(x, vert_offset);
    std::vector< std::vector< MatType > > ret;
    for (auto&& r : rows) ret.push_back(MatType::horzsplit(r, horz_offset));
    return ret;
  }
 
  template<typename MatType>
  std::vector< std::vector< MatType > >
  SparsityInterface<MatType>::blocksplit(const MatType& x,
      casadi_int vert_incr, casadi_int horz_incr) {
    casadi_assert(vert_incr>=1, "blocksplit: vert_incr (" + str(vert_incr) + ") must be positive");
    casadi_assert(horz_incr>=1, "blocksplit: horz_incr (" + str(horz_incr) + ") must be positive");
    casadi_int sz1 = x.size1();
    std::vector<casadi_int> offset1 = range(0, sz1, vert_incr);
    offset1.push_back(sz1);
    casadi_int sz2 = x.size2();
    std::vector<casadi_int> offset2 = range(0, sz2, horz_incr);
    offset2.push_back(sz2);
    return blocksplit(x, offset1, offset2);
  }
 
  template<typename MatType>
  std::vector<casadi_int>
  SparsityInterface<MatType>::offset(const std::vector< MatType > &v, bool vert) {
    std::vector<casadi_int> ret(v.size()+1);
    ret[0]=0;
    for (casadi_int i=0; i<v.size(); ++i) {
      ret[i+1] = ret[i] + (vert ? v[i].size1() : v[i].size2());
    }
    return ret;
  }
 
  template<typename MatType>
  MatType SparsityInterface<MatType>::veccat(const std::vector< MatType >& x) {
    std::vector< MatType > x_vec = x;
    for (typename std::vector< MatType >::iterator it=x_vec.begin();
         it!=x_vec.end(); ++it) {
      *it = vec(*it);
    }
    if (x_vec.empty()) {
      return MatType(0, 1);
    } else {
      return vertcat(x_vec);
    }
  }
 
  template<typename MatType>
  std::vector< MatType >
  SparsityInterface<MatType>::diagsplit(const MatType& x,
      const std::vector<casadi_int>& output_offset) {
    casadi_assert(x.is_square(), "diagsplit(x,incr)::input must be square but got "
                          + x.dim()  + ".");
    return MatType::diagsplit(x, output_offset, output_offset);
  }
 
  template<typename MatType>
  std::vector< MatType >
  SparsityInterface<MatType>::diagsplit(const MatType& x, casadi_int incr) {
    casadi_assert(incr>=1, "diagsplit(x,incr): incr (" + str(incr) + ") must be positive");
    casadi_assert(x.is_square(), "diagsplit(x,incr)::input must be square but got "
                          + x.dim()  + ".");
    std::vector<casadi_int> offset2 = range(0, x.size2(), incr);
    offset2.push_back(x.size2());
    return MatType::diagsplit(x, offset2);
  }
 
  template<typename MatType>
  std::vector< MatType >
  SparsityInterface<MatType>::diagsplit(const MatType& x, casadi_int incr1, casadi_int incr2) {
    casadi_assert(incr1>=1, "diagsplit: incr1 (" + str(incr1) + ") must be positive");
    casadi_assert(incr2>=1, "diagsplit: incr2 (" + str(incr2) + ") must be positive");
    std::vector<casadi_int> offset1 = range(0, x.size1(), incr1);
    offset1.push_back(x.size1());
    std::vector<casadi_int> offset2 = range(0, x.size2(), incr2);
    offset2.push_back(x.size2());
    return MatType::diagsplit(x, offset1, offset2);
  }
 
  template<typename MatType>
  MatType SparsityInterface<MatType>::mtimes(const std::vector<MatType> &args,
                                             const std::string& blas) {
    casadi_assert(!args.empty(),
                          "mul(std::vector<MatType> &args): "
                          "supplied list must not be empty.");
    MatType ret = args[0];
    for (casadi_int i=1; i<args.size(); ++i) ret = MatType::mtimes(ret, args[i], blas);
    return ret;
  }
 
  template<typename MatType>
  std::vector<MatType > SparsityInterface<MatType>::horzsplit(const MatType& x, casadi_int incr) {
    casadi_assert(incr>=1, "horzsplit(x,incr): incr (" + str(incr) + ") must be positive");
    casadi_int sz2 = x.size2();
    std::vector<casadi_int> offset2 = range(0, sz2, incr);
    offset2.push_back(sz2);
    return MatType::horzsplit(x, offset2);
  }
 
  template<typename MatType>
  std::vector<MatType > SparsityInterface<MatType>::vertsplit(const MatType& x, casadi_int incr) {
    casadi_assert(incr>=1, "vertsplit(x,incr): incr (" + str(incr) + ") must be positive");
    casadi_int sz1 = x.size1();
    std::vector<casadi_int> offset1 = range(0, sz1, incr);
    offset1.push_back(sz1);
    return MatType::vertsplit(x, offset1);
  }
  template<typename MatType>
  std::vector<MatType > SparsityInterface<MatType>::horzsplit_n(const MatType& x, casadi_int n) {
    casadi_assert(n>=0, "horzsplit_n(x,n): n (" + str(n) + ") must be non-negative");
    if (x.size2()==0) return std::vector<MatType>(n, x);
    casadi_assert(n>0 && x.size2() % n==0, "horzsplit_n(x,n): x.size2() (" + str(x.size2()) +
                     ") must be a multiple of n (" + str(n) + ")");
    return horzsplit(x, x.size2()/n);
  }
  template<typename MatType>
  std::vector<MatType > SparsityInterface<MatType>::vertsplit_n(const MatType& x, casadi_int n) {
    casadi_assert(n>=0, "vertsplit_n(x,n): n (" + str(n) + ") must be non-negative");
    if (x.size1()==0) return std::vector<MatType>(n, x);
    casadi_assert(n>0 && x.size1() % n==0, "vertsplit_n(x,n): x.size1() (" + str(x.size1()) +
                     ") must be a multiple of n (" + str(n) + ")");
    return vertsplit(x, x.size1()/n);
  }
 
#endif // SWIG
 
} // namespace casadi
 
#endif // CASADI_SPARSITY_INTERFACE_HPP