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