scilab/modules/ast/src/cpp/types/sparse.cpp

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