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