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