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