* Bug 15599 fixed: now degree of zero polynomial is -Inf
[scilab.git] / scilab / modules / ast / includes / types / sparse.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 __SPARSE_HH__
17 #define __SPARSE_HH__
18
19 //#include <Eigen/Sparse>
20 #include <complex>
21 #include "double.hxx"
22 #include "bool.hxx"
23 #include "keepForSparse.hxx"
24
25 #define SPARSE_CONST
26
27 namespace Eigen
28 {
29 template<typename _Scalar, int _Flags, typename _Index>  class SparseMatrix;
30 }
31
32 namespace types
33 {
34 /* Utility function to create a new var on the heap from another type
35  */
36 template<typename Dest, typename Arg>
37 Dest* create_new(Arg const&);
38
39 struct SparseBool;
40
41 /**
42    Sparse is a wrapper over Eigen sparse matrices templates for either double or std::complex<double> values.
43  */
44 struct EXTERN_AST Sparse : GenericType
45 {
46     virtual ~Sparse();
47     /* @param src: Double matrix to copy into a new sparse matrix
48     **/
49     Sparse(Double SPARSE_CONST& src);
50     /* @param src : Double matrix to copy into a new sparse matrix
51        @param idx : Double matrix to use as indexes to get values from the src
52     **/
53     Sparse(Double SPARSE_CONST& src, Double SPARSE_CONST& idx);
54     /* @param src : Double matrix to copy into a new sparse matrix
55        @param idx : Double matrix to use as indexes to get values from the src
56        @param dims : Double matrix containing the dimensions of the new matrix
57     **/
58     Sparse(Double SPARSE_CONST& src, Double SPARSE_CONST& idx, Double SPARSE_CONST& dims);
59     /*
60       @param rows : nb of rows of the new matrix
61       @param rows : nb of columns of the new matrix
62       @param cplx : if the new matrix contains complex numbers
63     **/
64     Sparse(int rows, int cols, bool cplx = false);
65     Sparse(Sparse const& o);
66     /* cf. adj2sp()
67       @param xadj : adjacency matrix for the new matrix
68       @param adjncy : adjacency matrix (row indexes) for the new matrix
69       @param src : data for the new matrix
70       @param r : nb of rows for the new matrix
71       @param c : nb of columns for the new matrix
72     **/
73     Sparse(Double SPARSE_CONST& xadj, Double SPARSE_CONST& adjncy, Double SPARSE_CONST& src, std::size_t r, std::size_t c);
74
75     //constructor to create a sparse from value extract to another ( save / load operation typically)
76     Sparse(int rows, int cols, int nonzeros, int* inner, int* outer, double* real, double* img);
77
78     bool isSparse()
79     {
80         return true;
81     }
82     void finalize();
83
84     /*data management member function defined for compatibility with the Double API*/
85     Sparse* set(int _iRows, int _iCols, double _dblReal, bool _bFinalize = true);
86     Sparse* set(int _iIndex, double _dblReal, bool _bFinalize = true)
87     {
88         return set(_iIndex % m_iRows, _iIndex / m_iRows, _dblReal, _bFinalize);
89     }
90
91     Sparse* set(int _iRows, int _iCols, std::complex<double> v, bool _bFinalize = true);
92     Sparse* set(int _iIndex, std::complex<double> v, bool _bFinalize = true)
93     {
94         return set(_iIndex % m_iRows, _iIndex / m_iRows, v, _bFinalize);
95     }
96     /*
97       set non zero values to 1.
98     **/
99     bool one_set();
100     /* get real value at coords (r,c)
101     **/
102     double getReal(int r, int c) const;
103     double getReal(int _iIndex) const
104     {
105         return getReal(_iIndex % m_iRows, _iIndex / m_iRows);
106     }
107
108     double* get();
109     double get(int r, int c) const;
110     double get(int _iIndex) const
111     {
112         return get(_iIndex % m_iRows, _iIndex / m_iRows);
113     }
114
115     std::complex<double>* getImg();
116     std::complex<double> getImg(int r, int c) const;
117     std::complex<double> getImg(int _iIndex) const
118     {
119         return getImg(_iIndex % m_iRows, _iIndex / m_iRows);
120     }
121
122     /* return true if matrix contains complex numbers, false otherwise.
123     **/
124     bool isComplex() const;
125     // overload of GenericType methode.
126     bool isComplex()
127     {
128         // force const to call isComplex const method.
129         const Sparse* sp = this;
130         return sp->isComplex();
131     }
132
133     inline bool isScalar()
134     {
135         return (getRows() == 1 && getCols() == 1);
136     }
137     /* clear all the values of the matrix to 0. (or 0.+0.i if complex)
138     **/
139     bool zero_set();
140
141     /*
142       Config management and GenericType methods overrides
143     */
144     void whoAmI() SPARSE_CONST;
145     bool isExtract() const;
146     Sparse* clone(void);
147     bool toString(std::wostringstream& ostr);
148
149     /* post condition: dimensions are at least _iNewRows, _iNewCols
150        preserving existing data.
151        If dimensions where already >=, this is a no-op.
152
153        @param _iNewRows new minimum nb of rows
154        @param _iNewCols new minimum nb of cols
155        @return true upon succes, false otherwise.
156      */
157     Sparse* resize(int _iNewRows, int _iNewCols);
158     /* post condition: new total size must be equal to the old size.
159                        Two dimensions maximum.
160
161        @param _iNewRows new nb of rows
162        @param _iNewCols new nb of cols
163        @param _piNewDims new nb of dimension
164        @param _iNewDims new size for each dimension
165        @return true upon succes, false otherwise.
166     */
167     Sparse* reshape(int* _piNewDims, int _iNewDims);
168     Sparse* reshape(int _iNewRows, int _iNewCols);
169     /*
170       insert _iSeqCount elements from _poSource at coords given by _piSeqCoord (max in _piMaxDim).
171       coords are considered 1D if _bAsVector, 2D otherwise.
172       @param _iSeqCount nb of elts to insert
173       @param _piSeqCoord dest coords
174       @param _poSource src
175       @param  _bAsVector if _piSeqCoord contains 1D coords.
176      */
177     Sparse* insert(typed_list* _pArgs, InternalType* _pSource);
178
179     GenericType* remove(typed_list* _pArgs);
180
181     GenericType* insertNew(typed_list* _pArgs);
182
183     /* append _poSource from coords _iRows, _iCols
184        @param _iRows row to append from
185        @param _iCols col to append from
186        @param _poSource src data to append
187      */
188     Sparse* append(int r, int c, types::Sparse SPARSE_CONST* src);
189
190     /*
191       extract a submatrix
192       @param _iSeqCount nb of elts to extract
193       @param _piSeqCoord src coords
194       @param _piMaxDim max coords
195       @param _piDimSize size of the extracted matrix
196       @param  _bAsVector if _piSeqCoord contains 1D coords.
197
198      */
199     GenericType* extract(typed_list* _pArgs);
200     Sparse* extract(int _iSeqCount, int* _piSeqCoord, int* _piMaxDim, int* _piDimSize, bool _bAsVector) SPARSE_CONST;
201     virtual bool invoke(typed_list & in, optional_list & /*opt*/, int /*_iRetCount*/, typed_list & out, const ast::Exp & e);
202     virtual bool isInvokable() const;
203     virtual bool hasInvokeOption() const;
204     virtual int getInvokeNbIn();
205     virtual int getInvokeNbOut();
206
207     /*
208        change the sign (inplace).
209      */
210     void opposite();
211
212     /*
213       compares with an other value for equality (same nb of elts, with same values)
214       TODO: should it handle other types ?
215      */
216     bool operator==(const InternalType& it) SPARSE_CONST;
217     /*
218       compares with an other value for inequality (same nb of elts, with same values)
219       TODO: should it handle other types ?
220      */
221     bool operator!=(const InternalType& it) SPARSE_CONST
222     {
223         return !(*this == it);
224     }
225
226
227     /* return type as string ( double, int, cell, list, ... )*/
228     virtual std::wstring         getTypeStr() const
229     {
230         return std::wstring(L"sparse");
231     }
232
233     /* return type as short string ( s, i, ce, l, ... ), as in overloading macros*/
234     virtual std::wstring         getShortTypeStr() const
235     {
236         return std::wstring(L"sp");
237     }
238
239     /* create a new sparse matrix containing the result of an addition
240        @param o other matrix to add
241        @return ptr to the new matrix, 0 in case of failure
242      */
243     Sparse* add(Sparse const& o) const;
244
245     /* create a new sparse matrix containing the result of an addition
246        @param d scalar to add
247        @return ptr to the new matrix, 0 in case of failure
248      */
249     Double* add(double d) const;
250
251     /* create a new sparse matrix containing the result of an addition
252        @param c complex  to add
253        @return ptr to the new matrix, 0 in case of failure
254      */
255     Double* add(std::complex<double> c) const;
256
257
258
259
260     /* create a new sparse matrix containing the result of a substraction
261        @param o other matrix to substract
262        @return ptr to the new matrix, 0 in case of failure
263      */
264     Sparse* substract(Sparse const& o) const;
265
266     /* create a new sparse matrix containing the result of an subtraction
267        @param d scalar to subtract
268        @return ptr to the new matrix, 0 in case of failure
269      */
270     Double* substract(double d) const;
271
272     /* create a new sparse matrix containing the result of an subtraction
273        @param c scalar to subtract
274        @return ptr to the new matrix, 0 in case of failure
275      */
276     Double* substract(std::complex<double> c) const;
277
278
279     /* create a new sparse matrix containing the result of a multiplication
280        @param o other matrix to substract
281        @return ptr to the new matrix, 0 in case of failure
282      */
283     Sparse* multiply(Sparse const& o) const;
284
285     /* create a new sparse matrix containing the result of an multiplication
286        @param s scalar to multiply by
287        @return ptr to the new matrix, 0 in case of failure
288      */
289     Sparse* multiply(double s) const;
290
291     /* create a new sparse matrix containing the result of an multiplication
292        @param c scalar to subtract
293        @return ptr to the new matrix, 0 in case of failure
294      */
295     Sparse* multiply(std::complex<double> c) const;
296
297     /* create a new matrix containing the result of an .*
298        @param o sparse matrix to .*
299        @return ptr to the new matrix, 0 in case of failure
300      */
301     Sparse* dotMultiply(Sparse SPARSE_CONST& o) const;
302
303     /* create a new matrix containing the result of an ./
304       @param o sparse matrix to ./
305       @return ptr to the new matrix, 0 in case of failure
306     */
307     Sparse* dotDivide(Sparse SPARSE_CONST& o) const;
308
309     bool neg(InternalType *& out);
310
311     bool transpose(InternalType *& out);
312     bool adjoint(InternalType *& out);
313     int newCholLLT(Sparse** permut, Sparse** factor) const;
314
315     /** create a new sparse matrix containing the non zero values set to 1.
316        equivalent but faster than calling one_set() on a new copy of the
317        current matrix.
318        @return ptr to the new matrix, 0 in case of failure
319      */
320     Sparse* newOnes() const;
321
322     Sparse* newReal() const;
323     /** @return the nb of non zero values.
324      */
325     std::size_t nonZeros() const;
326
327     /* @param i row of the current sparse matrix
328        @return the nb of non zero values in row i
329      */
330     std::size_t nonZeros(std::size_t i) const;
331
332     int* getNbItemByRow(int* _piNbItemByRows);
333     int* getColPos(int* _piColPos);
334     int* getInnerPtr(int* count);
335     int* getOuterPtr(int* count);
336
337
338     /**
339        "in-place" cast into a sparse matrix of comlpex values
340      */
341     void toComplex();
342
343     /* coefficient wise relational operator < between *this sparse matrix and an other.
344        Matrices must have the same dimensions except if one of them is of size (1,1)
345        (i.e. a scalar) : it is then treated as a constant matrix of thre required dimensions.
346
347        @param o other sparse matrix
348
349        @return ptr to a new Sparse matrix where each element is the result of the logical operator
350         '<' between the elements of *this and those of o.
351      */
352     SparseBool* newLessThan(Sparse &o);
353
354     /* coefficient wise relational operator != between *this sparse matrix and an other.
355        Matrices must have the same dimensions except if one of them is of size (1,1)
356        (i.e. a scalar) : it is then treated as a constant matrix of thre required dimensions.
357
358        @param o other sparse matrix
359
360        @return ptr to a new Sparse matrix where each element is the result of the logical operator
361         '!=' between the elements of *this and those of o.
362      */
363     SparseBool* newNotEqualTo(Sparse const&o) const;
364
365     /* coefficient wise relational operator <= between *this sparse matrix and an other.
366        Matrices must have the same dimensions except if one of them is of size (1,1)
367        (i.e. a scalar) : it is then treated as a constant matrix of thre required dimensions.
368
369        Do not use this function is possible as the result will be dense because
370        0. <= 0. is true, hence the result matrix will hold a non default value (i.e. true)
371        for each pair of default values (0.) of the sparse arguments !
372
373        @param o other sparse matrix
374
375        @return ptr to a new Sparse matrix where each element is the result of the logical operator
376         '<=' between the elements of *this and those of o.
377      */
378     SparseBool* newLessOrEqual(Sparse &o);
379
380     /* coefficient wise relational operator == between *this sparse matrix and an other.
381        Matrices must have the same dimensions except if one of them is of size (1,1)
382        (i.e. a scalar) : it is then treated as a constant matrix of thre required dimensions.
383
384        Do not use this function is possible as the result will be dense because
385        0. == 0. is true, hence the result matrix will hold a non default value (i.e. true)
386        for each pair of default values (0.) of the sparse arguments !
387
388        @param o other sparse matrix
389
390        @return ptr to a new Sparse matrix where each element is the result of the logical operator
391         '==' between the elements of *this and those of o.
392      */
393     SparseBool* newEqualTo(Sparse &o);
394
395     /**
396        output 1-base column numbers of the non zero elements
397        @param out : ptr used as an output iterator over double values
398        @return past-the-end output iterator after ouput is done
399      */
400     double* outputCols(double* out) const;
401
402     /**
403        output real and imaginary values of the non zero elements
404        @param outReal : ptr used as an output iterator over double values for real values
405        @param outImag : ptr used as an output iterator over double values for imaginary values if any
406        @return pair of past-the-end output iterators after ouput is done
407      */
408     std::pair<double*, double*> outputValues(double* outReal, double* outImag)const;
409
410     /**
411        ouput rows and afterwards columns of the non zero elements
412        @param out : ptr used as an output iterator over double values
413        @return past-the-end output iterators after ouput is done
414      */
415     int* outputRowCol(int* out)const;
416
417     /**
418        @param dest Double to be filled with values from the current sparse matrix.
419      */
420     void fill(Double& dest, int r = 0, int c = 0) SPARSE_CONST;
421
422
423     inline ScilabType  getType(void) SPARSE_CONST
424     {
425         return ScilabSparse;
426     }
427
428     inline ScilabId    getId(void) SPARSE_CONST
429     {
430         if (isComplex())
431         {
432             return IdSparseComplex;
433         }
434         return IdSparse;
435     }
436
437
438
439     SparseBool* newLesserThan(Sparse const&o);
440
441     typedef Eigen::SparseMatrix<double, 0x1, int>                   RealSparse_t;
442     typedef Eigen::SparseMatrix<std::complex<double>, 0x1, int>     CplxSparse_t;
443     /**
444        One and only one of the args should be 0.
445        @param realSp ptr to an Eigen sparse matrix of double values
446        @param cplxSp ptr to an Eigen sparse matrix of std::complex<double> elements
447      */
448     Sparse(RealSparse_t* realSp, CplxSparse_t* cplxSp);
449
450     RealSparse_t* matrixReal;
451     CplxSparse_t* matrixCplx;
452
453 protected :
454 private :
455
456     /** utility function used by constructors
457         @param rows : nb of rows
458         @param cols : nb of columns
459         @param src : Double matrix data source
460         @param : iterator (cf MatrixIterator.hxx) with indices
461         @param n : nb of elements to copy from data source.
462      */
463     template<typename DestIter>
464     void create(int rows, int cols, Double SPARSE_CONST& src, DestIter o, std::size_t n);
465     void create2(int rows, int cols, Double SPARSE_CONST& src, Double SPARSE_CONST& idx);
466
467     /** utility function used by insert functions conceptually : sp[destTrav]= src[srcTrav]
468         @param src : data source
469         @param SrcTrav : iterator over the data source
470         @param n : nb of elements to copy
471         @param sp : sparse destination matrix
472         @param destTrav : iterator over the data sink (i.e. sp)
473      */
474     template<typename Src, typename SrcTraversal, typename Sz, typename DestTraversal>
475     static bool copyToSparse(Src SPARSE_CONST& src, SrcTraversal srcTrav, Sz n, Sparse& sp, DestTraversal destTrav);
476
477     Sparse* insert(typed_list* _pArgs, Sparse* _pSource);
478 };
479
480 template<typename T>
481 void neg(const int r, const int c, const T * const in, Eigen::SparseMatrix<bool, 1, int> * const out);
482
483
484 /*
485   Implement sparse boolean matrix
486  */
487 struct EXTERN_AST SparseBool : GenericType
488 {
489     virtual ~SparseBool();
490     /* @param src: Bool matrix to copy into a new sparse matrix
491     **/
492     SparseBool(Bool SPARSE_CONST& src);
493     /* @param src : Bool matrix to copy into a new sparse matrix
494        @param idx : Double matrix to use as indexes to get values from the src
495     **/
496     SparseBool(Bool SPARSE_CONST& src, Double SPARSE_CONST& idx);
497     /* @param src : Bool matrix to copy into a new sparse matrix
498        @param idx : Double matrix to use as indexes to get values from the src
499        @param dims : Double matrix containing the dimensions of the new matrix
500     **/
501     SparseBool(Bool SPARSE_CONST& src, Double SPARSE_CONST& idx, Double SPARSE_CONST& dims);
502     /*
503        @param rows : nb of rows of the new matrix
504        @param rows : nb of columns of the new matrix
505     */
506     SparseBool(int rows, int cols);
507
508     SparseBool(SparseBool const& o);
509
510     //constructor to create a sparse from value extract to another ( save / load operation typically)
511     SparseBool(int rows, int cols, int trues, int* inner, int* outer);
512
513     bool isSparseBool()
514     {
515         return true;
516     }
517     void finalize();
518
519     bool toString(std::wostringstream& ostr);
520
521     /* Config management and GenericType methods overrides */
522     SparseBool* clone(void);
523
524     SparseBool* resize(int _iNewRows, int _iNewCols);
525     SparseBool* reshape(int* _piNewDims, int _iNewDims);
526     SparseBool* reshape(int _iNewRows, int _iNewCols);
527     SparseBool* insert(typed_list* _pArgs, InternalType* _pSource);
528     SparseBool* append(int _iRows, int _iCols, SparseBool SPARSE_CONST* _poSource);
529
530     GenericType* remove(typed_list* _pArgs);
531     GenericType* insertNew(typed_list* _pArgs);
532     GenericType* extract(typed_list* _pArgs);
533
534     SparseBool* extract(int _iSeqCount, int* _piSeqCoord, int* _piMaxDim, int* _piDimSize, bool _bAsVector) SPARSE_CONST;
535
536     virtual bool invoke(typed_list & in, optional_list &/*opt*/, int /*_iRetCount*/, typed_list & out, const ast::Exp & e);
537     virtual bool isInvokable() const;
538     virtual bool hasInvokeOption() const;
539     virtual int getInvokeNbIn();
540     virtual int getInvokeNbOut();
541
542     bool transpose(InternalType *& out);
543
544     /** @return the nb of non zero values.
545      */
546     std::size_t nbTrue() const;
547     /* @param i row of the current sparse matrix
548        @return the nb of non zero values in row i
549      */
550     std::size_t nbTrue(std::size_t i) const;
551
552     void setTrue(bool finalize = true);
553     void setFalse(bool finalize = true);
554
555     int* getNbItemByRow(int* _piNbItemByRows);
556     int* getColPos(int* _piColPos);
557     int* getInnerPtr(int* count);
558     int* getOuterPtr(int* count);
559     /**
560        output 1-base column numbers of the non zero elements
561        @param out : ptr used as an output iterator over double values
562        @return past-the-end output iterator after ouput is done
563      */
564
565     int* outputRowCol(int* out)const;
566
567     bool operator==(const InternalType& it) SPARSE_CONST;
568     bool operator!=(const InternalType& it) SPARSE_CONST;
569
570     /* return type as string ( double, int, cell, list, ... )*/
571     virtual std::wstring getTypeStr() const
572     {
573         return std::wstring(L"boolean sparse");
574     }
575     /* return type as short string ( s, i, ce, l, ... )*/
576     virtual std::wstring getShortTypeStr() const
577     {
578         return std::wstring(L"spb");
579     }
580
581     inline ScilabType getType(void) SPARSE_CONST
582     {
583         return ScilabSparseBool;
584     }
585
586     inline ScilabId getId(void) SPARSE_CONST
587     {
588         return IdSparseBool;
589     }
590
591     inline bool isScalar()
592     {
593         return (getRows() == 1 && getCols() == 1);
594     }
595
596     bool isTrue()
597     {
598         if (static_cast<int>(nbTrue()) == m_iSize)
599         {
600             return true;
601         }
602         return false;
603     }
604
605     bool neg(InternalType *& out)
606     {
607         SparseBool * _out = new SparseBool(getRows(), getCols());
608         types::neg(getRows(), getCols(), matrixBool, _out->matrixBool);
609         _out->finalize();
610         out = _out;
611         return true;
612     }
613
614     void whoAmI() SPARSE_CONST;
615
616     bool* get();
617     bool get(int r, int c) SPARSE_CONST;
618     bool get(int _iIndex) SPARSE_CONST
619     {
620         return get(_iIndex % m_iRows, _iIndex / m_iRows);
621     }
622
623     SparseBool* set(int r, int c, bool b, bool _bFinalize = true) SPARSE_CONST;
624     SparseBool* set(int _iIndex, bool b, bool _bFinalize = true) SPARSE_CONST
625     {
626         return set(_iIndex % m_iRows, _iIndex / m_iRows, b, _bFinalize);
627     }
628
629     void fill(Bool& dest, int r = 0, int c = 0) SPARSE_CONST;
630
631     Sparse* newOnes() const;
632     SparseBool* newNotEqualTo(SparseBool const&o) const;
633     SparseBool* newEqualTo(SparseBool& o);
634
635     SparseBool* newLogicalOr(SparseBool const&o) const;
636     SparseBool* newLogicalAnd(SparseBool const&o) const;
637
638     typedef Eigen::SparseMatrix<bool, 0x1, int> BoolSparse_t;
639     SparseBool(BoolSparse_t* o);
640     BoolSparse_t* matrixBool;
641
642 private:
643     void create2(int rows, int cols, Bool SPARSE_CONST& src, Double SPARSE_CONST& idx);
644     SparseBool* insert(typed_list* _pArgs, SparseBool* _pSource);
645 };
646
647 template<typename T>
648 struct SparseTraits
649 {
650     typedef types::Sparse type;
651 };
652 template<>
653 struct SparseTraits<types::Bool>
654 {
655     typedef types::SparseBool type;
656 };
657 template<>
658 struct SparseTraits<types::SparseBool>
659 {
660     typedef types::SparseBool type;
661 };
662
663 }
664
665 #endif /* !__SPARSE_HH__ */