00001
00015 #ifndef DLR_NUMERIC_UTILITIES_H
00016 #define DLR_NUMERIC_UTILITIES_H
00017
00018 #include <cmath>
00019
00020 #include <cstdlib>
00021
00022 #include <dlrCommon/physicalConstants.h>
00023 #include <dlrCommon/tags.h>
00024 #include <dlrNumeric/numericTraits.h>
00025 #include <dlrNumeric/array1D.h>
00026 #include <dlrNumeric/array2D.h>
00027 #include <dlrNumeric/array3D.h>
00028 #include <dlrNumeric/index2D.h>
00029 #include <dlrNumeric/solveQuadratic.h>
00030 #include <dlrNumeric/vector2D.h>
00031 #include <dlrNumeric/vector3D.h>
00032
00033
00034 namespace dlr {
00035
00036 namespace numeric {
00037
00038 namespace constants {
00039
00040 using dlr::common::constants::pi;
00041
00042 }
00043
00044
00057 template <class Type>
00058 Array1D<Type>
00059 abs(const Array1D<Type>& array0);
00060
00061
00074 template <class Type>
00075 Array2D<Type>
00076 abs(const Array2D<Type>& array0);
00077
00078
00088 template <class Type>
00089 bool
00090 allFalse(const Array1D<Type>& array0);
00091
00092
00102 template <class Type>
00103 bool
00104 allTrue(const Array1D<Type>& array0);
00105
00106
00116 template <class Type>
00117 inline bool
00118 anyFalse(const Array1D<Type>& array0);
00119
00120
00130 template <class Type>
00131 inline bool
00132 anyTrue(const Array1D<Type>& array0);
00133
00134
00145 template <class Type>
00146 inline size_t
00147 argmax(const Array1D<Type>& array0);
00148
00149
00163 template <class IterType>
00164 inline size_t
00165 argmax(IterType beginIter, IterType endIter);
00166
00167
00189 template <class Type, class Functor>
00190 inline size_t
00191 argmax(const Array1D<Type>& array0, Functor comparator);
00192
00193
00218 template <class IterType, class Functor>
00219 inline size_t
00220 argmax(IterType beginIter, IterType endIter, Functor comparator);
00221
00222
00232 template <class Type>
00233 inline Index2D
00234 argmax2D(const Array2D<Type>& array0);
00235
00236
00255 template <class Type, class Functor>
00256 inline Index2D
00257 argmax2D(const Array2D<Type>& array0, Functor comparator);
00258
00259
00270 template <class Type>
00271 inline size_t
00272 argmin(const Array1D<Type>& array0);
00273
00274
00292 template <class Type, class Functor>
00293 inline size_t
00294 argmin(const Array1D<Type>& array0, Functor comparator);
00295
00296
00306 template <class Type>
00307 inline Index2D
00308 argmin2D(const Array2D<Type>& array0);
00309
00310
00325 template <class Type, class Functor>
00326 inline Index2D
00327 argmin2D(const Array2D<Type>& array0, Functor comparator);
00328
00329
00343 template <class Type>
00344 Array1D<size_t>
00345 argsort(const Array1D<Type>& array0);
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367
00377 template <class Type>
00378 Array1D<Type>
00379 axisMax(const Array2D<Type>& array0,
00380 size_t axis) {return axisMaximum(array0, axis);}
00381
00398 template <class Type>
00399 inline Array1D<Type>
00400 axisMaximum(const Array2D<Type>& array0, size_t axis);
00401
00432 template <class Type, class Functor>
00433 Array1D<Type>
00434 axisMaximum(const Array2D<Type>& array0, size_t axis, Functor comparator);
00435
00445 template <class Type>
00446 inline Array1D<Type>
00447 axisMin(const Array2D<Type>& array0,
00448 size_t axis) {return axisMinimum(array0, axis);}
00449
00466 template <class Type>
00467 inline Array1D<Type>
00468 axisMinimum(const Array2D<Type>& array0, size_t axis);
00469
00500 template <class Type, class Functor>
00501 Array1D<Type>
00502 axisMinimum(const Array2D<Type>& array0, size_t axis, Functor comparator);
00503
00528 template <class Type>
00529 inline Array1D<typename NumericTraits<Type>::SumType>
00530 axisSum(const Array2D<Type>& array0, size_t axis);
00531
00554 template <class Type, class ResultType>
00555 inline Array1D<ResultType>
00556 axisSum(const Array2D<Type>& array0, size_t axis,
00557 type_tag<ResultType> resultTag);
00558
00598 template <class Type, class ResultType, class Functor>
00599 Array1D<ResultType>
00600 axisSum(const Array2D<Type>& array0, size_t axis,
00601 type_tag<ResultType> resultTag,
00602 const ResultType& initialValue, Functor adder);
00603
00631 template <class Type>
00632 Array2D<Type>
00633 columnIndices(size_t rows, size_t columns, type_tag<Type> typeTag);
00634
00666 template <class Type0, class Type1>
00667 inline Array1D<Type1>
00668 compress(const Array1D<Type0>& condition, const Array1D<Type1>& input);
00669
00696 template <class Type0, class Type1>
00697 Array1D<Type1>
00698 compress(const Array1D<Type0>& condition, const Array1D<Type1>& input,
00699 size_t numTrue);
00700
00701
00712 template <class Type>
00713 inline size_t
00714 count(const Array1D<Type>& array0);
00715
00729 inline Vector3D
00730 cross(const Vector3D& vector0, const Vector3D& vector1);
00731
00743 template <class Type>
00744 inline typename NumericTraits<Type>::ProductType
00745 dot(const Array1D<Type>& array0, const Array1D<Type>& array1);
00746
00761 template <class Type0, class Type1, class Type2>
00762 inline Type2
00763 dot(const Array1D<Type0>& array0, const Array1D<Type1>& array1,
00764 type_tag<Type2> typeTag);
00765
00775 inline double
00776 dot(const Vector2D& vector0, const Vector2D& vector1);
00777
00787 inline double
00788 dot(const Vector3D& vector0, const Vector3D& vector1);
00789
00811 template <class Type>
00812 Array2D<Type>
00813 equivalentMatrix(const Array1D<Type>& vector0, size_t rowsInMatrix);
00814
00815
00838 template <class Type>
00839 void
00840 getMeanAndCovariance(const Array2D<Type>& sampleArray,
00841 Array1D<double>& mean,
00842 Array2D<double>& covariance,
00843 size_t majorAxis=0);
00844
00845
00863 template <class Type>
00864 inline Array2D<Type>
00865 identity(size_t rows, size_t columns, type_tag<Type> typeTag);
00866
00867
00882 template <class Type>
00883 Array2D<Type>
00884 identity(size_t rows, size_t columns);
00885
00886
00896 template <class Type>
00897 Array1D<Type>
00898 ln(const Array1D<Type>& array0);
00899
00900
00910 template <class Type>
00911 Array2D<Type>
00912 ln(const Array2D<Type>& array0);
00913
00914
00935 template <class Type>
00936 Array2D<Type>
00937 logicalNot(const Array2D<Type>& array0);
00938
00948 template <class Type>
00949 inline Type
00950 magnitude(const Array1D<Type>& array0);
00951
00961 inline double
00962 magnitude(const Vector2D& vector0);
00963
00973 inline double
00974 magnitude(const Vector3D& vector0);
00975
00976
00987 template <class Type>
00988 inline typename NumericTraits<Type>::ProductType
00989 magnitudeSquared(const Array1D<Type>& array0);
00990
00991
01002 inline double
01003 magnitudeSquared(const Vector2D& vector0);
01004
01005
01016 inline double
01017 magnitudeSquared(const Vector3D& vector0);
01018
01019
01037 template <class Type>
01038 inline Array1D<typename NumericTraits<Type>::ProductType>
01039 matrixMultiply(const Array1D<Type>& vector0, const Array2D<Type>& matrix0);
01040
01058 template <class Type0, class Type1, class Type2>
01059 Array1D<Type2>
01060 matrixMultiply(const Array1D<Type0>& vector0, const Array2D<Type0>& matrix0,
01061 type_tag<Type2> typeTag);
01062
01079 template <class Type>
01080 inline Array1D<typename NumericTraits<Type>::ProductType>
01081 matrixMultiply(const Array2D<Type>& matrix0, const Array1D<Type>& vector0);
01082
01100 template <class Type0, class Type1, class Type2>
01101 Array1D<Type2>
01102 matrixMultiply(const Array2D<Type0>& matrix0, const Array1D<Type1>& vector0,
01103 type_tag<Type2> typeTag);
01104
01105
01120 template <class Type>
01121 inline Array2D<typename NumericTraits<Type>::ProductType>
01122 matrixMultiply(const Array2D<Type>& matrix0, const Array2D<Type>& matrix1);
01123
01141 template <class Type0, class Type1, class Type2>
01142 Array2D<Type2>
01143 matrixMultiply(const Array2D<Type0>& matrix0, const Array2D<Type1>& matrix1,
01144 type_tag<Type2> typeTag);
01145
01146
01147
01148
01149
01150
01151
01152
01153
01154
01155
01156
01157
01158
01159
01160
01161
01162
01163
01164
01174 template <class Type>
01175 inline Type
01176 maximum(const Array1D<Type>& array0);
01177
01197 template <class Type, class Functor>
01198 Type
01199 maximum(const Array1D<Type>& array0, Functor comparator);
01200
01201
01213 template <class Iterator, class Type>
01214 Type
01215 mean(const Iterator& beginIter, const Iterator& endIter);
01216
01217
01229 template <class Type>
01230 inline double
01231 mean(const Array1D<Type>& array0);
01232
01233
01250 template <class Type0, class Type1>
01251 inline Type1
01252 mean(const Array1D<Type0>& array0, type_tag<Type1> typeTag);
01253
01254
01255
01256
01257
01258
01259
01260
01261
01262
01263
01264
01265
01266
01267
01268
01269
01270
01271
01272
01282 template <class Type>
01283 inline Type
01284 minimum(const Array1D<Type>& array0);
01285
01305 template <class Type, class Functor>
01306 Type
01307 minimum(const Array1D<Type>& array0, Functor comparator);
01308
01309
01337 template <class FUNCTOR>
01338 double newtonRaphson(double startPoint, FUNCTOR objectiveFunction,
01339 double epsilon, size_t maxIterations);
01340
01362 template <class Type>
01363 inline double
01364 normalizedCorrelation(const Array1D<Type>& signal0,
01365 const Array1D<Type>& signal1);
01366
01367
01385 template <class Type, class Type2>
01386 Type2
01387 normalizedCorrelation(const Array1D<Type>& signal0,
01388 const Array1D<Type>& signal1,
01389 type_tag<Type2> typeTag);
01390
01403 template <class Type>
01404 inline Array1D<Type>
01405 ones(int size, type_tag<Type> typeTag);
01406
01422 template <class Type>
01423 inline Array2D<Type>
01424 ones(int rows, int columns, type_tag<Type> typeTag);
01425
01441 template <class Type>
01442 inline Array2D<typename NumericTraits<Type>::ProductType>
01443 outerProduct(const Array1D<Type>& array0, const Array1D<Type>& array1);
01444
01464 template <class Type0, class Type1, class Type2>
01465 Array2D<Type2>
01466 outerProduct(const Array1D<Type0>& array0, const Array1D<Type1>& array1,
01467 type_tag<Type2> typeTag);
01468
01508 template <class Type>
01509 Array1D<Type>
01510 range(Type start, Type stop, Type stride=1);
01511
01527 template <class Type>
01528 inline Array1D<Type>
01529 ravel(Array2D<Type>& inputArray) {return inputArray.ravel();}
01530
01542 template <class Type>
01543 inline const Array1D<Type>
01544 ravel(const Array2D<Type>& inputArray) {return inputArray.ravel();}
01545
01546
01558 template <class Type>
01559 inline Type
01560 rms(const Array1D<Type>& array0);
01561
01562
01579 template <class Type0, class Type1>
01580 Type1
01581 rms(const Array1D<Type0>& array0, type_tag<Type1> typeTag);
01582
01583
01610 template <class Type>
01611 Array2D<Type>
01612 rowIndices(size_t rows, size_t columns, type_tag<Type> typeTag);
01613
01614
01629 template <class Type0, class Type1>
01630 inline bool
01631 shapeMatch(const Array1D<Type0>& array0, const Array1D<Type1>& array1);
01632
01633
01648 template <class Type0, class Type1>
01649 inline bool
01650 shapeMatch(const Array2D<Type0>& array0, const Array2D<Type1>& array1);
01651
01652
01667 template <class Type0, class Type1>
01668 inline bool
01669 shapeMatch(const Array3D<Type0>& array0, const Array3D<Type1>& array1);
01670
01671
01680 template <class Type>
01681 inline Array2D<Type>
01682 skewSymmetric(const Array1D<Type>& vector0);
01683
01684
01743 template <class Type>
01744 void
01745 solveQuadratic(Type c0, Type c1, Type c2,
01746 Type& root0, Type& root1, bool& valid,
01747 bool checkValidity=true);
01748
01749
01767 template <class Type>
01768 inline double
01769 standardDeviation(const Array1D<Type>& array0);
01770
01794 template <class Type0, class Type1>
01795 inline Type1
01796 standardDeviation(const Array1D<Type0>& array0, type_tag<Type1> typeTag);
01797
01798
01808 template <class Type>
01809 inline typename NumericTraits<Type>::SumType
01810 sum(const Array1D<Type>& array0);
01811
01812
01826 template <class Type, class Type2>
01827 Type2
01828 sum(const Array1D<Type>& array0, type_tag<Type2> typeTag);
01829
01830
01851 template <class Type>
01852 inline typename NumericTraits<Type>::SumType
01853 sum(const Array2D<Type>& array0,
01854 const Index2D& upperLeftCorner,
01855 const Index2D& lowerRightCorner);
01856
01857
01882 template <class Type, class Type2>
01883 Type2
01884 sum(const Array2D<Type>& array0,
01885 const Index2D& upperLeftCorner,
01886 const Index2D& lowerRightCorner,
01887 type_tag<Type2> typeTag);
01888
01889
01914 template <class Type, class IntegralType>
01915 Array1D<Type>
01916 take(const Array1D<Type>& dataArray,
01917 const Array1D<IntegralType>& indexArray);
01918
01919
01947 template <class Type, class IntegralType>
01948 inline Array1D<Type>
01949 take(const Array2D<Type>& dataArray,
01950 const Array1D<IntegralType>& indexArray);
01951
01952
01991 template <class Type, class IntegralType>
01992 Array2D<Type>
01993 take(const Array2D<Type>& dataArray,
01994 const Array1D<IntegralType>& indexArray,
01995 unsigned int axis);
01996
01997
02015 template <class Type>
02016 inline double
02017 variance(const Array1D<Type>& array0);
02018
02041 template <class Type0, class Type1>
02042 inline Type1
02043 variance(const Array1D<Type0>& array0, type_tag<Type1> typeTag);
02044
02057 template <class Type>
02058 inline Array1D<Type>
02059 zeros(size_t size, type_tag<Type> typeTag);
02060
02076 template <class Type>
02077 inline Array2D<Type>
02078 zeros(size_t rows, size_t columns, type_tag<Type> typeTag);
02079
02098 template <class Type>
02099 inline Array3D<Type>
02100 zeros(size_t shape0, size_t shape1, size_t shape2, type_tag<Type> typeTag);
02101
02102 }
02103
02104 }
02105
02106
02107
02108
02109 namespace dlr {
02110
02111 using numeric::abs;
02112 using numeric::abs;
02113 using numeric::allFalse;
02114 using numeric::allTrue;
02115 using numeric::anyFalse;
02116 using numeric::anyTrue;
02117 using numeric::argmax;
02118 using numeric::argmin;
02119 using numeric::argsort;
02120 using numeric::axisMax;
02121 using numeric::axisMaximum;
02122 using numeric::axisMin;
02123 using numeric::axisMinimum;
02124 using numeric::axisSum;
02125 using numeric::columnIndices;
02126 using numeric::compress;
02127 using numeric::count;
02128 using numeric::cross;
02129 using numeric::dot;
02130 using numeric::equivalentMatrix;
02131 using numeric::getMeanAndCovariance;
02132 using numeric::identity;
02133 using numeric::ln;
02134 using numeric::logicalNot;
02135 using numeric::magnitude;
02136 using numeric::magnitudeSquared;
02137 using numeric::matrixMultiply;
02138 using numeric::maximum;
02139 using numeric::mean;
02140 using numeric::minimum;
02141 using numeric::newtonRaphson;
02142 using numeric::normalizedCorrelation;
02143 using numeric::ones;
02144 using numeric::outerProduct;
02145 using numeric::range;
02146 using numeric::ravel;
02147 using numeric::rms;
02148 using numeric::rowIndices;
02149 using numeric::shapeMatch;
02150 using numeric::skewSymmetric;
02151 using numeric::solveQuadratic;
02152 using numeric::standardDeviation;
02153 using numeric::sum;
02154 using numeric::variance;
02155 using numeric::zeros;
02156
02157 }
02158
02159
02160
02161
02162 #include <cmath>
02163 #include <sstream>
02164 #include <numeric>
02165
02166
02167 namespace dlr {
02168
02169 namespace numeric {
02170
02182 namespace privateCode {
02183
02189 template <class Type>
02190 struct absFunctor : public std::unary_function<Type, Type> {
02198 inline Type operator()(const Type& input) {
02199 DLR_THROW3(NotImplementedException,
02200 "absFunctor<Type>::operator()(const Type&)",
02201 "absFunctor must be specialized for each type.");
02202 return static_cast<Type>(0);
02203 }
02204 };
02205
02209 template <>
02210 struct absFunctor<long double>
02211 : public std::unary_function<long double, long double> {
02218 inline double operator()(const long double& input) {
02219
02220
02221 return (input > 0.0) ? input : -input;
02222 }
02223 };
02224
02225
02229 template <>
02230 struct absFunctor<double> : public std::unary_function<double, double> {
02237 inline double operator()(const double& input) {
02238 return std::fabs(input);
02239 }
02240 };
02241
02242
02246 template <>
02247 struct absFunctor<float> : public std::unary_function<float, float> {
02254 inline float operator()(const float& input) {
02255
02256
02257 return fabsf(input);
02258 }
02259 };
02260
02261
02265 template <>
02266 struct absFunctor<long int>
02267 : public std::unary_function<long int, long int> {
02274 inline long int operator()(const long int& input) {
02275 return std::labs(input);
02276 }
02277 };
02278
02279
02283 template <>
02284 struct absFunctor<long long int>
02285 : public std::unary_function<long long int, long long int> {
02292 inline long long int operator()(const long long int& input) {
02293
02294
02295 return input >= 0LL ? input : -input;
02296 }
02297 };
02298
02299
02303 template <>
02304 struct absFunctor<int> : public std::unary_function<int, int> {
02311 inline int operator()(const int& input) {
02312 return std::abs(input);
02313 }
02314 };
02315
02316
02317 }
02321
02322
02323
02324 template <class Type>
02325 Array1D<Type>
02326 abs(const Array1D<Type>& array0)
02327 {
02328 Array1D<Type> result(array0.size());
02329 std::transform(array0.begin(), array0.end(), result.begin(),
02330 privateCode::absFunctor<Type>());
02331 return result;
02332 }
02333
02334
02335
02336
02337 template <class Type>
02338 Array2D<Type>
02339 abs(const Array2D<Type>& array0)
02340 {
02341 Array2D<Type> result(array0.rows(), array0.columns());
02342 std::transform(array0.begin(), array0.end(), result.begin(),
02343 privateCode::absFunctor<Type>());
02344 return result;
02345 }
02346
02347
02348
02349
02350 template <class Type>
02351 bool
02352 allFalse(const Array1D<Type>& array0)
02353 {
02354 for(typename Array1D<Type>::const_iterator iter = array0.begin();
02355 iter != array0.end();
02356 ++iter) {
02357 if(*iter) {
02358 return false;
02359 }
02360 }
02361 return true;
02362 }
02363
02364
02365
02366
02367 template <class Type>
02368 bool
02369 allTrue(const Array1D<Type>& array0)
02370 {
02371 for(typename Array1D<Type>::const_iterator iter = array0.begin();
02372 iter != array0.end();
02373 ++iter) {
02374 if(!(*iter)) {
02375 return false;
02376 }
02377 }
02378 return true;
02379 }
02380
02381
02382
02383
02384 template <class Type>
02385 inline bool
02386 anyFalse(const Array1D<Type>& array0)
02387 {
02388 return !allTrue(array0);
02389 }
02390
02391
02392
02393
02394 template <class Type>
02395 inline bool
02396 anyTrue(const Array1D<Type>& array0)
02397 {
02398 return !allFalse(array0);
02399 }
02400
02401
02402
02403
02404 template <class Type>
02405 inline size_t
02406 argmax(const Array1D<Type>& array0)
02407 {
02408 return argmax(array0, std::less<Type>());
02409 }
02410
02411
02412
02413
02414 template <class IterType>
02415 inline size_t
02416 argmax(IterType beginIter, IterType endIter)
02417 {
02418 return argmax(beginIter, endIter, std::less<IterType>());
02419 }
02420
02421
02422
02423
02424 template <class Type, class Functor>
02425 inline size_t
02426 argmax(const Array1D<Type>& array0, Functor comparator)
02427 {
02428 return argmax(array0.begin(), array0.end(), comparator);
02429 }
02430
02431
02432
02433
02434 template <class IterType, class Functor>
02435 inline size_t
02436 argmax(IterType beginIter, IterType endIter, Functor comparator)
02437 {
02438 return static_cast<size_t>(
02439 std::max_element(beginIter, endIter, comparator) - beginIter);
02440 }
02441
02442
02443
02444
02445 template <class Type>
02446 inline Index2D
02447 argmax2D(const Array2D<Type>& array0)
02448 {
02449 return argmax2D(array0, std::less<Type>());
02450 }
02451
02452
02453
02454
02455
02456 template <class Type, class Functor>
02457 inline Index2D
02458 argmax2D(const Array2D<Type>& array0, Functor comparator)
02459 {
02460 size_t ravelIndex = argmax(array0.ravel(), comparator);
02461 size_t row = ravelIndex / array0.columns();
02462 size_t column = ravelIndex - row * array0.columns();
02463 return Index2D(row, column);
02464 }
02465
02466
02467
02468
02469 template <class Type>
02470 inline size_t
02471 argmin(const Array1D<Type>& array0)
02472 {
02473 return argmin(array0, std::less<Type>());
02474 }
02475
02476
02477
02478
02479 template <class Type, class Functor>
02480 inline size_t
02481 argmin(const Array1D<Type>& array0, Functor comparator)
02482 {
02483 const Type* minPtr = std::min_element(array0.begin(), array0.end(),
02484 comparator);
02485 return static_cast<size_t>(minPtr - array0.begin());
02486 }
02487
02488
02489
02490
02491 template <class Type>
02492 inline Index2D
02493 argmin2D(const Array2D<Type>& array0)
02494 {
02495 return argmin2D(array0, std::less<Type>());
02496 }
02497
02498
02499
02500
02501
02502 template <class Type, class Functor>
02503 inline Index2D
02504 argmin2D(const Array2D<Type>& array0, Functor comparator)
02505 {
02506 size_t ravelIndex = argmin(array0.ravel(), comparator);
02507 size_t row = ravelIndex / array0.columns();
02508 size_t column = ravelIndex - row * array0.columns();
02509 return Index2D(row, column);
02510 }
02511
02512
02513
02514
02515
02516
02517 template <class Type>
02518 Array1D<size_t>
02519 argsort(const Array1D<Type>& array0)
02520 {
02521 Array1D< std::pair<Type, size_t> > sortVector(array0.size());
02522 Array1D<size_t> resultVector(array0.size());
02523
02524 for(size_t index0 = 0; index0 < array0.size(); ++index0) {
02525 sortVector[index0] = std::pair<Type, size_t>(array0[index0], index0);
02526 }
02527 std::sort(sortVector.begin(), sortVector.end());
02528 std::transform(sortVector.begin(), sortVector.end(), resultVector.begin(),
02529 ExtractSecondFunctor<Type, size_t>());
02530 return resultVector;
02531 }
02532
02533
02534
02535
02536
02537
02538
02539
02540
02541
02542
02543
02544
02545
02546
02547
02548
02549
02550
02551
02552
02553
02554
02555
02556
02557
02558 template <class Type>
02559 inline Array1D<Type>
02560 axisMaximum(const Array2D<Type>& array0, size_t axis)
02561 {
02562 return axisMaximum(array0, axis, std::less<Type>());
02563 }
02564
02565
02566
02567
02568
02569 template <class Type, class Functor>
02570 Array1D<Type>
02571 axisMaximum(const Array2D<Type>& array0, size_t axis, Functor comparator)
02572 {
02573 Array1D<Type> result;
02574 switch(axis) {
02575 case 0:
02576 result.reinit(array0.columns());
02577 for(size_t column = 0; column < array0.columns(); ++column) {
02578 const Type* dataPtr = array0.data(0, column);
02579 Type columnMax = *dataPtr;
02580 size_t stride = array0.columns();
02581 for(size_t row = 0; row < array0.rows(); ++row) {
02582 if(!comparator(*dataPtr, columnMax)) {
02583 columnMax = *dataPtr;
02584 }
02585 dataPtr += stride;
02586 }
02587 result(column) = columnMax;
02588 }
02589 break;
02590 case 1:
02591 result.reinit(array0.rows());
02592 for(size_t row = 0; row < array0.rows(); ++row) {
02593 const Type* dataPtr = array0.data(row, 0);
02594 Type rowMax = *dataPtr;
02595 for(size_t column = 0; column < array0.columns(); ++column) {
02596 if(!comparator(*dataPtr, rowMax)) {
02597 rowMax = *dataPtr;
02598 }
02599 ++dataPtr;
02600 }
02601 result(row) = rowMax;
02602 }
02603 break;
02604 default:
02605 std::ostringstream message;
02606 message << "Axis " << axis << " is invalid for an Array2D.";
02607 DLR_THROW3(IndexException, "axisMaximum(const Array2D&, size_t, ...)",
02608 message.str().c_str());
02609 break;
02610 }
02611 return result;
02612 }
02613
02614
02615
02616
02617
02618 template <class Type>
02619 inline Array1D<Type>
02620 axisMinimum(const Array2D<Type>& array0, size_t axis)
02621 {
02622 return axisMinimum(array0, axis, std::less<Type>());
02623 }
02624
02625
02626
02627
02628
02629 template <class Type, class Functor>
02630 Array1D<Type>
02631 axisMinimum(const Array2D<Type>& array0, size_t axis, Functor comparator)
02632 {
02633 Array1D<Type> result;
02634 switch(axis) {
02635 case 0:
02636 result.reinit(array0.columns());
02637 for(size_t column = 0; column < array0.columns(); ++column) {
02638 const Type* dataPtr = array0.data(0, column);
02639 Type columnMax = *dataPtr;
02640 size_t stride = array0.columns();
02641 for(size_t row = 0; row < array0.rows(); ++row) {
02642 if(comparator(*dataPtr, columnMax)) {
02643 columnMax = *dataPtr;
02644 }
02645 dataPtr += stride;
02646 }
02647 result(column) = columnMax;
02648 }
02649 break;
02650 case 1:
02651 result.reinit(array0.rows());
02652 for(size_t row = 0; row < array0.rows(); ++row) {
02653 const Type* dataPtr = array0.data(row, 0);
02654 Type rowMax = *dataPtr;
02655 for(size_t column = 0; column < array0.columns(); ++column) {
02656 if(comparator(*dataPtr, rowMax)) {
02657 rowMax = *dataPtr;
02658 }
02659 ++dataPtr;
02660 }
02661 result(row) = rowMax;
02662 }
02663 break;
02664 default:
02665 std::ostringstream message;
02666 message << "Axis " << axis << " is invalid for an Array2D.";
02667 DLR_THROW3(IndexException, "axisMinimum(const Array2D&, size_t, ...)",
02668 message.str().c_str());
02669 break;
02670 }
02671 return result;
02672 }
02673
02674
02675
02676
02677 template <class Type>
02678 inline Array1D<typename NumericTraits<Type>::SumType>
02679 axisSum(const Array2D<Type>& array0, size_t axis)
02680 {
02681 return axisSum(array0, axis,
02682 type_tag<typename NumericTraits<Type>::SumType>());
02683 }
02684
02685
02686
02687 template <class Type, class ResultType>
02688 inline Array1D<ResultType>
02689 axisSum(const Array2D<Type>& array0, size_t axis, type_tag<ResultType>)
02690 {
02691 return axisSum(array0, axis, type_tag<ResultType>(),
02692 static_cast<ResultType>(0), std::plus<ResultType>());
02693 }
02694
02695
02696
02697 template <class Type, class ResultType, class Functor>
02698 Array1D<ResultType>
02699 axisSum(const Array2D<Type>& array0, size_t axis, type_tag<ResultType>,
02700 const ResultType& initialValue, Functor adder)
02701 {
02702 Array1D<ResultType> result;
02703 switch(axis) {
02704 case 0:
02705 result.reinit(array0.columns());
02706 for(size_t column = 0; column < array0.columns(); ++column) {
02707 ResultType columnSum = initialValue;
02708 const Type* dataPtr = array0.data(0, column);
02709 size_t stride = array0.columns();
02710 for(size_t row = 0; row < array0.rows(); ++row) {
02711 columnSum = adder(columnSum, *dataPtr);
02712 dataPtr += stride;
02713 }
02714 result(column) = columnSum;
02715 }
02716 break;
02717 case 1:
02718 result.reinit(array0.rows());
02719 for(size_t row = 0; row < array0.rows(); ++row) {
02720 ResultType rowSum = initialValue;
02721 const Type* dataPtr = array0.data(row, 0);
02722 for(size_t column = 0; column < array0.columns(); ++column) {
02723 rowSum = adder(rowSum, *dataPtr);
02724 ++dataPtr;
02725 }
02726 result(row) = rowSum;
02727 }
02728 break;
02729 default:
02730 std::ostringstream message;
02731 message << "Axis " << axis << " is invalid for an Array2D.";
02732 DLR_THROW3(IndexException, "axisSum(const Array2D&, size_t, ...)",
02733 message.str().c_str());
02734 break;
02735 }
02736 return result;
02737 }
02738
02739
02740
02741 template <class Type>
02742 Array2D<Type>
02743 columnIndices(size_t rows, size_t columns, type_tag<Type>)
02744 {
02745 Array2D<Type> returnValue(rows, columns);
02746 Array1D<Type> modelRow = range(static_cast<Type>(0),
02747 static_cast<Type>(columns),
02748 static_cast<Type>(1));
02749 for(size_t row=0; row < rows; ++row) {
02750 std::copy(modelRow.begin(), modelRow.end(), returnValue.data(row, 0));
02751 }
02752 return returnValue;
02753 }
02754
02755
02756
02757
02758 template <class Type0, class Type1>
02759 inline Array1D<Type1>
02760 compress(const Array1D<Type0>& condition,
02761 const Array1D<Type1>& input)
02762 {
02763 size_t numTrue = count(condition);
02764 return compress(condition, input, numTrue);
02765 }
02766
02767
02768
02769
02770 template <class Type0, class Type1>
02771 Array1D<Type1>
02772 compress(const Array1D<Type0>& condition,
02773 const Array1D<Type1>& input,
02774 size_t numTrue)
02775 {
02776 if(condition.size() != input.size()) {
02777 std::ostringstream message;
02778 message << "Condition and input arguments must have the same "
02779 << "size, but condition has size = " << condition.size()
02780 << ", while input has size = " << input.size() << ".";
02781 DLR_THROW3(ValueException,
02782 "compress(const Array1D&, const Array1D&, size_t)",
02783 message.str().c_str());
02784 }
02785 Array1D<Type1> result(numTrue);
02786
02787
02788
02789
02790
02791
02792
02793
02794
02795
02796
02797
02798
02799
02800
02801
02802
02803
02804 typename Array1D<Type1>::const_iterator first = input.begin();
02805 typename Array1D<Type1>::const_iterator last = input.end();
02806 typename Array1D<Type0>::const_iterator mask = condition.begin();
02807 typename Array1D<Type1>::iterator target = result.begin();
02808 while(first != last) {
02809 if(*mask) {
02810 *(target++) = *first;
02811 }
02812 ++mask;
02813 ++first;
02814 }
02815
02816 return result;
02817 }
02818
02819
02820
02821 template <class Type>
02822 inline size_t
02823 count(const Array1D<Type>& x)
02824 {
02825 return std::count_if(x.begin(), x.end(), StaticCastFunctor<Type, bool>());
02826 }
02827
02828
02829
02830 inline Vector3D
02831 cross(const Vector3D& vector0, const Vector3D& vector1) {
02832 return Vector3D((vector0.y() * vector1.z()) - (vector0.z() * vector1.y()),
02833 (vector0.z() * vector1.x()) - (vector0.x() * vector1.z()),
02834 (vector0.x() * vector1.y()) - (vector0.y() * vector1.x()));
02835 }
02836
02837
02838
02839
02840 template <class Type>
02841 inline typename NumericTraits<Type>::ProductType
02842 dot(const Array1D<Type>& x, const Array1D<Type>& y)
02843 {
02844 return dot(x, y, type_tag<typename NumericTraits<Type>::ProductType>());
02845 }
02846
02847
02848
02849
02850 template <class Type0, class Type1, class Type2>
02851 inline Type2
02852 dot(const Array1D<Type0>& x, const Array1D<Type1>& y, type_tag<Type2>)
02853 {
02854 if(x.size() != y.size()) {
02855 std::ostringstream message;
02856 message << "Input arguments must have matching sizes, but x.size() == "
02857 << x.size() << " and y.size() == " << y.size() << "."
02858 << std::endl;
02859 DLR_THROW3(ValueException, "dot()", message.str().c_str());
02860 }
02861 return std::inner_product(x.begin(), x.end(), y.begin(),
02862 static_cast<Type2>(0));
02863 }
02864
02865
02866 inline double
02867 dot(const Vector2D& vector0, const Vector2D& vector1)
02868 {
02869 return vector0.x() * vector1.x() + vector0.y() * vector1.y();
02870 }
02871
02872
02873 inline double
02874 dot(const Vector3D& vector0, const Vector3D& vector1)
02875 {
02876 return (vector0.x() * vector1.x() + vector0.y() * vector1.y()
02877 + vector0.z() * vector1.z());
02878 }
02879
02880
02881 template <class Type>
02882 Array2D<Type>
02883 equivalentMatrix(const Array1D<Type>& vector0, size_t rowsInMatrix)
02884 {
02885 Array2D<Type> XMatrix = zeros(rowsInMatrix, vector0.size() * rowsInMatrix,
02886 type_tag<Type>());
02887 for(size_t XRow = 0; XRow != rowsInMatrix; ++XRow) {
02888 for(size_t index0 = 0; index0 < vector0.length(); ++index0) {
02889 XMatrix(XRow, index0 * rowsInMatrix + XRow) = vector0(index0);
02890 }
02891 }
02892 return XMatrix;
02893 }
02894
02895
02896
02897
02898 template <class Type>
02899 Array2D<Type>
02900 identity(size_t rows, size_t columns, type_tag<Type>)
02901 {
02902 return identity<Type>(rows, columns);
02903 }
02904
02905
02906
02907
02908
02909 template <class Type>
02910 Array2D<Type>
02911 identity(size_t rows, size_t columns)
02912 {
02913 Array2D<Type> returnArray = zeros(rows, columns, type_tag<Type>());
02914 size_t rank = std::min(rows, columns);
02915 for(size_t index = 0; index < rank; ++index) {
02916 returnArray(index, index) = static_cast<Type>(1);
02917 }
02918 return returnArray;
02919 }
02920
02921
02922
02923
02924 template <class Type>
02925 Array1D<Type>
02926 ln(const Array1D<Type>& array0)
02927 {
02928 Array1D<Type> result(array0.size());
02929
02930 for(size_t index0 = 0; index0 < array0.size(); ++index0) {
02931 result[index0] = std::log(array0[index0]);
02932 }
02933 return result;
02934 }
02935
02936
02937
02938
02939
02940
02941
02942
02943
02944
02945
02946
02947
02948
02949
02950
02951
02952
02953
02954
02955
02956
02957
02958
02959
02960
02961
02962
02963
02964
02965 template <class Type>
02966 Array2D<Type>
02967 ln(const Array2D<Type>& array0)
02968 {
02969 Array2D<Type> result(array0.size());
02970
02971 for(size_t index0 = 0; index0 < array0.size(); ++index0) {
02972 result[index0] = std::log(array0[index0]);
02973 }
02974 return result;
02975 }
02976
02977
02978
02979
02980
02981 template <class Type>
02982 Array2D<Type>
02983 logicalNot(const Array2D<Type>& array0)
02984 {
02985 Array2D<Type> result(array0.rows(), array0.columns());
02986 std::transform(array0.begin(), array0.end(), result.begin(),
02987 unaryComposeFunctor(StaticCastFunctor<bool, Type>(),
02988 std::logical_not<Type>()));
02989 return result;
02990 }
02991
02992
02993
02994 template <class Type>
02995 inline Type
02996 magnitude(const Array1D<Type>& array0)
02997 {
02998 return static_cast<Type>(std::sqrt(magnitudeSquared(array0)));
02999 }
03000
03001
03002
03003 inline double
03004 magnitude(const Vector2D& vector0) {
03005 return std::sqrt(magnitudeSquared(vector0));
03006 }
03007
03008
03009
03010 inline double
03011 magnitude(const Vector3D& vector0) {
03012 return std::sqrt(magnitudeSquared(vector0));
03013 }
03014
03015
03016
03017
03018 template <class Type>
03019 inline typename NumericTraits<Type>::ProductType
03020 magnitudeSquared(const Array1D<Type>& array0)
03021 {
03022 typedef typename NumericTraits<Type>::ProductType ProductType;
03023 return std::inner_product(array0.begin(), array0.end(),
03024 array0.begin(), static_cast<ProductType>(0));
03025 }
03026
03027
03028
03029
03030 inline double
03031 magnitudeSquared(const Vector2D& vector0) {
03032 return (vector0.x() * vector0.x() + vector0.y() * vector0.y());
03033 }
03034
03035
03036
03037
03038 inline double
03039 magnitudeSquared(const Vector3D& vector0) {
03040 return (vector0.x() * vector0.x() + vector0.y() * vector0.y()
03041 + vector0.z() * vector0.z());
03042 }
03043
03044
03045
03046
03047
03048
03049
03050
03051 template <class Type>
03052 inline Array1D<typename NumericTraits<Type>::ProductType>
03053 matrixMultiply(const Array1D<Type>& vector0, const Array2D<Type>& matrix0)
03054 {
03055 return matrixMultiply(
03056 vector0, matrix0, type_tag<typename NumericTraits<Type>::ProductType>());
03057 }
03058
03059
03060
03061
03062
03063 template <class Type0, class Type1, class Type2>
03064 Array1D<Type2>
03065 matrixMultiply(const Array1D<Type0>& vector0, const Array2D<Type0>& matrix0,
03066 type_tag<Type2>)
03067 {
03068 if(vector0.size() != matrix0.rows()) {
03069 std::ostringstream message;
03070 message << "Can't left-multiply a "
03071 << matrix0.rows() << " x " << matrix0.columns()
03072 << " matrix by a " << vector0.size() << " element vector\n";
03073 DLR_THROW(ValueException, "matrixMultiply()", message.str().c_str());
03074 }
03075 Array1D<Type2> result = zeros(matrix0.columns(), type_tag<Type2>());
03076 size_t index = 0;
03077 for(size_t row = 0; row < matrix0.rows(); ++row) {
03078 for(size_t column = 0; column < matrix0.columns(); ++column) {
03079 result[column] += (static_cast<Type2>(vector0[row])
03080 * static_cast<Type2>(matrix0[index]));
03081 ++index;
03082 }
03083 }
03084 return result;
03085 }
03086
03087
03088
03089
03090
03091
03092 template <class Type>
03093 inline Array1D<typename NumericTraits<Type>::ProductType>
03094 matrixMultiply(const Array2D<Type>& matrix0, const Array1D<Type>& vector0)
03095 {
03096 return matrixMultiply(
03097 matrix0, vector0, type_tag<typename NumericTraits<Type>::ProductType>());
03098 }
03099
03100
03101
03102
03103
03104
03105 template <class Type0, class Type1, class Type2>
03106 Array1D<Type2>
03107 matrixMultiply(const Array2D<Type0>& matrix0, const Array1D<Type1>& vector0,
03108 type_tag<Type2>)
03109 {
03110 if(vector0.size() != matrix0.columns()) {
03111 std::ostringstream message;
03112 message << "matrixMultiply() -- can't right-multiply a "
03113 << matrix0.rows() << " x " << matrix0.columns()
03114 << " matrix by a " << vector0.size() << " element vector\n";
03115 DLR_THROW(ValueException, "matrixMultiply()", message.str().c_str());
03116 }
03117 Array1D<Type2> result(matrix0.rows());
03118 for(size_t row = 0; row < matrix0.rows(); ++row) {
03119 result[row] = dot(vector0, matrix0.row(row), type_tag<Type2>());
03120 }
03121 return result;
03122 }
03123
03124
03125
03126
03127 template <class Type>
03128 inline Array2D<typename NumericTraits<Type>::ProductType>
03129 matrixMultiply(const Array2D<Type>& matrix0, const Array2D<Type>& matrix1)
03130 {
03131 return matrixMultiply(
03132 matrix0, matrix1, type_tag<typename NumericTraits<Type>::ProductType>());
03133 }
03134
03135
03136
03137
03138
03139 template <class Type0, class Type1, class Type2>
03140 Array2D<Type2>
03141 matrixMultiply(const Array2D<Type0>& matrix0, const Array2D<Type1>& matrix1,
03142 type_tag<Type2>)
03143 {
03144 if(matrix1.rows() != matrix0.columns()) {
03145 std::ostringstream message;
03146 message << "matrixMultiply() -- can't left-multiply a "
03147 << matrix1.rows() << " x " << matrix1.columns()
03148 << " matrix by a "
03149 << matrix0.rows() << " x " << matrix0.columns()
03150 << " matrix\n";
03151 DLR_THROW(ValueException, "matrixMultiply()", message.str().c_str());
03152 }
03153 Array2D<Type2> result = zeros(matrix0.rows(), matrix1.columns(),
03154 type_tag<Type2>());
03155 for(size_t resultRow = 0; resultRow < result.rows(); ++resultRow) {
03156 for(size_t resultColumn = 0; resultColumn < result.columns();
03157 ++resultColumn) {
03158 for(size_t index = 0; index < matrix0.columns(); ++index) {
03159 result(resultRow, resultColumn) +=
03160 (static_cast<Type2>(matrix0(resultRow, index))
03161 * static_cast<Type2>(matrix1(index, resultColumn)));
03162 }
03163 }
03164 }
03165 return result;
03166 }
03167
03168
03169
03170 template <class Type>
03171 inline Type
03172 maximum(const Array1D<Type>& array0)
03173 {
03174 return maximum(array0, std::less<Type>());
03175 }
03176
03177
03178
03179
03180 template <class Type, class Functor>
03181 Type
03182 maximum(const Array1D<Type>& array0, Functor comparator)
03183 {
03184 if(array0.size() == 0) {
03185 DLR_THROW3(ValueException, "maximum()",
03186 "Can't find the maximum element of an empty array.");
03187 }
03188 return *std::max_element(array0.begin(), array0.end(), comparator);
03189 }
03190
03191
03192
03193
03194 template <class Iterator, class Type>
03195 Type
03196 mean(const Iterator& beginIter, const Iterator& endIter)
03197 {
03198 Type meanValue;
03199 size_t count = 0;
03200 while(beginIter != endIter) {
03201 meanValue += *beginIter;
03202 ++count;
03203 ++beginIter;
03204 }
03205 return meanValue / count;
03206 }
03207
03208
03209
03210
03211 template <class Type>
03212 inline double
03213 mean(const Array1D<Type>& array0)
03214 {
03215 if(array0.size() == 0) {
03216 DLR_THROW3(ValueException, "mean()",
03217 "Can't find the mean of an empty array.");
03218 }
03219 NumericTypeConversionFunctor<typename NumericTraits<Type>::SumType,
03220 Type> functor;
03221 return functor(
03222 mean(array0, type_tag<typename NumericTraits<Type>::SumType>()));
03223 }
03224
03225
03226
03227
03228
03229 template <class Type0, class Type1>
03230 inline Type1
03231 mean(const Array1D<Type0>& array0, type_tag<Type1>)
03232 {
03233 return sum(array0, type_tag<Type1>()) / static_cast<Type1>(array0.size());
03234 }
03235
03236
03237
03238
03239
03240 template <class Type>
03241 void
03242 getMeanAndCovariance(const Array2D<Type>& sampleArray,
03243 Array1D<double>& meanArray,
03244 Array2D<double>& covarianceArray,
03245 size_t majorAxis)
03246 {
03247
03248 if(sampleArray.size() == 0) {
03249 meanArray.reinit(0);
03250 covarianceArray.reinit(0, 0);
03251 return;
03252 }
03253
03254
03255 size_t dimensionality;
03256 size_t numberOfSamples;
03257 if(majorAxis == 0) {
03258 dimensionality = sampleArray.columns();
03259 numberOfSamples = sampleArray.rows();
03260 } else {
03261 dimensionality = sampleArray.rows();
03262 numberOfSamples = sampleArray.columns();
03263 }
03264
03265
03266
03267
03268
03269
03270
03271
03272
03273
03274
03275
03276
03277
03278
03279
03280
03281
03282 meanArray.reinit(dimensionality);
03283 meanArray = 0.0;
03284 covarianceArray.reinit(dimensionality, dimensionality);
03285 covarianceArray = 0.0;
03286 typename Array2D<Type>::const_iterator sampleIterator =
03287 sampleArray.begin();
03288
03289 if(majorAxis == 0) {
03290 for(size_t rowIndex = 0; rowIndex < sampleArray.rows(); ++rowIndex) {
03291 for(size_t columnIndex = 0; columnIndex < sampleArray.columns();
03292 ++columnIndex) {
03293
03294 meanArray[columnIndex] += *sampleIterator;
03295
03296
03297
03298 typename Array2D<Type>::const_iterator sampleIterator2 =
03299 sampleIterator;
03300 for(size_t covarianceIndex = columnIndex;
03301 covarianceIndex < dimensionality;
03302 ++covarianceIndex) {
03303
03304 covarianceArray(columnIndex, covarianceIndex) +=
03305 *sampleIterator * *sampleIterator2;
03306 ++sampleIterator2;
03307 }
03308 ++sampleIterator;
03309 }
03310 }
03311 } else {
03312 for(size_t rowIndex = 0; rowIndex < sampleArray.rows(); ++rowIndex) {
03313 for(size_t columnIndex = 0; columnIndex < sampleArray.columns();
03314 ++columnIndex) {
03315
03316 meanArray[rowIndex] += *sampleIterator;
03317
03318
03319
03320 typename Array2D<Type>::const_iterator sampleIterator2 =
03321 sampleIterator;
03322 for(size_t covarianceIndex = rowIndex;
03323 covarianceIndex < dimensionality;
03324 ++covarianceIndex) {
03325
03326 covarianceArray(rowIndex, covarianceIndex) +=
03327 *sampleIterator * *sampleIterator2;
03328 sampleIterator2 += numberOfSamples;
03329 }
03330 ++sampleIterator;
03331 }
03332 }
03333 }
03334
03335
03336 meanArray /= static_cast<double>(numberOfSamples);
03337
03338
03339 for(size_t index0 = 0; index0 < dimensionality; ++index0) {
03340 for(size_t index1 = index0; index1 < dimensionality; ++index1) {
03341
03342 covarianceArray(index0, index1) -=
03343 numberOfSamples * meanArray[index0] * meanArray[index1];
03344 covarianceArray(index0, index1) /= (numberOfSamples - 1);
03345
03346
03347 if(index0 != index1) {
03348 covarianceArray(index1, index0) = covarianceArray(index0, index1);
03349 }
03350 }
03351 }
03352
03353 }
03354
03355
03356
03357
03358 template <class Type>
03359 inline Type
03360 minimum(const Array1D<Type>& array0)
03361 {
03362 if(array0.size() == 0) {
03363 DLR_THROW3(ValueException, "minimum()",
03364 "Can't find the minimum element of an empty array.");
03365 }
03366 return minimum(array0, std::less<Type>());
03367 }
03368
03369
03370
03371
03372 template <class Type, class Functor>
03373 Type
03374 minimum(const Array1D<Type>& array0, Functor comparator)
03375 {
03376 return *std::min_element(array0.begin(), array0.end(), comparator);
03377 }
03378
03379
03380
03381
03382 template <class FUNCTOR>
03383 double newtonRaphson(double startPoint, FUNCTOR objectiveFunction,
03384 double epsilon, size_t maxIterations)
03385 {
03386 double argument = startPoint;
03387 for(size_t count = 0; count < maxIterations; ++count) {
03388 double result = objectiveFunction(argument);
03389 if(std::fabs(result) < epsilon) {
03390 return argument;
03391 }
03392 double resultPrime = objectiveFunction.derivative(argument);
03393 if(resultPrime == 0.0) {
03394 DLR_THROW(ValueException, "newtonRaphson()", "Zero derivative.");
03395 }
03396 argument = argument - (result / resultPrime);
03397 }
03398 DLR_THROW(ValueException, "newtonRaphson()",
03399 "Root finding failed to converge.");
03400 return 0.0;
03401 }
03402
03403
03404
03405 template <class Type>
03406 inline double
03407 normalizedCorrelation(const Array1D<Type>& signal0,
03408 const Array1D<Type>& signal1)
03409 {
03410 return normalizedCorrelation(signal0, signal1, type_tag<double>());
03411 }
03412
03413
03414
03415
03416 template <class Type, class Type2>
03417 Type2
03418 normalizedCorrelation(const Array1D<Type>& signal0,
03419 const Array1D<Type>& signal1,
03420 type_tag<Type2>)
03421 {
03422 if(signal0.size() != signal1.size()) {
03423 DLR_THROW(ValueException, "normalizedCorrelation()",
03424 "Input arrays must have the same size.");
03425 }
03426 Type2 oneOverN =
03427 static_cast<Type2>(1) / static_cast<Type2>(signal0.size());
03428 Type2 sum0 = static_cast<Type2>(0);
03429 Type2 sum1 = static_cast<Type2>(0);
03430 Type2 sum00 = static_cast<Type2>(0);
03431 Type2 sum01 = static_cast<Type2>(0);
03432 Type2 sum11 = static_cast<Type2>(0);
03433 typedef typename Array1D<Type>::const_iterator SignalIter;
03434 SignalIter iter0 = signal0.begin();
03435 SignalIter iter1 = signal1.begin();
03436 SignalIter end0 = signal0.end();
03437 while(iter0 != end0) {
03438 sum0 += *iter0;
03439 sum1 += *iter1;
03440 sum00 += (*iter0) * (*iter0);
03441 sum01 += (*iter0) * (*iter1);
03442 sum11 += (*iter1) * (*iter1);
03443 ++iter0;
03444 ++iter1;
03445 }
03446 Type2 numerator = sum01 - oneOverN * sum0 * sum1;
03447 Type2 denominator = std::sqrt((sum00 - oneOverN * sum0 * sum0)
03448 * (sum11 - oneOverN * sum1 * sum1));
03449 return numerator / denominator;
03450 }
03451
03452
03453
03454 template <class Type>
03455 inline Array1D<Type>
03456 ones(int size, type_tag<Type>)
03457 {
03458 Array1D<Type> result(size);
03459 result = 1;
03460 return result;
03461 }
03462
03463
03464
03465 template <class Type>
03466 inline Array2D<Type>
03467 ones(int rows, int columns, type_tag<Type>)
03468 {
03469 Array2D<Type> result(rows, columns);
03470 result = 1;
03471 return result;
03472 }
03473
03474
03475
03476 template <class Type>
03477 inline Array2D<typename NumericTraits<Type>::ProductType>
03478 outerProduct(const Array1D<Type>& x, const Array1D<Type>& y)
03479 {
03480 return outerProduct(
03481 x, y, type_tag<typename NumericTraits<Type>::ProductType>());
03482 }
03483
03484
03485
03486
03487 template <class Type0, class Type1, class Type2>
03488 Array2D<Type2>
03489 outerProduct(const Array1D<Type0>& x, const Array1D<Type1>& y,
03490 type_tag<Type2>)
03491 {
03492 Array2D<Type2> result(x.size(), y.size());
03493 typename Array2D<Type2>::iterator runningIterator = result.begin();
03494 typename Array1D<Type1>::const_iterator yBegin = y.begin();
03495 typename Array1D<Type1>::const_iterator yEnd = y.end();
03496 for(size_t index = 0; index < x.size(); ++index) {
03497 std::transform(yBegin, yEnd, runningIterator,
03498 std::bind2nd(std::multiplies<Type2>(), x[index]));
03499 runningIterator += y.size();
03500 }
03501 return result;
03502 }
03503
03504
03505
03506
03507 template <class Type>
03508 Array1D<Type>
03509 range(Type start, Type stop, Type stride)
03510 {
03511 if(stride == 0) {
03512 DLR_THROW3(ValueException, "range(Type, Type, Type)",
03513 "Argument \"stride\" must not be equal to 0.");
03514 }
03515 int length = static_cast<int>((stop - start) / stride);
03516 if(stride > 0) {
03517 if((length * stride) < (stop - start)) {
03518 ++length;
03519 }
03520 } else {
03521 if((length * stride) > (stop - start)) {
03522 ++length;
03523 }
03524 }
03525 Array1D<Type> result(length);
03526 for(int index = 0; index < length; ++index) {
03527 result(index) = start;
03528 start += stride;
03529 }
03530 return result;
03531 }
03532
03533
03534
03535
03536 template <class Type>
03537 inline Type
03538 rms(const Array1D<Type>& array0)
03539 {
03540 NumericTypeConversionFunctor<
03541 typename NumericTraits<Type>::SumType, Type> functor;
03542 return functor(rms(array0,
03543 type_tag<typename NumericTraits<Type>::ProductType>()));
03544 }
03545
03546
03547
03548
03549 template <class Type0, class Type1>
03550 Type1
03551 rms(const Array1D<Type0>& array0, type_tag<Type1>)
03552 {
03553 Type1 accumulator = 0;
03554 for(int index = 0; index < array0.size(); ++index) {
03555 Type1 element = static_cast<Type1>(array0[index]);
03556 accumulator += element * element;
03557 }
03558 return ::sqrt(accumulator / static_cast<Type1>(array0.size()));
03559 }
03560
03561
03562
03563 template <class Type>
03564 Array2D<Type>
03565 rowIndices(size_t rows, size_t columns, type_tag<Type>)
03566 {
03567 Array2D<Type> returnValue(rows, columns);
03568 for(size_t row=0; row < rows; ++row) {
03569 Type* rowPtr = returnValue.data(row, 0);
03570 std::fill(rowPtr, rowPtr + columns, static_cast<Type>(row));
03571 }
03572 return returnValue;
03573 }
03574
03575
03576
03577 template <class Type0, class Type1>
03578 inline bool
03579 shapeMatch(const Array1D<Type0>& array0, const Array1D<Type1>& array1)
03580 {
03581 return array0.size() == array1.size();
03582 }
03583
03584
03585
03586 template <class Type0, class Type1>
03587 inline bool
03588 shapeMatch(const Array2D<Type0>& array0, const Array2D<Type1>& array1)
03589 {
03590 return ((array0.rows() == array1.rows())
03591 && (array0.columns() == array1.columns()));
03592 }
03593
03594
03595
03596 template <class Type0, class Type1>
03597 inline bool
03598 shapeMatch(const Array3D<Type0>& array0, const Array3D<Type1>& array1)
03599 {
03600 return ((array0.shape0() == array1.shape0())
03601 && (array0.shape1() == array1.shape1())
03602 && (array0.shape2() == array1.shape2()));
03603 }
03604
03605
03606
03607 template <class Type>
03608 inline Array2D<Type>
03609 skewSymmetric(const Array1D<Type>& vector0)
03610 {
03611 if(vector0.size() != 3) {
03612 std::ostringstream message;
03613 message << "Argument must have exactly 3 elements, but instead has "
03614 << vector0.size() << "elements.";
03615 DLR_THROW(ValueException, "skewSymmetric()", message.str().c_str());
03616 }
03617 Array2D<Type> returnVal(3, 3);
03618 returnVal(0, 0) = 0;
03619 returnVal(0, 1) = -vector0(2);
03620 returnVal(0, 2) = vector0(1);
03621
03622 returnVal(1, 0) = vector0(2);
03623 returnVal(1, 1) = 0;
03624 returnVal(1, 2) = -vector0(0);
03625
03626 returnVal(2, 0) = -vector0(1);
03627 returnVal(2, 1) = vector0(0);
03628 returnVal(2, 2) = 0;
03629 return returnVal;
03630 }
03631
03632
03633
03634
03635 template <class Type>
03636 void
03637 solveQuadratic(Type c0, Type c1, Type c2,
03638 Type& root0, Type& root1, bool& valid, bool)
03639 {
03640 valid = solveQuadratic(c0, c1, c2, root0, root1);
03641 }
03642
03643
03644
03645
03646 template <class Type>
03647 inline double
03648 standardDeviation(const Array1D<Type>& array0)
03649 {
03650 return standardDeviation(array0, type_tag<double>());
03651 }
03652
03653
03654
03655
03656 template <class Type0, class Type1>
03657 inline Type1
03658 standardDeviation(const Array1D<Type0>& array0, type_tag<Type1>)
03659 {
03660 return ::sqrt(variance(array0, type_tag<Type1>()));
03661 }
03662
03663
03664
03665
03666 template <class Type>
03667 inline typename NumericTraits<Type>::SumType
03668 sum(const Array1D<Type>& array0)
03669 {
03670 return sum(array0, type_tag<typename NumericTraits<Type>::SumType>());
03671 }
03672
03673
03674
03675
03676 template <class Type, class Type2>
03677 Type2
03678 sum(const Array1D<Type>& array0, type_tag<Type2>)
03679 {
03680 return std::accumulate(array0.begin(), array0.end(),
03681 static_cast<Type2>(0));
03682 }
03683
03684
03685
03686
03687 template <class Type>
03688 inline typename NumericTraits<Type>::SumType
03689 sum(const Array2D<Type>& array0,
03690 const Index2D& upperLeftCorner,
03691 const Index2D& lowerRightCorner)
03692 {
03693 return sum(array0, upperLeftCorner, lowerRightCorner,
03694 type_tag<typename NumericTraits<Type>::SumType>());
03695 }
03696
03697
03698
03699
03700 template <class Type, class Type2>
03701 Type2
03702 sum(const Array2D<Type>& array0,
03703 const Index2D& upperLeftCorner,
03704 const Index2D& lowerRightCorner,
03705 type_tag<Type2>)
03706 {
03707 Type2 result = static_cast<Type2>(0);
03708 for(int row = upperLeftCorner.getRow();
03709 row < lowerRightCorner.getRow();
03710 ++row) {
03711 typename Array2D<Type>::const_iterator rowBegin = array0.rowBegin(row);
03712 result += std::accumulate(
03713 rowBegin + upperLeftCorner.getColumn(),
03714 rowBegin + lowerRightCorner.getColumn(),
03715 static_cast<Type2>(0));
03716 }
03717 return result;
03718 }
03719
03720
03721
03722
03723 template <class Type, class IntegralType>
03724 Array1D<Type>
03725 take(const Array1D<Type>& dataArray,
03726 const Array1D<IntegralType>& indexArray)
03727 {
03728 Array1D<Type> resultArray(indexArray.size());
03729 for(unsigned int ii = 0; ii < indexArray.size(); ++ii) {
03730 resultArray[ii] = dataArray[indexArray[ii]];
03731 }
03732 return resultArray;
03733 }
03734
03735
03736
03737
03738 template <class Type, class IntegralType>
03739 Array1D<Type>
03740 take(const Array2D<Type>& dataArray,
03741 const Array1D<IntegralType>& indexArray)
03742 {
03743 return take(dataArray.ravel(), indexArray);
03744 }
03745
03746
03747
03748
03749
03750 template <class Type, class IntegralType>
03751 Array2D<Type>
03752 take(const Array2D<Type>& dataArray,
03753 const Array1D<IntegralType>& indexArray,
03754 unsigned int axis)
03755 {
03756 Array2D<Type> resultArray;
03757 switch(axis) {
03758 case 0:
03759 resultArray.reinit(indexArray.size(), dataArray.columns());
03760 for(unsigned int ii = 0; ii < indexArray.size(); ++ii) {
03761 unsigned int jj = indexArray[ii];
03762 std::copy(dataArray.rowBegin(jj), dataArray.rowEnd(jj),
03763 resultArray.rowBegin(ii));
03764 }
03765 break;
03766 default:
03767 resultArray.reinit(dataArray.rows(), indexArray.size());
03768 for(unsigned int ii = 0; ii < dataArray.rows(); ++ii) {
03769 typename Array2D<Type>::iterator resultIterator =
03770 resultArray.rowBegin(ii);
03771 typename Array2D<Type>::const_iterator dataRowBegin =
03772 dataArray.rowBegin(ii);
03773 for(unsigned int jj = 0; jj < indexArray.size(); ++jj) {
03774 *resultIterator = *(dataRowBegin + indexArray[jj]);
03775 ++resultIterator;
03776 }
03777 }
03778 break;
03779 }
03780 return resultArray;
03781 }
03782
03783
03784
03785
03786 template <class Type>
03787 inline double
03788 variance(const Array1D<Type>& array0)
03789 {
03790 return variance(array0, type_tag<double>());
03791 }
03792
03793
03794
03795
03796 template <class Type0, class Type1>
03797 inline Type1
03798 variance(const Array1D<Type0>& array0, type_tag<Type1>)
03799 {
03800 Type1 meanValue = mean(array0, type_tag<Type1>());
03801 Type1 accumulator = 0;
03802 for(size_t index0 = 0; index0 < array0.size(); ++index0) {
03803 Type1 offset = static_cast<Type1>(array0[index0]) - meanValue;
03804 accumulator += offset * offset;
03805 }
03806 return accumulator / static_cast<Type1>(array0.size());
03807 }
03808
03809
03810
03811 template <class Type>
03812 inline Array1D<Type>
03813 zeros(size_t size, type_tag<Type>)
03814 {
03815 Array1D<Type> result(size);
03816 result = 0;
03817 return result;
03818 }
03819
03820
03821
03822 template <class Type>
03823 inline Array2D<Type>
03824 zeros(size_t rows, size_t columns, type_tag<Type>)
03825 {
03826 Array2D<Type> result(rows, columns);
03827 result = 0;
03828 return result;
03829 }
03830
03831
03832
03833
03834 template <class Type>
03835 inline Array3D<Type>
03836 zeros(size_t shape0, size_t shape1, size_t shape2, type_tag<Type>)
03837 {
03838 Array3D<Type> result(shape0, shape1, shape2);
03839 result = 0;
03840 return result;
03841 }
03842
03843 }
03844
03845 }
03846
03847 #endif