PageRenderTime 54ms CodeModel.GetById 23ms RepoModel.GetById 0ms app.codeStats 0ms

/_3rdParty/MathLib/Vector.h

https://bitbucket.org/byzhang/mldemos
C Header | 1087 lines | 840 code | 100 blank | 147 comment | 85 complexity | 2fa31eb6d553b10e802c6f552f0e7603 MD5 | raw file
Possible License(s): BSD-3-Clause
  1. /*
  2. * Copyright (C) 2010 Learning Algorithms and Systems Laboratory, EPFL, Switzerland
  3. * Author: Eric Sauser
  4. * email: eric.sauser@a3.epf.ch
  5. * website: lasa.epfl.ch
  6. *
  7. * Permission is granted to copy, distribute, and/or modify this program
  8. * under the terms of the GNU General Public License, version 2 or any
  9. * later version published by the Free Software Foundation.
  10. *
  11. * This program is distributed in the hope that it will be useful, but
  12. * WITHOUT ANY WARRANTY; without even the implied warranty of
  13. * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
  14. * Public License for more details
  15. */
  16. #ifndef VECTOR_H
  17. #define VECTOR_H
  18. #include "MathLibCommon.h"
  19. #include "Macros.h"
  20. #include <assert.h>
  21. #include <math.h>
  22. #include <iostream>
  23. #include <vector>
  24. using namespace std;
  25. #ifndef NULL
  26. #define NULL 0
  27. #endif
  28. #ifdef USE_MATHLIB_NAMESPACE
  29. namespace MathLib {
  30. #endif
  31. #ifdef USE_T_EXTENSIONS
  32. template<unsigned int ROW> class TVector;
  33. #endif
  34. typedef vector<unsigned int> IndicesVector;
  35. /**
  36. * \defgroup MathLib MathLib
  37. *
  38. * \brief A mathematical library for doing various things such as
  39. * matrix manipulation
  40. */
  41. /**
  42. * \class Vector
  43. *
  44. * \ingroup MathLib
  45. *
  46. * \brief The basic vector class
  47. *
  48. * This vector class can be used for doing various vector manipulation.
  49. * This should be combined with the Matrix class for doing almost anything
  50. * you ever dreamt of.
  51. */
  52. class Matrix;
  53. class Vector
  54. {
  55. friend class Matrix;
  56. #ifdef USE_T_EXTENSIONS
  57. template<unsigned int ROW> friend class TVector;
  58. #endif
  59. protected:
  60. static REALTYPE undef; ///< The default value returned when a vector out-of-bound is reached
  61. unsigned int row; ///< The size of the vector
  62. REALTYPE *_; ///< The data array.
  63. public:
  64. /// Create an empty vector
  65. inline Vector() {
  66. row = 0;
  67. _ = NULL;
  68. }
  69. /// The destructor
  70. inline virtual ~Vector(){
  71. Release();
  72. }
  73. /// The copy constructor
  74. inline Vector(const Vector &vector){
  75. row = 0;
  76. _ = NULL;
  77. Set(vector);
  78. }
  79. /**
  80. * \brief Create a sized vector
  81. * \param size The size
  82. * \param clear Tell if the vector sould be set to 0
  83. */
  84. inline Vector(unsigned int size, bool clear = true){
  85. row = 0;
  86. _ = NULL;
  87. Resize(size,false);
  88. if(clear) Zero();
  89. }
  90. /**
  91. * \brief Create a vector by copying a memory area
  92. * \param array The array of REALTYPE to be copied
  93. * \param size The size
  94. */
  95. inline Vector(const REALTYPE *array, unsigned int size){
  96. row = 0;
  97. _ = NULL;
  98. Set(array,size);
  99. }
  100. #ifdef USE_T_EXTENSIONS
  101. /// Copy Contructor of a template Vector (see TVector)
  102. template<unsigned int ROW> inline Vector(const TVector<ROW> &vector){
  103. row = 0;
  104. _ = NULL;
  105. Set(vector._,ROW);
  106. }
  107. #endif
  108. /// An assignment function
  109. inline Vector& Set(const Vector &vector){
  110. Resize(vector.row,false);
  111. memcpy(_,vector._,row*sizeof(REALTYPE));
  112. return *this;
  113. }
  114. /// An assignment function
  115. inline Vector& Set(const REALTYPE * array, unsigned int size){
  116. Resize(size,false);
  117. if((row)&&(array))
  118. memcpy(_,array,size*sizeof(REALTYPE));
  119. return *this;
  120. }
  121. /// Set all values to 0
  122. inline Vector& Zero(){
  123. if(row)
  124. memset(_,0,row*sizeof(REALTYPE));
  125. return *this;
  126. }
  127. /// Set all values to 1
  128. inline Vector& One(){
  129. REALTYPE *src = _;
  130. unsigned int len = row;
  131. while(len--)
  132. *(src++) = R_ONE;
  133. return *this;
  134. }
  135. /// Set random values uniformly distributed between 0 and 1
  136. inline Vector& Random(){
  137. REALTYPE *src = _;
  138. unsigned int len = row;
  139. while(len--){
  140. *(src++) = RND(R_ONE);
  141. }
  142. return *this;
  143. }
  144. /// Get the size of the vector
  145. inline unsigned int Size() const{
  146. return row;
  147. }
  148. /// Get the data pointer
  149. inline REALTYPE *Array() const{
  150. return _;
  151. }
  152. /**
  153. * \brief Get the reference to a vector value
  154. * \param row The row number (starting from 0)
  155. * \return The corresponding value
  156. */
  157. inline REALTYPE& operator[] (const unsigned int row){
  158. if(row<this->row)
  159. return _[row];
  160. return undef;
  161. }
  162. inline REALTYPE& operator() (const unsigned int row){
  163. if(row<this->row)
  164. return _[row];
  165. return undef;
  166. }
  167. inline REALTYPE& Ref(const unsigned int row){
  168. if(row<this->row)
  169. return _[row];
  170. return undef;
  171. }
  172. inline REALTYPE& RefNoCheck(const unsigned int row){
  173. return _[row];
  174. }
  175. /**
  176. * \brief Get a vector value
  177. * \param row The row number (starting from 0)
  178. * \return The corresponding value
  179. */
  180. inline REALTYPE At(const unsigned int row) const {
  181. if(row<this->row)
  182. return _[row];
  183. return undef;
  184. }
  185. inline REALTYPE AtNoCheck(const unsigned int row) const {
  186. return _[row];
  187. }
  188. inline Vector& operator = (const Vector &vector){
  189. return Set(vector);
  190. }
  191. inline Vector& operator = (const REALTYPE val){
  192. REALTYPE *src = _;
  193. unsigned int len = row;
  194. while(len--)
  195. *(src++) = val;
  196. return *this;
  197. }
  198. /**
  199. * \brief The - operator
  200. * \return A new vector
  201. */
  202. inline Vector operator - () const{
  203. Vector result;
  204. REALTYPE *src = _;
  205. REALTYPE *dst = result._;
  206. unsigned int len = row;
  207. while(len--)
  208. *(dst++) = -(*(src++));
  209. return result;
  210. }
  211. inline Vector& SMinus()
  212. {
  213. REALTYPE *src = _;
  214. unsigned int len = row;
  215. while(len--){
  216. *(src) = -(*(src));
  217. src++;
  218. }
  219. return *this;
  220. }
  221. inline Vector& operator += (const Vector &vector){
  222. REALTYPE *src = _;
  223. REALTYPE *dst = vector._;
  224. unsigned int len = (row<=vector.row?row:vector.row);
  225. while(len--)
  226. *(src++) += *(dst++);
  227. return *this;
  228. }
  229. inline Vector& operator -= (const Vector &vector){
  230. REALTYPE *src = _;
  231. REALTYPE *dst = vector._;
  232. unsigned int len = (row<=vector.row?row:vector.row);
  233. while(len--)
  234. *(src++) -= *(dst++);
  235. return *this;
  236. }
  237. /// Element-wise multiplication
  238. inline Vector& operator ^= (const Vector &vector){
  239. REALTYPE *src = _;
  240. REALTYPE *dst = vector._;
  241. unsigned int len = (row<=vector.row?row:vector.row);
  242. while(len--)
  243. *(src++) *= *(dst++);
  244. return *this;
  245. }
  246. /// Element-wise division
  247. inline Vector& operator /= (const Vector &vector){
  248. REALTYPE *src = _;
  249. REALTYPE *dst = vector._;
  250. unsigned int len = (row<=vector.row?row:vector.row);
  251. while(len--)
  252. *(src++) /= *(dst++);
  253. return *this;
  254. }
  255. inline Vector& operator += (REALTYPE scalar){
  256. REALTYPE *src = _;
  257. unsigned int len = row;
  258. while(len--) *(src++) += scalar;
  259. return *this;
  260. }
  261. inline Vector& operator -= (REALTYPE scalar){
  262. REALTYPE *src = _;
  263. unsigned int len = row;
  264. while(len--) *(src++) -= scalar;
  265. return *this;
  266. }
  267. inline Vector& operator *= (REALTYPE scalar){
  268. REALTYPE *src = _;
  269. unsigned int len = row;
  270. while(len--) *(src++) *= scalar;
  271. return *this;
  272. }
  273. inline Vector& operator /= (REALTYPE scalar){
  274. scalar = R_ONE / scalar;
  275. REALTYPE *src = _;
  276. unsigned int len = row;
  277. while(len--) *(src++) *= scalar;
  278. return *this;
  279. }
  280. inline Vector operator + (const Vector &vector) const{
  281. Vector result;
  282. return Add(vector,result);
  283. }
  284. inline Vector operator - (const Vector &vector) const{
  285. Vector result;
  286. return Sub(vector,result);
  287. }
  288. /// Element-wise multiplication
  289. inline Vector operator ^ (const Vector &vector) const{
  290. Vector result;
  291. return PMult(vector,result);
  292. }
  293. /// Element-wise division
  294. inline Vector operator / (const Vector &vector) const{
  295. Vector result;
  296. return PDiv(vector,result);
  297. }
  298. inline REALTYPE operator * (const Vector &vector) const{
  299. return Dot(vector);
  300. }
  301. inline Vector operator + (REALTYPE scalar) const{
  302. Vector result;
  303. return Add(scalar,result);
  304. }
  305. inline Vector operator - (REALTYPE scalar) const{
  306. Vector result;
  307. return Sub(scalar,result);
  308. }
  309. inline Vector operator * (REALTYPE scalar) const{
  310. Vector result;
  311. return Mult(scalar,result);
  312. }
  313. inline Vector operator / (REALTYPE scalar) const{
  314. Vector result;
  315. return Div(scalar,result);
  316. }
  317. inline bool operator == (const Vector& vector) const{
  318. if(row!=vector.row) return false;
  319. REALTYPE *src = _;
  320. REALTYPE *dst = vector._;
  321. unsigned int len = row;
  322. while(len--)
  323. if(*(src++) != *(dst++)) return false;
  324. return true;
  325. }
  326. inline bool operator != (const Vector& vector) const{
  327. return !(*this == vector);
  328. }
  329. inline Vector& ScaleAddTo(REALTYPE scale, Vector &result) const
  330. {
  331. REALTYPE *src = _;
  332. REALTYPE *dst = result._;
  333. unsigned int len = (row<=result.row?row:result.row);
  334. while(len--)
  335. *(dst++) += *(src++) * scale;
  336. return result;
  337. }
  338. /**
  339. * \brief Sum two vector in a faster way than using the + operator
  340. * \param vector The second vector to be summed up
  341. * \param result The result
  342. * \return The result vector
  343. */
  344. inline Vector& Add(const Vector &vector, Vector& result) const{
  345. result.Resize(row,false);
  346. REALTYPE *src0 = _;
  347. REALTYPE *src = vector._;
  348. REALTYPE *dst = result._;
  349. if(row<=vector.row){
  350. unsigned int len = row;
  351. while(len--)
  352. *(dst++) = *(src0++) + *(src++);
  353. }else{
  354. unsigned int len = vector.row;
  355. while(len--)
  356. *(dst++) = *(src0++) + *(src++);
  357. len = row - vector.row;
  358. while(len--)
  359. *(dst++) = *(src0++);
  360. }
  361. return result;
  362. }
  363. /**
  364. * \brief Substract two vector in a faster way than using the + operator
  365. * \param vector The substracting vector
  366. * \param result The result
  367. * \return The result vector
  368. */
  369. inline Vector& Sub(const Vector &vector, Vector& result) const{
  370. result.Resize(row,false);
  371. REALTYPE *src0 = _;
  372. REALTYPE *src = vector._;
  373. REALTYPE *dst = result._;
  374. if(row<=vector.row){
  375. unsigned int len = row;
  376. while(len--)
  377. *(dst++) = *(src0++) - *(src++);
  378. }else{
  379. unsigned int len = vector.row;
  380. while(len--)
  381. *(dst++) = *(src0++) - *(src++);
  382. len = row - vector.row;
  383. while(len--)
  384. *(dst++) = *(src0++);
  385. }
  386. return result;
  387. }
  388. /// Element-wise multiplication
  389. inline Vector& PMult(const Vector &vector, Vector& result) const{
  390. result.Resize(row,false);
  391. REALTYPE *src0 = _;
  392. REALTYPE *src = vector._;
  393. REALTYPE *dst = result._;
  394. if(row<=vector.row){
  395. unsigned int len = row;
  396. while(len--)
  397. *(dst++) = *(src0++) * *(src++);
  398. }else{
  399. unsigned int len = vector.row;
  400. while(len--)
  401. *(dst++) = *(src0++) * *(src++);
  402. len = row - vector.row;
  403. while(len--)
  404. *(dst++) = *(src0++);
  405. }
  406. return result;
  407. }
  408. /// Element-wise division
  409. inline Vector& PDiv(const Vector &vector, Vector& result) const{
  410. result.Resize(row,false);
  411. REALTYPE *src0 = _;
  412. REALTYPE *src = vector._;
  413. REALTYPE *dst = result._;
  414. if(row<=vector.row){
  415. unsigned int len = row;
  416. while(len--)
  417. *(dst++) = *(src0++) / *(src++);
  418. }else{
  419. unsigned int len = vector.row;
  420. while(len--)
  421. *(dst++) = *(src0++) / *(src++);
  422. len = row - vector.row;
  423. while(len--)
  424. *(dst++) = *(src0++);
  425. }
  426. return result;
  427. }
  428. /// Scalar summation
  429. inline Vector& Add(REALTYPE scalar, Vector& result) const{
  430. result.Resize(row,false);
  431. REALTYPE *src = _;
  432. REALTYPE *dst = result._;
  433. unsigned int len = row;
  434. while(len--)
  435. *(dst++) = *(src++) + scalar;
  436. return result;
  437. }
  438. /// Scalar substraction
  439. inline Vector& Sub(REALTYPE scalar, Vector& result) const{
  440. result.Resize(row,false);
  441. REALTYPE *src = _;
  442. REALTYPE *dst = result._;
  443. unsigned int len = row;
  444. while(len--)
  445. *(dst++) = *(src++) - scalar;
  446. return result;
  447. }
  448. /// Scalar multiplication
  449. inline Vector& Mult(REALTYPE scalar, Vector& result) const{
  450. result.Resize(row,false);
  451. REALTYPE *src = _;
  452. REALTYPE *dst = result._;
  453. unsigned int len = row;
  454. while(len--)
  455. *(dst++) = *(src++) * scalar;
  456. return result;
  457. }
  458. /// Scalar division
  459. inline Vector& Div(REALTYPE scalar, Vector& result) const{
  460. scalar = R_ONE/scalar;
  461. result.Resize(row,false);
  462. REALTYPE *src = _;
  463. REALTYPE *dst = result._;
  464. unsigned int len = row;
  465. while(len--)
  466. *(dst++) = *(src++) * scalar;
  467. return result;
  468. }
  469. Matrix& MultTranspose(const Vector & vec, Matrix& result);
  470. Matrix MultTranspose(const Vector & vec);
  471. /// Sum up all elements of the vector
  472. inline REALTYPE Sum() const {
  473. REALTYPE result = R_ZERO;
  474. REALTYPE *src = _;
  475. unsigned int len = row;
  476. while(len--)
  477. result += *(src++);
  478. return result;
  479. }
  480. /// The norm of the vector
  481. inline REALTYPE Norm() const {
  482. #ifdef MATHLIB_USE_DOUBLE_AS_REAL
  483. return sqrt(Norm2());
  484. #else
  485. return sqrtf(Norm2());
  486. #endif
  487. }
  488. /// The squared norm of the vector
  489. inline REALTYPE Norm2() const {
  490. REALTYPE result = R_ZERO;
  491. REALTYPE *src = _;
  492. unsigned int len = row;
  493. while(len--){
  494. result += *(src) * (*(src));
  495. src++;
  496. }
  497. return result;
  498. }
  499. /// Normalize the vector to 1
  500. inline void Normalize(){
  501. REALTYPE norm = Norm();
  502. if(norm>EPSILON){
  503. REALTYPE scalar = R_ONE / Norm();
  504. (*this) *= scalar;
  505. }else{
  506. Zero();
  507. }
  508. }
  509. /// The distance between two vectors
  510. inline REALTYPE Distance(const Vector &vector) const{
  511. return (*this-vector).Norm();
  512. }
  513. /// The squared distance between two vectors
  514. inline REALTYPE Distance2(const Vector &vector) const{
  515. return (*this-vector).Norm2();
  516. }
  517. /// The dot product with another vector
  518. inline REALTYPE Dot(const Vector &vector) const{
  519. unsigned int k = (row<=vector.row?row:vector.row);
  520. REALTYPE result = R_ZERO;
  521. REALTYPE *src = _;
  522. REALTYPE *dst = vector._;
  523. while(k--)
  524. result += (*(src++)) * (*(dst++));
  525. return result;
  526. }
  527. /// The maximum value of the vector
  528. inline REALTYPE Max(unsigned int *index=NULL){
  529. if(row==0){
  530. if(index) *index = 0;
  531. return R_ZERO;
  532. }
  533. REALTYPE *src =_;
  534. REALTYPE res =*(src++);
  535. if(index){
  536. *index = 0;
  537. for(unsigned int i=1;i<row;i++){
  538. if(*src>res){
  539. res = *src;
  540. *index = i;
  541. }
  542. src++;
  543. }
  544. }else{
  545. for(unsigned int i=1;i<row;i++){
  546. if(*src>res)
  547. res = *src;
  548. src++;
  549. }
  550. }
  551. return res;
  552. }
  553. /// The minimum value of the vector
  554. inline REALTYPE Min(unsigned int *index=NULL){
  555. if(row==0){
  556. if(index) *index = 0;
  557. return R_ZERO;
  558. }
  559. REALTYPE *src =_;
  560. REALTYPE res =*(src++);
  561. if(index){
  562. *index = 0;
  563. for(unsigned int i=1;i<row;i++){
  564. if(*src<res){
  565. res = *src;
  566. *index = i;
  567. }
  568. src++;
  569. }
  570. }else{
  571. for(unsigned int i=1;i<row;i++){
  572. if(*src<res)
  573. res = *src;
  574. src++;
  575. }
  576. }
  577. return res;
  578. }
  579. /// The index of the maximum value of the vector
  580. inline int MaxId(){
  581. unsigned int res;
  582. Max(&res);
  583. return res;
  584. }
  585. /// The index of the minimum value of the vector
  586. inline int MinId(){
  587. unsigned int res;
  588. Min(&res);
  589. return res;
  590. }
  591. /// The maximum value of two vector
  592. inline Vector& Max(const Vector& v0,const Vector& v1){
  593. const unsigned int kmin = (v0.row<v1.row?v0.row:v1.row);
  594. const unsigned int kmax = (v0.row<v1.row?v1.row:v0.row);
  595. Resize(kmax,false);
  596. for(unsigned int i=0;i<kmin;i++){
  597. _[i] = (v0._[i]>v1._[i]?v0._[i]:v1._[i]);
  598. }
  599. if(kmin<kmax){
  600. if(v0.row<v1.row){
  601. for(unsigned int i=kmin;i<kmax;i++){
  602. _[i] = v1._[i];
  603. }
  604. }else{
  605. for(unsigned int i=kmin;i<kmax;i++){
  606. _[i] = v0._[i];
  607. }
  608. }
  609. }
  610. return *this;
  611. }
  612. /// The minimum value of two vector
  613. inline Vector& Min(const Vector& v0,const Vector& v1){
  614. const unsigned int kmin = (v0.row<v1.row?v0.row:v1.row);
  615. const unsigned int kmax = (v0.row<v1.row?v1.row:v0.row);
  616. Resize(kmax,false);
  617. for(unsigned int i=0;i<kmin;i++){
  618. _[i] = (v0._[i]<v1._[i]?v0._[i]:v1._[i]);
  619. }
  620. if(kmin<kmax){
  621. if(v0.row<v1.row){
  622. for(unsigned int i=kmin;i<kmax;i++){
  623. _[i] = v1._[i];
  624. }
  625. }else{
  626. for(unsigned int i=kmin;i<kmax;i++){
  627. _[i] = v0._[i];
  628. }
  629. }
  630. }
  631. return *this;
  632. }
  633. /// Return the absolute value of the vector
  634. inline Vector Abs(){
  635. Vector result;
  636. return Abs(result);
  637. }
  638. /// Set the result to the absolute value of the vector
  639. inline Vector& Abs(Vector &result) const{
  640. result.Resize(row,false);
  641. REALTYPE *src = _;
  642. REALTYPE *dst = result._;
  643. unsigned int len = row;
  644. while(len--)
  645. *(dst++) = fabs(*(src++));
  646. return result;
  647. }
  648. /// Absolute value of the vector
  649. inline Vector& SAbs() {
  650. REALTYPE *src = _;
  651. unsigned int len = row;
  652. while(len--){
  653. *src = fabs(*src);
  654. src++;
  655. }
  656. return *this;
  657. }
  658. inline Vector& Sort(IndicesVector * indices=NULL){
  659. if(indices){
  660. indices->resize(row);
  661. for(unsigned int i=0;i<row;i++) indices->at(i)=i;
  662. }
  663. REALTYPE cmax;
  664. unsigned int maxId;
  665. for(unsigned int i=0;i<row-1;i++){
  666. cmax = _[i];
  667. maxId = i;
  668. for(unsigned int j=i+1;j<row;j++){
  669. if(cmax<_[j]){
  670. cmax = _[j];
  671. maxId = j;
  672. }
  673. }
  674. if(maxId!=i){
  675. REALTYPE tmp = _[i];
  676. _[i] = _[maxId];
  677. _[maxId] = tmp;
  678. if(indices){
  679. unsigned int idx = indices->at(i);
  680. indices->at(i) = indices->at(maxId);
  681. indices->at(maxId) = idx;
  682. }
  683. }
  684. }
  685. return *this;
  686. }
  687. inline Vector& AbsSort(IndicesVector * indices=NULL){
  688. if(indices){
  689. indices->resize(row);
  690. for(unsigned int i=0;i<row;i++) indices->at(i)=i;
  691. }
  692. REALTYPE cmax;
  693. unsigned int maxId;
  694. for(unsigned int i=0;i<row-1;i++){
  695. cmax = fabs(_[i]);
  696. maxId = i;
  697. for(unsigned int j=i+1;j<row;j++){
  698. if(cmax<fabs(_[j])){
  699. cmax = fabs(_[j]);
  700. maxId = j;
  701. }
  702. }
  703. if(maxId!=i){
  704. REALTYPE tmp = _[i];
  705. _[i] = _[maxId];
  706. _[maxId] = tmp;
  707. if(indices){
  708. unsigned int idx = indices->at(i);
  709. indices->at(i) = indices->at(maxId);
  710. indices->at(maxId) = idx;
  711. }
  712. }
  713. }
  714. return *this;
  715. }
  716. /**
  717. * \brief Set the value of another vector into the current vector
  718. * \param startPos The index starting from which the data of the passed vector will be copied
  719. * \param vector The input vector
  720. */
  721. inline Vector& SetSubVector(unsigned int startPos, const Vector &vector)
  722. {
  723. if(startPos<row){
  724. const unsigned int k = (row-startPos<=vector.row?row-startPos:vector.row);
  725. for (unsigned int i = 0; i < k; i++){
  726. _[startPos+i] = vector._[i];
  727. }
  728. }
  729. return *this;
  730. }
  731. /**
  732. * \brief Get a vector containing some of the vector values
  733. * \param startPos The starting index
  734. * \param len The length of data to be copied
  735. */
  736. inline Vector GetSubVector(unsigned int startPos, unsigned int len) const
  737. {
  738. Vector result(len,false);
  739. return GetSubVector(startPos,len,result);
  740. }
  741. /**
  742. * \brief Get a vector containing some of the vector values
  743. * \param startPos The starting index
  744. * \param len The length of data to be copied
  745. * \param result The output vector
  746. */
  747. inline Vector& GetSubVector(unsigned int startPos, unsigned int len, Vector &result) const
  748. {
  749. result.Resize(len,false);
  750. if(startPos<row){
  751. const unsigned int k = (row-startPos<=len?row-startPos:len);
  752. for (unsigned int i = 0; i < k; i++){
  753. result[i] = _[startPos+i];
  754. }
  755. for (unsigned int i = k; i < len; i++){
  756. result[i] = R_ZERO;
  757. }
  758. }else{
  759. result.Zero();
  760. }
  761. return result;
  762. }
  763. inline Vector& GetSubVector(const Vector& ids, Vector &result) const
  764. {
  765. int s = ids.Size();
  766. IndicesVector id;
  767. for(int i=0;i<s;i++) id.push_back(int(ROUND(ids.At(i))));
  768. return GetSubVector(id,result);
  769. }
  770. /**
  771. * \brief Get a vector containing the vector values a given indices
  772. * \param ids A vecotr containing the indices
  773. * \param result The output vector
  774. */
  775. inline Vector& GetSubVector(const IndicesVector& ids, Vector &result) const
  776. {
  777. const unsigned int k=ids.size();
  778. result.Resize(k);
  779. for(unsigned int i=0;i<k;i++){
  780. const unsigned int g = ids[i];
  781. if(g<row){
  782. result._[i] = _[g];
  783. }else{
  784. result._[i] = R_ZERO;
  785. }
  786. }
  787. return result;
  788. }
  789. /**
  790. * \brief Set at given indices a given vector values
  791. * \param ids The indices
  792. * \param result The vector
  793. */
  794. inline Vector& SetSubVector(const IndicesVector& ids, const Vector &source)
  795. {
  796. const unsigned int j=ids.size();
  797. const unsigned int k= (source.row<j?source.row:j);
  798. for(unsigned int i=0;i<k;i++){
  799. const unsigned int g = ids[i];
  800. if(g<row){
  801. _[g] = source._[i];
  802. }
  803. }
  804. return *this;
  805. }
  806. inline Vector& InsertSubVector(unsigned int start,
  807. const Vector& vector,
  808. unsigned int vectorStart, unsigned int vectorLength){
  809. if(vectorStart >= vector.row) return *this;
  810. if(start >= row) return *this;
  811. if(vectorStart+vectorLength > vector.row) vectorLength = vector.row-vectorStart;
  812. if(start+vectorLength > row) vectorLength = row-start;
  813. unsigned int rowOffset = start;
  814. unsigned int vectorRowOffset = vectorStart;
  815. for(unsigned int j=0;j<vectorLength;j++){
  816. _[rowOffset] = vector._[vectorRowOffset];
  817. rowOffset++;
  818. vectorRowOffset++;
  819. }
  820. return *this;
  821. }
  822. /// Shift the value to the right (rightmost value is set to the left)
  823. inline Vector& ShiftRight(){
  824. if(row>1){
  825. REALTYPE *src = _ + row-1;
  826. REALTYPE zero = *src;
  827. unsigned int len = row-1;
  828. while(len--){
  829. *src = *(src-1);
  830. src--;
  831. }
  832. *src = zero;
  833. }
  834. return *this;
  835. }
  836. /// Shift the value toi the left (leftmost value is set to the right)
  837. inline Vector& ShiftLeft(){
  838. if(row>1){
  839. REALTYPE *src = _;
  840. REALTYPE zero = *src;
  841. unsigned int len = row-1;
  842. while(len--){
  843. *src = *(src+1);
  844. src++;
  845. }
  846. *src = zero;
  847. }
  848. return *this;
  849. }
  850. inline Vector& Trunc(const Vector& min, const Vector& max){
  851. unsigned int k = MIN(min.row,max.row);
  852. unsigned int len = MIN(row,k);
  853. REALTYPE *dst = _;
  854. REALTYPE *v1 = min._;
  855. REALTYPE *v2 = max._;
  856. while(len--){
  857. *dst = TRUNC((*dst),(*(v1)),(*(v2)));
  858. v1++;
  859. v2++;
  860. dst++;
  861. }
  862. if(row>k){
  863. unsigned int len = row-k;
  864. while(len--)
  865. *(dst++) = R_ZERO;
  866. }
  867. return *this;
  868. }
  869. inline Vector& Trunc(const REALTYPE min, const REALTYPE max){
  870. REALTYPE *dst = _;
  871. unsigned int len = row;
  872. while(len--){
  873. *dst = TRUNC((*dst),min,max);
  874. dst++;
  875. }
  876. return *this;
  877. }
  878. /// Print the vector to stdout
  879. void Print() const;
  880. void Print(string name) const;
  881. /*{
  882. std::cout << "Vector " <<row<<std::endl;;
  883. for (unsigned int i = 0; i < row; i++)
  884. std::cout << _[i] <<" ";
  885. std::cout << "\n";
  886. }*/
  887. /// Print the vector to stdout
  888. friend std::ostream & operator<<(std::ostream& out, const Vector & a){
  889. PRINT_BEGIN(out);
  890. for (unsigned int i = 0; i < a.Size(); i++){
  891. out.width(PRINT_WIDTH);
  892. out<< a.AtNoCheck(i)<<" " ;
  893. }
  894. PRINT_END(out);
  895. return out;
  896. }
  897. int PrintToString(char * str, int maxsize=0){
  898. int cIndex = 0;
  899. str[0] = 0;
  900. if(maxsize<=0){
  901. for (unsigned int i = 0; i < row; i++){
  902. cIndex += sprintf(str+cIndex,"%1.12f ",_[i]);
  903. }
  904. }else{
  905. for (unsigned int i = 0; i < row; i++){
  906. int nb = snprintf(str+cIndex,maxsize-cIndex,"%1.12f ",_[i]);
  907. if(nb>=maxsize-cIndex) break;
  908. cIndex += nb;
  909. }
  910. }
  911. return cIndex;
  912. }
  913. protected:
  914. inline void Release(){
  915. if(_!=NULL) delete [] _;
  916. row = 0;
  917. _ = NULL;
  918. }
  919. public:
  920. /**
  921. * \brief Resize the vector
  922. * \param size The new size
  923. * \param copy Should keep the original data or just resize.
  924. */
  925. inline virtual Vector& Resize(unsigned int size, bool copy = true){
  926. if(row!=size){
  927. if(size){
  928. REALTYPE *arr = new REALTYPE[size];
  929. if(copy){
  930. unsigned int len = (row<size?row:size);
  931. memcpy(arr,_,len*sizeof(REALTYPE));
  932. if(len<size){
  933. memset(arr+len,0,(size-len)*sizeof(REALTYPE));
  934. }
  935. }
  936. if(_!=NULL) delete [] _;
  937. _ = arr;
  938. row = size;
  939. }else{
  940. Release();
  941. }
  942. }
  943. return *this;
  944. }
  945. };
  946. class SharedVector : public Vector
  947. {
  948. protected:
  949. unsigned int maxMemorySize;
  950. public:
  951. inline SharedVector():Vector(){}
  952. inline SharedVector(REALTYPE *array, unsigned int rowSize):Vector(){
  953. SetSharedPtr(array,rowSize);
  954. }
  955. inline SharedVector(const Vector & vector):Vector(){
  956. SetSharedPtr(vector.Array(),vector.Size());
  957. }
  958. #ifdef USE_T_EXTENSIONS
  959. /// Copy Contructor of a template Vector (see TVector)
  960. template<unsigned int ROW> inline SharedVector(const TVector<ROW> &vector):Vector(){
  961. SetSharedPtr(vector.Array(),ROW);
  962. }
  963. #endif
  964. inline virtual ~SharedVector(){
  965. Release();
  966. }
  967. inline SharedVector& SetSharedPtr(REALTYPE *array, unsigned int size){
  968. row = size;
  969. maxMemorySize = size;
  970. _ = array;
  971. return *this;
  972. }
  973. protected:
  974. inline virtual void Release(){
  975. _ = NULL;
  976. }
  977. public:
  978. inline virtual SharedVector& Resize(unsigned int rowSize, bool copy = true){
  979. assert(rowSize<=maxMemorySize);
  980. row = rowSize;
  981. return *this;
  982. }
  983. };
  984. #ifdef USE_MATHLIB_NAMESPACE
  985. }
  986. #endif
  987. #endif