a4d2295a7d854113d3c22568c3f0a5583ac79771
[scilab.git] / scilab / modules / ast / includes / types / sparseOp.hxx
1 /*
2  *  Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3  *  Copyright (C) 2010-2010 - DIGITEO - Bernard Hugueney
4  *
5  * Copyright (C) 2012 - 2016 - Scilab Enterprises
6  *
7  * This file is hereby licensed under the terms of the GNU GPL v2.0,
8  * pursuant to article 5.3.4 of the CeCILL v.2.1.
9  * This file was originally licensed under the terms of the CeCILL v2.1,
10  * and continues to be available under such terms.
11  * For more information, see the COPYING file which you should have received
12  * along with this program.
13  *
14  */
15
16 #ifndef __SPARSEOP_HH__
17 #define __SPARSEOP_HH__
18
19 #include "sparse.hxx"
20 #include "keepForSparse.hxx"
21
22
23 /*
24  * Logical operators implementation for types::Sparse and types::SparseBool
25  *
26  * Many operations result in sparse matrices with a lot of true elements
27  * they are slow and take a lot of memory to implement.
28  *
29  * Mixed mode operators between real and complex values are commented out
30  * until Eigen::Sparse provides a .coeff() member function to cast<> retults.
31  *
32  * complex ordering is selected through conditional compilation to show
33  * how SEP 21 could be implemented.
34  *
35  * pruning [c|sh]ould be tied to relative precision as defined in elem_common.h
36  *
37  */
38
39 namespace std
40 {
41 /* std namespace does not implement complex ordering: we use the lexicographical order   */
42 typedef double (*cplxFun_t)(std::complex<double> const&);
43 #ifndef SCI_CPLX_ORDER_SEP21
44 /* on abs() and arg() according to current practise */
45 cplxFun_t const mostSignificant (&std::abs);
46 cplxFun_t const leastSignificant(&std::arg);
47 #else
48 /*  on real and imaginary parts as per SEP 21 recommendation.  */
49 cplxFun_t const mostSignificant (&std::real);
50 cplxFun_t const leastSignificant(&std::imag);
51 #endif
52
53 template<> struct less<std::complex<double> >: std::binary_function<std::complex<double>, std::complex<double>, bool>
54 {
55     bool operator()(std::complex<double> const& op1, std::complex<double> const& op2) const
56     {
57         if (mostSignificant(op1) < mostSignificant(op2))
58         {
59             return true;
60         }
61         if (mostSignificant(op2) < mostSignificant(op1))
62         {
63             return false;
64         }
65         return leastSignificant(op1) < leastSignificant(op2);
66     }
67 };
68 template<> struct greater<std::complex<double> >: std::binary_function<std::complex<double>, std::complex<double>, bool>
69 {
70     bool operator()(std::complex<double> const& op1, std::complex<double> const& op2) const
71     {
72         if (mostSignificant(op1) > mostSignificant(op2))
73         {
74             return true;
75         }
76         if (mostSignificant(op2) > mostSignificant(op1))
77         {
78             return false;
79         }
80         return leastSignificant(op1) > leastSignificant(op2);
81     }
82 };
83
84 template<> struct less_equal<std::complex<double> >: std::binary_function<std::complex<double>, std::complex<double>, bool>
85 {
86     bool operator()(std::complex<double> const& op1, std::complex<double> const& op2) const
87     {
88         return !std::greater<std::complex<double> >()(op2, op1);
89     }
90 };
91 template<> struct greater_equal<std::complex<double> >: std::binary_function<std::complex<double>, std::complex<double>, bool>
92 {
93     bool operator()(std::complex<double> const& op1, std::complex<double> const& op2) const
94     {
95         return !std::less<std::complex<double> >()(op2, op1);
96     }
97 };
98 }
99
100 namespace
101 {
102 /*
103  * Operations on sparse matrices can require full traversal of the sparse args if the default (0 or false) value
104  * of one arg does not set the result to the default result value.
105  * For example, '==' does require full traversal because "0. == 0" is true (and the default SparseBool value is false)
106  * However, '!=' does not require it because "0. != 0." is false so the default value of the sparse result is good.
107  * we dispatch according the operators thanks to the OperatorTraits struct (as is usually done for example to
108  * dispatch amongst std:: algorithms implementations according to the iterator_category in std::iterator_traits<>
109  */
110 struct FullTraversal {};
111 struct UnionTraversal {};
112
113 struct WithFullTraversal
114 {
115     typedef FullTraversal traversal;
116 };
117
118 struct WithUnionTraversal
119 {
120     typedef UnionTraversal traversal;
121 };
122
123 template<typename T> struct OperatorTraits : WithFullTraversal
124 {
125 };
126
127
128 template<typename Scalar> struct OperatorTraits<std::less<Scalar> > : WithUnionTraversal
129 {
130 };
131 template<typename Scalar> struct OperatorTraits<std::greater<Scalar> > : WithUnionTraversal
132 {
133 };
134 template<typename Scalar> struct OperatorTraits<std::not_equal_to<Scalar> > : WithUnionTraversal
135 {
136 };
137 template<typename Scalar> struct OperatorTraits<std::less_equal<Scalar> > : WithFullTraversal
138 {
139 };
140 template<typename Scalar> struct OperatorTraits<std::greater_equal<Scalar> > : WithFullTraversal
141 {
142 };
143 template<typename Scalar> struct OperatorTraits<std::equal_to<Scalar> > : WithFullTraversal
144 {
145 };
146
147 // would be even better with intersectionTraversal if Eigen exposed such traversal
148 template<typename Scalar> struct OperatorTraits<std::logical_and<Scalar> > : WithUnionTraversal
149 {
150 };
151
152 template<typename Scalar> struct OperatorTraits<std::logical_or<Scalar> > : WithUnionTraversal
153 {
154 };
155
156 /*
157  * Special case when the second operand is a Scalar
158  * @param op1 Eigen first operand
159  * @param op2 scalar second operand
160  * @param op binary operation to perform
161  * @return ptr to the new Eigen::Sparse result containing the cwise binary op with the scalar.
162  */
163 template<typename Sp1, typename Scalar2, typename Op>
164 Eigen::SparseMatrix<typename Op::result_type, Eigen::RowMajor>* scalarOp(Eigen::EigenBase<Sp1> const& op1, Scalar2 op2, Op op)
165 {
166     typedef typename Eigen::internal::traits<Sp1>::Scalar Scalar1;
167     typedef typename Op::result_type result_scalar;
168     typedef Eigen::SparseMatrix<result_scalar, Eigen::RowMajor> result_t;
169     result_t* res;
170     if (op(Scalar1(), op2) == result_scalar())
171     {
172         res = new result_t(op1.derived().unaryExpr(std::bind2nd(op, op2)));
173         res->prune(&keepForSparse<result_scalar>);
174     }
175     else
176     {
177         /* it is not possible to  perform full traversal (including virtual 0 elts) of the sparse
178          matrix hence the need for a dense tmp :( */
179         // TODO remove dense temp when Eigen provides sparse full traversal API
180         Eigen::Matrix<Scalar1, Eigen::Dynamic, Eigen::Dynamic> tmp(op1);
181         res = new result_t(tmp.unaryExpr(std::bind2nd(op, op2)).sparseView());
182     }
183     return res;
184 }
185
186 /*
187  * Special case when the first operand is a Scalar
188  * @param op1 scalar first operand
189  * @param op2 Eigen second operand
190  * @param op binary operation to perform
191  * @return ptr to the new Eigen::Sparse result containing the cwise binary op with the scalar.
192  */
193 template<typename Scalar1, typename Sp2, typename Op>
194 Eigen::SparseMatrix<typename Op::result_type, Eigen::RowMajor>* scalarOp(Scalar1 op1, Eigen::EigenBase<Sp2> const& op2, Op op)
195 {
196     typedef typename Eigen::internal::traits<Sp2>::Scalar Scalar2;
197     typedef typename Op::result_type result_scalar;
198     typedef Eigen::SparseMatrix<result_scalar, Eigen::RowMajor> result_t;
199     result_t* res;
200     if (op(op1, Scalar2()) == result_scalar())
201     {
202         res = new result_t(op2.derived().unaryExpr(std::bind1st(op, op1)));
203         res->prune(&keepForSparse<result_scalar>);
204     }
205     else
206     {
207         // TODO remove dense temp when Eigen provides sparse full traversal API
208         Eigen::Matrix<Scalar2, Eigen::Dynamic, Eigen::Dynamic> tmp(op2);
209         res = new result_t(tmp.unaryExpr(std::bind1st(op, op1)).sparseView());
210     }
211     return res;
212 }
213 /*
214  * Common case between two Eigen matrices. Performs the dispatch to scalar specialized versions
215  * or according to the traversal.
216  *
217  * @param op1 Eigen first operand
218  * @param op2 Eigen second operand
219  * @param op binary operation to perform
220  * @return ptr to the new Eigen::Sparse result containing the cwise binary op with the two matrices.
221  */
222 template< typename Sp1, typename Sp2 , typename Op>
223 Eigen::SparseMatrix<typename Op::result_type, Eigen::RowMajor>* cwiseOp(Eigen::EigenBase<Sp1> const& op1, Eigen::EigenBase<Sp2> const& op2, Op op)
224 {
225     if (op1.rows() == 1 && op1.cols() == 1)
226     {
227         return scalarOp(op1.derived().coeff(0, 0), op2, op);
228     }
229     if (op2.rows() == 1 && op2.cols() == 1)
230     {
231         return scalarOp(op1, op2.derived().coeff(0, 0), op);
232     }
233     return cwiseOp(op1, op2, op, typename OperatorTraits<Op>::traversal());
234 }
235 /*
236  * Performs the cwise binary operator between two non scalar matrices, with full traversal.
237  * SLOW and using memory, but you asked for it by using inapropriate operators.
238  *
239  * @param op1 Eigen first operand
240  * @param op2 Eigen second operand
241  * @param op binary operation to perform
242  * @return ptr to the new Eigen::Sparse result containing the cwise binary op between the two matrices
243  */
244 template< typename Sp1, typename Sp2 , typename Op>
245 Eigen::SparseMatrix<typename Op::result_type, Eigen::RowMajor>* cwiseOp(Eigen::EigenBase<Sp1> const& op1, Eigen::EigenBase<Sp2> const& op2, Op op, FullTraversal /*unused*/)
246 {
247     typedef typename Op::result_type result_scalar;
248     typedef Eigen::SparseMatrix<result_scalar, Eigen::RowMajor> result_t;
249     typedef typename Eigen::internal::traits<Sp1>::Scalar Scalar1;
250     typedef typename Eigen::internal::traits<Sp2>::Scalar Scalar2;
251     // TODO remove dense temp when Eigen provides sparse full traversal API
252     Eigen::Matrix<Scalar1, Eigen::Dynamic, Eigen::Dynamic> tmp1(op1);
253     Eigen::Matrix<Scalar2, Eigen::Dynamic, Eigen::Dynamic> tmp2(op2);
254
255     return new result_t(tmp1.binaryExpr(tmp2, op).sparseView());
256 }
257 /*
258  * Performs the cwise binary operator between two non scalar matrices, with full traversal.
259  * SLOW and using memory, but you asked for it by using inapropriate operators.
260  *
261  * @param op1 Eigen first operand
262  * @param op2 Eigen second operand
263  * @param op binary operation to perform
264  * @return ptr to the new Eigen::Sparse result containing the cwise binary op with the scalar.
265  */
266 template< typename Sp1, typename Sp2 , typename Op>
267 Eigen::SparseMatrix<typename Op::result_type, Eigen::RowMajor>* cwiseOp(Eigen::EigenBase<Sp1> const& op1, Eigen::EigenBase<Sp2> const& op2, Op op, UnionTraversal /*unused*/)
268 {
269     typedef typename Op::result_type result_scalar;
270     typedef Eigen::SparseMatrix<result_scalar, Eigen::RowMajor> result_t;
271     result_t* res(new result_t(op1.derived().binaryExpr(op2.derived(), op)));
272     res->prune(&keepForSparse<result_scalar>);
273     return res;
274 }
275
276 /*
277  * Performs the cwise binary operator between two types::Sparse arguments
278  * dispatching according to whether they are complex or not.
279  *
280  * @param op1 types::Sparse first operand
281  * @param op2 types::Sparse second operand
282  * @param op binary operation to perform
283  * @return ptr to the new types::SparseBool result containing the cwise binary op with the scalar.
284  */
285 template<template <typename ElementType> class Op>
286 types::SparseBool* cwiseOp(types::Sparse const& op1, types::Sparse const& op2)
287 {
288     types::SparseBool::BoolSparse_t* res(0);
289     if (op1.isComplex())
290     {
291         if (op2.isComplex())
292         {
293             res = cwiseOp(*op1.matrixCplx, *op2.matrixCplx, Op<std::complex<double> >());
294         }
295         else
296         {
297             types::Sparse temp(op2);
298             temp.toComplex();
299             res = cwiseOp(*op1.matrixCplx, *temp.matrixCplx, Op<std::complex<double> >());
300         }
301     }
302     else
303     {
304
305         if (op2.isComplex())
306         {
307             types::Sparse temp(op1);
308             temp.toComplex();
309             res = cwiseOp(*temp.matrixCplx, *op2.matrixCplx, Op<std::complex<double> >());
310         }
311         else
312         {
313             res = cwiseOp(*op1.matrixReal, *op2.matrixReal, Op<double>());
314         }
315     }
316     return new types::SparseBool(res);
317 }
318
319 template<template <typename ElementType> class Op>
320 types::SparseBool* cwiseOp(types::SparseBool const& op1, types::SparseBool const& op2)
321 {
322     return new types::SparseBool(cwiseOp(*op1.matrixBool, *op2.matrixBool, Op<bool>()));
323 }
324
325 }
326 #endif