Changeset 1027 for trunk/yat/utility/VectorMutable.cc
 Timestamp:
 Feb 2, 2008, 10:29:29 PM (14 years ago)
 File:

 1 copied
Legend:
 Unmodified
 Added
 Removed

trunk/yat/utility/VectorMutable.cc
r1026 r1027 25 25 */ 26 26 27 #include "Vector Base.h"27 #include "VectorMutable.h" 28 28 #include "matrix.h" 29 29 #include "utility.h" … … 43 43 44 44 45 Vector Base::VectorBase(void)46 : vec_(NULL), const_vec_(NULL)45 VectorMutable::VectorMutable(void) 46 : VectorBase(NULL), vec_(NULL) 47 47 { 48 48 } 49 49 50 50 51 Vector Base::VectorBase(gsl_vector* v)52 : vec_(v), const_vec_(v)51 VectorMutable::VectorMutable(gsl_vector* v) 52 : VectorBase(v), vec_(v) 53 53 { 54 54 } 55 55 56 56 57 Vector Base::VectorBase(const gsl_vector* v)58 : vec_(NULL), const_vec_(v)57 VectorMutable::VectorMutable(const gsl_vector* v) 58 : VectorBase(v), vec_(NULL) 59 59 { 60 60 } 61 61 62 62 63 Vector Base::~VectorBase(void)63 VectorMutable::~VectorMutable(void) 64 64 { 65 65 } 66 66 67 67 68 void Vector Base::all(double value)68 void VectorMutable::all(double value) 69 69 { 70 70 assert(vec_); … … 73 73 74 74 75 VectorBase::const_iterator VectorBase::begin(void) const 76 { 77 return const_iterator(*this, 0); 78 } 79 80 81 VectorBase::iterator VectorBase::begin(void) 75 VectorMutable::iterator VectorMutable::begin(void) 82 76 { 83 77 return iterator(*this, 0); … … 85 79 86 80 87 void Vector Base::div(const VectorBase& other)81 void VectorMutable::div(const VectorBase& other) 88 82 { 89 83 assert(vec_); 90 84 int status=gsl_vector_div(vec_,other.gsl_vector_p()); 91 85 if (status) 92 throw utility::GSL_error(std::string("Vector Base::div",status));86 throw utility::GSL_error(std::string("VectorMutable::div",status)); 93 87 } 94 88 95 89 96 VectorBase::const_iterator VectorBase::end(void) const 97 { 98 return const_iterator(*this, size()); 99 } 100 101 102 VectorBase::iterator VectorBase::end(void) 90 VectorMutable::iterator VectorMutable::end(void) 103 91 { 104 92 return iterator(*this, size()); … … 106 94 107 95 108 bool VectorBase::equal(const VectorBase& other, const double d) const 109 { 110 if (this==&other) 111 return true; 112 if (size()!=other.size()) 113 return false; 114 // if gsl error handler disabled, out of bounds index will not 115 // abort the program. 116 for (size_t i=0; i<size(); ++i) 117 if (fabs( (*this)(i)other(i) ) > d  118 std::isnan((*this)(i))  std::isnan(other(i)) ) 119 return false; 120 return true; 121 } 122 123 124 const gsl_vector* VectorBase::gsl_vector_p(void) const 125 { 126 return const_vec_; 127 } 128 129 130 gsl_vector* VectorBase::gsl_vector_p(void) 96 gsl_vector* VectorMutable::gsl_vector_p(void) 131 97 { 132 98 return vec_; … … 134 100 135 101 136 void Vector Base::mul(const VectorBase& other)102 void VectorMutable::mul(const VectorBase& other) 137 103 { 138 104 assert(vec_); 139 105 int status=gsl_vector_mul(vec_,other.gsl_vector_p()); 140 106 if (status) 141 throw utility::GSL_error(std::string("Vector Base::div",status));107 throw utility::GSL_error(std::string("VectorMutable::div",status)); 142 108 } 143 109 144 110 145 size_t VectorBase::size(void) const111 double& VectorMutable::operator()(size_t i) 146 112 { 147 if (!const_vec_) 148 return 0; 149 return const_vec_>size; 150 } 151 152 153 const double& VectorBase::operator()(size_t i) const 154 { 155 const double* d=gsl_vector_const_ptr(const_vec_, i); 113 double* d=gsl_vector_ptr(vec_, i); 156 114 if (!d) 157 throw utility::GSL_error("Vector Base::operator()",GSL_EINVAL);115 throw utility::GSL_error("VectorMutable::operator()",GSL_EINVAL); 158 116 return *d; 159 117 } 160 118 161 119 162 double& VectorBase::operator()(size_t i) 163 { 164 double* d=gsl_vector_ptr(vec_, i); 165 if (!d) 166 throw utility::GSL_error("VectorBase::operator()",GSL_EINVAL); 167 return *d; 168 } 169 170 171 double& VectorBase::operator[](size_t i) 120 double& VectorMutable::operator[](size_t i) 172 121 { 173 122 return this>operator()(i); … … 175 124 176 125 177 const double& VectorBase::operator[](size_t i) const 178 { 179 return this>operator()(i); 180 } 181 182 183 bool VectorBase::operator==(const VectorBase& other) const 184 { 185 return equal(other); 186 } 187 188 189 bool VectorBase::operator!=(const VectorBase& other) const 190 { 191 return !equal(other); 192 } 193 194 195 double VectorBase::operator*( const VectorBase &other ) const 196 { 197 double res = 0.0;; 198 for ( size_t i = 0; i < size(); ++i ) 199 res += other(i) * (*this)(i); 200 return res; 201 } 202 203 204 const VectorBase& VectorBase::operator+=(double d) 126 const VectorMutable& VectorMutable::operator+=(double d) 205 127 { 206 128 assert(vec_); … … 210 132 211 133 212 const Vector Base& VectorBase::operator=(const VectorBase& other)134 const VectorMutable& VectorMutable::operator=(const VectorBase& other) 213 135 { 214 136 assert(vec_); 215 137 int status=gsl_vector_sub(vec_, other.gsl_vector_p()); 216 138 if (status) 217 throw utility::GSL_error(std::string("Vector Base::sub", status));139 throw utility::GSL_error(std::string("VectorMutable::sub", status)); 218 140 return *this; 219 141 } 220 142 221 143 222 const Vector Base& VectorBase::operator=(const double d)144 const VectorMutable& VectorMutable::operator=(const double d) 223 145 { 224 146 assert(vec_); … … 228 150 229 151 230 const Vector Base& VectorBase::operator*=(double d)152 const VectorMutable& VectorMutable::operator*=(double d) 231 153 { 232 154 assert(vec_); … … 236 158 237 159 238 bool isnull(const VectorBase& v)160 void shuffle(VectorMutable& invec) 239 161 { 240 r eturn gsl_vector_isnull(v.gsl_vector_p());162 random::random_shuffle(invec.begin(), invec.end()); 241 163 } 242 164 243 165 244 double max(const VectorBase& v)166 void sort(VectorMutable& invec) 245 167 { 246 return gsl_vector_max(v.gsl_vector_p());168 std::sort(invec.begin(), invec.end()); 247 169 } 248 170 249 171 250 size_t max_index(const VectorBase& v)172 VectorMutable::operator proxy() 251 173 { 252 return gsl_vector_max_index(v.gsl_vector_p()); 253 } 254 255 256 double min(const VectorBase& v) 257 { 258 return gsl_vector_min(v.gsl_vector_p()); 259 } 260 261 262 size_t min_index(const VectorBase& v) 263 { 264 return gsl_vector_min_index(v.gsl_vector_p()); 265 } 266 267 268 bool nan(const VectorBase& templat, vector& flag) 269 { 270 size_t vsize(templat.size()); 271 flag = vector(vsize, 1.0); 272 bool nan=false; 273 for (size_t i=0; i<vsize; i++) 274 if (std::isnan(templat(i))) { 275 flag(i)=0; 276 nan=true; 277 } 278 return nan; 279 } 280 281 282 void shuffle(VectorBase& invec) 283 { 284 random::DiscreteUniform rnd; 285 std::random_shuffle(invec.begin(), invec.end(), rnd); 286 } 287 288 289 void sort_index(std::vector<size_t>& sort_index, const VectorBase& invec) 290 { 291 assert(invec.gsl_vector_p()); 292 gsl_permutation* p = gsl_permutation_alloc(invec.size()); 293 int status=gsl_sort_vector_index(p,invec.gsl_vector_p()); 294 if (status) { 295 gsl_permutation_free(p); 296 throw utility::GSL_error(std::string("sort_index(vector&,const VectorBase&)",status)); 297 } 298 sort_index=std::vector<size_t>(p>data,p>data+p>size); 299 gsl_permutation_free(p); 300 } 301 302 303 void sort_smallest_index(std::vector<size_t>& sort_index, size_t k, 304 const VectorBase& invec) 305 { 306 assert(invec.gsl_vector_p()); 307 assert(k<=invec.size()); 308 sort_index.resize(k); 309 gsl_sort_vector_smallest_index(&sort_index[0],k,invec.gsl_vector_p()); 310 } 311 312 void sort_largest_index(std::vector<size_t>& sort_index, size_t k, 313 const VectorBase& invec) 314 { 315 assert(invec.gsl_vector_p()); 316 assert(k<=invec.size()); 317 sort_index.resize(k); 318 gsl_sort_vector_largest_index(&sort_index[0],k,invec.gsl_vector_p()); 319 } 320 321 322 double sum(const VectorBase& v) 323 { 324 double sum = 0; 325 size_t vsize=v.size(); 326 for (size_t i=0; i<vsize; ++i) 327 sum += v(i); 328 return sum; 329 } 330 331 332 std::ostream& operator<<(std::ostream& s, const VectorBase& a) 333 { 334 s.setf(std::ios::dec); 335 s.precision(12); 336 for (size_t j = 0; j < a.size(); ++j) { 337 s << a(j); 338 if ( (j+1)<a.size() ) 339 s << s.fill(); 340 } 341 return s; 174 proxy p; 175 p.vec_ = vec_; 176 vec_ = NULL; // proxy takes ownership and delivers to its receiver 177 return p; 342 178 } 343 179
Note: See TracChangeset
for help on using the changeset viewer.