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