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