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