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  static std::vector<MatType > horzsplit(const MatType& x, casadi_int incr);
74  static std::vector<MatType > vertsplit(const MatType& x, casadi_int incr);
75  static std::vector<MatType > horzsplit_n(const MatType& x, casadi_int n);
76  static std::vector<MatType > vertsplit_n(const MatType& x, casadi_int n);
77  static MatType repmat(const MatType &A, const std::pair<casadi_int, casadi_int>& rc) {
78  return MatType::repmat(A, rc.first, rc.second);
79  }
81 
97  inline friend MatType horzcat(const std::vector<MatType> &v) {
98  return MatType::horzcat(v);
99  }
100 
108  inline friend MatType vertcat(const std::vector<MatType> &v) {
109  return MatType::vertcat(v);
110  }
111 
120  inline friend std::vector<MatType >
121  horzsplit(const MatType &x, const std::vector<casadi_int>& offset) {
122  return MatType::horzsplit(x, offset);
123  }
124 
134  inline friend std::vector<MatType > horzsplit(const MatType& x, casadi_int incr=1) {
135  return MatType::horzsplit(x, incr);
136  }
137 
149  inline friend std::vector<MatType > horzsplit_n(const MatType& x, casadi_int n) {
150  return MatType::horzsplit_n(x, n);
151  }
152 
161  inline friend std::vector<MatType >
162  vertsplit(const MatType& x, const std::vector<casadi_int>& offset) {
163  return MatType::vertsplit(x, offset);
164  }
165 
169  inline friend std::vector<casadi_int > offset(const std::vector<MatType> &v, bool vert=true) {
170  return MatType::offset(v, vert);
171  }
172 
204  inline friend std::vector<MatType > vertsplit(const MatType &x, casadi_int incr=1) {
205  return MatType::vertsplit(x, incr);
206  }
207 
219  inline friend std::vector<MatType > vertsplit_n(const MatType& x, casadi_int n) {
220  return MatType::vertsplit_n(x, n);
221  }
222 
226  inline friend MatType blockcat(const std::vector< std::vector<MatType > > &v) {
227  return MatType::blockcat(v);
228  }
229 
233  inline friend MatType
234  blockcat(const MatType &A, const MatType &B, const MatType &C, const MatType &D) {
235  return vertcat(horzcat(A, B), horzcat(C, D));
236  }
237 
246  inline friend std::vector< std::vector< MatType > >
247  blocksplit(const MatType& x,
248  const std::vector<casadi_int>& vert_offset,
249  const std::vector<casadi_int>& horz_offset) {
250  return MatType::blocksplit(x, vert_offset, horz_offset);
251  }
252 
261  inline friend std::vector< std::vector< MatType > >
262  blocksplit(const MatType& x, casadi_int vert_incr=1, casadi_int horz_incr=1) {
263  return MatType::blocksplit(x, vert_incr, horz_incr);
264  }
265 
269  inline friend MatType diagcat(const std::vector<MatType> &A) {
270  return MatType::diagcat(A);
271  }
272 
283  friend std::vector< MatType >
284  diagsplit(const MatType& x,
285  const std::vector<casadi_int>& output_offset1,
286  const std::vector<casadi_int>& output_offset2) {
287  return MatType::diagsplit(x, output_offset1, output_offset2);
288  }
289 
298  inline friend std::vector< MatType >
299  diagsplit(const MatType& x, const std::vector<casadi_int>& output_offset) {
300  return MatType::diagsplit(x, output_offset);
301  }
302 
310  inline friend std::vector< MatType >
311  diagsplit(const MatType& x, casadi_int incr=1) {
312  return MatType::diagsplit(x, incr);
313  }
314 
323  inline friend std::vector< MatType >
324  diagsplit(const MatType& x, casadi_int incr1, casadi_int incr2) {
325  return MatType::diagsplit(x, incr1, incr2);
326  }
327 
331  inline friend MatType veccat(const std::vector< MatType >& x) {
332  return MatType::veccat(x);
333  }
334 
338  inline friend MatType mtimes(const MatType &x, const MatType &y) {
339  return MatType::mtimes(x, y);
340  }
341 
345  inline friend MatType mtimes(const std::vector<MatType> &args) {
346  return MatType::mtimes(args);
347  }
348 
357  inline friend MatType
358  mac(const MatType &x, const MatType &y, const MatType &z) {
359  return MatType::mac(x, y, z);
360  }
361 
365  inline friend MatType transpose(const MatType& X) {
366  return X.T();
367  }
368 
386  inline friend MatType vec(const MatType& x) {
387  return MatType::vec(x);
388  }
389 
393  inline friend MatType reshape(const MatType& x, casadi_int nrow, casadi_int ncol) {
394  return MatType::reshape(x, nrow, ncol);
395  }
396 
400  inline friend MatType reshape(const MatType& x, std::pair<casadi_int, casadi_int> rc) {
401  return MatType::reshape(x, rc.first, rc.second);
402  }
403 
407  inline friend MatType reshape(const MatType& x, const Sparsity& sp) {
408  return MatType::reshape(x, sp);
409  }
410 
414  inline friend MatType sparsity_cast(const MatType& x, const Sparsity& sp) {
415  return MatType::sparsity_cast(x, sp);
416  }
417 
421  inline friend casadi_int sprank(const MatType& x) {
422  return MatType::sprank(x);
423  }
424 
428  inline friend casadi_int norm_0_mul(const MatType &x, const MatType &y) {
429  return MatType::norm_0_mul(x, y);
430  }
431 
435  inline friend MatType triu(const MatType& x, bool includeDiagonal=true) {
436  return MatType::triu(x, includeDiagonal);
437  }
438 
442  inline friend MatType tril(const MatType& x, bool includeDiagonal=true) {
443  return MatType::tril(x, includeDiagonal);
444  }
445 
451  inline friend MatType kron(const MatType& a, const MatType& b) {
452  return MatType::kron(a, b);
453  }
454 
458  inline friend MatType repmat(const MatType &A, casadi_int n, casadi_int m=1) {
459  return MatType::repmat(A, n, m);
460  }
461 
465  inline friend MatType repmat(const MatType &A, const std::pair<casadi_int, casadi_int>& rc) {
466  return MatType::repmat(A, rc);
467  }
468 
472  inline friend MatType horzcat(const MatType &x, const MatType &y) {
473  return horzcat(std::vector<MatType>{x, y});
474  }
475 
479  inline friend MatType horzcat(const MatType &x, const MatType &y, const MatType &z) {
480  return horzcat(std::vector<MatType>{x, y, z});
481  }
482 
486  inline friend MatType horzcat(const MatType &x, const MatType &y, const MatType &z,
487  const MatType &w) {
488  return horzcat(std::vector<MatType>{x, y, z, w});
489  }
490 
494  inline friend MatType horzcat(const MatType &x, const MatType &y, const MatType &z,
495  const MatType &w, const MatType &v) {
496  return horzcat(std::vector<MatType>{x, y, z, w, v});
497  }
498 
502  inline friend MatType horzcat(const MatType &x, const MatType &y, const MatType &z,
503  const MatType &w, const MatType &v, const MatType &u) {
504  return horzcat(std::vector<MatType>{x, y, z, w, v, u});
505  }
506 
510  inline friend MatType vertcat(const MatType &x, const MatType &y) {
511  return vertcat(std::vector<MatType>{x, y});
512  }
513 
517  inline friend MatType vertcat(const MatType &x, const MatType &y, const MatType &z) {
518  return vertcat(std::vector<MatType>{x, y, z});
519  }
520 
524  inline friend MatType vertcat(const MatType &x, const MatType &y, const MatType &z,
525  const MatType &w) {
526  return vertcat(std::vector<MatType>{x, y, z, w});
527  }
528 
532  inline friend MatType vertcat(const MatType &x, const MatType &y, const MatType &z,
533  const MatType &w, const MatType &v) {
534  return vertcat(std::vector<MatType>{x, y, z, w, v});
535  }
536 
540  inline friend MatType vertcat(const MatType &x, const MatType &y, const MatType &z,
541  const MatType &w, const MatType &v, const MatType &u) {
542  return vertcat(std::vector<MatType>{x, y, z, w, v, u});
543  }
544 
548  inline friend MatType diagcat(const MatType &x, const MatType &y) {
549  return diagcat(std::vector<MatType>{x, y});
550  }
551 
555  inline friend MatType diagcat(const MatType &x, const MatType &y, const MatType &z) {
556  return diagcat(std::vector<MatType>{x, y, z});
557  }
558 
562  inline friend MatType diagcat(const MatType &x, const MatType &y, const MatType &z,
563  const MatType &w) {
564  return diagcat(std::vector<MatType>{x, y, z, w});
565  }
566 
570  inline friend MatType diagcat(const MatType &x, const MatType &y, const MatType &z,
571  const MatType &w, const MatType &v) {
572  return diagcat(std::vector<MatType>{x, y, z, w, v});
573  }
574 
578  inline friend MatType diagcat(const MatType &x, const MatType &y, const MatType &z,
579  const MatType &w, const MatType &v, const MatType &u) {
580  return diagcat(std::vector<MatType>{x, y, z, w, v, u});
581  }
582 
586  inline friend MatType sum1(const MatType &x) { return MatType::sum1(x);}
587 
591  inline friend MatType sum2(const MatType &x) { return MatType::sum2(x);}
592 
596  inline friend MatType sum(const MatType &x) {
597  if (x.is_vector()) {
598  if (x.is_column()) {
599  return MatType::sum1(x);
600  } else {
601  return MatType::sum2(x);
602  }
603  }
604  if (x.size1()>x.size2()) {
605  return MatType::sum2(MatType::sum1(x));
606  } else {
607  return MatType::sum1(MatType::sum2(x));
608  }
609  }
610 
612  };
613 #endif // SWIG
614 
615 #ifndef SWIG
616  template<typename MatType>
617  MatType SparsityInterface<MatType>::vec(const MatType& x) {
618  if (x.size2()==1) return x;
619  return reshape(x, x.numel(), 1);
620  }
621 
622  template<typename MatType>
623  MatType SparsityInterface<MatType>::repmat(const MatType& x, casadi_int n, casadi_int m) {
624  if (n==1 && m==1) return x;
625  MatType allrows = vertcat(std::vector<MatType>(n, x));
626  if (n==0) allrows = MatType(0, x.size2());
627  MatType ret = horzcat(std::vector<MatType>(m, allrows));
628  if (m==0) ret = MatType(allrows.size1(), 0);
629  return ret;
630  }
631 
632  template<typename MatType>
633  std::vector< std::vector< MatType > >
634  SparsityInterface<MatType>::blocksplit(const MatType& x,
635  const std::vector<casadi_int>& vert_offset,
636  const std::vector<casadi_int>& horz_offset) {
637  std::vector<MatType> rows = MatType::vertsplit(x, vert_offset);
638  std::vector< std::vector< MatType > > ret;
639  for (auto&& r : rows) ret.push_back(MatType::horzsplit(r, horz_offset));
640  return ret;
641  }
642 
643  template<typename MatType>
644  std::vector< std::vector< MatType > >
645  SparsityInterface<MatType>::blocksplit(const MatType& x,
646  casadi_int vert_incr, casadi_int horz_incr) {
647  casadi_assert_dev(horz_incr>=1);
648  casadi_assert_dev(vert_incr>=1);
649  casadi_int sz1 = x.size1();
650  std::vector<casadi_int> offset1 = range(0, sz1, vert_incr);
651  offset1.push_back(sz1);
652  casadi_int sz2 = x.size2();
653  std::vector<casadi_int> offset2 = range(0, sz2, horz_incr);
654  offset2.push_back(sz2);
655  return blocksplit(x, offset1, offset2);
656  }
657 
658  template<typename MatType>
659  std::vector<casadi_int>
660  SparsityInterface<MatType>::offset(const std::vector< MatType > &v, bool vert) {
661  std::vector<casadi_int> ret(v.size()+1);
662  ret[0]=0;
663  for (casadi_int i=0; i<v.size(); ++i) {
664  ret[i+1] = ret[i] + (vert ? v[i].size1() : v[i].size2());
665  }
666  return ret;
667  }
668 
669  template<typename MatType>
670  MatType SparsityInterface<MatType>::veccat(const std::vector< MatType >& x) {
671  std::vector< MatType > x_vec = x;
672  for (typename std::vector< MatType >::iterator it=x_vec.begin();
673  it!=x_vec.end(); ++it) {
674  *it = vec(*it);
675  }
676  if (x_vec.empty()) {
677  return MatType(0, 1);
678  } else {
679  return vertcat(x_vec);
680  }
681  }
682 
683  template<typename MatType>
684  std::vector< MatType >
685  SparsityInterface<MatType>::diagsplit(const MatType& x,
686  const std::vector<casadi_int>& output_offset) {
687  casadi_assert(x.is_square(), "diagsplit(x,incr)::input must be square but got "
688  + x.dim() + ".");
689  return MatType::diagsplit(x, output_offset, output_offset);
690  }
691 
692  template<typename MatType>
693  std::vector< MatType >
694  SparsityInterface<MatType>::diagsplit(const MatType& x, casadi_int incr) {
695  casadi_assert_dev(incr>=1);
696  casadi_assert(x.is_square(), "diagsplit(x,incr)::input must be square but got "
697  + x.dim() + ".");
698  std::vector<casadi_int> offset2 = range(0, x.size2(), incr);
699  offset2.push_back(x.size2());
700  return MatType::diagsplit(x, offset2);
701  }
702 
703  template<typename MatType>
704  std::vector< MatType >
705  SparsityInterface<MatType>::diagsplit(const MatType& x, casadi_int incr1, casadi_int incr2) {
706  casadi_assert_dev(incr1>=1);
707  casadi_assert_dev(incr2>=1);
708  std::vector<casadi_int> offset1 = range(0, x.size1(), incr1);
709  offset1.push_back(x.size1());
710  std::vector<casadi_int> offset2 = range(0, x.size2(), incr2);
711  offset2.push_back(x.size2());
712  return MatType::diagsplit(x, offset1, offset2);
713  }
714 
715  template<typename MatType>
716  MatType SparsityInterface<MatType>::mtimes(const std::vector<MatType> &args) {
717  casadi_assert(!args.empty(),
718  "mul(std::vector<MatType> &args): "
719  "supplied list must not be empty.");
720  MatType ret = args[0];
721  for (casadi_int i=1; i<args.size(); ++i) ret = MatType::mtimes(ret, args[i]);
722  return ret;
723  }
724 
725  template<typename MatType>
726  std::vector<MatType > SparsityInterface<MatType>::horzsplit(const MatType& x, casadi_int incr) {
727  casadi_assert_dev(incr>=1);
728  casadi_int sz2 = x.size2();
729  std::vector<casadi_int> offset2 = range(0, sz2, incr);
730  offset2.push_back(sz2);
731  return MatType::horzsplit(x, offset2);
732  }
733 
734  template<typename MatType>
735  std::vector<MatType > SparsityInterface<MatType>::vertsplit(const MatType& x, casadi_int incr) {
736  casadi_assert_dev(incr>=1);
737  casadi_int sz1 = x.size1();
738  std::vector<casadi_int> offset1 = range(0, sz1, incr);
739  offset1.push_back(sz1);
740  return MatType::vertsplit(x, offset1);
741  }
742  template<typename MatType>
743  std::vector<MatType > SparsityInterface<MatType>::horzsplit_n(const MatType& x, casadi_int n) {
744  casadi_assert(n>=0, "horzsplit_n(x,n): n (" + str(n) + ") must be non-negative");
745  if (x.size2()==0) return std::vector<MatType>(n, x);
746  casadi_assert(x.size2() % n==0, "horzsplit_n(x,n): x.size2() (" + str(x.size2()) +
747  ") must be a multiple of n (" + str(n) + ")");
748  return horzsplit(x, x.size2()/n);
749  }
750  template<typename MatType>
751  std::vector<MatType > SparsityInterface<MatType>::vertsplit_n(const MatType& x, casadi_int n) {
752  casadi_assert(n>=0, "vertsplit_n(x,n): n (" + str(n) + ") must be non-negative");
753  if (x.size1()==0) return std::vector<MatType>(n, x);
754  casadi_assert(x.size1() % n==0, "vertsplit(x,n): x.size1() (" + str(x.size1()) +
755  ") must be a multiple of n (" + str(n) + ")");
756  return vertsplit(x, x.size1()/n);
757  }
758 
759 #endif // SWIG
760 
761 } // namespace casadi
762 
763 #endif // CASADI_SPARSITY_INTERFACE_HPP
casadi
The casadi namespace.
Definition: archiver.hpp:32