flex fixed about bool token
[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 const&o) const
2239 {
2240     return cwiseOp<std::less>(*this, o);
2241 }
2242
2243 SparseBool* Sparse::newGreaterThan(Sparse const&o) const
2244 {
2245     return cwiseOp<std::greater>(*this, o);
2246 }
2247
2248 SparseBool* Sparse::newNotEqualTo(Sparse const&o) const
2249 {
2250     return cwiseOp<std::not_equal_to>(*this, o);
2251 }
2252
2253 SparseBool* Sparse::newLessOrEqual(Sparse const&o) const
2254 {
2255     return cwiseOp<std::less_equal>(*this, o);
2256 }
2257
2258 SparseBool* Sparse::newGreaterOrEqual(Sparse const&o) const
2259 {
2260     return cwiseOp<std::greater_equal>(*this, o);
2261 }
2262
2263 SparseBool* Sparse::newEqualTo(Sparse const&o) const
2264 {
2265     return cwiseOp<std::equal_to>(*this, o);
2266 }
2267
2268 bool Sparse::reshape(int* _piDims, int _iDims)
2269 {
2270     bool bOk = false;
2271     int iCols = 1;
2272
2273     if (_iDims == 2)
2274     {
2275         iCols = _piDims[1];
2276     }
2277
2278     if (_iDims <= 2)
2279     {
2280         bOk = reshape(_piDims[0], iCols);
2281     }
2282
2283     return bOk;
2284 }
2285
2286 bool Sparse::reshape(int _iNewRows, int _iNewCols)
2287 {
2288     if (_iNewRows * _iNewCols != getRows() * getCols())
2289     {
2290         return false;
2291     }
2292
2293     bool res = false;
2294     try
2295     {
2296         if (matrixReal)
2297         {
2298             //item count
2299             size_t iNonZeros = nonZeros();
2300             RealSparse_t *newReal = new RealSparse_t(_iNewRows, _iNewCols);
2301             newReal->reserve((int)iNonZeros);
2302
2303             //coords
2304             int* pRows = new int[iNonZeros * 2];
2305             outputRowCol(pRows);
2306             int* pCols = pRows + iNonZeros;
2307
2308             //values
2309             double* pNonZeroR = new double[iNonZeros];
2310             double* pNonZeroI = new double[iNonZeros];
2311             outputValues(pNonZeroR, pNonZeroI);
2312
2313             typedef Eigen::Triplet<double> triplet;
2314             std::vector<triplet> tripletList;
2315
2316             for (size_t i = 0 ; i < iNonZeros ; i++)
2317             {
2318                 int iCurrentPos = ((int)pCols[i] - 1) * getRows() + ((int)pRows[i] - 1);
2319                 tripletList.push_back(triplet((int)(iCurrentPos % _iNewRows), (int)(iCurrentPos / _iNewRows), pNonZeroR[i]));
2320             }
2321
2322             newReal->setFromTriplets(tripletList.begin(), tripletList.end());
2323
2324             delete matrixReal;
2325             matrixReal = newReal;
2326             delete[] pRows;
2327             delete[] pNonZeroR;
2328             delete[] pNonZeroI;
2329         }
2330         else
2331         {
2332             //item count
2333             size_t iNonZeros = nonZeros();
2334             CplxSparse_t *newCplx = new CplxSparse_t(_iNewRows, _iNewCols);
2335             newCplx->reserve((int)iNonZeros);
2336
2337             //coords
2338             int* pRows = new int[iNonZeros * 2];
2339             outputRowCol(pRows);
2340             int* pCols = pRows + iNonZeros;
2341
2342             //values
2343             double* pNonZeroR = new double[iNonZeros];
2344             double* pNonZeroI = new double[iNonZeros];
2345             outputValues(pNonZeroR, pNonZeroI);
2346
2347             typedef Eigen::Triplet<std::complex<double> > triplet;
2348             std::vector<triplet> tripletList;
2349
2350             for (size_t i = 0 ; i < iNonZeros ; i++)
2351             {
2352                 int iCurrentPos = ((int)pCols[i] - 1) * getRows() + ((int)pRows[i] - 1);
2353                 tripletList.push_back(triplet((int)(iCurrentPos % _iNewRows), (int)(iCurrentPos / _iNewRows), std::complex<double>(pNonZeroR[i], pNonZeroI[i])));
2354             }
2355
2356             newCplx->setFromTriplets(tripletList.begin(), tripletList.end());
2357
2358             delete matrixCplx;
2359             matrixCplx = newCplx;
2360             delete[] pRows;
2361             delete[] pNonZeroR;
2362             delete[] pNonZeroI;
2363         }
2364
2365         m_iRows = _iNewRows;
2366         m_iCols = _iNewCols;
2367         m_iSize = _iNewRows * _iNewCols;
2368
2369         m_iDims = 2;
2370         m_piDims[0] = m_iRows;
2371         m_piDims[1] = m_iCols;
2372
2373         finalize();
2374
2375         res = true;
2376     }
2377     catch (...)
2378     {
2379         res = false;
2380     }
2381     return res;
2382 }
2383
2384 //    SparseBool* SparseBool::new
2385
2386 SparseBool::SparseBool(Bool SPARSE_CONST& src)
2387 {
2388     //compute idx
2389     int size = src.getSize();
2390     int row = src.getRows();
2391     Double* idx = new Double(src.getSize(), 2);
2392     double* p = idx->get();
2393     for (int i = 0; i < size; ++i)
2394     {
2395         p[i] = (double)(i % row) + 1;
2396         p[i + size] = (double)(i / row) + 1;
2397     }
2398     create2(src.getRows(), src.getCols(), src, *idx);
2399     idx->killMe();
2400 #ifndef NDEBUG
2401     Inspector::addItem(this);
2402 #endif
2403 }
2404 /* @param src : Bool matrix to copy into a new sparse matrix
2405 @param idx : Double matrix to use as indexes to get values from the src
2406 **/
2407 SparseBool::SparseBool(Bool SPARSE_CONST& src, Double SPARSE_CONST& idx)
2408 {
2409     int idxrow = idx.getRows();
2410     int rows = static_cast<int>(*std::max_element(idx.get(), idx.get() + idxrow));
2411     int cols = static_cast<int>(*std::max_element(idx.get() + idxrow, idx.get() + idxrow * 2));
2412     create2(rows, cols, src, idx);
2413 #ifndef NDEBUG
2414     Inspector::addItem(this);
2415 #endif
2416 }
2417
2418 /* @param src : Bool matrix to copy into a new sparse matrix
2419 @param idx : Double matrix to use as indexes to get values from the src
2420 @param dims : Double matrix containing the dimensions of the new matrix
2421 **/
2422 SparseBool::SparseBool(Bool SPARSE_CONST& src, Double SPARSE_CONST& idx, Double SPARSE_CONST& dims)
2423 {
2424     create2(static_cast<int>(dims.get(0)), static_cast<int>(dims.get(1)), src, idx);
2425 #ifndef NDEBUG
2426     Inspector::addItem(this);
2427 #endif
2428 }
2429
2430 SparseBool::SparseBool(int _iRows, int _iCols) : matrixBool(new BoolSparse_t(_iRows, _iCols))
2431 {
2432     m_iRows = _iRows;
2433     m_iCols = _iCols;
2434     m_iSize = _iRows * _iCols;
2435     m_iDims = 2;
2436     m_piDims[0] = _iRows;
2437     m_piDims[1] = _iCols;
2438 #ifndef NDEBUG
2439     Inspector::addItem(this);
2440 #endif
2441 }
2442
2443 SparseBool::SparseBool(SparseBool const& src) : matrixBool(new BoolSparse_t(*src.matrixBool))
2444 {
2445     m_iDims = 2;
2446     m_iRows = const_cast<SparseBool*>(&src)->getRows();
2447     m_iCols = const_cast<SparseBool*>(&src)->getCols();
2448     m_iSize = m_iRows * m_iCols;
2449     m_piDims[0] = m_iRows;
2450     m_piDims[1] = m_iCols;
2451 #ifndef NDEBUG
2452     Inspector::addItem(this);
2453 #endif
2454 }
2455
2456 SparseBool::SparseBool(BoolSparse_t* src) : matrixBool(src)
2457 {
2458     m_iRows = src->rows();
2459     m_iCols = src->cols();
2460     m_iSize = m_iRows * m_iCols;
2461     m_iDims = 2;
2462     m_piDims[0] = m_iRows;
2463     m_piDims[1] = m_iCols;
2464 #ifndef NDEBUG
2465     Inspector::addItem(this);
2466 #endif
2467 }
2468
2469 SparseBool::SparseBool(int rows, int cols, int trues, int* inner, int* outer)
2470 {
2471     int* out = nullptr;
2472     int* in = nullptr;
2473
2474     matrixBool = new BoolSparse_t(rows, cols);
2475     matrixBool->reserve((int)trues);
2476     out = matrixBool->outerIndexPtr();
2477     in = matrixBool->innerIndexPtr();
2478
2479     //update outerIndexPtr
2480     memcpy(out, outer, sizeof(int) * (rows + 1));
2481     //update innerIndexPtr
2482     memcpy(in, inner, sizeof(int) * trues);
2483
2484     bool* data = matrixBool->valuePtr();
2485     for (int i = 0; i < trues; ++i)
2486     {
2487         data[i] = true;
2488     }
2489
2490     m_iCols = cols;
2491     m_iRows = rows;
2492     m_iSize = cols * rows;
2493     m_iDims = 2;
2494     m_piDims[0] = m_iRows;
2495     m_piDims[1] = m_iCols;
2496
2497     matrixBool->resizeNonZeros(trues);
2498 }
2499
2500 void SparseBool::create2(int rows, int cols, Bool SPARSE_CONST& src, Double SPARSE_CONST& idx)
2501 {
2502     int nnz = src.getSize();
2503     double* i = idx.get();
2504     double* j = i + idx.getRows();
2505     int* val = src.get();
2506
2507     typedef Eigen::Triplet<bool> T;
2508     std::vector<T> tripletList;
2509     tripletList.reserve((int)nnz);
2510
2511     for (int k = 0; k < nnz; ++k)
2512     {
2513         tripletList.push_back(T(static_cast<int>(i[k]) - 1, static_cast<int>(j[k]) - 1, val[k] == 1));
2514     }
2515
2516     matrixBool = new BoolSparse_t(rows, cols);
2517     matrixBool->setFromTriplets(tripletList.begin(), tripletList.end());
2518
2519     m_iRows = matrixBool->rows();
2520     m_iCols = matrixBool->cols();
2521     m_iSize = cols * rows;
2522     m_iDims = 2;
2523     m_piDims[0] = m_iRows;
2524     m_piDims[1] = m_iCols;
2525     finalize();
2526 }
2527
2528 SparseBool::~SparseBool()
2529 {
2530     delete matrixBool;
2531 #ifndef NDEBUG
2532     Inspector::removeItem(this);
2533 #endif
2534 }
2535
2536 bool SparseBool::toString(std::wostringstream& ostr) const
2537 {
2538     ostr << ::toString(*matrixBool, 0);
2539     return true;
2540 }
2541
2542 void SparseBool::whoAmI() SPARSE_CONST
2543 {
2544     std::cout << "types::SparseBool";
2545 }
2546
2547 SparseBool* SparseBool::clone(void) const
2548 {
2549     return new SparseBool(*this);
2550 }
2551
2552 bool SparseBool::resize(int _iNewRows, int _iNewCols)
2553 {
2554     if (_iNewRows <= getRows() && _iNewCols <= getCols())
2555     {
2556         //nothing to do: hence we do NOT fail
2557         return true;
2558     }
2559
2560     bool res = false;
2561     try
2562     {
2563         //item count
2564         size_t iNonZeros = nbTrue();
2565
2566         BoolSparse_t *newBool = new BoolSparse_t(_iNewRows, _iNewCols);
2567         newBool->reserve((int)iNonZeros);
2568
2569         //coords
2570         int* pRows = new int[iNonZeros * 2];
2571         outputRowCol(pRows);
2572         int* pCols = pRows + iNonZeros;
2573
2574         typedef Eigen::Triplet<bool> triplet;
2575         std::vector<triplet> tripletList;
2576
2577         for (size_t i = 0 ; i < iNonZeros ; i++)
2578         {
2579             tripletList.push_back(triplet((int)pRows[i] - 1, (int)pCols[i] - 1, true));
2580         }
2581
2582         newBool->setFromTriplets(tripletList.begin(), tripletList.end());
2583
2584         delete matrixBool;
2585         matrixBool = newBool;
2586         delete[] pRows;
2587
2588         m_iRows = _iNewRows;
2589         m_iCols = _iNewCols;
2590         m_iSize = _iNewRows * _iNewCols;
2591         m_piDims[0] = m_iRows;
2592         m_piDims[1] = m_iCols;
2593
2594         res = true;
2595
2596     }
2597     catch (...)
2598     {
2599         res = false;
2600     }
2601     return res;
2602 }
2603
2604 SparseBool* SparseBool::insert(typed_list* _pArgs, SparseBool* _pSource)
2605 {
2606     bool bNeedToResize  = false;
2607     int iDims           = (int)_pArgs->size();
2608     if (iDims > 2)
2609     {
2610         //sparse are only in 2 dims
2611         return NULL;
2612     }
2613
2614     typed_list pArg;
2615
2616     int piMaxDim[2];
2617     int piCountDim[2];
2618
2619     //on case of resize
2620     int iNewRows    = 0;
2621     int iNewCols    = 0;
2622
2623     //evaluate each argument and replace by appropriate value and compute the count of combinations
2624     int iSeqCount = checkIndexesArguments(this, _pArgs, &pArg, piMaxDim, piCountDim);
2625     if (iSeqCount == 0)
2626     {
2627         //free pArg content
2628         cleanIndexesArguments(_pArgs, &pArg);
2629         return this;
2630     }
2631
2632     if (iDims < 2)
2633     {
2634         //see as vector
2635         if (getRows() == 1 || getCols() == 1)
2636         {
2637             //vector or scalar
2638             if (getSize() < piMaxDim[0])
2639             {
2640                 bNeedToResize = true;
2641
2642                 //need to enlarge sparse dimensions
2643                 if (getCols() == 1 || getSize() == 0)
2644                 {
2645                     //column vector
2646                     iNewRows    = piMaxDim[0];
2647                     iNewCols    = 1;
2648                 }
2649                 else if (getRows() == 1)
2650                 {
2651                     //row vector
2652                     iNewRows    = 1;
2653                     iNewCols    = piMaxDim[0];
2654                 }
2655             }
2656         }
2657         else if (getSize() < piMaxDim[0])
2658         {
2659             //free pArg content
2660             cleanIndexesArguments(_pArgs, &pArg);
2661             //out of range
2662             return NULL;
2663         }
2664     }
2665     else
2666     {
2667         if (piMaxDim[0] > getRows() || piMaxDim[1] > getCols())
2668         {
2669             bNeedToResize = true;
2670             iNewRows = std::max(getRows(), piMaxDim[0]);
2671             iNewCols = std::max(getCols(), piMaxDim[1]);
2672         }
2673     }
2674
2675     //check number of insertion
2676     if (_pSource->isScalar() == false && _pSource->getSize() != iSeqCount)
2677     {
2678         //free pArg content
2679         cleanIndexesArguments(_pArgs, &pArg);
2680         return NULL;
2681     }
2682
2683     //now you are sure to be able to insert values
2684     if (bNeedToResize)
2685     {
2686         if (resize(iNewRows, iNewCols) == false)
2687         {
2688             //free pArg content
2689             cleanIndexesArguments(_pArgs, &pArg);
2690             return NULL;
2691         }
2692     }
2693
2694     if (iDims == 1)
2695     {
2696         double* pIdx = pArg[0]->getAs<Double>()->get();
2697         for (int i = 0 ; i < iSeqCount ; i++)
2698         {
2699             int iRow = static_cast<int>(pIdx[i] - 1) % getRows();
2700             int iCol = static_cast<int>(pIdx[i] - 1) / getRows();
2701
2702             if (_pSource->isScalar())
2703             {
2704                 set(iRow, iCol, _pSource->get(0, 0), false);
2705             }
2706             else
2707             {
2708                 int iRowOrig = i % _pSource->getRows();
2709                 int iColOrig = i / _pSource->getRows();
2710                 set(iRow, iCol, _pSource->get(iRowOrig, iColOrig), false);
2711             }
2712         }
2713     }
2714     else
2715     {
2716         double* pIdxRow = pArg[0]->getAs<Double>()->get();
2717         int iRowSize    = pArg[0]->getAs<Double>()->getSize();
2718         double* pIdxCol = pArg[1]->getAs<Double>()->get();
2719
2720         for (int i = 0 ; i < iSeqCount ; i++)
2721         {
2722             if (_pSource->isScalar())
2723             {
2724                 set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, _pSource->get(0, 0), false);
2725             }
2726             else
2727             {
2728                 int iRowOrig = i % _pSource->getRows();
2729                 int iColOrig = i / _pSource->getRows();
2730                 set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, _pSource->get(iRowOrig, iColOrig), false);
2731             }
2732         }
2733     }
2734
2735     finalize();
2736
2737     //free pArg content
2738     cleanIndexesArguments(_pArgs, &pArg);
2739
2740     return this;
2741 }
2742
2743 SparseBool* SparseBool::insert(typed_list* _pArgs, InternalType* _pSource)
2744 {
2745     bool bNeedToResize  = false;
2746     int iDims           = (int)_pArgs->size();
2747     if (iDims > 2)
2748     {
2749         //sparse are only in 2 dims
2750         return NULL;
2751     }
2752
2753     typed_list pArg;
2754
2755     int piMaxDim[2];
2756     int piCountDim[2];
2757
2758     //on case of resize
2759     int iNewRows    = 0;
2760     int iNewCols    = 0;
2761     Bool* pSource = _pSource->getAs<Bool>();
2762
2763     //evaluate each argument and replace by appropriate value and compute the count of combinations
2764     int iSeqCount = checkIndexesArguments(this, _pArgs, &pArg, piMaxDim, piCountDim);
2765     if (iSeqCount == 0)
2766     {
2767         //free pArg content
2768         cleanIndexesArguments(_pArgs, &pArg);
2769         return this;
2770     }
2771
2772     if (iDims < 2)
2773     {
2774         //see as vector
2775         if (getRows() == 1 || getCols() == 1)
2776         {
2777             //vector or scalar
2778             bNeedToResize = true;
2779             if (getSize() < piMaxDim[0])
2780             {
2781                 //need to enlarge sparse dimensions
2782                 if (getCols() == 1 || getSize() == 0)
2783                 {
2784                     //column vector
2785                     iNewRows    = piMaxDim[0];
2786                     iNewCols    = 1;
2787                 }
2788                 else if (getRows() == 1)
2789                 {
2790                     //row vector
2791                     iNewRows    = 1;
2792                     iNewCols    = piMaxDim[0];
2793                 }
2794             }
2795         }
2796         else if (getSize() < piMaxDim[0])
2797         {
2798             //free pArg content
2799             cleanIndexesArguments(_pArgs, &pArg);
2800             //out of range
2801             return NULL;
2802         }
2803     }
2804     else
2805     {
2806         if (piMaxDim[0] > getRows() || piMaxDim[1] > getCols())
2807         {
2808             bNeedToResize = true;
2809             iNewRows = std::max(getRows(), piMaxDim[0]);
2810             iNewCols = std::max(getCols(), piMaxDim[1]);
2811         }
2812     }
2813
2814     //check number of insertion
2815     if (pSource->isScalar() == false && pSource->getSize() != iSeqCount)
2816     {
2817         //free pArg content
2818         cleanIndexesArguments(_pArgs, &pArg);
2819         return NULL;
2820     }
2821
2822     //now you are sure to be able to insert values
2823     if (bNeedToResize)
2824     {
2825         if (resize(iNewRows, iNewCols) == false)
2826         {
2827             //free pArg content
2828             cleanIndexesArguments(_pArgs, &pArg);
2829             return NULL;
2830         }
2831     }
2832
2833     if (iDims == 1)
2834     {
2835         double* pIdx = pArg[0]->getAs<Double>()->get();
2836         for (int i = 0 ; i < iSeqCount ; i++)
2837         {
2838             int iRow = static_cast<int>(pIdx[i] - 1) % getRows();
2839             int iCol = static_cast<int>(pIdx[i] - 1) / getRows();
2840             if (pSource->isScalar())
2841             {
2842                 set(iRow, iCol, pSource->get(0) != 0, false);
2843             }
2844             else
2845             {
2846                 set(iRow, iCol, pSource->get(i) != 0, false);
2847             }
2848         }
2849     }
2850     else
2851     {
2852         double* pIdxRow = pArg[0]->getAs<Double>()->get();
2853         int iRowSize    = pArg[0]->getAs<Double>()->getSize();
2854         double* pIdxCol = pArg[1]->getAs<Double>()->get();
2855
2856         for (int i = 0 ; i < iSeqCount ; i++)
2857         {
2858             if (pSource->isScalar())
2859             {
2860                 set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, pSource->get(0) != 0, false);
2861             }
2862             else
2863             {
2864                 int iRowOrig = i % pSource->getRows();
2865                 int iColOrig = i / pSource->getRows();
2866
2867                 set((int)pIdxRow[i % iRowSize] - 1, (int)pIdxCol[i / iRowSize] - 1, pSource->get(iRowOrig, iColOrig) != 0, false);
2868             }
2869         }
2870     }
2871
2872     finalize();
2873
2874     //free pArg content
2875     cleanIndexesArguments(_pArgs, &pArg);
2876     return this;
2877 }
2878
2879 SparseBool* SparseBool::remove(typed_list* _pArgs)
2880 {
2881     SparseBool* pOut = NULL;
2882     int iDims = (int)_pArgs->size();
2883     if (iDims > 2)
2884     {
2885         //sparse are only in 2 dims
2886         return NULL;
2887     }
2888
2889     typed_list pArg;
2890
2891     int piMaxDim[2];
2892     int piCountDim[2];
2893
2894     //evaluate each argument and replace by appropriate value and compute the count of combinations
2895     int iSeqCount = checkIndexesArguments(this, _pArgs, &pArg, piMaxDim, piCountDim);
2896     if (iSeqCount == 0)
2897     {
2898         //free pArg content
2899         cleanIndexesArguments(_pArgs, &pArg);
2900         return this;
2901     }
2902
2903     bool* pbFull = new bool[iDims];
2904     //coord must represent all values on a dimension
2905     for (int i = 0 ; i < iDims ; i++)
2906     {
2907         pbFull[i]       = false;
2908         int iDimToCheck = getVarMaxDim(i, iDims);
2909         int iIndexSize  = pArg[i]->getAs<GenericType>()->getSize();
2910
2911         //we can have index more than once
2912         if (iIndexSize >= iDimToCheck)
2913         {
2914             //size is good, now check datas
2915             double* pIndexes = getDoubleArrayFromDouble(pArg[i]);
2916             for (int j = 0 ; j < iDimToCheck ; j++)
2917             {
2918                 bool bFind = false;
2919                 for (int k = 0 ; k < iIndexSize ; k++)
2920                 {
2921                     if ((int)pIndexes[k] == j + 1)
2922                     {
2923                         bFind = true;
2924                         break;
2925                     }
2926                 }
2927                 pbFull[i]  = bFind;
2928             }
2929         }
2930     }
2931
2932     //only one dims can be not full/entire
2933     bool bNotEntire = false;
2934     int iNotEntire  = 0;
2935     bool bTooMuchNotEntire = false;
2936     for (int i = 0 ; i < iDims ; i++)
2937     {
2938         if (pbFull[i] == false)
2939         {
2940             if (bNotEntire == false)
2941             {
2942                 bNotEntire = true;
2943                 iNotEntire = i;
2944             }
2945             else
2946             {
2947                 bTooMuchNotEntire = true;
2948                 break;
2949             }
2950         }
2951     }
2952
2953     if (bTooMuchNotEntire == true)
2954     {
2955         //free pArg content
2956         cleanIndexesArguments(_pArgs, &pArg);
2957         return NULL;
2958     }
2959
2960     delete[] pbFull;
2961
2962     //find index to keep
2963     int iNotEntireSize          = pArg[iNotEntire]->getAs<GenericType>()->getSize();
2964     double* piNotEntireIndex    = getDoubleArrayFromDouble(pArg[iNotEntire]);
2965     int iKeepSize               = getVarMaxDim(iNotEntire, iDims);
2966     bool* pbKeep                = new bool[iKeepSize];
2967
2968     //fill pbKeep with true value
2969     for (int i = 0 ; i < iKeepSize ; i++)
2970     {
2971         pbKeep[i] = true;
2972     }
2973
2974     for (int i = 0 ; i < iNotEntireSize ; i++)
2975     {
2976         int idx = (int)piNotEntireIndex[i] - 1;
2977
2978         //don't care of value out of bounds
2979         if (idx < iKeepSize)
2980         {
2981             pbKeep[idx] = false;
2982         }
2983     }
2984
2985     int iNewDimSize = 0;
2986     for (int i = 0 ; i < iKeepSize ; i++)
2987     {
2988         if (pbKeep[i] == true)
2989         {
2990             iNewDimSize++;
2991         }
2992     }
2993     delete[] pbKeep;
2994
2995     int* piNewDims = new int[iDims];
2996     for (int i = 0 ; i < iDims ; i++)
2997     {
2998         if (i == iNotEntire)
2999         {
3000             piNewDims[i] = iNewDimSize;
3001         }
3002         else
3003         {
3004             piNewDims[i] = getVarMaxDim(i, iDims);
3005         }
3006     }
3007
3008     //remove last dimension if are == 1
3009     int iOrigDims = iDims;
3010     for (int i = (iDims - 1) ; i >= 2 ; i--)
3011     {
3012         if (piNewDims[i] == 1)
3013         {
3014             iDims--;
3015         }
3016         else
3017         {
3018             break;
3019         }
3020     }
3021
3022     if (iDims == 1)
3023     {
3024         if (iNewDimSize == 0)
3025         {
3026             //free pArg content
3027             cleanIndexesArguments(_pArgs, &pArg);
3028             return new SparseBool(0, 0);
3029         }
3030         else
3031         {
3032             //two cases, depends of original matrix/vector
3033             if ((*_pArgs)[0]->isColon() == false && m_iDims == 2 && m_piDims[0] == 1 && m_piDims[1] != 1)
3034             {
3035                 //special case for row vector
3036                 pOut = new SparseBool(1, iNewDimSize);
3037                 //in this case we have to care of 2nd dimension
3038                 //iNotEntire = 1;
3039             }
3040             else
3041             {
3042                 pOut = new SparseBool(iNewDimSize, 1);
3043             }
3044         }
3045     }
3046     else
3047     {
3048         pOut = new SparseBool(piNewDims[0], piNewDims[0]);
3049     }
3050
3051     delete[] piNewDims;
3052     //find a way to copy existing data to new variable ...
3053     int iNewPos = 0;
3054     int* piIndexes = new int[iOrigDims];
3055     int* piViewDims = new int[iOrigDims];
3056     for (int i = 0 ; i < iOrigDims ; i++)
3057     {
3058         piViewDims[i] = getVarMaxDim(i, iOrigDims);
3059     }
3060
3061     for (int i = 0 ; i < getSize() ; i++)
3062     {
3063         bool bByPass = false;
3064         getIndexesWithDims(i, piIndexes, piViewDims, iOrigDims);
3065
3066         //check if piIndexes use removed indexes
3067         for (int j = 0 ; j < iNotEntireSize ; j++)
3068         {
3069             if ((piNotEntireIndex[j] - 1) == piIndexes[iNotEntire])
3070             {
3071                 //by pass this value
3072                 bByPass = true;
3073                 break;
3074             }
3075         }
3076
3077         if (bByPass == false)
3078         {
3079             //compute new index
3080             pOut->set(iNewPos, get(i));
3081             iNewPos++;
3082         }
3083     }
3084
3085     //free allocated data
3086     for (int i = 0 ; i < iDims ; i++)
3087     {
3088         if (pArg[i] != (*_pArgs)[i])
3089         {
3090             delete pArg[i];
3091         }
3092     }
3093
3094     delete[] piIndexes;
3095     delete[] piViewDims;
3096
3097     //free pArg content
3098     cleanIndexesArguments(_pArgs, &pArg);
3099
3100     return pOut;
3101 }
3102
3103 bool SparseBool::append(int r, int c, SparseBool SPARSE_CONST* src)
3104 {
3105     doAppend(*src, r, c, *matrixBool);
3106     finalize();
3107     return true;
3108 }
3109
3110 InternalType* SparseBool::insertNew(typed_list* _pArgs, InternalType* _pSource)
3111 {
3112     typed_list pArg;
3113     InternalType *pOut  = NULL;
3114     SparseBool* pSource = _pSource->getAs<SparseBool>();
3115
3116     int iDims           = (int)_pArgs->size();
3117     int* piMaxDim       = new int[iDims];
3118     int* piCountDim     = new int[iDims];
3119     bool bUndefine      = false;
3120
3121     //evaluate each argument and replace by appropriate value and compute the count of combinations
3122     int iSeqCount = checkIndexesArguments(NULL, _pArgs, &pArg, piMaxDim, piCountDim);
3123     if (iSeqCount == 0)
3124     {
3125         //free pArg content
3126         cleanIndexesArguments(_pArgs, &pArg);
3127         return createEmptyDouble();
3128     }
3129
3130     if (iSeqCount < 0)
3131     {
3132         iSeqCount = -iSeqCount;
3133         bUndefine = true;
3134     }
3135
3136     if (bUndefine)
3137     {
3138         //manage : and $ in creation by insertion
3139         int iSource         = 0;
3140         int *piSourceDims   = pSource->getDimsArray();
3141
3142         for (int i = 0 ; i < iDims ; i++)
3143         {
3144             if (pArg[i] == NULL)
3145             {
3146                 //undefine value
3147                 if (pSource->isScalar())
3148                 {
3149                     piMaxDim[i]     = 1;
3150                     piCountDim[i]   = 1;
3151                 }
3152                 else
3153                 {
3154                     piMaxDim[i]     = piSourceDims[iSource];
3155                     piCountDim[i]   = piSourceDims[iSource];
3156                 }
3157                 iSource++;
3158                 //replace pArg value by the new one
3159                 pArg[i] = createDoubleVector(piMaxDim[i]);
3160             }
3161             //else
3162             //{
3163             //    piMaxDim[i] = piCountDim[i];
3164             //}
3165         }
3166     }
3167
3168     //remove last dimension at size 1
3169     //remove last dimension if are == 1
3170     for (int i = (iDims - 1) ; i >= 2 ; i--)
3171     {
3172         if (piMaxDim[i] == 1)
3173         {
3174             iDims--;
3175             pArg.pop_back();
3176         }
3177         else
3178         {
3179             break;
3180         }
3181     }
3182
3183     if (checkArgValidity(pArg) == false)
3184     {
3185         //free pArg content
3186         cleanIndexesArguments(_pArgs, &pArg);
3187         //contain bad index, like <= 0, ...
3188         return NULL;
3189     }
3190
3191     if (iDims == 1)
3192     {
3193         if (pSource->getCols() == 1)
3194         {
3195             pOut = new SparseBool(piCountDim[0], 1);
3196         }
3197         else
3198         {
3199             //rows == 1
3200             pOut = new SparseBool(1, piCountDim[0]);
3201         }
3202     }
3203     else
3204     {
3205         pOut = new SparseBool(piMaxDim[0], piMaxDim[1]);
3206     }
3207
3208     //fill with null item
3209     SparseBool* pSpOut = pOut->getAs<SparseBool>();
3210
3211     //insert values in new matrix
3212     InternalType* pOut2 = pSpOut->insert(&pArg, pSource);
3213     if (pOut != pOut2)
3214     {
3215         delete pOut;
3216     }
3217
3218     //free pArg content
3219     cleanIndexesArguments(_pArgs, &pArg);
3220
3221     return pOut2;
3222 }
3223
3224 SparseBool* SparseBool::extract(int nbCoords, int SPARSE_CONST* coords, int SPARSE_CONST* maxCoords, int SPARSE_CONST* resSize, bool asVector) SPARSE_CONST
3225 {
3226     if ( (asVector && maxCoords[0] > getSize()) ||
3227     (asVector == false && maxCoords[0] > getRows()) ||
3228     (asVector == false && maxCoords[1] > getCols()))
3229     {
3230         return 0;
3231     }
3232
3233     SparseBool * pSp (0);
3234     if (asVector)
3235     {
3236         pSp = (getRows() == 1) ?  new SparseBool(1, resSize[0]) : new SparseBool(resSize[0], 1);
3237         mycopy_n(makeMatrixIterator<bool>(*this,  Coords<true>(coords, getRows())), nbCoords
3238         , makeMatrixIterator<bool>(*(pSp->matrixBool), RowWiseFullIterator(pSp->getRows(), pSp->getCols())));
3239     }
3240     else
3241     {
3242         pSp = new SparseBool(resSize[0], resSize[1]);
3243         mycopy_n(makeMatrixIterator<bool>(*this,  Coords<false>(coords, getRows())), nbCoords
3244         , makeMatrixIterator<bool>(*(pSp->matrixBool), RowWiseFullIterator(pSp->getRows(), pSp->getCols())));
3245
3246     }
3247     return pSp;
3248 }
3249
3250 /*
3251 * create a new SparseBool of dims according to resSize and fill it from currentSparseBool (along coords)
3252 */
3253 InternalType* SparseBool::extract(typed_list* _pArgs)
3254 {
3255     SparseBool* pOut    = NULL;
3256     int iDims           = (int)_pArgs->size();
3257     typed_list pArg;
3258
3259     int* piMaxDim       = new int[iDims];
3260     int* piCountDim     = new int[iDims];
3261
3262     //evaluate each argument and replace by appropriate value and compute the count of combinations
3263     int iSeqCount = checkIndexesArguments(this, _pArgs, &pArg, piMaxDim, piCountDim);
3264     if (iSeqCount == 0)
3265     {
3266         //free pArg content
3267         cleanIndexesArguments(_pArgs, &pArg);
3268         if (_pArgs->size() == 0)
3269         {
3270             delete[] piMaxDim;
3271             delete[] piCountDim;
3272             //a()
3273             return this;
3274         }
3275         else
3276         {
3277             delete[] piMaxDim;
3278             delete[] piCountDim;
3279             //a([])
3280             return Double::Empty();
3281         }
3282     }
3283
3284     if (iDims < 2)
3285     {
3286         // Check that we stay inside the input size.
3287         if (piMaxDim[0] <= getSize())
3288         {
3289             int iNewRows = 0;
3290             int iNewCols = 0;
3291
3292             if (getRows() == 1 && getCols() != 1 && (*_pArgs)[0]->isColon() == false)
3293             {
3294                 //special case for row vector
3295                 iNewRows = 1;
3296                 iNewCols = piCountDim[0];
3297             }
3298             else
3299             {
3300                 iNewRows = piCountDim[0];
3301                 iNewCols = 1;
3302             }
3303
3304             pOut = new SparseBool(iNewRows, iNewCols);
3305             double* pIdx = pArg[0]->getAs<Double>()->get();
3306             // Write in output all elements extract from input.
3307             for (int i = 0 ; i < iSeqCount ; i++)
3308             {
3309                 if (pIdx[i] < 1)
3310                 {
3311                     delete pOut;
3312                     pOut = NULL;
3313                     break;
3314                 }
3315                 int iRowRead = static_cast<int>(pIdx[i] - 1) % getRows();
3316                 int iColRead = static_cast<int>(pIdx[i] - 1) / getRows();
3317
3318                 int iRowWrite = static_cast<int>(i) % iNewRows;
3319                 int iColWrite = static_cast<int>(i) / iNewRows;
3320
3321                 bool bValue = get(iRowRead, iColRead);
3322                 if (bValue)
3323                 {
3324                     //only non zero values
3325                     pOut->set(iRowWrite, iColWrite, true, false);
3326                 }
3327             }
3328         }
3329         else
3330         {
3331             delete[] piMaxDim;
3332             delete[] piCountDim;
3333             //free pArg content
3334             cleanIndexesArguments(_pArgs, &pArg);
3335             return NULL;
3336         }
3337     }
3338     else
3339     {
3340         // Check that we stay inside the input size.
3341         if (piMaxDim[0] <= getRows() && piMaxDim[1] <= getCols())
3342         {
3343             double* pIdxRow = pArg[0]->getAs<Double>()->get();
3344             double* pIdxCol = pArg[1]->getAs<Double>()->get();
3345
3346             int iNewRows = pArg[0]->getAs<Double>()->getSize();
3347             int iNewCols = pArg[1]->getAs<Double>()->getSize();
3348
3349             pOut = new SparseBool(iNewRows, iNewCols);
3350
3351             int iPos = 0;
3352             // Write in output all elements extract from input.
3353             for (int iRow = 0 ; iRow < iNewRows ; iRow++)
3354             {
3355                 for (int iCol = 0 ; iCol < iNewCols ; iCol++)
3356                 {
3357                     if ((pIdxRow[iRow] < 1) || (pIdxCol[iCol] < 1))
3358                     {
3359                         delete pOut;
3360                         pOut = NULL;
3361                         break;
3362                     }
3363                     bool bValue = get((int)pIdxRow[iRow] - 1, (int)pIdxCol[iCol] - 1);
3364                     if (bValue)
3365                     {
3366                         //only non zero values
3367                         pOut->set(iRow, iCol, true, false);
3368                     }
3369                     iPos++;
3370                 }
3371             }
3372         }
3373         else
3374         {
3375             delete[] piMaxDim;
3376             delete[] piCountDim;
3377             //free pArg content
3378             cleanIndexesArguments(_pArgs, &pArg);
3379             return NULL;
3380         }
3381     }
3382
3383     finalize();
3384
3385     delete[] piMaxDim;
3386     delete[] piCountDim;
3387     //free pArg content
3388     cleanIndexesArguments(_pArgs, &pArg);
3389
3390     return pOut;
3391 }
3392
3393 std::size_t SparseBool::nbTrue() const
3394 {
3395     return  matrixBool->nonZeros() ;
3396 }
3397 std::size_t SparseBool::nbTrue(std::size_t r) const
3398 {
3399     int* piIndex = matrixBool->outerIndexPtr();
3400     return piIndex[r + 1] - piIndex[r];
3401 }
3402
3403
3404 void SparseBool::setTrue(bool finalize)
3405 {
3406     int rows = getRows();
3407     int cols = getCols();
3408
3409     typedef Eigen::Triplet<bool> triplet;
3410     std::vector<triplet> tripletList;
3411
3412     for (int i = 0; i < rows; ++i)
3413     {
3414         for (int j = 0; j < cols; ++j)
3415         {
3416             tripletList.push_back(triplet(i, j, true));
3417         }
3418     }
3419
3420     matrixBool->setFromTriplets(tripletList.begin(), tripletList.end());
3421
3422     if (finalize)
3423     {
3424         matrixBool->finalize();
3425     }
3426 }
3427
3428 void SparseBool::setFalse(bool finalize)
3429 {
3430     int rows = getRows();
3431     int cols = getCols();
3432
3433     typedef Eigen::Triplet<bool> triplet;
3434     std::vector<triplet> tripletList;
3435
3436     for (int i = 0; i < rows; ++i)
3437     {
3438         for (int j = 0; j < cols; ++j)
3439         {
3440             tripletList.push_back(triplet(i, j, false));
3441         }
3442     }
3443
3444     matrixBool->setFromTriplets(tripletList.begin(), tripletList.end());
3445
3446     if (finalize)
3447     {
3448         matrixBool->finalize();
3449     }
3450 }
3451
3452 int* SparseBool::getNbItemByRow(int* _piNbItemByRows)
3453 {
3454     int* piNbItemByRows = new int[getRows() + 1];
3455     mycopy_n(matrixBool->outerIndexPtr(), getRows() + 1, piNbItemByRows);
3456
3457     for (int i = 0 ; i < getRows() ; i++)
3458     {
3459         _piNbItemByRows[i] = piNbItemByRows[i + 1] - piNbItemByRows[i];
3460     }
3461
3462     delete[] piNbItemByRows;
3463     return _piNbItemByRows;
3464 }
3465
3466 int* SparseBool::getColPos(int* _piColPos)
3467 {
3468     mycopy_n(matrixBool->innerIndexPtr(), nbTrue(), _piColPos);
3469     for (int i = 0; i < nbTrue(); i++)
3470     {
3471         _piColPos[i]++;
3472     }
3473
3474     return _piColPos;
3475 }
3476
3477 int* SparseBool::outputRowCol(int* out)const
3478 {
3479     return sparseTransform(*matrixBool, sparseTransform(*matrixBool, out, GetRow<BoolSparse_t>()), GetCol<BoolSparse_t>());
3480 }
3481
3482 int* SparseBool::getInnerPtr(int* count)
3483 {
3484     *count = matrixBool->innerSize();
3485     return matrixBool->innerIndexPtr();
3486 }
3487
3488 int* SparseBool::getOuterPtr(int* count)
3489 {
3490     *count = matrixBool->outerSize();
3491     return matrixBool->outerIndexPtr();
3492 }
3493
3494
3495 bool SparseBool::operator==(const InternalType& it) SPARSE_CONST
3496 {
3497     SparseBool* otherSparse = const_cast<SparseBool*>(dynamic_cast<SparseBool const*>(&it));/* types::GenericType is not const-correct :( */
3498     return (otherSparse
3499     && (otherSparse->getRows() == getRows())
3500     && (otherSparse->getCols() == getCols())
3501     && equal(*matrixBool, *otherSparse->matrixBool));
3502 }
3503
3504 bool SparseBool::operator!=(const InternalType& it) SPARSE_CONST
3505 {
3506     return !(*this == it);
3507 }
3508
3509 void SparseBool::finalize()
3510 {
3511     matrixBool->prune(&keepForSparse<bool>);
3512     matrixBool->finalize();
3513 }
3514
3515 bool SparseBool::get(int r, int c) SPARSE_CONST
3516 {
3517     return matrixBool->coeff(r, c);
3518 }
3519
3520 bool SparseBool::set(int _iRows, int _iCols, bool _bVal, bool _bFinalize) SPARSE_CONST
3521 {
3522     matrixBool->coeffRef(_iRows, _iCols) = _bVal;
3523
3524     if (_bFinalize)
3525     {
3526         finalize();
3527     }
3528
3529     return true;
3530 }
3531
3532 void SparseBool::fill(Bool& dest, int r, int c) SPARSE_CONST
3533 {
3534     mycopy_n(makeMatrixIterator<bool >(*matrixBool, RowWiseFullIterator(getRows(), getCols())), getSize()
3535     , makeMatrixIterator<bool >(dest, RowWiseFullIterator(dest.getRows(), dest.getCols(), r, c)));
3536 }
3537
3538 Sparse* SparseBool::newOnes() const
3539 {
3540     return new Sparse(new types::Sparse::RealSparse_t(matrixBool->cast<double>()), 0);
3541 }
3542
3543 SparseBool* SparseBool::newNotEqualTo(SparseBool const&o) const
3544 {
3545     return cwiseOp<std::not_equal_to>(*this, o);
3546 }
3547
3548 SparseBool* SparseBool::newEqualTo(SparseBool const&o) const
3549 {
3550     return cwiseOp<std::equal_to>(*this, o);
3551 }
3552
3553 SparseBool* SparseBool::newLogicalOr(SparseBool const&o) const
3554 {
3555     return cwiseOp<std::logical_or>(*this, o);
3556 }
3557
3558 SparseBool* SparseBool::newLogicalAnd(SparseBool const&o) const
3559 {
3560     return cwiseOp<std::logical_and>(*this, o);
3561 }
3562
3563 bool SparseBool::reshape(int* _piDims, int _iDims)
3564 {
3565     bool bOk = false;
3566     int iCols = 1;
3567
3568     if (_iDims == 2)
3569     {
3570         iCols = _piDims[1];
3571     }
3572
3573     if (_iDims <= 2)
3574     {
3575         bOk = reshape(_piDims[0], iCols);
3576     }
3577
3578     return bOk;
3579 }
3580
3581 bool SparseBool::reshape(int _iNewRows, int _iNewCols)
3582 {
3583     if (_iNewRows * _iNewCols != getRows() * getCols())
3584     {
3585         return false;
3586     }
3587
3588     bool res = false;
3589     try
3590     {
3591         //item count
3592         size_t iNonZeros = matrixBool->nonZeros();
3593         BoolSparse_t *newBool = new BoolSparse_t(_iNewRows, _iNewCols);
3594         newBool->reserve((int)iNonZeros);
3595
3596         //coords
3597         int* pRows = new int[iNonZeros * 2];
3598         outputRowCol(pRows);
3599         int* pCols = pRows + iNonZeros;
3600
3601         typedef Eigen::Triplet<bool> triplet;
3602         std::vector<triplet> tripletList;
3603
3604         for (size_t i = 0 ; i < iNonZeros ; i++)
3605         {
3606             int iCurrentPos = ((int)pCols[i] - 1) * getRows() + ((int)pRows[i] - 1);
3607             tripletList.push_back(triplet((int)(iCurrentPos % _iNewRows), (int)(iCurrentPos / _iNewRows), true));
3608         }
3609
3610         newBool->setFromTriplets(tripletList.begin(), tripletList.end());
3611
3612         delete matrixBool;
3613         matrixBool = newBool;
3614         delete[] pRows;
3615
3616         m_iRows = _iNewRows;
3617         m_iCols = _iNewCols;
3618         m_iSize = _iNewRows * _iNewCols;
3619
3620         m_iDims = 2;
3621         m_piDims[0] = m_iRows;
3622         m_piDims[1] = m_iCols;
3623
3624         finalize();
3625
3626         res = true;
3627     }
3628     catch (...)
3629     {
3630         res = false;
3631     }
3632     return res;
3633 }
3634
3635 bool SparseBool::transpose(InternalType *& out)
3636 {
3637     out = new SparseBool(new BoolSparse_t(matrixBool->transpose()));
3638     return true;
3639 }
3640
3641 template<typename T>
3642 void neg(const int r, const int c, const T * const in, Eigen::SparseMatrix<bool, 1> * const out)
3643 {
3644     for (int i = 0; i < r; i++)
3645     {
3646         for (int j = 0; j < c; j++)
3647         {
3648             out->coeffRef(i, j) = !in->coeff(i, j);
3649         }
3650     }
3651
3652     out->prune(&keepForSparse<bool>);
3653     out->finalize();
3654 }
3655
3656
3657 }