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