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