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

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