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