Merge remote-tracking branch 'origin/master' into jit
[scilab.git] / scilab / modules / ast / src / cpp / types / sparse.cpp
1 /*
2 *  Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 *  Copyright (C) 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 #include <sstream>
14 #include <math.h>
15 #include <Eigen/Sparse>
16 #include <complex>
17 #include <iterator>
18 #include <algorithm>
19
20 #include <Eigen/Core>
21 #include <Eigen/IterativeLinearSolvers>
22 #include <Eigen/SparseCholesky>
23
24 #include "sparse.hxx"
25 #include "types.hxx"
26 #include "tostring_common.hxx"
27 #include "double.hxx"
28 #include "matrixiterator.hxx"
29 #include "types_subtraction.hxx"
30 #include "types_addition.hxx"
31 #include "types_multiplication.hxx"
32 #include "configvariable.hxx"
33 #include "scilabWrite.hxx"
34 #include "exp.hxx"
35
36 #include "sparseOp.hxx"
37
38 extern "C"
39 {
40 #include "elem_common.h"
41 }
42 namespace
43 {
44
45 /* used for debuging output
46 */
47 template<typename Os, typename In, typename Sz> Os& writeData(wchar_t const* title, In beg, Sz n, Os& os)
48 {
49     os << title;
50     /* TODO: use tostring_common (with a kind of std::boolalpha for boolean output)
51     */
52     mycopy_n(beg, n, std::ostream_iterator<typename std::iterator_traits<In>::value_type, char>(os, L" "));
53     os << std::endl;
54     return os;
55 }
56
57 struct Printer
58 {
59     Printer (int precision) : p(precision)
60     {
61     }
62     template<typename T>
63     std::wstring emptyName( /* */) const
64     {
65         return L" zero";
66     }
67
68     template<typename T>
69     std::wstring operator()(T const& t) const
70     {
71         //never call ?
72         std::wostringstream ostr;
73         ostr.precision(p);
74         ostr << t;
75         return ostr.str();
76     }
77     int p;
78 };
79
80 template<>
81 std::wstring Printer::operator()(bool const& b) const
82 {
83     if (b)
84     {
85         return L"T";
86     }
87     else
88     {
89         return L"F";
90     }
91 }
92
93 template<>
94 std::wstring Printer::operator()(double const& d) const
95 {
96     std::wostringstream ostr;
97     DoubleFormat df;
98     getDoubleFormat(d, &df);
99     addDoubleValue(&ostr, d, &df);
100     return ostr.str();
101 }
102
103 template<>
104 std::wstring Printer::operator()(std::complex<double > const& c) const
105 {
106     std::wostringstream ostr;
107     int iLen = 0;
108     DoubleFormat dfR, dfI;
109     getComplexFormat(c.real(), c.imag(), &iLen, &dfR, &dfI);
110     addDoubleComplexValue(&ostr, c.real(), c.imag(), iLen, &dfR, &dfI);
111     return ostr.str();
112 }
113
114 template<>
115 std::wstring Printer::emptyName<bool>() const
116 {
117     return L"False";
118 }
119
120
121 template<typename T> std::wstring toString(T const& m, int precision)
122 {
123     std::wostringstream ostr;
124
125     int iWidthRows  = 0;
126     int iWidthCols  = 0;
127     getSignedIntFormat(m.rows(), &iWidthRows);
128     getSignedIntFormat(m.cols(), &iWidthCols);
129
130     ostr << L"(" ;
131     addUnsignedIntValue<unsigned long long>(&ostr, m.rows(), iWidthRows);
132     ostr << ",";
133     addUnsignedIntValue<unsigned long long>(&ostr, m.cols(), iWidthCols);
134     ostr << L")";
135
136     Printer p(precision);
137     if (!m.nonZeros())
138     {
139         ostr << ( p.emptyName<typename Eigen::internal::traits<T>::Scalar>());
140     }
141     ostr << L" sparse matrix\n\n";
142
143     auto * pIColPos      = m.innerIndexPtr();
144     auto * pINbItemByRow = m.outerIndexPtr();
145
146     int iPos = 0;
147
148     for (size_t j = 1 ; j < m.rows() + 1 ; j++)
149     {
150         for (size_t i = pINbItemByRow[j - 1] ; i < pINbItemByRow[j] ; i++)
151         {
152             ostr << L"(";
153             addUnsignedIntValue<unsigned long long>(&ostr, (int)j, iWidthRows);
154             ostr << L",";
155             addUnsignedIntValue<unsigned long long>(&ostr, pIColPos[iPos] + 1, iWidthCols);
156             ostr << L")\t" << p(m.valuePtr()[iPos]) << std::endl;
157
158             iPos++;
159         }
160     }
161
162     return ostr.str();
163 }
164
165 /** utility function to compare two Eigen::Sparse matrices to equality
166 */
167 template<typename T> bool equal(T const& s1, T const& s2)
168 {
169     bool res(true);
170     // only compares elts when both inner iterators are "defined", so we assert that we compared all the non zero values
171     // i.e. the inner iterators where defined for the same values
172     std::size_t nbElts(0);
173
174     for (int k = 0; res && k != s1.outerSize(); ++k)
175     {
176         for (typename T::InnerIterator it1(s1, k), it2(s2, k); res && it1 && it2 ; ++it1, ++it2, ++nbElts)
177         {
178             res = (it1.value() == it2.value()
179                    && it1.row() == it2.row()
180                    && it1.col() == it2.col());
181         }
182     }
183     return res && (nbElts == s1.nonZeros()) && (nbElts == s2.nonZeros());
184 }
185 /**
186 utility function to set non zero values of an Eigen::Sparse matrix to a fixed values
187 @param s : sparse matrix to modify
188 @param v : value to set (default to 1.)
189 */
190 template<typename T> bool setNonZero(T& s, typename Eigen::internal::traits<T>::Scalar v = 1.)
191 {
192     for (auto j = 0; j < s.outerSize(); ++j)
193     {
194         for (typename T::InnerIterator it(s, j); it; ++it)
195         {
196             it.valueRef() = v;
197         }
198     }
199     return true;
200 }
201
202
203
204 template<typename Src, typename Sp>
205 void doAppend(Src SPARSE_CONST& src, int r, int c, Sp& dest)
206 {
207     typedef typename Eigen::internal::traits<Sp>::Scalar data_t;
208     mycopy_n(makeMatrixIterator<data_t>(src, makeNonZerosIterator(src)), nonZeros(src)
209              , makeMatrixIterator<data_t>(dest, makeTranslatedIterator(makeNonZerosIterator(src), Coords2D(r, c))));
210 }
211
212 template<typename Scalar1, typename Scalar2>
213 void doAppend(Eigen::SparseMatrix<Scalar1, Eigen::RowMajor> SPARSE_CONST& src, int r, int c, Eigen::SparseMatrix<Scalar2, Eigen::RowMajor>& dest)
214 {
215     typedef typename Eigen::SparseMatrix<Scalar1, Eigen::RowMajor>::InnerIterator srcIt_t;
216     for (std::size_t k = 0; k != src.outerSize(); ++k)
217     {
218         for (srcIt_t it(src, (int)k); it; ++it)
219         {
220             dest.insert( it.row() + r, it.col() + c) =  it.value();
221         }
222     }
223 }
224 /*
225 Sp is an Eigen::SparseMatrix
226 */
227 template<typename Sp, typename M>
228 void cwiseInPlaceProduct(Sp& sp, M SPARSE_CONST& m)
229 {
230     // should be a transform_n() over makeNonZerosIterator(src)
231     for (std::size_t k = 0; k != sp.outerSize(); ++k)
232     {
233         for (typename Sp::InnerIterator it(sp, k); it; ++it)
234         {
235             it.valueRef() *= get<typename Eigen::internal::traits<Sp>::Scalar >(m, it.row(), it.col());
236         }
237     }
238
239 }
240 }
241 namespace types
242 {
243
244 template<typename T, typename Arg>
245 T* create_new(Arg const& a)
246 {
247     return 0;
248 }
249
250 template<>
251 Double* create_new(double const& d)
252 {
253     Double* res(new Double(1, 1, false));
254     res->set(0, 0, d);
255     return res;
256 }
257
258 template<>
259 Double* create_new(std::complex<double>const& c)
260 {
261     Double* res(new Double(1, 1, true));
262     res->set(0, 0, c.real());
263     res->setImg(0, 0, c.imag());
264     return res;
265 }
266
267 template<>
268 Double* create_new(Sparse const& s)
269 {
270     Sparse& cs(const_cast<Sparse&>(s)); // inherited member functions are not const-correct
271     Double* res(new Double(cs.getRows(), cs.getCols(), cs.isComplex()));
272     const_cast<Sparse&>(s).fill(*res);
273     return res;
274 }
275
276
277 Sparse::~Sparse()
278 {
279     delete matrixReal;
280     delete matrixCplx;
281 #ifndef NDEBUG
282     Inspector::removeItem(this);
283 #endif
284 }
285
286 Sparse::Sparse(Sparse const& src)
287     : matrixReal(src.matrixReal ? new RealSparse_t(*src.matrixReal) : 0)
288     , matrixCplx(src.matrixCplx ? new CplxSparse_t(*src.matrixCplx) : 0)
289
290 {
291     m_iRows = const_cast<Sparse*>(&src)->getRows();
292     m_iCols = const_cast<Sparse*>(&src)->getCols();
293     m_iSize = m_iRows * m_iCols;
294     m_iDims = 2;
295     m_piDims[0] = m_iRows;
296     m_piDims[1] = m_iCols;
297 #ifndef NDEBUG
298     Inspector::addItem(this);
299 #endif
300 }
301
302 Sparse::Sparse(int _iRows, int _iCols, bool cplx)
303     : matrixReal(cplx ? 0 : new RealSparse_t(_iRows, _iCols))
304     , matrixCplx(cplx ? new CplxSparse_t(_iRows, _iCols) : 0)
305 {
306     m_iRows = _iRows;
307     m_iCols = _iCols;
308     m_iSize = _iRows * _iCols;
309     m_iDims = 2;
310     m_piDims[0] = _iRows;
311     m_piDims[1] = _iCols;
312 #ifndef NDEBUG
313     Inspector::addItem(this);
314 #endif
315 }
316
317 Sparse::Sparse(Double SPARSE_CONST& src)
318 {
319     //compute idx
320     int size = src.getSize();
321     int row = src.getRows();
322     Double* idx = new Double(src.getSize(), 2);
323     double* p = idx->get();
324     for (int i = 0; i < size; ++i)
325     {
326         p[i]        = (double)(i % row) + 1;
327         p[i + size] = (double)(i / row) + 1;
328     }
329     create2(src.getRows(), src.getCols(), src, *idx);
330     idx->killMe();
331 #ifndef NDEBUG
332     Inspector::addItem(this);
333 #endif
334 }
335
336 Sparse::Sparse(Double SPARSE_CONST& src, Double SPARSE_CONST& idx)
337 {
338     int idxrow = idx.getRows();
339     int rows = static_cast<int>(*std::max_element(idx.get(), idx.get() + idxrow));
340     int cols = static_cast<int>(*std::max_element(idx.get() + idxrow, idx.get() + idxrow * 2));
341
342     create2(rows, cols, src, idx);
343 #ifndef NDEBUG
344     Inspector::removeItem(this);
345 #endif
346 }
347
348 Sparse::Sparse(Double SPARSE_CONST& src, Double SPARSE_CONST& idx, Double SPARSE_CONST& dims)
349 {
350     create2(static_cast<int>(dims.get(0)), static_cast<int>(dims.get(1)), src, idx);
351 #ifndef NDEBUG
352     Inspector::addItem(this);
353 #endif
354 }
355
356 Sparse::Sparse(RealSparse_t* realSp, CplxSparse_t* cplxSp):  matrixReal(realSp), matrixCplx(cplxSp)
357 {
358     if (realSp)
359     {
360         m_iCols = realSp->cols();
361         m_iRows = realSp->rows();
362     }
363     else
364     {
365         m_iCols = cplxSp->cols();
366         m_iRows = cplxSp->rows();
367     }
368     m_iSize = m_iCols * m_iRows;
369     m_iDims = 2;
370     m_piDims[0] = m_iRows;
371     m_piDims[1] = m_iCols;
372
373     finalize();
374 #ifndef NDEBUG
375     Inspector::addItem(this);
376 #endif
377 }
378
379 Sparse::Sparse(Double SPARSE_CONST& xadj, Double SPARSE_CONST& adjncy, Double SPARSE_CONST& src, std::size_t r, std::size_t c)
380 {
381     Adjacency a(xadj.get(), adjncy.get());
382     create(static_cast<int>(r), static_cast<int>(c), src, makeIteratorFromVar(a), src.getSize());
383 #ifndef NDEBUG
384     Inspector::addItem(this);
385 #endif
386 }
387
388 Sparse::Sparse(int rows, int cols, int nonzeros, int* inner, int* outer, double* real, double* img)
389 {
390     int* out = nullptr;
391     int* in = nullptr;
392
393     if (img)
394     {
395         matrixCplx = new CplxSparse_t(rows, cols);
396         matrixCplx->reserve((int)nonzeros);
397         out = matrixCplx->outerIndexPtr();
398         in = matrixCplx->innerIndexPtr();
399         matrixReal = nullptr;
400     }
401     else
402     {
403         matrixReal = new RealSparse_t(rows, cols);
404         matrixReal->reserve((int)nonzeros);
405         out = matrixReal->outerIndexPtr();
406         in = matrixReal->innerIndexPtr();
407         matrixCplx = nullptr;
408     }
409
410     //update outerIndexPtr
411     memcpy(out, outer, sizeof(int) * (rows + 1));
412     //update innerIndexPtr
413     memcpy(in, inner, sizeof(int) * nonzeros);
414
415     if (img)
416     {
417         std::complex<double>* data = matrixCplx->valuePtr();
418         for (int i = 0; i < nonzeros; ++i)
419         {
420             data[i] = std::complex<double>(real[i], img[i]);
421         }
422     }
423     else
424     {
425         double* data = matrixReal->valuePtr();
426         for (int i = 0; i < nonzeros; ++i)
427         {
428             data[i] = real[i];
429         }
430
431     }
432
433     m_iCols = cols;
434     m_iRows = rows;
435     m_iSize = cols * rows;
436     m_iDims = 2;
437     m_piDims[0] = m_iRows;
438     m_piDims[1] = m_iCols;
439
440     matrixCplx ? matrixCplx->resizeNonZeros(nonzeros) : matrixReal->resizeNonZeros(nonzeros);
441     //finalize();
442 }
443
444 template<typename DestIter>
445 void Sparse::create(int rows, int cols, Double SPARSE_CONST& src, DestIter o, std::size_t n)
446 {
447     m_iCols = cols;
448     m_iRows = rows;
449     m_iSize = cols * rows;
450     m_iDims = 2;
451     m_piDims[0] = m_iRows;
452     m_piDims[1] = m_iCols;
453
454     if (src.isComplex())
455     {
456         matrixReal = 0;
457         matrixCplx = new CplxSparse_t(rows, cols);
458         matrixCplx->reserve((int)n);
459         mycopy_n(makeMatrixIterator<std::complex<double> >(src, RowWiseFullIterator(src.getRows(), src.getCols())), n, makeMatrixIterator<std::complex<double> >(*matrixCplx, o));
460     }
461     else
462     {
463         matrixReal = new RealSparse_t(rows, cols);
464         matrixReal->reserve((int)n);
465         matrixCplx = 0;
466         mycopy_n(makeMatrixIterator<double >(src, RowWiseFullIterator(src.getRows(), src.getCols())), n
467                  , makeMatrixIterator<double>(*matrixReal, o));
468     }
469     finalize();
470 }
471
472 void Sparse::create2(int rows, int cols, Double SPARSE_CONST& src, Double SPARSE_CONST& idx)
473 {
474     int nnz = src.getSize();
475     double* i = idx.get();
476     double* j = i + idx.getRows();
477     double* valR = src.get();
478
479     if (src.isComplex())
480     {
481         matrixReal = 0;
482
483         typedef Eigen::Triplet<std::complex<double> > T;
484         std::vector<T> tripletList;
485         tripletList.reserve((int)nnz);
486
487         double* valI = src.getImg();
488
489         for (int k = 0; k < nnz; ++k)
490         {
491             tripletList.push_back(T(static_cast<int>(i[k]) - 1, static_cast<int>(j[k]) - 1, std::complex<double>(valR[k], valI[k])));
492         }
493
494         matrixCplx = new CplxSparse_t(rows, cols);
495         matrixCplx->setFromTriplets(tripletList.begin(), tripletList.end());
496         m_iRows = matrixCplx->rows();
497         m_iCols = matrixCplx->cols();
498     }
499     else
500     {
501         matrixCplx = 0;
502
503         typedef Eigen::Triplet<double> T;
504         std::vector<T> tripletList;
505         tripletList.reserve((int)nnz);
506
507         for (int k = 0; k < nnz; ++k)
508         {
509             tripletList.push_back(T(static_cast<int>(i[k]) - 1, static_cast<int>(j[k]) - 1, valR[k]));
510         }
511
512         matrixReal = new RealSparse_t(rows, cols);
513         matrixReal->setFromTriplets(tripletList.begin(), tripletList.end());
514
515         m_iRows = matrixReal->rows();
516         m_iCols = matrixReal->cols();
517     }
518
519     m_iSize = m_iCols * m_iRows;
520     m_iDims = 2;
521     m_piDims[0] = m_iRows;
522     m_piDims[1] = m_iCols;
523     finalize();
524 }
525
526 void Sparse::fill(Double& dest, int r, int c) SPARSE_CONST
527 {
528     Sparse & cthis(const_cast<Sparse&>(*this));
529     if (isComplex())
530     {
531         mycopy_n(makeMatrixIterator<std::complex<double> >(*matrixCplx, RowWiseFullIterator(cthis.getRows(), cthis.getCols())), cthis.getSize()
532         , makeMatrixIterator<std::complex<double> >(dest, RowWiseFullIterator(dest.getRows(), dest.getCols(), r, c)));
533     }
534     else
535     {
536         mycopy_n( makeMatrixIterator<double>(*matrixReal,  RowWiseFullIterator(cthis.getRows(), cthis.getCols())), cthis.getSize()
537         , makeMatrixIterator<double >(dest, RowWiseFullIterator(dest.getRows(), dest.getCols(), r, c)));
538     }
539 }
540
541 bool Sparse::set(int _iRows, int _iCols, std::complex<double> v, bool _bFinalize)
542 {
543     if (_iRows >= getRows() || _iCols >= getCols())
544     {
545         return false;
546     }
547
548     if (matrixReal)
549     {
550         matrixReal->coeffRef(_iRows, _iCols) = v.real();
551     }
552     else
553     {
554         matrixCplx->coeffRef(_iRows, _iCols) = v;
555     }
556
557     if (_bFinalize)
558     {
559         finalize();
560     }
561     return true;
562 }
563
564 bool Sparse::set(int _iRows, int _iCols, double _dblReal, bool _bFinalize)
565 {
566     if (_iRows >= getRows() || _iCols >= getCols())
567     {
568         return false;
569     }
570
571     if (matrixReal)
572     {
573         matrixReal->coeffRef(_iRows, _iCols) = _dblReal;
574     }
575     else
576     {
577         matrixCplx->coeffRef(_iRows, _iCols) = std::complex<double>(_dblReal, 0);
578     }
579
580
581     if (_bFinalize)
582     {
583         finalize();
584     }
585     return true;
586 }
587
588 void Sparse::finalize()
589 {
590     if (isComplex())
591     {
592         matrixCplx->prune(&keepForSparse<std::complex<double> >);
593         matrixCplx->finalize();
594     }
595     else
596     {
597         matrixReal->prune(&keepForSparse<double>);
598         matrixReal->finalize();
599     }
600
601 }
602
603 bool Sparse::neg(InternalType *& out)
604 {
605     SparseBool * _out = new SparseBool(getRows(), getCols());
606     types::neg(getRows(), getCols(), matrixReal, _out->matrixBool);
607     out = _out;
608
609     return true;
610 }
611
612
613 bool Sparse::isComplex() const
614 {
615     return static_cast<bool>(matrixCplx != NULL);
616 }
617
618 // TODO: should have both a bounds checking and a non-checking interface to elt access
619 double* Sparse::get()
620 {
621     if (isComplex() == false)
622     {
623         return matrixReal->valuePtr();
624     }
625
626     return nullptr;
627 }
628
629 double  Sparse::get(int _iRows, int _iCols) const
630 {
631     return getReal(_iRows, _iCols);
632 }
633
634 double Sparse::getReal(int _iRows, int _iCols) const
635 {
636     double res = 0;
637     if (matrixReal)
638     {
639         res = matrixReal->coeff(_iRows, _iCols);
640     }
641     else
642     {
643         res = matrixCplx->coeff(_iRows, _iCols).real();
644     }
645     return res;
646 }
647
648 std::complex<double>* Sparse::getImg()
649 {
650     if (isComplex())
651     {
652         return matrixCplx->valuePtr();
653     }
654
655     return nullptr;
656 }
657
658 std::complex<double> Sparse::getImg(int _iRows, int _iCols) const
659 {
660     std::complex<double> res;
661     if (matrixCplx)
662     {
663         res = matrixCplx->coeff(_iRows, _iCols);
664     }
665     else
666     {
667         res = std::complex<double>(matrixReal->coeff(_iRows, _iCols), 0.);
668     }
669
670     return res;
671 }
672
673 void Sparse::whoAmI() SPARSE_CONST
674 {
675     std::cout << "types::Sparse";
676 }
677
678 Sparse* Sparse::clone(void) const
679 {
680     return new Sparse(*this);
681 }
682
683 bool Sparse::zero_set()
684 {
685     if (matrixReal)
686     {
687         matrixReal->setZero();
688     }
689     else
690     {
691         matrixCplx->setZero();
692     }
693
694     return true;
695 }
696
697 // TODO: handle precision and line length
698 bool Sparse::toString(std::wostringstream& ostr) const
699 {
700     int iPrecision = ConfigVariable::getFormatSize();
701     std::wstring res;
702     if (matrixReal)
703     {
704         res = ::toString(*matrixReal, iPrecision);
705     }
706     else
707     {
708         res = ::toString(*matrixCplx, iPrecision);
709     }
710
711     ostr << res;
712     return true;
713 }
714
715 bool Sparse::resize(int _iNewRows, int _iNewCols)
716 {
717     if (_iNewRows <= getRows() && _iNewCols <= getCols())
718     {
719         //nothing to do: hence we do NOT fail
720         return true;
721     }
722
723     bool res = false;
724     try
725     {
726         if (matrixReal)
727         {
728             //item count
729             size_t iNonZeros = nonZeros();
730             RealSparse_t *newReal = new RealSparse_t(_iNewRows, _iNewCols);
731             newReal->reserve((int)iNonZeros);
732
733
734             //coords
735             int* pRows = new int[iNonZeros * 2];
736             outputRowCol(pRows);
737             int* pCols = pRows + iNonZeros;
738
739             //values
740             double* pNonZeroR = new double[iNonZeros];
741             double* pNonZeroI = new double[iNonZeros];
742             outputValues(pNonZeroR, pNonZeroI);
743
744             typedef Eigen::Triplet<double> triplet;
745             std::vector<triplet> tripletList;
746
747             for (size_t i = 0 ; i < iNonZeros ; i++)
748             {
749                 tripletList.push_back(triplet((int)pRows[i] - 1, (int)pCols[i] - 1, pNonZeroR[i]));
750             }
751
752             newReal->setFromTriplets(tripletList.begin(), tripletList.end());
753
754             delete matrixReal;
755             matrixReal = newReal;
756             delete[] pRows;
757             delete[] pNonZeroR;
758             delete[] pNonZeroI;
759         }
760         else
761         {
762             //item count
763             size_t iNonZeros = nonZeros();
764             CplxSparse_t *newCplx = new CplxSparse_t(_iNewRows, _iNewCols);
765             newCplx->reserve((int)iNonZeros);
766
767             //coords
768             int* pRows = new int[iNonZeros * 2];
769             outputRowCol(pRows);
770             int* pCols = pRows + iNonZeros;
771
772             //values
773             double* pNonZeroR = new double[iNonZeros];
774             double* pNonZeroI = new double[iNonZeros];
775             outputValues(pNonZeroR, pNonZeroI);
776
777             typedef Eigen::Triplet<std::complex<double> > triplet;
778             std::vector<triplet> tripletList;
779
780             for (size_t i = 0 ; i < iNonZeros ; i++)
781             {
782                 tripletList.push_back(triplet((int)pRows[i] - 1, (int)pCols[i] - 1, std::complex<double>(pNonZeroR[i], pNonZeroI[i])));
783             }
784
785             newCplx->setFromTriplets(tripletList.begin(), tripletList.end());
786
787
788             delete matrixCplx;
789             matrixCplx = newCplx;
790             delete[] pRows;
791             delete[] pNonZeroR;
792             delete[] pNonZeroI;
793         }
794
795         m_iRows = _iNewRows;
796         m_iCols = _iNewCols;
797         m_iSize = _iNewRows * _iNewCols;
798         m_piDims[0] = m_iRows;
799         m_piDims[1] = m_iCols;
800
801         res = true;
802     }
803     catch (...)
804     {
805         res = false;
806     }
807     return res;
808 }
809 // TODO decide if a complex matrix with 0 imag can be == to a real matrix
810 // not true for dense (cf double.cpp)
811 bool Sparse::operator==(const InternalType& it) SPARSE_CONST
812 {
813     Sparse* otherSparse = const_cast<Sparse*>(dynamic_cast<Sparse const*>(&it));/* types::GenericType is not const-correct :( */
814     Sparse & cthis (const_cast<Sparse&>(*this));
815
816     if (otherSparse == NULL)
817     {
818         return false;
819     }
820
821     if (otherSparse->getRows() != cthis.getRows())
822     {
823         return false;
824     }
825
826     if (otherSparse->getCols() != cthis.getCols())
827     {
828         return false;
829     }
830
831     if (otherSparse->isComplex() != isComplex())
832     {
833         return false;
834     }
835
836     if (isComplex())
837     {
838         return equal(*matrixCplx, *otherSparse->matrixCplx);
839     }
840     else
841     {
842         return equal(*matrixReal, *otherSparse->matrixReal);
843     }
844 }
845
846 bool Sparse::one_set()
847 {
848     if (isComplex())
849     {
850         return setNonZero(*matrixCplx);
851     }
852     else
853     {
854         return setNonZero(*matrixReal);
855     }
856 }
857
858 void Sparse::toComplex()
859 {
860     if (!isComplex())
861     {
862         try
863         {
864             matrixCplx = new CplxSparse_t(matrixReal->cast<std::complex<double> >());
865             delete matrixReal;
866             matrixReal = NULL;
867         }
868         catch (...)
869         {
870             delete matrixCplx;
871             matrixCplx = NULL;
872             throw;
873         }
874     }
875 }
876
877 InternalType* Sparse::insertNew(typed_list* _pArgs, InternalType* _pSource)
878 {
879     typed_list pArg;
880     InternalType *pOut  = NULL;
881     Sparse* pSource = _pSource->getAs<Sparse>();
882
883     int iDims           = (int)_pArgs->size();
884     int* piMaxDim       = new int[iDims];
885     int* piCountDim     = new int[iDims];
886     bool bComplex       = pSource->isComplex();
887     bool bUndefine      = false;
888
889     //evaluate each argument and replace by appropriate value and compute the count of combinations
890     int iSeqCount = checkIndexesArguments(NULL, _pArgs, &pArg, piMaxDim, piCountDim);
891     if (iSeqCount == 0)
892     {
893         //free pArg content
894         cleanIndexesArguments(_pArgs, &pArg);
895         return createEmptyDouble();
896     }
897
898     if (iSeqCount < 0)
899     {
900         iSeqCount = -iSeqCount;
901         bUndefine = true;
902     }
903
904     if (bUndefine)
905     {
906         //manage : and $ in creation by insertion
907         int iSource         = 0;
908         int *piSourceDims   = pSource->getDimsArray();
909
910         for (int i = 0 ; i < iDims ; i++)
911         {
912             if (pArg[i] == NULL)
913             {
914                 //undefine value
915                 if (pSource->isScalar())
916                 {
917                     piMaxDim[i]     = 1;
918                     piCountDim[i]   = 1;
919                 }
920                 else
921                 {
922                     piMaxDim[i]     = piSourceDims[iSource];
923                     piCountDim[i]   = piSourceDims[iSource];
924                 }
925                 iSource++;
926                 //replace pArg value by the new one
927                 pArg[i] = createDoubleVector(piMaxDim[i]);
928             }
929             //else
930             //{
931             //    piMaxDim[i] = piCountDim[i];
932             //}
933         }
934     }
935
936     //remove last dimension at size 1
937     //remove last dimension if are == 1
938     for (int i = (iDims - 1) ; i >= 2 ; i--)
939     {
940         if (piMaxDim[i] == 1)
941         {
942             iDims--;
943             pArg.pop_back();
944         }
945         else
946         {
947             break;
948         }
949     }
950
951     if (checkArgValidity(pArg) == false)
952     {
953         //free pArg content
954         cleanIndexesArguments(_pArgs, &pArg);
955         //contain bad index, like <= 0, ...
956         return NULL;
957     }
958
959     if (iDims == 1)
960     {
961         if (pSource->getCols() == 1)
962         {
963             pOut = new Sparse(piCountDim[0], 1, bComplex);
964         }
965         else
966         {
967             //rows == 1
968             pOut = new Sparse(1, piCountDim[0], bComplex);
969         }
970     }
971     else
972     {
973         pOut = new Sparse(piMaxDim[0], piMaxDim[1], bComplex);
974         //pOut = pSource->createEmpty(iDims, piMaxDim, bComplex);
975     }
976
977     //fill with null item
978     Sparse* pSpOut = pOut->getAs<Sparse>();
979
980     //insert values in new matrix
981     InternalType* pOut2 = pSpOut->insert(&pArg, pSource);
982     if (pOut != pOut2)
983     {
984         delete pOut;
985     }
986
987     //free pArg content
988     cleanIndexesArguments(_pArgs, &pArg);
989
990     return pOut2;
991 }
992
993 Sparse* Sparse::insert(typed_list* _pArgs, InternalType* _pSource)
994 {
995     bool bNeedToResize  = false;
996     int iDims           = (int)_pArgs->size();
997     if (iDims > 2)
998     {
999         //sparse are only in 2 dims
1000         return NULL;
1001     }
1002
1003     typed_list pArg;
1004
1005     int piMaxDim[2];
1006     int piCountDim[2];
1007
1008     //on case of resize
1009     int iNewRows    = 0;
1010     int iNewCols    = 0;
1011     Double* pSource = _pSource->getAs<Double>();
1012
1013     //evaluate each argument and replace by appropriate value and compute the count of combinations
1014     int iSeqCount = checkIndexesArguments(this, _pArgs, &pArg, piMaxDim, piCountDim);
1015     if (iSeqCount == 0)
1016     {
1017         //free pArg content
1018         cleanIndexesArguments(_pArgs, &pArg);
1019         return this;
1020     }
1021
1022     if (iDims < 2)
1023     {
1024         //see as vector
1025         if (getRows() == 1 || getCols() == 1)
1026         {
1027             //vector or scalar
1028             if (getSize() < piMaxDim[0])
1029             {
1030                 bNeedToResize = true;
1031
1032                 //need to enlarge sparse dimensions
1033                 if (getCols() == 1 || getSize() == 0)
1034                 {
1035                     //column vector
1036                     iNewRows    = piMaxDim[0];
1037                     iNewCols    = 1;
1038                 }
1039                 else if (getRows() == 1)
1040                 {
1041                     //row vector
1042                     iNewRows    = 1;
1043                     iNewCols    = piMaxDim[0];
1044                 }
1045             }
1046         }
1047         else if (getSize() < piMaxDim[0])
1048         {
1049             //free pArg content
1050             cleanIndexesArguments(_pArgs, &pArg);
1051             //out of range
1052             return NULL;
1053         }
1054     }
1055     else
1056     {
1057         if (piMaxDim[0] > getRows() || piMaxDim[1] > getCols())
1058         {
1059             bNeedToResize = true;
1060             iNewRows = std::max(getRows(), piMaxDim[0]);
1061             iNewCols = std::max(getCols(), piMaxDim[1]);
1062         }
1063     }
1064
1065     //check number of insertion
1066     if (pSource->isScalar() == false && pSource->getSize() != iSeqCount)
1067     {
1068         //free pArg content
1069         cleanIndexesArguments(_pArgs, &pArg);
1070         return NULL;
1071     }
1072
1073     //now you are sure to be able to insert values
1074     if (bNeedToResize)
1075     {
1076         if (resize(iNewRows, iNewCols) == false)
1077         {
1078             //free pArg content
1079             cleanIndexesArguments(_pArgs, &pArg);
1080             return NULL;
1081         }
1082     }
1083
1084     //update complexity
1085     if (pSource->isComplex() && isComplex() == false)
1086     {
1087         toComplex();
1088     }
1089
1090
1091     if (iDims == 1)
1092     {
1093         double* pIdx = pArg[0]->getAs<Double>()->get();
1094         for (int i = 0 ; i < iSeqCount ; i++)
1095         {
1096             int iRow = static_cast<int>(pIdx[i] - 1) % getRows();
1097             int iCol = static_cast<int>(pIdx[i] - 1) / getRows();
1098             if (pSource->isScalar())
1099             {
1100                 if (pSource->isComplex())
1101                 {
1102                     set(iRow, iCol, std::complex<double>(pSource->get(0), pSource->getImg(0)), false);
1103                 }
1104                 else
1105                 {
1106                     set(iRow, iCol, pSource->get(0), false);
1107                 }
1108             }
1109             else
1110             {
1111                 if (pSource->isComplex())
1112                 {
1113                     set(iRow, iCol, std::complex<double>(pSource->get(i), pSource->getImg(i)), false);
1114                 }
1115                 else
1116                 {
1117                     set(iRow, iCol, pSource->get(i), false);
1118                 }
1119             }
1120         }
1121     }
1122     else
1123     {
1124         double* pIdxRow = pArg[0]->getAs<Double>()->get();
1125         int iRowSize    = pArg[0]->getAs<Double>()->getSize();
1126         double* pIdxCol = pArg[1]->getAs<Double>()->get();
1127
1128         for (int i = 0 ; i < iSeqCount ; i++)
1129         {
1130             if (pSource->isScalar())
1131             {
1132                 if (pSource->isComplex())
1133                 {
1134                     set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, std::complex<double>(pSource->get(0), pSource->getImg(0)), false);
1135                 }
1136                 else
1137                 {
1138                     set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, pSource->get(0), false);
1139                 }
1140             }
1141             else
1142             {
1143                 int iRowOrig = i % pSource->getRows();
1144                 int iColOrig = i / pSource->getRows();
1145
1146                 if (pSource->isComplex())
1147                 {
1148                     set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, std::complex<double>(pSource->get(iRowOrig, iColOrig), pSource->getImg(iRowOrig, iColOrig)), false);
1149                 }
1150                 else
1151                 {
1152                     set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, pSource->get(iRowOrig, iColOrig), false);
1153                 }
1154             }
1155         }
1156     }
1157
1158     finalize();
1159
1160     //free pArg content
1161     cleanIndexesArguments(_pArgs, &pArg);
1162
1163     return this;
1164 }
1165
1166 Sparse* Sparse::insert(typed_list* _pArgs, Sparse* _pSource)
1167 {
1168     bool bNeedToResize  = false;
1169     int iDims           = (int)_pArgs->size();
1170     if (iDims > 2)
1171     {
1172         //sparse are only in 2 dims
1173         return NULL;
1174     }
1175
1176     typed_list pArg;
1177
1178     int piMaxDim[2];
1179     int piCountDim[2];
1180
1181     //on case of resize
1182     int iNewRows    = 0;
1183     int iNewCols    = 0;
1184
1185     //evaluate each argument and replace by appropriate value and compute the count of combinations
1186     int iSeqCount = checkIndexesArguments(this, _pArgs, &pArg, piMaxDim, piCountDim);
1187     if (iSeqCount == 0)
1188     {
1189         //free pArg content
1190         cleanIndexesArguments(_pArgs, &pArg);
1191         return this;
1192     }
1193
1194     if (iDims < 2)
1195     {
1196         //see as vector
1197         if (getRows() == 1 || getCols() == 1)
1198         {
1199             //vector or scalar
1200             bNeedToResize = true;
1201             if (getSize() < piMaxDim[0])
1202             {
1203                 //need to enlarge sparse dimensions
1204                 if (getCols() == 1 || getSize() == 0)
1205                 {
1206                     //column vector
1207                     iNewRows    = piMaxDim[0];
1208                     iNewCols    = 1;
1209                 }
1210                 else if (getRows() == 1)
1211                 {
1212                     //row vector
1213                     iNewRows    = 1;
1214                     iNewCols    = piMaxDim[0];
1215                 }
1216             }
1217         }
1218         else if (getSize() < piMaxDim[0])
1219         {
1220             //free pArg content
1221             cleanIndexesArguments(_pArgs, &pArg);
1222             //out of range
1223             return NULL;
1224         }
1225     }
1226     else
1227     {
1228         if (piMaxDim[0] > getRows() || piMaxDim[1] > getCols())
1229         {
1230             bNeedToResize = true;
1231             iNewRows = std::max(getRows(), piMaxDim[0]);
1232             iNewCols = std::max(getCols(), piMaxDim[1]);
1233         }
1234     }
1235
1236     //check number of insertion
1237     if (_pSource->isScalar() == false && _pSource->getSize() != iSeqCount)
1238     {
1239         //free pArg content
1240         cleanIndexesArguments(_pArgs, &pArg);
1241         return NULL;
1242     }
1243
1244     //now you are sure to be able to insert values
1245     if (bNeedToResize)
1246     {
1247         if (resize(iNewRows, iNewCols) == false)
1248         {
1249             //free pArg content
1250             cleanIndexesArguments(_pArgs, &pArg);
1251             return NULL;
1252         }
1253     }
1254
1255     //update complexity
1256     if (_pSource->isComplex() && isComplex() == false)
1257     {
1258         toComplex();
1259     }
1260
1261     if (iDims == 1)
1262     {
1263         double* pIdx = pArg[0]->getAs<Double>()->get();
1264         for (int i = 0 ; i < iSeqCount ; i++)
1265         {
1266             int iRow = static_cast<int>(pIdx[i] - 1) % getRows();
1267             int iCol = static_cast<int>(pIdx[i] - 1) / getRows();
1268
1269             if (_pSource->isScalar())
1270             {
1271                 if (_pSource->isComplex())
1272                 {
1273                     set(iRow, iCol, _pSource->getImg(0, 0), false);
1274                 }
1275                 else
1276                 {
1277                     set(iRow, iCol, _pSource->get(0, 0), false);
1278                 }
1279             }
1280             else
1281             {
1282                 int iRowOrig = i % _pSource->getRows();
1283                 int iColOrig = i / _pSource->getRows();
1284                 if (_pSource->isComplex())
1285                 {
1286                     set(iRow, iCol, _pSource->getImg(iRowOrig, iColOrig), false);
1287                 }
1288                 else
1289                 {
1290                     set(iRow, iCol, _pSource->get(iRowOrig, iColOrig), false);
1291                 }
1292             }
1293         }
1294     }
1295     else
1296     {
1297         double* pIdxRow = pArg[0]->getAs<Double>()->get();
1298         int iRowSize    = pArg[0]->getAs<Double>()->getSize();
1299         double* pIdxCol = pArg[1]->getAs<Double>()->get();
1300
1301         for (int i = 0 ; i < iSeqCount ; i++)
1302         {
1303             if (_pSource->isScalar())
1304             {
1305                 if (_pSource->isComplex())
1306                 {
1307                     set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, _pSource->getImg(0, 0), false);
1308                 }
1309                 else
1310                 {
1311                     set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, _pSource->get(0, 0), false);
1312                 }
1313             }
1314             else
1315             {
1316                 int iRowOrig = i % _pSource->getRows();
1317                 int iColOrig = i / _pSource->getRows();
1318                 if (_pSource->isComplex())
1319                 {
1320                     set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, _pSource->getImg(iRowOrig, iColOrig), false);
1321                 }
1322                 else
1323                 {
1324                     set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, _pSource->get(iRowOrig, iColOrig), false);
1325                 }
1326             }
1327         }
1328     }
1329
1330     finalize();
1331
1332     //free pArg content
1333     cleanIndexesArguments(_pArgs, &pArg);
1334
1335     return this;
1336 }
1337
1338 Sparse* Sparse::remove(typed_list* _pArgs)
1339 {
1340     Sparse* pOut = NULL;
1341     int iDims = (int)_pArgs->size();
1342     if (iDims > 2)
1343     {
1344         //sparse are only in 2 dims
1345         return NULL;
1346     }
1347
1348     typed_list pArg;
1349
1350     int piMaxDim[2];
1351     int piCountDim[2];
1352
1353     //evaluate each argument and replace by appropriate value and compute the count of combinations
1354     int iSeqCount = checkIndexesArguments(this, _pArgs, &pArg, piMaxDim, piCountDim);
1355     if (iSeqCount == 0)
1356     {
1357         //free pArg content
1358         cleanIndexesArguments(_pArgs, &pArg);
1359         return this;
1360     }
1361
1362     bool* pbFull = new bool[iDims];
1363     //coord must represent all values on a dimension
1364     for (int i = 0 ; i < iDims ; i++)
1365     {
1366         pbFull[i]       = false;
1367         int iDimToCheck = getVarMaxDim(i, iDims);
1368         int iIndexSize  = pArg[i]->getAs<GenericType>()->getSize();
1369
1370         //we can have index more than once
1371         if (iIndexSize >= iDimToCheck)
1372         {
1373             //size is good, now check datas
1374             double* pIndexes = getDoubleArrayFromDouble(pArg[i]);
1375             for (int j = 0 ; j < iDimToCheck ; j++)
1376             {
1377                 bool bFind = false;
1378                 for (int k = 0 ; k < iIndexSize ; k++)
1379                 {
1380                     if ((int)pIndexes[k] == j + 1)
1381                     {
1382                         bFind = true;
1383                         break;
1384                     }
1385                 }
1386                 pbFull[i]  = bFind;
1387             }
1388         }
1389     }
1390
1391     //only one dims can be not full/entire
1392     bool bNotEntire = false;
1393     int iNotEntire  = 0;
1394     bool bTooMuchNotEntire = false;
1395     for (int i = 0 ; i < iDims ; i++)
1396     {
1397         if (pbFull[i] == false)
1398         {
1399             if (bNotEntire == false)
1400             {
1401                 bNotEntire = true;
1402                 iNotEntire = i;
1403             }
1404             else
1405             {
1406                 bTooMuchNotEntire = true;
1407                 break;
1408             }
1409         }
1410     }
1411
1412     if (bTooMuchNotEntire == true)
1413     {
1414         //free pArg content
1415         cleanIndexesArguments(_pArgs, &pArg);
1416         return NULL;
1417     }
1418
1419     delete[] pbFull;
1420
1421     //find index to keep
1422     int iNotEntireSize          = pArg[iNotEntire]->getAs<GenericType>()->getSize();
1423     double* piNotEntireIndex    = getDoubleArrayFromDouble(pArg[iNotEntire]);
1424     int iKeepSize               = getVarMaxDim(iNotEntire, iDims);
1425     bool* pbKeep                = new bool[iKeepSize];
1426
1427     //fill pbKeep with true value
1428     for (int i = 0 ; i < iKeepSize ; i++)
1429     {
1430         pbKeep[i] = true;
1431     }
1432
1433     for (int i = 0 ; i < iNotEntireSize ; i++)
1434     {
1435         int idx = (int)piNotEntireIndex[i] - 1;
1436
1437         //don't care of value out of bounds
1438         if (idx < iKeepSize)
1439         {
1440             pbKeep[idx] = false;
1441         }
1442     }
1443
1444     int iNewDimSize = 0;
1445     for (int i = 0 ; i < iKeepSize ; i++)
1446     {
1447         if (pbKeep[i] == true)
1448         {
1449             iNewDimSize++;
1450         }
1451     }
1452     delete[] pbKeep;
1453
1454     int* piNewDims = new int[iDims];
1455     for (int i = 0 ; i < iDims ; i++)
1456     {
1457         if (i == iNotEntire)
1458         {
1459             piNewDims[i] = iNewDimSize;
1460         }
1461         else
1462         {
1463             piNewDims[i] = getVarMaxDim(i, iDims);
1464         }
1465     }
1466
1467     //remove last dimension if are == 1
1468     int iOrigDims = iDims;
1469     for (int i = (iDims - 1) ; i >= 2 ; i--)
1470     {
1471         if (piNewDims[i] == 1)
1472         {
1473             iDims--;
1474         }
1475         else
1476         {
1477             break;
1478         }
1479     }
1480
1481     if (iDims == 1)
1482     {
1483         if (iNewDimSize == 0)
1484         {
1485             delete[] piNewDims;
1486             //free pArg content
1487             cleanIndexesArguments(_pArgs, &pArg);
1488             return new Sparse(0, 0);
1489         }
1490         else
1491         {
1492             //two cases, depends of original matrix/vector
1493             if ((*_pArgs)[0]->isColon() == false && m_iDims == 2 && m_piDims[0] == 1 && m_piDims[1] != 1)
1494             {
1495                 //special case for row vector
1496                 pOut = new Sparse(1, iNewDimSize, isComplex());
1497                 //in this case we have to care of 2nd dimension
1498                 //iNotEntire = 1;
1499             }
1500             else
1501             {
1502                 pOut = new Sparse(iNewDimSize, 1, isComplex());
1503             }
1504         }
1505     }
1506     else
1507     {
1508         pOut = new Sparse(piNewDims[0], piNewDims[0], isComplex());
1509     }
1510
1511     delete[] piNewDims;
1512     //find a way to copy existing data to new variable ...
1513     int iNewPos = 0;
1514     int* piIndexes = new int[iOrigDims];
1515     int* piViewDims = new int[iOrigDims];
1516     for (int i = 0 ; i < iOrigDims ; i++)
1517     {
1518         piViewDims[i] = getVarMaxDim(i, iOrigDims);
1519     }
1520
1521     for (int i = 0 ; i < getSize() ; i++)
1522     {
1523         bool bByPass = false;
1524         getIndexesWithDims(i, piIndexes, piViewDims, iOrigDims);
1525
1526         //check if piIndexes use removed indexes
1527         for (int j = 0 ; j < iNotEntireSize ; j++)
1528         {
1529             if ((piNotEntireIndex[j] - 1) == piIndexes[iNotEntire])
1530             {
1531                 //by pass this value
1532                 bByPass = true;
1533                 break;
1534             }
1535         }
1536
1537         if (bByPass == false)
1538         {
1539             //compute new index
1540             if (isComplex())
1541             {
1542                 pOut->set(iNewPos, getImg(i));
1543             }
1544             else
1545             {
1546                 pOut->set(iNewPos, get(i));
1547             }
1548             iNewPos++;
1549         }
1550     }
1551
1552     //free allocated data
1553     for (int i = 0 ; i < iDims ; i++)
1554     {
1555         if (pArg[i] != (*_pArgs)[i])
1556         {
1557             delete pArg[i];
1558         }
1559     }
1560
1561     delete[] piIndexes;
1562     delete[] piViewDims;
1563
1564     //free pArg content
1565     cleanIndexesArguments(_pArgs, &pArg);
1566
1567     return pOut;
1568 }
1569
1570 bool Sparse::append(int r, int c, types::Sparse SPARSE_CONST* src)
1571 {
1572     //        std::wcerr << L"to a sparse of size"<<getRows() << L","<<getCols() << L" should append @"<<r << L","<<c<< "a sparse:"<< src->toString(32,80)<<std::endl;
1573     if (src->isComplex())
1574     {
1575         toComplex();
1576     }
1577     if (isComplex())
1578     {
1579         if (src->isComplex())
1580         {
1581             doAppend(*(src->matrixCplx), r, c, *matrixCplx);
1582         }
1583         else
1584         {
1585             doAppend(*(src->matrixReal), r, c, *matrixCplx);
1586         }
1587     }
1588     else
1589     {
1590         doAppend(*(src->matrixReal), r, c, *matrixReal);
1591     }
1592
1593     finalize();
1594
1595     return true; // realloc is meaningless for sparse matrices
1596 }
1597
1598 /*
1599 * create a new Sparse of dims according to resSize and fill it from currentSparse (along coords)
1600 */
1601 InternalType* Sparse::extract(typed_list* _pArgs)
1602 {
1603     Sparse* pOut        = NULL;
1604     int iDims           = (int)_pArgs->size();
1605     typed_list pArg;
1606
1607     int* piMaxDim       = new int[iDims];
1608     int* piCountDim     = new int[iDims];
1609
1610     //evaluate each argument and replace by appropriate value and compute the count of combinations
1611     int iSeqCount = checkIndexesArguments(this, _pArgs, &pArg, piMaxDim, piCountDim);
1612     if (iSeqCount == 0)
1613     {
1614         //free pArg content
1615         cleanIndexesArguments(_pArgs, &pArg);
1616         if (_pArgs->size() == 0)
1617         {
1618             //a()
1619             delete[] piMaxDim;
1620             delete[] piCountDim;
1621             //free pArg content
1622             cleanIndexesArguments(_pArgs, &pArg);
1623             return this;
1624         }
1625         else
1626         {
1627             //a([])
1628             delete[] piMaxDim;
1629             delete[] piCountDim;
1630             //free pArg content
1631             cleanIndexesArguments(_pArgs, &pArg);
1632             return Double::Empty();
1633         }
1634     }
1635
1636     if (iDims < 2)
1637     {
1638         if (piMaxDim[0] <= getSize())
1639         {
1640             int iNewRows = 0;
1641             int iNewCols = 0;
1642
1643             if (getRows() == 1 && getCols() != 1 && (*_pArgs)[0]->isColon() == false)
1644             {
1645                 //special case for row vector
1646                 iNewRows = 1;
1647                 iNewCols = piCountDim[0];
1648             }
1649             else
1650             {
1651                 iNewRows = piCountDim[0];
1652                 iNewCols = 1;
1653             }
1654
1655             pOut = new Sparse(iNewRows, iNewCols, isComplex());
1656             double* pIdx = pArg[0]->getAs<Double>()->get();
1657             for (int i = 0 ; i < iSeqCount ; i++)
1658             {
1659                 if (pIdx[i] < 1)
1660                 {
1661                     delete pOut;
1662                     pOut = NULL;
1663                     break;
1664                 }
1665                 int iRowRead = static_cast<int>(pIdx[i] - 1) % getRows();
1666                 int iColRead = static_cast<int>(pIdx[i] - 1) / getRows();
1667
1668                 int iRowWrite = static_cast<int>(i) % iNewRows;
1669                 int iColWrite = static_cast<int>(i) / iNewRows;
1670                 if (isComplex())
1671                 {
1672                     std::complex<double> dbl = getImg(iRowRead, iColRead);
1673                     if (dbl.real() != 0 || dbl.imag() != 0)
1674                     {
1675                         //only non zero values
1676                         pOut->set(iRowWrite, iColWrite, dbl, false);
1677                     }
1678                 }
1679                 else
1680                 {
1681                     double dbl = get(iRowRead, iColRead);
1682                     if (dbl != 0)
1683                     {
1684                         //only non zero values
1685                         pOut->set(iRowWrite, iColWrite, dbl, false);
1686                     }
1687                 }
1688             }
1689         }
1690         else
1691         {
1692             delete[] piMaxDim;
1693             delete[] piCountDim;
1694             //free pArg content
1695             cleanIndexesArguments(_pArgs, &pArg);
1696             return NULL;
1697         }
1698     }
1699     else
1700     {
1701         if (piMaxDim[0] <= getRows() && piMaxDim[1] <= getCols())
1702         {
1703             double* pIdxRow = pArg[0]->getAs<Double>()->get();
1704             double* pIdxCol = pArg[1]->getAs<Double>()->get();
1705
1706             int iNewRows = pArg[0]->getAs<Double>()->getSize();
1707             int iNewCols = pArg[1]->getAs<Double>()->getSize();
1708
1709             pOut = new Sparse(iNewRows, iNewCols, isComplex());
1710
1711             int iPos = 0;
1712             for (int iRow = 0 ; iRow < iNewRows ; iRow++)
1713             {
1714                 for (int iCol = 0 ; iCol < iNewCols ; iCol++)
1715                 {
1716                     if ((pIdxRow[iRow] < 1) || (pIdxCol[iCol] < 1))
1717                     {
1718                         delete pOut;
1719                         pOut = NULL;
1720                         break;
1721                     }
1722                     if (isComplex())
1723                     {
1724                         std::complex<double> dbl = getImg((int)pIdxRow[iRow] - 1, (int)pIdxCol[iCol] - 1);
1725                         if (dbl.real() != 0 || dbl.imag() != 0)
1726                         {
1727                             //only non zero values
1728                             pOut->set(iRow, iCol, dbl, false);
1729                         }
1730                     }
1731                     else
1732                     {
1733                         double dbl = get((int)pIdxRow[iRow] - 1, (int)pIdxCol[iCol] - 1);
1734                         if (dbl != 0)
1735                         {
1736                             //only non zero values
1737                             pOut->set(iRow, iCol, dbl, false);
1738                         }
1739                     }
1740                     iPos++;
1741                 }
1742             }
1743         }
1744         else
1745         {
1746             delete[] piMaxDim;
1747             delete[] piCountDim;
1748             //free pArg content
1749             cleanIndexesArguments(_pArgs, &pArg);
1750             return NULL;
1751         }
1752     }
1753
1754     pOut->finalize();
1755
1756     delete[] piMaxDim;
1757     delete[] piCountDim;
1758     //free pArg content
1759     cleanIndexesArguments(_pArgs, &pArg);
1760
1761     return pOut;
1762 }
1763
1764 Sparse* Sparse::extract(int nbCoords, int SPARSE_CONST* coords, int SPARSE_CONST* maxCoords, int SPARSE_CONST* resSize, bool asVector) SPARSE_CONST
1765 {
1766     if ( (asVector && maxCoords[0] > getSize()) ||
1767     (asVector == false && maxCoords[0] > getRows()) ||
1768     (asVector == false && maxCoords[1] > getCols()))
1769     {
1770         return 0;
1771     }
1772
1773     bool const cplx(isComplex());
1774     Sparse * pSp (0);
1775     if (asVector)
1776     {
1777         pSp = (getRows() == 1) ?  new Sparse(1, resSize[0], cplx) : new Sparse(resSize[0], 1, cplx);
1778     }
1779     else
1780     {
1781         pSp = new Sparse(resSize[0], resSize[1], cplx);
1782     }
1783     //        std::cerr<<"extracted sparse:"<<pSp->getRows()<<", "<<pSp->getCols()<<"seqCount="<<nbCoords<<"maxDim="<<maxCoords[0] <<","<< maxCoords[1]<<std::endl;
1784     if (! (asVector
1785     ? copyToSparse(*this,  Coords<true>(coords, getRows()), nbCoords
1786     , *pSp, RowWiseFullIterator(pSp->getRows(), pSp->getCols()))
1787     : copyToSparse(*this,  Coords<false>(coords), nbCoords
1788     , *pSp, RowWiseFullIterator(pSp->getRows(), pSp->getCols()))))
1789     {
1790         delete pSp;
1791         pSp = NULL;
1792     }
1793     return pSp;
1794 }
1795
1796 bool Sparse::invoke(typed_list & in, optional_list & /*opt*/, int /*_iRetCount*/, typed_list & out, const ast::Exp & e)
1797 {
1798     if (in.size() == 0)
1799     {
1800         out.push_back(this);
1801     }
1802     else
1803     {
1804         InternalType * _out = extract(&in);
1805         if (!_out)
1806         {
1807             std::wostringstream os;
1808             os << _W("Invalid index.\n");
1809             throw ast::InternalError(os.str(), 999, e.getLocation());
1810         }
1811         out.push_back(_out);
1812     }
1813
1814     return true;
1815 }
1816
1817
1818 bool Sparse::isInvokable() const
1819 {
1820     return true;
1821 }
1822
1823 bool Sparse::hasInvokeOption() const
1824 {
1825     return false;
1826 }
1827
1828 int Sparse::getInvokeNbIn()
1829 {
1830     return -1;
1831 }
1832
1833 int Sparse::getInvokeNbOut()
1834 {
1835     return 1;
1836 }
1837
1838 /*
1839 coords are Scilab 1-based
1840 extract std::make_pair(coords, asVector), rowIter
1841 */
1842 template<typename Src, typename SrcTraversal, typename Sz, typename DestTraversal>
1843 bool Sparse::copyToSparse(Src SPARSE_CONST& src, SrcTraversal srcTrav, Sz n, Sparse& sp, DestTraversal destTrav)
1844 {
1845     if (!(src.isComplex() || sp.isComplex()))
1846     {
1847         mycopy_n(makeMatrixIterator<double>(src, srcTrav), n
1848                  , makeMatrixIterator<double>(*sp.matrixReal, destTrav));
1849     }
1850     else
1851     {
1852         sp.toComplex();
1853         mycopy_n(makeMatrixIterator<std::complex<double> >(src, srcTrav), n
1854                  , makeMatrixIterator<std::complex<double> >(*sp.matrixCplx, destTrav));
1855     }
1856
1857     sp.finalize();
1858     return true;
1859 }
1860
1861 // GenericType because we might return a Double* for scalar operand
1862 Sparse* Sparse::add(Sparse const& o) const
1863 {
1864     RealSparse_t* realSp(0);
1865     CplxSparse_t* cplxSp(0);
1866     if (isComplex() == false && o.isComplex() == false)
1867     {
1868         //R + R -> R
1869         realSp = new RealSparse_t(*matrixReal + * (o.matrixReal));
1870     }
1871     else if (isComplex() == false && o.isComplex() == true)
1872     {
1873         cplxSp = new CplxSparse_t(matrixReal->cast<std::complex<double> >() + * (o.matrixCplx));
1874     }
1875     else if (isComplex() == true && o.isComplex() == false)
1876     {
1877         cplxSp = new CplxSparse_t(*matrixCplx + o.matrixReal->cast<std::complex<double> >());
1878     }
1879     else if (isComplex() == true && o.isComplex() == true)
1880     {
1881         cplxSp = new CplxSparse_t(*matrixCplx + * (o.matrixCplx));
1882     }
1883
1884     return new Sparse(realSp, cplxSp);
1885 }
1886
1887 Sparse* Sparse::substract(Sparse const& o) const
1888 {
1889     RealSparse_t* realSp(0);
1890     CplxSparse_t* cplxSp(0);
1891     if (isComplex() == false && o.isComplex() == false)
1892     {
1893         //R - R -> R
1894         realSp = new RealSparse_t(*matrixReal - * (o.matrixReal));
1895     }
1896     else if (isComplex() == false && o.isComplex() == true)
1897     {
1898         //R - C -> C
1899         cplxSp = new CplxSparse_t(matrixReal->cast<std::complex<double> >() - * (o.matrixCplx));
1900     }
1901     else if (isComplex() == true && o.isComplex() == false)
1902     {
1903         //C - R -> C
1904         cplxSp = new CplxSparse_t(*matrixCplx - o.matrixReal->cast<std::complex<double> >());
1905     }
1906     else if (isComplex() == true && o.isComplex() == true)
1907     {
1908         //C - C -> C
1909         cplxSp = new CplxSparse_t(*matrixCplx - * (o.matrixCplx));
1910     }
1911
1912     return new Sparse(realSp, cplxSp);
1913 }
1914
1915 Sparse* Sparse::multiply(double s) const
1916 {
1917     return new Sparse( isComplex() ? 0 : new RealSparse_t((*matrixReal)*s)
1918                        , isComplex() ? new CplxSparse_t((*matrixCplx)*s) : 0);
1919 }
1920
1921 Sparse* Sparse::multiply(std::complex<double> s) const
1922 {
1923     return new Sparse( 0
1924                        , isComplex() ? new CplxSparse_t((*matrixCplx) * s) : new CplxSparse_t((*matrixReal) * s));
1925 }
1926
1927 Sparse* Sparse::multiply(Sparse const& o) const
1928 {
1929     RealSparse_t* realSp(0);
1930     CplxSparse_t* cplxSp(0);
1931
1932     if (isComplex() == false && o.isComplex() == false)
1933     {
1934         realSp = new RealSparse_t(*matrixReal **(o.matrixReal));
1935     }
1936     else if (isComplex() == false && o.isComplex() == true)
1937     {
1938         cplxSp = new CplxSparse_t(matrixReal->cast<std::complex<double> >() **(o.matrixCplx));
1939     }
1940     else if (isComplex() == true && o.isComplex() == false)
1941     {
1942         cplxSp = new CplxSparse_t(*matrixCplx * o.matrixReal->cast<std::complex<double> >());
1943     }
1944     else if (isComplex() == true && o.isComplex() == true)
1945     {
1946         cplxSp = new CplxSparse_t(*matrixCplx **(o.matrixCplx));
1947     }
1948
1949     return new Sparse(realSp, cplxSp);
1950 }
1951
1952 Sparse* Sparse::dotMultiply(Sparse SPARSE_CONST& o) const
1953 {
1954     RealSparse_t* realSp(0);
1955     CplxSparse_t* cplxSp(0);
1956     if (isComplex() == false && o.isComplex() == false)
1957     {
1958         realSp = new RealSparse_t(matrixReal->cwiseProduct(*(o.matrixReal)));
1959     }
1960     else if (isComplex() == false && o.isComplex() == true)
1961     {
1962         cplxSp = new CplxSparse_t(matrixReal->cast<std::complex<double> >().cwiseProduct( *(o.matrixCplx)));
1963     }
1964     else if (isComplex() == true && o.isComplex() == false)
1965     {
1966         cplxSp = new CplxSparse_t(matrixCplx->cwiseProduct(o.matrixReal->cast<std::complex<double> >()));
1967     }
1968     else if (isComplex() == true && o.isComplex() == true)
1969     {
1970         cplxSp = new CplxSparse_t(matrixCplx->cwiseProduct(*(o.matrixCplx)));
1971     }
1972
1973     return new Sparse(realSp, cplxSp);
1974 }
1975
1976 Sparse* Sparse::dotDivide(Sparse SPARSE_CONST& o) const
1977 {
1978     RealSparse_t* realSp(0);
1979     CplxSparse_t* cplxSp(0);
1980     if (isComplex() == false && o.isComplex() == false)
1981     {
1982         realSp = new RealSparse_t(matrixReal->cwiseQuotient(*(o.matrixReal)));
1983     }
1984     else if (isComplex() == false && o.isComplex() == true)
1985     {
1986         cplxSp = new CplxSparse_t(matrixReal->cast<std::complex<double> >().cwiseQuotient( *(o.matrixCplx)));
1987     }
1988     else if (isComplex() == true && o.isComplex() == false)
1989     {
1990         cplxSp = new CplxSparse_t(matrixCplx->cwiseQuotient(o.matrixReal->cast<std::complex<double> >()));
1991     }
1992     else if (isComplex() == true && o.isComplex() == true)
1993     {
1994         cplxSp = new CplxSparse_t(matrixCplx->cwiseQuotient(*(o.matrixCplx)));
1995     }
1996
1997     return new Sparse(realSp, cplxSp);
1998 }
1999
2000 int Sparse::newCholLLT(Sparse** _SpPermut, Sparse** _SpFactor) const
2001 {
2002     typedef Eigen::SparseMatrix<double, Eigen::ColMajor> RealSparseCol_t;
2003     RealSparseCol_t spColMajor = RealSparseCol_t((const RealSparse_t&) * matrixReal);
2004
2005     // Constructs and performs the LLT factorization of sparse
2006     Eigen::SimplicialLLT<RealSparseCol_t> pLLT(spColMajor);
2007     int iInfo = pLLT.info();
2008     if (iInfo != Eigen::Success)
2009     {
2010         *_SpFactor = NULL;
2011         *_SpPermut = NULL;
2012         return iInfo;
2013     }
2014
2015     // Get the lower matrix of factorization.
2016     // The new RealSparse_t will be setted in Sparse without copy.
2017     *_SpFactor = new Sparse(new RealSparse_t(pLLT.matrixL()), NULL);
2018
2019     // Get the permutation matrix.
2020     Eigen::PermutationMatrix<Eigen::Dynamic, Eigen::Dynamic, int> p = pLLT.permutationP();
2021     *_SpPermut = new Sparse(p.rows(), p.cols());
2022     for (int i = 0; i < p.rows(); i++)
2023     {
2024         (*_SpPermut)->set(i, p.indices()[i], 1, false);
2025     }
2026
2027     (*_SpPermut)->finalize();
2028
2029     return iInfo;
2030 }
2031
2032 bool Sparse::transpose(InternalType *& out)
2033 {
2034     out = new Sparse(matrixReal ? new RealSparse_t(matrixReal->transpose()) : 0, matrixCplx ? new CplxSparse_t(matrixCplx->transpose()) : 0);
2035     return true;
2036 }
2037
2038 bool Sparse::adjoint(InternalType *& out)
2039 {
2040     out = new Sparse(matrixReal ? new RealSparse_t(matrixReal->adjoint()) : 0, matrixCplx ? new CplxSparse_t(matrixCplx->adjoint()) : 0);
2041     return true;
2042 }
2043
2044 struct BoolCast
2045 {
2046     BoolCast(std::complex<double> const& c): b(c.real() || c.imag()) {}
2047     operator bool () const
2048     {
2049         return b;
2050     }
2051     operator double() const
2052     {
2053         return b ? 1. : 0.;
2054     }
2055     bool b;
2056 };
2057 Sparse* Sparse::newOnes() const
2058 {
2059     // result is never cplx
2060     return new Sparse( matrixReal
2061                        ? new RealSparse_t(matrixReal->cast<bool>().cast<double>())
2062                        : new RealSparse_t(matrixCplx->cast<BoolCast>().cast<double>())
2063                        , 0);
2064 }
2065
2066 struct RealCast
2067 {
2068     RealCast(std::complex<double> const& c): b(c.real()) {}
2069     operator bool () const
2070     {
2071         return b != 0;
2072     }
2073     operator double() const
2074     {
2075         return b;
2076     }
2077     double b;
2078 };
2079 Sparse* Sparse::newReal() const
2080 {
2081     return new Sparse( matrixReal
2082                        ? matrixReal
2083                        : new RealSparse_t(matrixCplx->cast<RealCast>().cast<double>())
2084                        , 0);
2085 }
2086
2087 std::size_t Sparse::nonZeros() const
2088 {
2089     if (isComplex())
2090     {
2091         return matrixCplx->nonZeros();
2092     }
2093     else
2094     {
2095         return matrixReal->nonZeros();
2096     }
2097 }
2098 std::size_t Sparse::nonZeros(std::size_t r) const
2099 {
2100     std::size_t res;
2101     if (matrixReal)
2102     {
2103         int* piIndex = matrixReal->outerIndexPtr();
2104         res = piIndex[r + 1] - piIndex[r];
2105     }
2106     else
2107     {
2108         int* piIndex = matrixCplx->outerIndexPtr();
2109         res = piIndex[r + 1] - piIndex[r];
2110     }
2111
2112     return res;
2113 }
2114
2115 int* Sparse::getNbItemByRow(int* _piNbItemByRows)
2116 {
2117     int* piNbItemByCols = new int[getRows() + 1];
2118     if (isComplex())
2119     {
2120         mycopy_n(matrixCplx->outerIndexPtr(), getRows() + 1, piNbItemByCols);
2121     }
2122     else
2123     {
2124         mycopy_n(matrixReal->outerIndexPtr(), getRows() + 1, piNbItemByCols);
2125     }
2126
2127     for (int i = 0 ; i < getRows() ; i++)
2128     {
2129         _piNbItemByRows[i] = piNbItemByCols[i + 1] - piNbItemByCols[i];
2130     }
2131
2132     delete[] piNbItemByCols;
2133     return _piNbItemByRows;
2134 }
2135
2136 int* Sparse::getColPos(int* _piColPos)
2137 {
2138     if (isComplex())
2139     {
2140         mycopy_n(matrixCplx->innerIndexPtr(), nonZeros(), _piColPos);
2141     }
2142     else
2143     {
2144         mycopy_n(matrixReal->innerIndexPtr(), nonZeros(), _piColPos);
2145     }
2146
2147     for (int i = 0; i < nonZeros(); i++)
2148     {
2149         _piColPos[i]++;
2150     }
2151
2152     return _piColPos;
2153 }
2154
2155 int* Sparse::getInnerPtr(int* count)
2156 {
2157     int* ret = nullptr;
2158     if (isComplex())
2159     {
2160         ret = matrixCplx->innerIndexPtr();
2161         *count = matrixCplx->innerSize();
2162     }
2163     else
2164     {
2165         ret = matrixReal->innerIndexPtr();
2166         *count = matrixReal->innerSize();
2167     }
2168
2169     return ret;
2170 }
2171
2172 int* Sparse::getOuterPtr(int* count)
2173 {
2174     int* ret = nullptr;
2175     if (isComplex())
2176     {
2177         ret = matrixCplx->outerIndexPtr();
2178         *count = matrixCplx->outerSize();
2179     }
2180     else
2181     {
2182         ret = matrixReal->outerIndexPtr();
2183         *count = matrixReal->outerSize();
2184     }
2185
2186     return ret;
2187 }
2188
2189 namespace
2190 {
2191 template<typename S> struct GetReal: std::unary_function<typename S::InnerIterator, double>
2192 {
2193     double operator()(typename S::InnerIterator it) const
2194     {
2195         return it.value();
2196     }
2197 };
2198 template<> struct GetReal< Eigen::SparseMatrix<std::complex<double >, Eigen::RowMajor > >
2199         : std::unary_function<Sparse::CplxSparse_t::InnerIterator, double>
2200 {
2201     double operator()( Sparse::CplxSparse_t::InnerIterator it) const
2202     {
2203         return it.value().real();
2204     }
2205 };
2206 template<typename S> struct GetImag: std::unary_function<typename S::InnerIterator, double>
2207 {
2208     double operator()(typename S::InnerIterator it) const
2209     {
2210         return it.value().imag();
2211     }
2212 };
2213 template<typename S> struct GetRow: std::unary_function<typename S::InnerIterator, int>
2214 {
2215     int operator()(typename S::InnerIterator it) const
2216     {
2217         return it.row() + 1;
2218     }
2219 };
2220 template<typename S> struct GetCol: std::unary_function<typename S::InnerIterator, int>
2221 {
2222     int operator()(typename S::InnerIterator it) const
2223     {
2224         return it.col() + 1;
2225     }
2226 };
2227
2228 template<typename S, typename Out, typename F> Out sparseTransform(S& s, Out o, F f)
2229 {
2230     for (std::size_t k(0); k < s.outerSize(); ++k)
2231     {
2232         for (typename S::InnerIterator it(s, (int)k); it; ++it, ++o)
2233         {
2234             *o = f(it);
2235         }
2236     }
2237     return o;
2238 }
2239 }
2240
2241 std::pair<double*, double*> Sparse::outputValues(double* outReal, double* outImag)const
2242 {
2243     return matrixReal
2244            ? std::make_pair(sparseTransform(*matrixReal, outReal, GetReal<RealSparse_t>()), outImag)
2245            : std::make_pair(sparseTransform(*matrixCplx, outReal, GetReal<CplxSparse_t>())
2246                             , sparseTransform(*matrixCplx, outImag, GetImag<CplxSparse_t>()));
2247 }
2248
2249 int* Sparse::outputRowCol(int* out)const
2250 {
2251     return matrixReal
2252            ? sparseTransform(*matrixReal, sparseTransform(*matrixReal, out, GetRow<RealSparse_t>()), GetCol<RealSparse_t>())
2253            : sparseTransform(*matrixCplx, sparseTransform(*matrixCplx, out, GetRow<CplxSparse_t>()), GetCol<CplxSparse_t>());
2254 }
2255 double* Sparse::outputCols(double* out) const
2256 {
2257     if (isComplex())
2258     {
2259         mycopy_n(matrixCplx->innerIndexPtr(), nonZeros(), out);
2260     }
2261     else
2262     {
2263         mycopy_n(matrixReal->innerIndexPtr(), nonZeros(), out);
2264     }
2265
2266     return std::transform(out, out, out, std::bind2nd(std::plus<double>(), 1));
2267
2268 }
2269
2270 void Sparse::opposite(void)
2271 {
2272     if (isComplex())
2273     {
2274         std::complex<double>* data = matrixCplx->valuePtr();
2275         std::transform(data, data + matrixCplx->nonZeros(), data, std::negate<std::complex<double> >());
2276     }
2277     else
2278     {
2279         double* data = matrixReal->valuePtr();
2280         std::transform(data, data + matrixReal->nonZeros(), data, std::negate<double>());
2281     }
2282 }
2283
2284 SparseBool* Sparse::newLessThan(Sparse &o)
2285 {
2286     //only real values !
2287
2288     //return cwiseOp<std::less>(*this, o);
2289     int rowL = getRows();
2290     int colL = getCols();
2291
2292     int rowR = o.getRows();
2293     int colR = o.getCols();
2294     int row = std::max(rowL, rowR);
2295     int col = std::max(colL, colR);
2296
2297     //create a boolean sparse matrix with dims of sparses
2298     types::SparseBool* ret = new types::SparseBool(row, col);
2299     if (isScalar() && o.isScalar())
2300     {
2301         double l = get(0, 0);
2302         double r = o.get(0, 0);
2303         ret->set(0, 0, l < r, false);
2304     }
2305     else if (isScalar())
2306     {
2307         int nnzR = static_cast<int>(o.nonZeros());
2308         std::vector<int> rowcolR(nnzR * 2, 0);
2309         o.outputRowCol(rowcolR.data());
2310
2311         //compare all items of R with R[0]
2312         double l = get(0, 0);
2313         if (l < 0)
2314         {
2315             //set true
2316             ret->setTrue(false);
2317         }
2318
2319         for (int i = 0; i < nnzR; ++i)
2320         {
2321             double r = o.get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2322             ret->set(rowcolR[i] - 1, rowcolR[i + nnzR] - 1, l < r, false);
2323         }
2324     }
2325     else if (o.isScalar())
2326     {
2327         int nnzL = static_cast<int>(nonZeros());
2328         std::vector<int> rowcolL(nnzL * 2, 0);
2329         outputRowCol(rowcolL.data());
2330
2331         double r = o.get(0, 0);
2332         if (r >= 0)
2333         {
2334             ret->setTrue(true);
2335         }
2336
2337         for (int i = 0; i < nnzL; ++i)
2338         {
2339             double l = get(rowcolL[i] - 1, rowcolL[i + nnzL] - 1);
2340             ret->set(rowcolL[i] - 1, rowcolL[i + nnzL] - 1, l < r, false);
2341         }
2342     }
2343     else
2344     {
2345         int nnzR = static_cast<int>(o.nonZeros());
2346         std::vector<int> rowcolR(nnzR * 2, 0);
2347         o.outputRowCol(rowcolR.data());
2348         int nnzL = static_cast<int>(nonZeros());
2349         std::vector<int> rowcolL(nnzL * 2, 0);
2350         outputRowCol(rowcolL.data());
2351         //set all values to %t
2352         ret->setFalse(false);
2353
2354         for (int i = 0; i < nnzL; ++i)
2355         {
2356             double l = get(rowcolL[i] - 1, rowcolL[i + nnzL] - 1);
2357             ret->set(rowcolL[i] - 1, rowcolL[i + nnzL] - 1, l < 0, false);
2358         }
2359         ret->finalize();
2360
2361         //set _pR[i] == _pL[i] for each _pR values
2362         for (int i = 0; i < nnzR; ++i)
2363         {
2364             //get l and r following non zeros value of R
2365             double l = get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2366             double r = o.get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2367             //set value following non zeros value of R
2368             ret->set(rowcolR[i] - 1, rowcolR[i + nnzR] - 1, l < r, false);
2369         }
2370     }
2371
2372     ret->finalize();
2373     return ret;
2374 }
2375
2376 SparseBool* Sparse::newNotEqualTo(Sparse const&o) const
2377 {
2378     return cwiseOp<std::not_equal_to>(*this, o);
2379 }
2380
2381 SparseBool* Sparse::newLessOrEqual(Sparse &o)
2382 {
2383     //only real values !
2384
2385     //return cwiseOp<std::less>(*this, o);
2386     int rowL = getRows();
2387     int colL = getCols();
2388
2389     int rowR = o.getRows();
2390     int colR = o.getCols();
2391     int row = std::max(rowL, rowR);
2392     int col = std::max(colL, colR);
2393
2394     //create a boolean sparse matrix with dims of sparses
2395     types::SparseBool* ret = new types::SparseBool(row, col);
2396     if (isScalar() && o.isScalar())
2397     {
2398         double l = get(0, 0);
2399         double r = o.get(0, 0);
2400         ret->set(0, 0, l <= r, false);
2401     }
2402     else if (isScalar())
2403     {
2404         int nnzR = static_cast<int>(o.nonZeros());
2405         std::vector<int> rowcolR(nnzR * 2, 0);
2406         o.outputRowCol(rowcolR.data());
2407
2408         //compare all items of R with R[0]
2409         double l = get(0, 0);
2410         if (l <= 0)
2411         {
2412             //set true
2413             ret->setTrue(false);
2414         }
2415
2416         for (int i = 0; i < nnzR; ++i)
2417         {
2418             double r = o.get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2419             ret->set(rowcolR[i] - 1, rowcolR[i + nnzR] - 1, l <= r, false);
2420         }
2421     }
2422     else if (o.isScalar())
2423     {
2424         int nnzL = static_cast<int>(nonZeros());
2425         std::vector<int> rowcolL(nnzL * 2, 0);
2426         outputRowCol(rowcolL.data());
2427
2428         double r = o.get(0, 0);
2429         if (r > 0)
2430         {
2431             ret->setTrue(true);
2432         }
2433
2434         for (int i = 0; i < nnzL; ++i)
2435         {
2436             double l = get(rowcolL[i] - 1, rowcolL[i + nnzL] - 1);
2437             ret->set(rowcolL[i] - 1, rowcolL[i + nnzL] - 1, l <= r, false);
2438         }
2439     }
2440     else
2441     {
2442         int nnzR = static_cast<int>(o.nonZeros());
2443         std::vector<int> rowcolR(nnzR * 2, 0);
2444         o.outputRowCol(rowcolR.data());
2445         int nnzL = static_cast<int>(nonZeros());
2446         std::vector<int> rowcolL(nnzL * 2, 0);
2447         outputRowCol(rowcolL.data());
2448         //set all values to %t
2449         ret->setTrue(false);
2450
2451         for (int i = 0; i < nnzL; ++i)
2452         {
2453             double l = get(rowcolL[i] - 1, rowcolL[i + nnzL] - 1);
2454             ret->set(rowcolL[i] - 1, rowcolL[i + nnzL] - 1, l <= 0, false);
2455         }
2456         ret->finalize();
2457
2458         //set _pR[i] == _pL[i] for each _pR values
2459         for (int i = 0; i < nnzR; ++i)
2460         {
2461             //get l and r following non zeros value of R
2462             double l = get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2463             double r = o.get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2464             //set value following non zeros value of R
2465             ret->set(rowcolR[i] - 1, rowcolR[i + nnzR] - 1, l <= r, false);
2466         }
2467     }
2468
2469     ret->finalize();
2470     return ret;
2471 }
2472
2473 SparseBool* Sparse::newEqualTo(Sparse &o)
2474 {
2475     int rowL = getRows();
2476     int colL = getCols();
2477
2478     int rowR = o.getRows();
2479     int colR = o.getCols();
2480     int row = std::max(rowL, rowR);
2481     int col = std::max(colL, colR);
2482
2483     //create a boolean sparse matrix with dims of sparses
2484     types::SparseBool* ret = new types::SparseBool(row, col);
2485     if (isScalar() && o.isScalar())
2486     {
2487         if (isComplex() || o.isComplex())
2488         {
2489             std::complex<double> l = getImg(0, 0);
2490             std::complex<double> r = o.getImg(0, 0);
2491             ret->set(0, 0, l == r, false);
2492         }
2493         else
2494         {
2495             double l = get(0, 0);
2496             double r = o.get(0, 0);
2497             ret->set(0, 0, l == r, false);
2498         }
2499     }
2500     else if (isScalar())
2501     {
2502         int nnzR = static_cast<int>(o.nonZeros());
2503         std::vector<int> rowcolR(nnzR * 2, 0);
2504         o.outputRowCol(rowcolR.data());
2505
2506         //compare all items of R with R[0]
2507         if (isComplex() || o.isComplex())
2508         {
2509             std::complex<double> l = getImg(0, 0);
2510             for (int i = 0; i < nnzR; ++i)
2511             {
2512                 std::complex<double> r = o.getImg(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2513                 ret->set(rowcolR[i] - 1, rowcolR[i + nnzR] - 1, l == r, false);
2514             }
2515         }
2516         else
2517         {
2518             double l = get(0, 0);
2519             for (int i = 0; i < nnzR; ++i)
2520             {
2521                 double r = o.get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2522                 ret->set(rowcolR[i] - 1, rowcolR[i + nnzR] - 1, l == r, false);
2523             }
2524         }
2525     }
2526     else if (o.isScalar())
2527     {
2528         int nnzL = static_cast<int>(nonZeros());
2529         std::vector<int> rowcolL(nnzL * 2, 0);
2530         outputRowCol(rowcolL.data());
2531
2532         if (isComplex() || o.isComplex())
2533         {
2534             std::complex<double> r = o.getImg(0, 0);
2535             for (int i = 0; i < nnzL; ++i)
2536             {
2537                 std::complex<double> l = getImg(rowcolL[i] - 1, rowcolL[i + nnzL] - 1);
2538                 ret->set(rowcolL[i] - 1, rowcolL[i + nnzL] - 1, l == r, false);
2539             }
2540         }
2541         else
2542         {
2543             double r = get(0, 0);
2544             for (int i = 0; i < nnzL; ++i)
2545             {
2546                 double l = get(rowcolL[i] - 1, rowcolL[i + nnzL] - 1);
2547                 ret->set(rowcolL[i] - 1, rowcolL[i + nnzL] - 1, l == r, false);
2548             }
2549         }
2550     }
2551     else
2552     {
2553         int nnzR = static_cast<int>(o.nonZeros());
2554         std::vector<int> rowcolR(nnzR * 2, 0);
2555         o.outputRowCol(rowcolR.data());
2556         int nnzL = static_cast<int>(nonZeros());
2557         std::vector<int> rowcolL(nnzL * 2, 0);
2558         outputRowCol(rowcolL.data());
2559         //set all values to %t
2560         ret->setTrue(false);
2561         //set %f in each pL values
2562         for (int i = 0; i < nnzL; ++i)
2563         {
2564             ret->set(rowcolL[i] - 1, rowcolL[i + nnzL] - 1, false, false);
2565         }
2566         ret->finalize();
2567
2568         //set _pR[i] == _pL[i] for each _pR values
2569         if (isComplex() || o.isComplex())
2570         {
2571             for (int i = 0; i < nnzR; ++i)
2572             {
2573                 //get l and r following non zeros value of R
2574                 std::complex<double> l = getImg(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2575                 std::complex<double> r = o.getImg(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2576                 //set value following non zeros value of R
2577                 ret->set(rowcolR[i] - 1, rowcolR[i + nnzR] - 1, l == r, false);
2578             }
2579         }
2580         else
2581         {
2582             for (int i = 0; i < nnzR; ++i)
2583             {
2584                 //get l and r following non zeros value of R
2585                 double l = get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2586                 double r = o.get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
2587                 //set value following non zeros value of R
2588                 ret->set(rowcolR[i] - 1, rowcolR[i + nnzR] - 1, l == r, false);
2589             }
2590         }
2591     }
2592
2593     ret->finalize();
2594     return ret;
2595 }
2596
2597 bool Sparse::reshape(int* _piDims, int _iDims)
2598 {
2599     bool bOk = false;
2600     int iCols = 1;
2601
2602     if (_iDims == 2)
2603     {
2604         iCols = _piDims[1];
2605     }
2606
2607     if (_iDims <= 2)
2608     {
2609         bOk = reshape(_piDims[0], iCols);
2610     }
2611
2612     return bOk;
2613 }
2614
2615 bool Sparse::reshape(int _iNewRows, int _iNewCols)
2616 {
2617     if (_iNewRows * _iNewCols != getRows() * getCols())
2618     {
2619         return false;
2620     }
2621
2622     bool res = false;
2623     try
2624     {
2625         if (matrixReal)
2626         {
2627             //item count
2628             size_t iNonZeros = nonZeros();
2629             RealSparse_t *newReal = new RealSparse_t(_iNewRows, _iNewCols);
2630             newReal->reserve((int)iNonZeros);
2631
2632             //coords
2633             int* pRows = new int[iNonZeros * 2];
2634             outputRowCol(pRows);
2635             int* pCols = pRows + iNonZeros;
2636
2637             //values
2638             double* pNonZeroR = new double[iNonZeros];
2639             double* pNonZeroI = new double[iNonZeros];
2640             outputValues(pNonZeroR, pNonZeroI);
2641
2642             typedef Eigen::Triplet<double> triplet;
2643             std::vector<triplet> tripletList;
2644
2645             for (size_t i = 0 ; i < iNonZeros ; i++)
2646             {
2647                 int iCurrentPos = ((int)pCols[i] - 1) * getRows() + ((int)pRows[i] - 1);
2648                 tripletList.push_back(triplet((int)(iCurrentPos % _iNewRows), (int)(iCurrentPos / _iNewRows), pNonZeroR[i]));
2649             }
2650
2651             newReal->setFromTriplets(tripletList.begin(), tripletList.end());
2652
2653             delete matrixReal;
2654             matrixReal = newReal;
2655             delete[] pRows;
2656             delete[] pNonZeroR;
2657             delete[] pNonZeroI;
2658         }
2659         else
2660         {
2661             //item count
2662             size_t iNonZeros = nonZeros();
2663             CplxSparse_t *newCplx = new CplxSparse_t(_iNewRows, _iNewCols);
2664             newCplx->reserve((int)iNonZeros);
2665
2666             //coords
2667             int* pRows = new int[iNonZeros * 2];
2668             outputRowCol(pRows);
2669             int* pCols = pRows + iNonZeros;
2670
2671             //values
2672             double* pNonZeroR = new double[iNonZeros];
2673             double* pNonZeroI = new double[iNonZeros];
2674             outputValues(pNonZeroR, pNonZeroI);
2675
2676             typedef Eigen::Triplet<std::complex<double> > triplet;
2677             std::vector<triplet> tripletList;
2678
2679             for (size_t i = 0 ; i < iNonZeros ; i++)
2680             {
2681                 int iCurrentPos = ((int)pCols[i] - 1) * getRows() + ((int)pRows[i] - 1);
2682                 tripletList.push_back(triplet((int)(iCurrentPos % _iNewRows), (int)(iCurrentPos / _iNewRows), std::complex<double>(pNonZeroR[i], pNonZeroI[i])));
2683             }
2684
2685             newCplx->setFromTriplets(tripletList.begin(), tripletList.end());
2686
2687             delete matrixCplx;
2688             matrixCplx = newCplx;
2689             delete[] pRows;
2690             delete[] pNonZeroR;
2691             delete[] pNonZeroI;
2692         }
2693
2694         m_iRows = _iNewRows;
2695         m_iCols = _iNewCols;
2696         m_iSize = _iNewRows * _iNewCols;
2697
2698         m_iDims = 2;
2699         m_piDims[0] = m_iRows;
2700         m_piDims[1] = m_iCols;
2701
2702         finalize();
2703
2704         res = true;
2705     }
2706     catch (...)
2707     {
2708         res = false;
2709     }
2710     return res;
2711 }
2712
2713 //    SparseBool* SparseBool::new
2714
2715 SparseBool::SparseBool(Bool SPARSE_CONST& src)
2716 {
2717     //compute idx
2718     int size = src.getSize();
2719     int row = src.getRows();
2720     Double* idx = new Double(src.getSize(), 2);
2721     double* p = idx->get();
2722     for (int i = 0; i < size; ++i)
2723     {
2724         p[i] = (double)(i % row) + 1;
2725         p[i + size] = (double)(i / row) + 1;
2726     }
2727     create2(src.getRows(), src.getCols(), src, *idx);
2728     idx->killMe();
2729 #ifndef NDEBUG
2730     Inspector::addItem(this);
2731 #endif
2732 }
2733 /* @param src : Bool matrix to copy into a new sparse matrix
2734 @param idx : Double matrix to use as indexes to get values from the src
2735 **/
2736 SparseBool::SparseBool(Bool SPARSE_CONST& src, Double SPARSE_CONST& idx)
2737 {
2738     int idxrow = idx.getRows();
2739     int rows = static_cast<int>(*std::max_element(idx.get(), idx.get() + idxrow));
2740     int cols = static_cast<int>(*std::max_element(idx.get() + idxrow, idx.get() + idxrow * 2));
2741     create2(rows, cols, src, idx);
2742 #ifndef NDEBUG
2743     Inspector::addItem(this);
2744 #endif
2745 }
2746
2747 /* @param src : Bool matrix to copy into a new sparse matrix
2748 @param idx : Double matrix to use as indexes to get values from the src
2749 @param dims : Double matrix containing the dimensions of the new matrix
2750 **/
2751 SparseBool::SparseBool(Bool SPARSE_CONST& src, Double SPARSE_CONST& idx, Double SPARSE_CONST& dims)
2752 {
2753     create2(static_cast<int>(dims.get(0)), static_cast<int>(dims.get(1)), src, idx);
2754 #ifndef NDEBUG
2755     Inspector::addItem(this);
2756 #endif
2757 }
2758
2759 SparseBool::SparseBool(int _iRows, int _iCols) : matrixBool(new BoolSparse_t(_iRows, _iCols))
2760 {
2761     m_iRows = _iRows;
2762     m_iCols = _iCols;
2763     m_iSize = _iRows * _iCols;
2764     m_iDims = 2;
2765     m_piDims[0] = _iRows;
2766     m_piDims[1] = _iCols;
2767 #ifndef NDEBUG
2768     Inspector::addItem(this);
2769 #endif
2770 }
2771
2772 SparseBool::SparseBool(SparseBool const& src) : matrixBool(new BoolSparse_t(*src.matrixBool))
2773 {
2774     m_iDims = 2;
2775     m_iRows = const_cast<SparseBool*>(&src)->getRows();
2776     m_iCols = const_cast<SparseBool*>(&src)->getCols();
2777     m_iSize = m_iRows * m_iCols;
2778     m_piDims[0] = m_iRows;
2779     m_piDims[1] = m_iCols;
2780 #ifndef NDEBUG
2781     Inspector::addItem(this);
2782 #endif
2783 }
2784
2785 SparseBool::SparseBool(BoolSparse_t* src) : matrixBool(src)
2786 {
2787     m_iRows = src->rows();
2788     m_iCols = src->cols();
2789     m_iSize = m_iRows * m_iCols;
2790     m_iDims = 2;
2791     m_piDims[0] = m_iRows;
2792     m_piDims[1] = m_iCols;
2793 #ifndef NDEBUG
2794     Inspector::addItem(this);
2795 #endif
2796 }
2797
2798 SparseBool::SparseBool(int rows, int cols, int trues, int* inner, int* outer)
2799 {
2800     int* out = nullptr;
2801     int* in = nullptr;
2802
2803     matrixBool = new BoolSparse_t(rows, cols);
2804     matrixBool->reserve((int)trues);
2805     out = matrixBool->outerIndexPtr();
2806     in = matrixBool->innerIndexPtr();
2807
2808     //update outerIndexPtr
2809     memcpy(out, outer, sizeof(int) * (rows + 1));
2810     //update innerIndexPtr
2811     memcpy(in, inner, sizeof(int) * trues);
2812
2813     bool* data = matrixBool->valuePtr();
2814     for (int i = 0; i < trues; ++i)
2815     {
2816         data[i] = true;
2817     }
2818
2819     m_iCols = cols;
2820     m_iRows = rows;
2821     m_iSize = cols * rows;
2822     m_iDims = 2;
2823     m_piDims[0] = m_iRows;
2824     m_piDims[1] = m_iCols;
2825
2826     matrixBool->resizeNonZeros(trues);
2827 }
2828
2829 void SparseBool::create2(int rows, int cols, Bool SPARSE_CONST& src, Double SPARSE_CONST& idx)
2830 {
2831     int nnz = src.getSize();
2832     double* i = idx.get();
2833     double* j = i + idx.getRows();
2834     int* val = src.get();
2835
2836     typedef Eigen::Triplet<bool> T;
2837     std::vector<T> tripletList;
2838     tripletList.reserve((int)nnz);
2839
2840     for (int k = 0; k < nnz; ++k)
2841     {
2842         tripletList.push_back(T(static_cast<int>(i[k]) - 1, static_cast<int>(j[k]) - 1, val[k] == 1));
2843     }
2844
2845     matrixBool = new BoolSparse_t(rows, cols);
2846     matrixBool->setFromTriplets(tripletList.begin(), tripletList.end());
2847
2848     m_iRows = matrixBool->rows();
2849     m_iCols = matrixBool->cols();
2850     m_iSize = cols * rows;
2851     m_iDims = 2;
2852     m_piDims[0] = m_iRows;
2853     m_piDims[1] = m_iCols;
2854     finalize();
2855 }
2856
2857 SparseBool::~SparseBool()
2858 {
2859     delete matrixBool;
2860 #ifndef NDEBUG
2861     Inspector::removeItem(this);
2862 #endif
2863 }
2864
2865 bool SparseBool::toString(std::wostringstream& ostr) const
2866 {
2867     ostr << ::toString(*matrixBool, 0);
2868     return true;
2869 }
2870
2871 void SparseBool::whoAmI() SPARSE_CONST
2872 {
2873     std::cout << "types::SparseBool";
2874 }
2875
2876 SparseBool* SparseBool::clone(void) const
2877 {
2878     return new SparseBool(*this);
2879 }
2880
2881 bool SparseBool::resize(int _iNewRows, int _iNewCols)
2882 {
2883     if (_iNewRows <= getRows() && _iNewCols <= getCols())
2884     {
2885         //nothing to do: hence we do NOT fail
2886         return true;
2887     }
2888
2889     bool res = false;
2890     try
2891     {
2892         //item count
2893         size_t iNonZeros = nbTrue();
2894
2895         BoolSparse_t *newBool = new BoolSparse_t(_iNewRows, _iNewCols);
2896         newBool->reserve((int)iNonZeros);
2897
2898         //coords
2899         int* pRows = new int[iNonZeros * 2];
2900         outputRowCol(pRows);
2901         int* pCols = pRows + iNonZeros;
2902
2903         typedef Eigen::Triplet<bool> triplet;
2904         std::vector<triplet> tripletList;
2905
2906         for (size_t i = 0 ; i < iNonZeros ; i++)
2907         {
2908             tripletList.push_back(triplet((int)pRows[i] - 1, (int)pCols[i] - 1, true));
2909         }
2910
2911         newBool->setFromTriplets(tripletList.begin(), tripletList.end());
2912
2913         delete matrixBool;
2914         matrixBool = newBool;
2915         delete[] pRows;
2916
2917         m_iRows = _iNewRows;
2918         m_iCols = _iNewCols;
2919         m_iSize = _iNewRows * _iNewCols;
2920         m_piDims[0] = m_iRows;
2921         m_piDims[1] = m_iCols;
2922
2923         res = true;
2924
2925     }
2926     catch (...)
2927     {
2928         res = false;
2929     }
2930     return res;
2931 }
2932
2933 SparseBool* SparseBool::insert(typed_list* _pArgs, SparseBool* _pSource)
2934 {
2935     bool bNeedToResize  = false;
2936     int iDims           = (int)_pArgs->size();
2937     if (iDims > 2)
2938     {
2939         //sparse are only in 2 dims
2940         return NULL;
2941     }
2942
2943     typed_list pArg;
2944
2945     int piMaxDim[2];
2946     int piCountDim[2];
2947
2948     //on case of resize
2949     int iNewRows    = 0;
2950     int iNewCols    = 0;
2951
2952     //evaluate each argument and replace by appropriate value and compute the count of combinations
2953     int iSeqCount = checkIndexesArguments(this, _pArgs, &pArg, piMaxDim, piCountDim);
2954     if (iSeqCount == 0)
2955     {
2956         //free pArg content
2957         cleanIndexesArguments(_pArgs, &pArg);
2958         return this;
2959     }
2960
2961     if (iDims < 2)
2962     {
2963         //see as vector
2964         if (getRows() == 1 || getCols() == 1)
2965         {
2966             //vector or scalar
2967             if (getSize() < piMaxDim[0])
2968             {
2969                 bNeedToResize = true;
2970
2971                 //need to enlarge sparse dimensions
2972                 if (getCols() == 1 || getSize() == 0)
2973                 {
2974                     //column vector
2975                     iNewRows    = piMaxDim[0];
2976                     iNewCols    = 1;
2977                 }
2978                 else if (getRows() == 1)
2979                 {
2980                     //row vector
2981                     iNewRows    = 1;
2982                     iNewCols    = piMaxDim[0];
2983                 }
2984             }
2985         }
2986         else if (getSize() < piMaxDim[0])
2987         {
2988             //free pArg content
2989             cleanIndexesArguments(_pArgs, &pArg);
2990             //out of range
2991             return NULL;
2992         }
2993     }
2994     else
2995     {
2996         if (piMaxDim[0] > getRows() || piMaxDim[1] > getCols())
2997         {
2998             bNeedToResize = true;
2999             iNewRows = std::max(getRows(), piMaxDim[0]);
3000             iNewCols = std::max(getCols(), piMaxDim[1]);
3001         }
3002     }
3003
3004     //check number of insertion
3005     if (_pSource->isScalar() == false && _pSource->getSize() != iSeqCount)
3006     {
3007         //free pArg content
3008         cleanIndexesArguments(_pArgs, &pArg);
3009         return NULL;
3010     }
3011
3012     //now you are sure to be able to insert values
3013     if (bNeedToResize)
3014     {
3015         if (resize(iNewRows, iNewCols) == false)
3016         {
3017             //free pArg content
3018             cleanIndexesArguments(_pArgs, &pArg);
3019             return NULL;
3020         }
3021     }
3022
3023     if (iDims == 1)
3024     {
3025         double* pIdx = pArg[0]->getAs<Double>()->get();
3026         for (int i = 0 ; i < iSeqCount ; i++)
3027         {
3028             int iRow = static_cast<int>(pIdx[i] - 1) % getRows();
3029             int iCol = static_cast<int>(pIdx[i] - 1) / getRows();
3030
3031             if (_pSource->isScalar())
3032             {
3033                 set(iRow, iCol, _pSource->get(0, 0), false);
3034             }
3035             else
3036             {
3037                 int iRowOrig = i % _pSource->getRows();
3038                 int iColOrig = i / _pSource->getRows();
3039                 set(iRow, iCol, _pSource->get(iRowOrig, iColOrig), false);
3040             }
3041         }
3042     }
3043     else
3044     {
3045         double* pIdxRow = pArg[0]->getAs<Double>()->get();
3046         int iRowSize    = pArg[0]->getAs<Double>()->getSize();
3047         double* pIdxCol = pArg[1]->getAs<Double>()->get();
3048
3049         for (int i = 0 ; i < iSeqCount ; i++)
3050         {
3051             if (_pSource->isScalar())
3052             {
3053                 set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, _pSource->get(0, 0), false);
3054             }
3055             else
3056             {
3057                 int iRowOrig = i % _pSource->getRows();
3058                 int iColOrig = i / _pSource->getRows();
3059                 set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, _pSource->get(iRowOrig, iColOrig), false);
3060             }
3061         }
3062     }
3063
3064     finalize();
3065
3066     //free pArg content
3067     cleanIndexesArguments(_pArgs, &pArg);
3068
3069     return this;
3070 }
3071
3072 SparseBool* SparseBool::insert(typed_list* _pArgs, InternalType* _pSource)
3073 {
3074     bool bNeedToResize  = false;
3075     int iDims           = (int)_pArgs->size();
3076     if (iDims > 2)
3077     {
3078         //sparse are only in 2 dims
3079         return NULL;
3080     }
3081
3082     typed_list pArg;
3083
3084     int piMaxDim[2];
3085     int piCountDim[2];
3086
3087     //on case of resize
3088     int iNewRows    = 0;
3089     int iNewCols    = 0;
3090     Bool* pSource = _pSource->getAs<Bool>();
3091
3092     //evaluate each argument and replace by appropriate value and compute the count of combinations
3093     int iSeqCount = checkIndexesArguments(this, _pArgs, &pArg, piMaxDim, piCountDim);
3094     if (iSeqCount == 0)
3095     {
3096         //free pArg content
3097         cleanIndexesArguments(_pArgs, &pArg);
3098         return this;
3099     }
3100
3101     if (iDims < 2)
3102     {
3103         //see as vector
3104         if (getRows() == 1 || getCols() == 1)
3105         {
3106             //vector or scalar
3107             bNeedToResize = true;
3108             if (getSize() < piMaxDim[0])
3109             {
3110                 //need to enlarge sparse dimensions
3111                 if (getCols() == 1 || getSize() == 0)
3112                 {
3113                     //column vector
3114                     iNewRows    = piMaxDim[0];
3115                     iNewCols    = 1;
3116                 }
3117                 else if (getRows() == 1)
3118                 {
3119                     //row vector
3120                     iNewRows    = 1;
3121                     iNewCols    = piMaxDim[0];
3122                 }
3123             }
3124         }
3125         else if (getSize() < piMaxDim[0])
3126         {
3127             //free pArg content
3128             cleanIndexesArguments(_pArgs, &pArg);
3129             //out of range
3130             return NULL;
3131         }
3132     }
3133     else
3134     {
3135         if (piMaxDim[0] > getRows() || piMaxDim[1] > getCols())
3136         {
3137             bNeedToResize = true;
3138             iNewRows = std::max(getRows(), piMaxDim[0]);
3139             iNewCols = std::max(getCols(), piMaxDim[1]);
3140         }
3141     }
3142
3143     //check number of insertion
3144     if (pSource->isScalar() == false && pSource->getSize() != iSeqCount)
3145     {
3146         //free pArg content
3147         cleanIndexesArguments(_pArgs, &pArg);
3148         return NULL;
3149     }
3150
3151     //now you are sure to be able to insert values
3152     if (bNeedToResize)
3153     {
3154         if (resize(iNewRows, iNewCols) == false)
3155         {
3156             //free pArg content
3157             cleanIndexesArguments(_pArgs, &pArg);
3158             return NULL;
3159         }
3160     }
3161
3162     if (iDims == 1)
3163     {
3164         double* pIdx = pArg[0]->getAs<Double>()->get();
3165         for (int i = 0 ; i < iSeqCount ; i++)
3166         {
3167             int iRow = static_cast<int>(pIdx[i] - 1) % getRows();
3168             int iCol = static_cast<int>(pIdx[i] - 1) / getRows();
3169             if (pSource->isScalar())
3170             {
3171                 set(iRow, iCol, pSource->get(0) != 0, false);
3172             }
3173             else
3174             {
3175                 set(iRow, iCol, pSource->get(i) != 0, false);
3176             }
3177         }
3178     }
3179     else
3180     {
3181         double* pIdxRow = pArg[0]->getAs<Double>()->get();
3182         int iRowSize    = pArg[0]->getAs<Double>()->getSize();
3183         double* pIdxCol = pArg[1]->getAs<Double>()->get();
3184
3185         for (int i = 0 ; i < iSeqCount ; i++)
3186         {
3187             if (pSource->isScalar())
3188             {
3189                 set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, pSource->get(0) != 0, false);
3190             }
3191             else
3192             {
3193                 int iRowOrig = i % pSource->getRows();
3194                 int iColOrig = i / pSource->getRows();
3195
3196                 set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, pSource->get(iRowOrig, iColOrig) != 0, false);
3197             }
3198         }
3199     }
3200
3201     finalize();
3202
3203     //free pArg content
3204     cleanIndexesArguments(_pArgs, &pArg);
3205     return this;
3206 }
3207
3208 SparseBool* SparseBool::remove(typed_list* _pArgs)
3209 {
3210     SparseBool* pOut = NULL;
3211     int iDims = (int)_pArgs->size();
3212     if (iDims > 2)
3213     {
3214         //sparse are only in 2 dims
3215         return NULL;
3216     }
3217
3218     typed_list pArg;
3219
3220     int piMaxDim[2];
3221     int piCountDim[2];
3222
3223     //evaluate each argument and replace by appropriate value and compute the count of combinations
3224     int iSeqCount = checkIndexesArguments(this, _pArgs, &pArg, piMaxDim, piCountDim);
3225     if (iSeqCount == 0)
3226     {
3227         //free pArg content
3228         cleanIndexesArguments(_pArgs, &pArg);
3229         return this;
3230     }
3231
3232     bool* pbFull = new bool[iDims];
3233     //coord must represent all values on a dimension
3234     for (int i = 0 ; i < iDims ; i++)
3235     {
3236         pbFull[i]       = false;
3237         int iDimToCheck = getVarMaxDim(i, iDims);
3238         int iIndexSize  = pArg[i]->getAs<GenericType>()->getSize();
3239
3240         //we can have index more than once
3241         if (iIndexSize >= iDimToCheck)
3242         {
3243             //size is good, now check datas
3244             double* pIndexes = getDoubleArrayFromDouble(pArg[i]);
3245             for (int j = 0 ; j < iDimToCheck ; j++)
3246             {
3247                 bool bFind = false;
3248                 for (int k = 0 ; k < iIndexSize ; k++)
3249                 {
3250                     if ((int)pIndexes[k] == j + 1)
3251                     {
3252                         bFind = true;
3253                         break;
3254                     }
3255                 }
3256                 pbFull[i]  = bFind;
3257             }
3258         }
3259     }
3260
3261     //only one dims can be not full/entire
3262     bool bNotEntire = false;
3263     int iNotEntire  = 0;
3264     bool bTooMuchNotEntire = false;
3265     for (int i = 0 ; i < iDims ; i++)
3266     {
3267         if (pbFull[i] == false)
3268         {
3269             if (bNotEntire == false)
3270             {
3271                 bNotEntire = true;
3272                 iNotEntire = i;
3273             }
3274             else
3275             {
3276                 bTooMuchNotEntire = true;
3277                 break;
3278             }
3279         }
3280     }
3281
3282     if (bTooMuchNotEntire == true)
3283     {
3284         //free pArg content
3285         cleanIndexesArguments(_pArgs, &pArg);
3286         return NULL;
3287     }
3288
3289     delete[] pbFull;
3290
3291     //find index to keep
3292     int iNotEntireSize          = pArg[iNotEntire]->getAs<GenericType>()->getSize();
3293     double* piNotEntireIndex    = getDoubleArrayFromDouble(pArg[iNotEntire]);
3294     int iKeepSize               = getVarMaxDim(iNotEntire, iDims);
3295     bool* pbKeep                = new bool[iKeepSize];
3296
3297     //fill pbKeep with true value
3298     for (int i = 0 ; i < iKeepSize ; i++)
3299     {
3300         pbKeep[i] = true;
3301     }
3302
3303     for (int i = 0 ; i < iNotEntireSize ; i++)
3304     {
3305         int idx = (int)piNotEntireIndex[i] - 1;
3306
3307         //don't care of value out of bounds
3308         if (idx < iKeepSize)
3309         {
3310             pbKeep[idx] = false;
3311         }
3312     }
3313
3314     int iNewDimSize = 0;
3315     for (int i = 0 ; i < iKeepSize ; i++)
3316     {
3317         if (pbKeep[i] == true)
3318         {
3319             iNewDimSize++;
3320         }
3321     }
3322     delete[] pbKeep;
3323
3324     int* piNewDims = new int[iDims];
3325     for (int i = 0 ; i < iDims ; i++)
3326     {
3327         if (i == iNotEntire)
3328         {
3329             piNewDims[i] = iNewDimSize;
3330         }
3331         else
3332         {
3333             piNewDims[i] = getVarMaxDim(i, iDims);
3334         }
3335     }
3336
3337     //remove last dimension if are == 1
3338     int iOrigDims = iDims;
3339     for (int i = (iDims - 1) ; i >= 2 ; i--)
3340     {
3341         if (piNewDims[i] == 1)
3342         {
3343             iDims--;
3344         }
3345         else
3346         {
3347             break;
3348         }
3349     }
3350
3351     if (iDims == 1)
3352     {
3353         if (iNewDimSize == 0)
3354         {
3355             //free pArg content
3356             cleanIndexesArguments(_pArgs, &pArg);
3357             return new SparseBool(0, 0);
3358         }
3359         else
3360         {
3361             //two cases, depends of original matrix/vector
3362             if ((*_pArgs)[0]->isColon() == false && m_iDims == 2 && m_piDims[0] == 1 && m_piDims[1] != 1)
3363             {
3364                 //special case for row vector
3365                 pOut = new SparseBool(1, iNewDimSize);
3366                 //in this case we have to care of 2nd dimension
3367                 //iNotEntire = 1;
3368             }
3369             else
3370             {
3371                 pOut = new SparseBool(iNewDimSize, 1);
3372             }
3373         }
3374     }
3375     else
3376     {
3377         pOut = new SparseBool(piNewDims[0], piNewDims[0]);
3378     }
3379
3380     delete[] piNewDims;
3381     //find a way to copy existing data to new variable ...
3382     int iNewPos = 0;
3383     int* piIndexes = new int[iOrigDims];
3384     int* piViewDims = new int[iOrigDims];
3385     for (int i = 0 ; i < iOrigDims ; i++)
3386     {
3387         piViewDims[i] = getVarMaxDim(i, iOrigDims);
3388     }
3389
3390     for (int i = 0 ; i < getSize() ; i++)
3391     {
3392         bool bByPass = false;
3393         getIndexesWithDims(i, piIndexes, piViewDims, iOrigDims);
3394
3395         //check if piIndexes use removed indexes
3396         for (int j = 0 ; j < iNotEntireSize ; j++)
3397         {
3398             if ((piNotEntireIndex[j] - 1) == piIndexes[iNotEntire])
3399             {
3400                 //by pass this value
3401                 bByPass = true;
3402                 break;
3403             }
3404         }
3405
3406         if (bByPass == false)
3407         {
3408             //compute new index
3409             pOut->set(iNewPos, get(i));
3410             iNewPos++;
3411         }
3412     }
3413
3414     //free allocated data
3415     for (int i = 0 ; i < iDims ; i++)
3416     {
3417         if (pArg[i] != (*_pArgs)[i])
3418         {
3419             delete pArg[i];
3420         }
3421     }
3422
3423     delete[] piIndexes;
3424     delete[] piViewDims;
3425
3426     //free pArg content
3427     cleanIndexesArguments(_pArgs, &pArg);
3428
3429     return pOut;
3430 }
3431
3432 bool SparseBool::append(int r, int c, SparseBool SPARSE_CONST* src)
3433 {
3434     doAppend(*src, r, c, *matrixBool);
3435     finalize();
3436     return true;
3437 }
3438
3439 InternalType* SparseBool::insertNew(typed_list* _pArgs, InternalType* _pSource)
3440 {
3441     typed_list pArg;
3442     InternalType *pOut  = NULL;
3443     SparseBool* pSource = _pSource->getAs<SparseBool>();
3444
3445     int iDims           = (int)_pArgs->size();
3446     int* piMaxDim       = new int[iDims];
3447     int* piCountDim     = new int[iDims];
3448     bool bUndefine      = false;
3449
3450     //evaluate each argument and replace by appropriate value and compute the count of combinations
3451     int iSeqCount = checkIndexesArguments(NULL, _pArgs, &pArg, piMaxDim, piCountDim);
3452     if (iSeqCount == 0)
3453     {
3454         //free pArg content
3455         cleanIndexesArguments(_pArgs, &pArg);
3456         return createEmptyDouble();
3457     }
3458
3459     if (iSeqCount < 0)
3460     {
3461         iSeqCount = -iSeqCount;
3462         bUndefine = true;
3463     }
3464
3465     if (bUndefine)
3466     {
3467         //manage : and $ in creation by insertion
3468         int iSource         = 0;
3469         int *piSourceDims   = pSource->getDimsArray();
3470
3471         for (int i = 0 ; i < iDims ; i++)
3472         {
3473             if (pArg[i] == NULL)
3474             {
3475                 //undefine value
3476                 if (pSource->isScalar())
3477                 {
3478                     piMaxDim[i]     = 1;
3479                     piCountDim[i]   = 1;
3480                 }
3481                 else
3482                 {
3483                     piMaxDim[i]     = piSourceDims[iSource];
3484                     piCountDim[i]   = piSourceDims[iSource];
3485                 }
3486                 iSource++;
3487                 //replace pArg value by the new one
3488                 pArg[i] = createDoubleVector(piMaxDim[i]);
3489             }
3490             //else
3491             //{
3492             //    piMaxDim[i] = piCountDim[i];
3493             //}
3494         }
3495     }
3496
3497     //remove last dimension at size 1
3498     //remove last dimension if are == 1
3499     for (int i = (iDims - 1) ; i >= 2 ; i--)
3500     {
3501         if (piMaxDim[i] == 1)
3502         {
3503             iDims--;
3504             pArg.pop_back();
3505         }
3506         else
3507         {
3508             break;
3509         }
3510     }
3511
3512     if (checkArgValidity(pArg) == false)
3513     {
3514         //free pArg content
3515         cleanIndexesArguments(_pArgs, &pArg);
3516         //contain bad index, like <= 0, ...
3517         return NULL;
3518     }
3519
3520     if (iDims == 1)
3521     {
3522         if (pSource->getCols() == 1)
3523         {
3524             pOut = new SparseBool(piCountDim[0], 1);
3525         }
3526         else
3527         {
3528             //rows == 1
3529             pOut = new SparseBool(1, piCountDim[0]);
3530         }
3531     }
3532     else
3533     {
3534         pOut = new SparseBool(piMaxDim[0], piMaxDim[1]);
3535     }
3536
3537     //fill with null item
3538     SparseBool* pSpOut = pOut->getAs<SparseBool>();
3539
3540     //insert values in new matrix
3541     InternalType* pOut2 = pSpOut->insert(&pArg, pSource);
3542     if (pOut != pOut2)
3543     {
3544         delete pOut;
3545     }
3546
3547     //free pArg content
3548     cleanIndexesArguments(_pArgs, &pArg);
3549
3550     return pOut2;
3551 }
3552
3553 SparseBool* SparseBool::extract(int nbCoords, int SPARSE_CONST* coords, int SPARSE_CONST* maxCoords, int SPARSE_CONST* resSize, bool asVector) SPARSE_CONST
3554 {
3555     if ( (asVector && maxCoords[0] > getSize()) ||
3556     (asVector == false && maxCoords[0] > getRows()) ||
3557     (asVector == false && maxCoords[1] > getCols()))
3558     {
3559         return 0;
3560     }
3561
3562     SparseBool * pSp (0);
3563     if (asVector)
3564     {
3565         pSp = (getRows() == 1) ?  new SparseBool(1, resSize[0]) : new SparseBool(resSize[0], 1);
3566         mycopy_n(makeMatrixIterator<bool>(*this,  Coords<true>(coords, getRows())), nbCoords
3567         , makeMatrixIterator<bool>(*(pSp->matrixBool), RowWiseFullIterator(pSp->getRows(), pSp->getCols())));
3568     }
3569     else
3570     {
3571         pSp = new SparseBool(resSize[0], resSize[1]);
3572         mycopy_n(makeMatrixIterator<bool>(*this,  Coords<false>(coords, getRows())), nbCoords
3573         , makeMatrixIterator<bool>(*(pSp->matrixBool), RowWiseFullIterator(pSp->getRows(), pSp->getCols())));
3574
3575     }
3576     return pSp;
3577 }
3578
3579 /*
3580 * create a new SparseBool of dims according to resSize and fill it from currentSparseBool (along coords)
3581 */
3582 InternalType* SparseBool::extract(typed_list* _pArgs)
3583 {
3584     SparseBool* pOut    = NULL;
3585     int iDims           = (int)_pArgs->size();
3586     typed_list pArg;
3587
3588     int* piMaxDim       = new int[iDims];
3589     int* piCountDim     = new int[iDims];
3590
3591     //evaluate each argument and replace by appropriate value and compute the count of combinations
3592     int iSeqCount = checkIndexesArguments(this, _pArgs, &pArg, piMaxDim, piCountDim);
3593     if (iSeqCount == 0)
3594     {
3595         //free pArg content
3596         cleanIndexesArguments(_pArgs, &pArg);
3597         if (_pArgs->size() == 0)
3598         {
3599             delete[] piMaxDim;
3600             delete[] piCountDim;
3601             //a()
3602             return this;
3603         }
3604         else
3605         {
3606             delete[] piMaxDim;
3607             delete[] piCountDim;
3608             //a([])
3609             return Double::Empty();
3610         }
3611     }
3612
3613     if (iDims < 2)
3614     {
3615         // Check that we stay inside the input size.
3616         if (piMaxDim[0] <= getSize())
3617         {
3618             int iNewRows = 0;
3619             int iNewCols = 0;
3620
3621             if (getRows() == 1 && getCols() != 1 && (*_pArgs)[0]->isColon() == false)
3622             {
3623                 //special case for row vector
3624                 iNewRows = 1;
3625                 iNewCols = piCountDim[0];
3626             }
3627             else
3628             {
3629                 iNewRows = piCountDim[0];
3630                 iNewCols = 1;
3631             }
3632
3633             pOut = new SparseBool(iNewRows, iNewCols);
3634             double* pIdx = pArg[0]->getAs<Double>()->get();
3635             // Write in output all elements extract from input.
3636             for (int i = 0 ; i < iSeqCount ; i++)
3637             {
3638                 if (pIdx[i] < 1)
3639                 {
3640                     delete pOut;
3641                     pOut = NULL;
3642                     break;
3643                 }
3644                 int iRowRead = static_cast<int>(pIdx[i] - 1) % getRows();
3645                 int iColRead = static_cast<int>(pIdx[i] - 1) / getRows();
3646
3647                 int iRowWrite = static_cast<int>(i) % iNewRows;
3648                 int iColWrite = static_cast<int>(i) / iNewRows;
3649
3650                 bool bValue = get(iRowRead, iColRead);
3651                 if (bValue)
3652                 {
3653                     //only non zero values
3654                     pOut->set(iRowWrite, iColWrite, true, false);
3655                 }
3656             }
3657         }
3658         else
3659         {
3660             delete[] piMaxDim;
3661             delete[] piCountDim;
3662             //free pArg content
3663             cleanIndexesArguments(_pArgs, &pArg);
3664             return NULL;
3665         }
3666     }
3667     else
3668     {
3669         // Check that we stay inside the input size.
3670         if (piMaxDim[0] <= getRows() && piMaxDim[1] <= getCols())
3671         {
3672             double* pIdxRow = pArg[0]->getAs<Double>()->get();
3673             double* pIdxCol = pArg[1]->getAs<Double>()->get();
3674
3675             int iNewRows = pArg[0]->getAs<Double>()->getSize();
3676             int iNewCols = pArg[1]->getAs<Double>()->getSize();
3677
3678             pOut = new SparseBool(iNewRows, iNewCols);
3679
3680             int iPos = 0;
3681             // Write in output all elements extract from input.
3682             for (int iRow = 0 ; iRow < iNewRows ; iRow++)
3683             {
3684                 for (int iCol = 0 ; iCol < iNewCols ; iCol++)
3685                 {
3686                     if ((pIdxRow[iRow] < 1) || (pIdxCol[iCol] < 1))
3687                     {
3688                         delete pOut;
3689                         pOut = NULL;
3690                         break;
3691                     }
3692                     bool bValue = get((int)pIdxRow[iRow] - 1, (int)pIdxCol[iCol] - 1);
3693                     if (bValue)
3694                     {
3695                         //only non zero values
3696                         pOut->set(iRow, iCol, true, false);
3697                     }
3698                     iPos++;
3699                 }
3700             }
3701         }
3702         else
3703         {
3704             delete[] piMaxDim;
3705             delete[] piCountDim;
3706             //free pArg content
3707             cleanIndexesArguments(_pArgs, &pArg);
3708             return NULL;
3709         }
3710     }
3711
3712     finalize();
3713
3714     delete[] piMaxDim;
3715     delete[] piCountDim;
3716     //free pArg content
3717     cleanIndexesArguments(_pArgs, &pArg);
3718
3719     return pOut;
3720 }
3721
3722 bool SparseBool::invoke(typed_list & in, optional_list &/*opt*/, int /*_iRetCount*/, typed_list & out, const ast::Exp & e)
3723 {
3724     if (in.size() == 0)
3725     {
3726         out.push_back(this);
3727     }
3728     else
3729     {
3730         InternalType * _out = extract(&in);
3731         if (!_out)
3732         {
3733             std::wostringstream os;
3734             os << _W("Invalid index.\n");
3735             throw ast::InternalError(os.str(), 999, e.getLocation());
3736         }
3737         out.push_back(_out);
3738     }
3739
3740     return true;
3741 }
3742
3743 bool SparseBool::isInvokable() const
3744 {
3745     return true;
3746 }
3747
3748 bool SparseBool::hasInvokeOption() const
3749 {
3750     return false;
3751 }
3752
3753 int SparseBool::getInvokeNbIn()
3754 {
3755     return -1;
3756 }
3757
3758 int SparseBool::getInvokeNbOut()
3759 {
3760     return 1;
3761 }
3762
3763 std::size_t SparseBool::nbTrue() const
3764 {
3765     return  matrixBool->nonZeros() ;
3766 }
3767 std::size_t SparseBool::nbTrue(std::size_t r) const
3768 {
3769     int* piIndex = matrixBool->outerIndexPtr();
3770     return piIndex[r + 1] - piIndex[r];
3771 }
3772
3773
3774 void SparseBool::setTrue(bool finalize)
3775 {
3776     int rows = getRows();
3777     int cols = getCols();
3778
3779     typedef Eigen::Triplet<bool> triplet;
3780     std::vector<triplet> tripletList;
3781
3782     for (int i = 0; i < rows; ++i)
3783     {
3784         for (int j = 0; j < cols; ++j)
3785         {
3786             tripletList.push_back(triplet(i, j, true));
3787         }
3788     }
3789
3790     matrixBool->setFromTriplets(tripletList.begin(), tripletList.end());
3791
3792     if (finalize)
3793     {
3794         matrixBool->finalize();
3795     }
3796 }
3797
3798 void SparseBool::setFalse(bool finalize)
3799 {
3800     int rows = getRows();
3801     int cols = getCols();
3802
3803     typedef Eigen::Triplet<bool> triplet;
3804     std::vector<triplet> tripletList;
3805
3806     for (int i = 0; i < rows; ++i)
3807     {
3808         for (int j = 0; j < cols; ++j)
3809         {
3810             tripletList.push_back(triplet(i, j, false));
3811         }
3812     }
3813
3814     matrixBool->setFromTriplets(tripletList.begin(), tripletList.end());
3815
3816     if (finalize)
3817     {
3818         matrixBool->finalize();
3819     }
3820 }
3821
3822 int* SparseBool::getNbItemByRow(int* _piNbItemByRows)
3823 {
3824     int* piNbItemByRows = new int[getRows() + 1];
3825     mycopy_n(matrixBool->outerIndexPtr(), getRows() + 1, piNbItemByRows);
3826
3827     for (int i = 0 ; i < getRows() ; i++)
3828     {
3829         _piNbItemByRows[i] = piNbItemByRows[i + 1] - piNbItemByRows[i];
3830     }
3831
3832     delete[] piNbItemByRows;
3833     return _piNbItemByRows;
3834 }
3835
3836 int* SparseBool::getColPos(int* _piColPos)
3837 {
3838     mycopy_n(matrixBool->innerIndexPtr(), nbTrue(), _piColPos);
3839     for (int i = 0; i < nbTrue(); i++)
3840     {
3841         _piColPos[i]++;
3842     }
3843
3844     return _piColPos;
3845 }
3846
3847 int* SparseBool::outputRowCol(int* out)const
3848 {
3849     return sparseTransform(*matrixBool, sparseTransform(*matrixBool, out, GetRow<BoolSparse_t>()), GetCol<BoolSparse_t>());
3850 }
3851
3852 int* SparseBool::getInnerPtr(int* count)
3853 {
3854     *count = matrixBool->innerSize();
3855     return matrixBool->innerIndexPtr();
3856 }
3857
3858 int* SparseBool::getOuterPtr(int* count)
3859 {
3860     *count = matrixBool->outerSize();
3861     return matrixBool->outerIndexPtr();
3862 }
3863
3864
3865 bool SparseBool::operator==(const InternalType& it) SPARSE_CONST
3866 {
3867     SparseBool* otherSparse = const_cast<SparseBool*>(dynamic_cast<SparseBool const*>(&it));/* types::GenericType is not const-correct :( */
3868     return (otherSparse
3869     && (otherSparse->getRows() == getRows())
3870     && (otherSparse->getCols() == getCols())
3871     && equal(*matrixBool, *otherSparse->matrixBool));
3872 }
3873
3874 bool SparseBool::operator!=(const InternalType& it) SPARSE_CONST
3875 {
3876     return !(*this == it);
3877 }
3878
3879 void SparseBool::finalize()
3880 {
3881     matrixBool->prune(&keepForSparse<bool>);
3882     matrixBool->finalize();
3883 }
3884
3885 bool SparseBool::get(int r, int c) SPARSE_CONST
3886 {
3887     return matrixBool->coeff(r, c);
3888 }
3889
3890 bool SparseBool::set(int _iRows, int _iCols, bool _bVal, bool _bFinalize) SPARSE_CONST
3891 {
3892     matrixBool->coeffRef(_iRows, _iCols) = _bVal;
3893
3894     if (_bFinalize)
3895     {
3896         finalize();
3897     }
3898
3899     return true;
3900 }
3901
3902 void SparseBool::fill(Bool& dest, int r, int c) SPARSE_CONST
3903 {
3904     mycopy_n(makeMatrixIterator<bool >(*matrixBool, RowWiseFullIterator(getRows(), getCols())), getSize()
3905     , makeMatrixIterator<bool >(dest, RowWiseFullIterator(dest.getRows(), dest.getCols(), r, c)));
3906 }
3907
3908 Sparse* SparseBool::newOnes() const
3909 {
3910     return new Sparse(new types::Sparse::RealSparse_t(matrixBool->cast<double>()), 0);
3911 }
3912
3913 SparseBool* SparseBool::newNotEqualTo(SparseBool const&o) const
3914 {
3915     return cwiseOp<std::not_equal_to>(*this, o);
3916 }
3917
3918 SparseBool* SparseBool::newEqualTo(SparseBool& o)
3919 {
3920     int rowL = getRows();
3921     int colL = getCols();
3922
3923     int rowR = o.getRows();
3924     int colR = o.getCols();
3925     int row = std::max(rowL, rowR);
3926     int col = std::max(colL, colR);
3927
3928     //create a boolean sparse matrix with dims of sparses
3929     types::SparseBool* ret = new types::SparseBool(row, col);
3930
3931     if (isScalar() && o.isScalar())
3932     {
3933         bool l = get(0, 0);
3934         bool r = o.get(0, 0);
3935         ret->set(0, 0, l == r, false);
3936     }
3937     else if (isScalar())
3938     {
3939         int nnzR = static_cast<int>(o.nbTrue());
3940         std::vector<int> rowcolR(nnzR * 2, 0);
3941         o.outputRowCol(rowcolR.data());
3942
3943         //compare all items of R with R[0]
3944         bool l = get(0, 0);
3945         for (int i = 0; i < nnzR; ++i)
3946         {
3947             bool r = o.get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
3948             ret->set(rowcolR[i] - 1, rowcolR[i + nnzR] - 1, l == r, false);
3949         }
3950     }
3951     else if (o.isScalar())
3952     {
3953         int nnzL = static_cast<int>(nbTrue());
3954         std::vector<int> rowcolL(nnzL * 2, 0);
3955         outputRowCol(rowcolL.data());
3956
3957         bool r = get(0, 0);
3958         for (int i = 0; i < nnzL; ++i)
3959         {
3960             bool l = get(rowcolL[i] - 1, rowcolL[i + nnzL] - 1);
3961             ret->set(rowcolL[i] - 1, rowcolL[i + nnzL] - 1, l == r, false);
3962         }
3963     }
3964     else
3965     {
3966         int nnzR = static_cast<int>(o.nbTrue());
3967         std::vector<int> rowcolR(nnzR * 2, 0);
3968         o.outputRowCol(rowcolR.data());
3969         int nnzL = static_cast<int>(nbTrue());
3970         std::vector<int> rowcolL(nnzL * 2, 0);
3971         outputRowCol(rowcolL.data());
3972         //set all values to %t
3973         ret->setTrue(false);
3974         //set %f in each pL values
3975         for (int i = 0; i < nnzL; ++i)
3976         {
3977             ret->set(rowcolL[i] - 1, rowcolL[i + nnzL] - 1, false, false);
3978         }
3979         ret->finalize();
3980
3981         //set _pR[i] == _pL[i] for each _pR values
3982         for (int i = 0; i < nnzR; ++i)
3983         {
3984             //get l and r following non zeros value of R
3985             bool l = get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
3986             bool r = o.get(rowcolR[i] - 1, rowcolR[i + nnzR] - 1);
3987             //set value following non zeros value of R
3988             ret->set(rowcolR[i] - 1, rowcolR[i + nnzR] - 1, l == r, false);
3989         }
3990     }
3991
3992     ret->finalize();
3993     return ret;
3994 }
3995
3996 SparseBool* SparseBool::newLogicalOr(SparseBool const&o) const
3997 {
3998     return cwiseOp<std::logical_or>(*this, o);
3999 }
4000
4001 SparseBool* SparseBool::newLogicalAnd(SparseBool const&o) const
4002 {
4003     return cwiseOp<std::logical_and>(*this, o);
4004 }
4005
4006 bool SparseBool::reshape(int* _piDims, int _iDims)
4007 {
4008     bool bOk = false;
4009     int iCols = 1;
4010
4011     if (_iDims == 2)
4012     {
4013         iCols = _piDims[1];
4014     }
4015
4016     if (_iDims <= 2)
4017     {
4018         bOk = reshape(_piDims[0], iCols);
4019     }
4020
4021     return bOk;
4022 }
4023
4024 bool SparseBool::reshape(int _iNewRows, int _iNewCols)
4025 {
4026     if (_iNewRows * _iNewCols != getRows() * getCols())
4027     {
4028         return false;
4029     }
4030
4031     bool res = false;
4032     try
4033     {
4034         //item count
4035         size_t iNonZeros = matrixBool->nonZeros();
4036         BoolSparse_t *newBool = new BoolSparse_t(_iNewRows, _iNewCols);
4037         newBool->reserve((int)iNonZeros);
4038
4039         //coords
4040         int* pRows = new int[iNonZeros * 2];
4041         outputRowCol(pRows);
4042         int* pCols = pRows + iNonZeros;
4043
4044         typedef Eigen::Triplet<bool> triplet;
4045         std::vector<triplet> tripletList;
4046
4047         for (size_t i = 0 ; i < iNonZeros ; i++)
4048         {
4049             int iCurrentPos = ((int)pCols[i] - 1) * getRows() + ((int)pRows[i] - 1);
4050             tripletList.push_back(triplet((int)(iCurrentPos % _iNewRows), (int)(iCurrentPos / _iNewRows), true));
4051         }
4052
4053         newBool->setFromTriplets(tripletList.begin(), tripletList.end());
4054
4055         delete matrixBool;
4056         matrixBool = newBool;
4057         delete[] pRows;
4058
4059         m_iRows = _iNewRows;
4060         m_iCols = _iNewCols;
4061         m_iSize = _iNewRows * _iNewCols;
4062
4063         m_iDims = 2;
4064         m_piDims[0] = m_iRows;
4065         m_piDims[1] = m_iCols;
4066
4067         finalize();
4068
4069         res = true;
4070     }
4071     catch (...)
4072     {
4073         res = false;
4074     }
4075     return res;
4076 }
4077
4078 bool SparseBool::transpose(InternalType *& out)
4079 {
4080     out = new SparseBool(new BoolSparse_t(matrixBool->transpose()));
4081     return true;
4082 }
4083
4084 template<typename T>
4085 void neg(const int r, const int c, const T * const in, Eigen::SparseMatrix<bool, 1> * const out)
4086 {
4087     for (int i = 0; i < r; i++)
4088     {
4089         for (int j = 0; j < c; j++)
4090         {
4091             out->coeffRef(i, j) = !in->coeff(i, j);
4092         }
4093     }
4094
4095     out->prune(&keepForSparse<bool>);
4096     out->finalize();
4097 }
4098
4099
4100 }