sparsity_interface.hpp
1 /*
2  * This file is part of CasADi.
3  *
4  * CasADi -- A symbolic framework for dynamic optimization.
5  * Copyright (C) 2010-2023 Joel Andersson, Joris Gillis, Moritz Diehl,
6  * KU Leuven. All rights reserved.
7  * Copyright (C) 2011-2014 Greg Horn
8  *
9  * CasADi is free software; you can redistribute it and/or
10  * modify it under the terms of the GNU Lesser General Public
11  * License as published by the Free Software Foundation; either
12  * version 3 of the License, or (at your option) any later version.
13  *
14  * CasADi is distributed in the hope that it will be useful,
15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17  * Lesser General Public License for more details.
18  *
19  * You should have received a copy of the GNU Lesser General Public
20  * License along with CasADi; if not, write to the Free Software
21  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
22  *
23  */
24 
25 
26 #ifndef CASADI_SPARSITY_INTERFACE_HPP
27 #define CASADI_SPARSITY_INTERFACE_HPP
28 
29 #include "casadi_misc.hpp"
30 
31 namespace casadi {
37  struct CASADI_EXPORT SparsityInterfaceCommon {};
38 
39 #ifndef SWIG
50  template<typename MatType>
51  class SparsityInterface : public SparsityInterfaceCommon {
52  protected:
53  // Helper functions
54  inline const MatType& self() const { return static_cast<const MatType&>(*this); }
55  inline MatType& self() { return static_cast<MatType&>(*this); }
56  public:
57 
59  static std::vector< std::vector< MatType > >
60  blocksplit(const MatType& x, const std::vector<casadi_int>& vert_offset,
61  const std::vector<casadi_int>& horz_offset);
62  static std::vector< std::vector< MatType > >
63  blocksplit(const MatType& x, casadi_int vert_incr, casadi_int horz_incr);
64  static MatType veccat(const std::vector< MatType >& x);
65  static MatType vec(const MatType& x);
66  static MatType repmat(const MatType& x, casadi_int n, casadi_int m=1);
67  static std::vector<casadi_int> offset(const std::vector< MatType > &v, bool vert=true);
68  static std::vector< MatType > diagsplit(const MatType& x,
69  const std::vector<casadi_int>& output_offset);
70  static std::vector< MatType > diagsplit(const MatType& x, casadi_int incr);
71  static std::vector< MatType > diagsplit(const MatType& x, casadi_int incr1, casadi_int incr2);
72  static MatType mtimes(const std::vector<MatType> &args,
73  const std::string& blas = "reference");
74  static std::vector<MatType > horzsplit(const MatType& x, casadi_int incr);
75  static std::vector<MatType > vertsplit(const MatType& x, casadi_int incr);
76  static std::vector<MatType > horzsplit_n(const MatType& x, casadi_int n);
77  static std::vector<MatType > vertsplit_n(const MatType& x, casadi_int n);
78  static MatType repmat(const MatType &A, const std::pair<casadi_int, casadi_int>& rc) {
79  return MatType::repmat(A, rc.first, rc.second);
80  }
82 
98  inline friend MatType horzcat(const std::vector<MatType> &v) {
99  return MatType::horzcat(v);
100  }
101 
109  inline friend MatType vertcat(const std::vector<MatType> &v) {
110  return MatType::vertcat(v);
111  }
112 
121  inline friend std::vector<MatType >
122  horzsplit(const MatType &x, const std::vector<casadi_int>& offset) {
123  return MatType::horzsplit(x, offset);
124  }
125 
135  inline friend std::vector<MatType > horzsplit(const MatType& x, casadi_int incr=1) {
136  return MatType::horzsplit(x, incr);
137  }
138 
150  inline friend std::vector<MatType > horzsplit_n(const MatType& x, casadi_int n) {
151  return MatType::horzsplit_n(x, n);
152  }
153 
162  inline friend std::vector<MatType >
163  vertsplit(const MatType& x, const std::vector<casadi_int>& offset) {
164  return MatType::vertsplit(x, offset);
165  }
166 
170  inline friend std::vector<casadi_int > offset(const std::vector<MatType> &v, bool vert=true) {
171  return MatType::offset(v, vert);
172  }
173 
205  inline friend std::vector<MatType > vertsplit(const MatType &x, casadi_int incr=1) {
206  return MatType::vertsplit(x, incr);
207  }
208 
220  inline friend std::vector<MatType > vertsplit_n(const MatType& x, casadi_int n) {
221  return MatType::vertsplit_n(x, n);
222  }
223 
227  inline friend MatType blockcat(const std::vector< std::vector<MatType > > &v) {
228  return MatType::blockcat(v);
229  }
230 
234  inline friend MatType
235  blockcat(const MatType &A, const MatType &B, const MatType &C, const MatType &D) {
236  return vertcat(horzcat(A, B), horzcat(C, D));
237  }
238 
247  inline friend std::vector< std::vector< MatType > >
248  blocksplit(const MatType& x,
249  const std::vector<casadi_int>& vert_offset,
250  const std::vector<casadi_int>& horz_offset) {
251  return MatType::blocksplit(x, vert_offset, horz_offset);
252  }
253 
262  inline friend std::vector< std::vector< MatType > >
263  blocksplit(const MatType& x, casadi_int vert_incr=1, casadi_int horz_incr=1) {
264  return MatType::blocksplit(x, vert_incr, horz_incr);
265  }
266 
270  inline friend MatType diagcat(const std::vector<MatType> &A) {
271  return MatType::diagcat(A);
272  }
273 
284  friend std::vector< MatType >
285  diagsplit(const MatType& x,
286  const std::vector<casadi_int>& output_offset1,
287  const std::vector<casadi_int>& output_offset2) {
288  return MatType::diagsplit(x, output_offset1, output_offset2);
289  }
290 
299  inline friend std::vector< MatType >
300  diagsplit(const MatType& x, const std::vector<casadi_int>& output_offset) {
301  return MatType::diagsplit(x, output_offset);
302  }
303 
311  inline friend std::vector< MatType >
312  diagsplit(const MatType& x, casadi_int incr=1) {
313  return MatType::diagsplit(x, incr);
314  }
315 
324  inline friend std::vector< MatType >
325  diagsplit(const MatType& x, casadi_int incr1, casadi_int incr2) {
326  return MatType::diagsplit(x, incr1, incr2);
327  }
328 
332  inline friend MatType veccat(const std::vector< MatType >& x) {
333  return MatType::veccat(x);
334  }
335 
356  inline friend MatType mtimes(const MatType &x, const MatType &y,
357  const std::string& blas = "reference") {
358  return MatType::mtimes(x, y, blas);
359  }
360 
366  inline friend MatType mtimes(const std::vector<MatType> &args,
367  const std::string& blas = "reference") {
368  return MatType::mtimes(args, blas);
369  }
370 
380  inline friend MatType
381  mac(const MatType &x, const MatType &y, const MatType &z,
382  const std::string& blas = "reference") {
383  return MatType::mac(x, y, z, blas);
384  }
385 
389  inline friend MatType transpose(const MatType& X) {
390  return X.T();
391  }
392 
410  inline friend MatType vec(const MatType& x) {
411  return MatType::vec(x);
412  }
413 
417  inline friend MatType reshape(const MatType& x, casadi_int nrow, casadi_int ncol) {
418  return MatType::reshape(x, nrow, ncol);
419  }
420 
424  inline friend MatType reshape(const MatType& x, std::pair<casadi_int, casadi_int> rc) {
425  return MatType::reshape(x, rc.first, rc.second);
426  }
427 
431  inline friend MatType reshape(const MatType& x, const Sparsity& sp) {
432  return MatType::reshape(x, sp);
433  }
434 
438  inline friend MatType sparsity_cast(const MatType& x, const Sparsity& sp) {
439  return MatType::sparsity_cast(x, sp);
440  }
441 
445  inline friend casadi_int sprank(const MatType& x) {
446  return MatType::sprank(x);
447  }
448 
452  inline friend casadi_int norm_0_mul(const MatType &x, const MatType &y) {
453  return MatType::norm_0_mul(x, y);
454  }
455 
459  inline friend MatType triu(const MatType& x, bool includeDiagonal=true) {
460  return MatType::triu(x, includeDiagonal);
461  }
462 
466  inline friend MatType tril(const MatType& x, bool includeDiagonal=true) {
467  return MatType::tril(x, includeDiagonal);
468  }
469 
475  inline friend MatType kron(const MatType& a, const MatType& b) {
476  return MatType::kron(a, b);
477  }
478 
482  inline friend MatType repmat(const MatType &A, casadi_int n, casadi_int m=1) {
483  return MatType::repmat(A, n, m);
484  }
485 
489  inline friend MatType repmat(const MatType &A, const std::pair<casadi_int, casadi_int>& rc) {
490  return MatType::repmat(A, rc);
491  }
492 
496  inline friend MatType horzcat(const MatType &x, const MatType &y) {
497  return horzcat(std::vector<MatType>{x, y});
498  }
499 
503  inline friend MatType horzcat(const MatType &x, const MatType &y, const MatType &z) {
504  return horzcat(std::vector<MatType>{x, y, z});
505  }
506 
510  inline friend MatType horzcat(const MatType &x, const MatType &y, const MatType &z,
511  const MatType &w) {
512  return horzcat(std::vector<MatType>{x, y, z, w});
513  }
514 
518  inline friend MatType horzcat(const MatType &x, const MatType &y, const MatType &z,
519  const MatType &w, const MatType &v) {
520  return horzcat(std::vector<MatType>{x, y, z, w, v});
521  }
522 
526  inline friend MatType horzcat(const MatType &x, const MatType &y, const MatType &z,
527  const MatType &w, const MatType &v, const MatType &u) {
528  return horzcat(std::vector<MatType>{x, y, z, w, v, u});
529  }
530 
534  inline friend MatType vertcat(const MatType &x, const MatType &y) {
535  return vertcat(std::vector<MatType>{x, y});
536  }
537 
541  inline friend MatType vertcat(const MatType &x, const MatType &y, const MatType &z) {
542  return vertcat(std::vector<MatType>{x, y, z});
543  }
544 
548  inline friend MatType vertcat(const MatType &x, const MatType &y, const MatType &z,
549  const MatType &w) {
550  return vertcat(std::vector<MatType>{x, y, z, w});
551  }
552 
556  inline friend MatType vertcat(const MatType &x, const MatType &y, const MatType &z,
557  const MatType &w, const MatType &v) {
558  return vertcat(std::vector<MatType>{x, y, z, w, v});
559  }
560 
564  inline friend MatType vertcat(const MatType &x, const MatType &y, const MatType &z,
565  const MatType &w, const MatType &v, const MatType &u) {
566  return vertcat(std::vector<MatType>{x, y, z, w, v, u});
567  }
568 
572  inline friend MatType diagcat(const MatType &x, const MatType &y) {
573  return diagcat(std::vector<MatType>{x, y});
574  }
575 
579  inline friend MatType diagcat(const MatType &x, const MatType &y, const MatType &z) {
580  return diagcat(std::vector<MatType>{x, y, z});
581  }
582 
586  inline friend MatType diagcat(const MatType &x, const MatType &y, const MatType &z,
587  const MatType &w) {
588  return diagcat(std::vector<MatType>{x, y, z, w});
589  }
590 
594  inline friend MatType diagcat(const MatType &x, const MatType &y, const MatType &z,
595  const MatType &w, const MatType &v) {
596  return diagcat(std::vector<MatType>{x, y, z, w, v});
597  }
598 
602  inline friend MatType diagcat(const MatType &x, const MatType &y, const MatType &z,
603  const MatType &w, const MatType &v, const MatType &u) {
604  return diagcat(std::vector<MatType>{x, y, z, w, v, u});
605  }
606 
610  inline friend MatType sum1(const MatType &x) { return MatType::sum1(x);}
611 
615  inline friend MatType sum2(const MatType &x) { return MatType::sum2(x);}
616 
620  inline friend MatType sum(const MatType &x) {
621  if (x.is_vector()) {
622  if (x.is_column()) {
623  return MatType::sum1(x);
624  } else {
625  return MatType::sum2(x);
626  }
627  }
628  if (x.size1()>x.size2()) {
629  return MatType::sum2(MatType::sum1(x));
630  } else {
631  return MatType::sum1(MatType::sum2(x));
632  }
633  }
634 
636  };
637 #endif // SWIG
638 
639 #ifndef SWIG
640  template<typename MatType>
641  MatType SparsityInterface<MatType>::vec(const MatType& x) {
642  if (x.size2()==1) return x;
643  return reshape(x, x.numel(), 1);
644  }
645 
646  template<typename MatType>
647  MatType SparsityInterface<MatType>::repmat(const MatType& x, casadi_int n, casadi_int m) {
648  if (n==1 && m==1) return x;
649  MatType allrows = vertcat(std::vector<MatType>(n, x));
650  if (n==0) allrows = MatType(0, x.size2());
651  MatType ret = horzcat(std::vector<MatType>(m, allrows));
652  if (m==0) ret = MatType(allrows.size1(), 0);
653  return ret;
654  }
655 
656  template<typename MatType>
657  std::vector< std::vector< MatType > >
658  SparsityInterface<MatType>::blocksplit(const MatType& x,
659  const std::vector<casadi_int>& vert_offset,
660  const std::vector<casadi_int>& horz_offset) {
661  std::vector<MatType> rows = MatType::vertsplit(x, vert_offset);
662  std::vector< std::vector< MatType > > ret;
663  for (auto&& r : rows) ret.push_back(MatType::horzsplit(r, horz_offset));
664  return ret;
665  }
666 
667  template<typename MatType>
668  std::vector< std::vector< MatType > >
669  SparsityInterface<MatType>::blocksplit(const MatType& x,
670  casadi_int vert_incr, casadi_int horz_incr) {
671  casadi_assert_dev(horz_incr>=1);
672  casadi_assert_dev(vert_incr>=1);
673  casadi_int sz1 = x.size1();
674  std::vector<casadi_int> offset1 = range(0, sz1, vert_incr);
675  offset1.push_back(sz1);
676  casadi_int sz2 = x.size2();
677  std::vector<casadi_int> offset2 = range(0, sz2, horz_incr);
678  offset2.push_back(sz2);
679  return blocksplit(x, offset1, offset2);
680  }
681 
682  template<typename MatType>
683  std::vector<casadi_int>
684  SparsityInterface<MatType>::offset(const std::vector< MatType > &v, bool vert) {
685  std::vector<casadi_int> ret(v.size()+1);
686  ret[0]=0;
687  for (casadi_int i=0; i<v.size(); ++i) {
688  ret[i+1] = ret[i] + (vert ? v[i].size1() : v[i].size2());
689  }
690  return ret;
691  }
692 
693  template<typename MatType>
694  MatType SparsityInterface<MatType>::veccat(const std::vector< MatType >& x) {
695  std::vector< MatType > x_vec = x;
696  for (typename std::vector< MatType >::iterator it=x_vec.begin();
697  it!=x_vec.end(); ++it) {
698  *it = vec(*it);
699  }
700  if (x_vec.empty()) {
701  return MatType(0, 1);
702  } else {
703  return vertcat(x_vec);
704  }
705  }
706 
707  template<typename MatType>
708  std::vector< MatType >
709  SparsityInterface<MatType>::diagsplit(const MatType& x,
710  const std::vector<casadi_int>& output_offset) {
711  casadi_assert(x.is_square(), "diagsplit(x,incr)::input must be square but got "
712  + x.dim() + ".");
713  return MatType::diagsplit(x, output_offset, output_offset);
714  }
715 
716  template<typename MatType>
717  std::vector< MatType >
718  SparsityInterface<MatType>::diagsplit(const MatType& x, casadi_int incr) {
719  casadi_assert_dev(incr>=1);
720  casadi_assert(x.is_square(), "diagsplit(x,incr)::input must be square but got "
721  + x.dim() + ".");
722  std::vector<casadi_int> offset2 = range(0, x.size2(), incr);
723  offset2.push_back(x.size2());
724  return MatType::diagsplit(x, offset2);
725  }
726 
727  template<typename MatType>
728  std::vector< MatType >
729  SparsityInterface<MatType>::diagsplit(const MatType& x, casadi_int incr1, casadi_int incr2) {
730  casadi_assert_dev(incr1>=1);
731  casadi_assert_dev(incr2>=1);
732  std::vector<casadi_int> offset1 = range(0, x.size1(), incr1);
733  offset1.push_back(x.size1());
734  std::vector<casadi_int> offset2 = range(0, x.size2(), incr2);
735  offset2.push_back(x.size2());
736  return MatType::diagsplit(x, offset1, offset2);
737  }
738 
739  template<typename MatType>
740  MatType SparsityInterface<MatType>::mtimes(const std::vector<MatType> &args,
741  const std::string& blas) {
742  casadi_assert(!args.empty(),
743  "mul(std::vector<MatType> &args): "
744  "supplied list must not be empty.");
745  MatType ret = args[0];
746  for (casadi_int i=1; i<args.size(); ++i) ret = MatType::mtimes(ret, args[i], blas);
747  return ret;
748  }
749 
750  template<typename MatType>
751  std::vector<MatType > SparsityInterface<MatType>::horzsplit(const MatType& x, casadi_int incr) {
752  casadi_assert_dev(incr>=1);
753  casadi_int sz2 = x.size2();
754  std::vector<casadi_int> offset2 = range(0, sz2, incr);
755  offset2.push_back(sz2);
756  return MatType::horzsplit(x, offset2);
757  }
758 
759  template<typename MatType>
760  std::vector<MatType > SparsityInterface<MatType>::vertsplit(const MatType& x, casadi_int incr) {
761  casadi_assert_dev(incr>=1);
762  casadi_int sz1 = x.size1();
763  std::vector<casadi_int> offset1 = range(0, sz1, incr);
764  offset1.push_back(sz1);
765  return MatType::vertsplit(x, offset1);
766  }
767  template<typename MatType>
768  std::vector<MatType > SparsityInterface<MatType>::horzsplit_n(const MatType& x, casadi_int n) {
769  casadi_assert(n>=0, "horzsplit_n(x,n): n (" + str(n) + ") must be non-negative");
770  if (x.size2()==0) return std::vector<MatType>(n, x);
771  casadi_assert(x.size2() % n==0, "horzsplit_n(x,n): x.size2() (" + str(x.size2()) +
772  ") must be a multiple of n (" + str(n) + ")");
773  return horzsplit(x, x.size2()/n);
774  }
775  template<typename MatType>
776  std::vector<MatType > SparsityInterface<MatType>::vertsplit_n(const MatType& x, casadi_int n) {
777  casadi_assert(n>=0, "vertsplit_n(x,n): n (" + str(n) + ") must be non-negative");
778  if (x.size1()==0) return std::vector<MatType>(n, x);
779  casadi_assert(x.size1() % n==0, "vertsplit(x,n): x.size1() (" + str(x.size1()) +
780  ") must be a multiple of n (" + str(n) + ")");
781  return vertsplit(x, x.size1()/n);
782  }
783 
784 #endif // SWIG
785 
786 } // namespace casadi
787 
788 #endif // CASADI_SPARSITY_INTERFACE_HPP
casadi
The casadi namespace.
Definition: archiver.hpp:32