utilities.h

Go to the documentation of this file.
00001 
00015 #ifndef DLR_NUMERIC_UTILITIES_H
00016 #define DLR_NUMERIC_UTILITIES_H
00017 
00018 #include <cmath>
00019 // #include <iostream>
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     } // nameapace contants
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 //    * This function returns an array of indices, result, so that the
00350 //    * sequence (array0[result[0]], array0[result[1]],
00351 //    * array0[result[2]], ...) is sorted from smallest to largest using
00352 //    * the supplied comparison operator.
00353 //    *
00354 //    * @param array0 The of this array will be compared to each other in
00355 //    * order to establish the correct sequence of indices.
00356 //    *
00357 //    * @param comparator This argument is a functor which supports the
00358 //    * std::binary_function<Type, Type, bool> interface, and returns
00359 //    * true if its first argument is smaller than its second argument.
00360 //    *
00361 //    * @return An array of indices as described above.
00362 //    */
00363 //   template <class Type, class Functor>
00364 //   Array1D<size_t>
00365 //   argsort(const Array1D<Type>& array0, Functor comparator);
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     // The next function is commented out because some environments use
01148     // #define directives for min() and max(), which don't obey
01149     // namespaces, and makes the code not compile if we define our own
01150     // min/max.
01151     // 
01152     // 
01153     // /** 
01154     // * This function is an alias for the function maximum(const
01155     // * Array1D&) below.
01156     // * 
01157     // * @param array0 See documentation for maximum(const Array1D&).
01158     // *
01159     // * @return See documentation for maximum(const Array1D&).
01160     // */
01161     // template <class Type>
01162     // inline Type
01163     // max(const Array1D<Type>& array0) {return maximum(array0);}
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     // The next function is commented out because some environments use
01256     // #define directives for min() and max(), which don't obey
01257     // namespaces, and makes the code not compile if we define our own
01258     // min/max.
01259     // 
01260     // 
01261     // /** 
01262     // * This function is an alias for the function minimum(const
01263     // * Array1D&) below.
01264     // * 
01265     // * @param array0 See documentation for minimum(const Array1D&).
01266     // *
01267     // * @return See documentation for minimum(const Array1D&).
01268     // */
01269     // template <class Type>
01270     // inline Type
01271     // min(const Array1D<Type>& array0) {return minimum(array0);}
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   } // namespace numeric
02103 
02104 } // namespace dlr
02105 
02106 
02107 /* ======= Declarations to maintain compatibility with legacy code. ======= */
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 } // namespace dlr
02158 
02159 
02160 /* ================ Definitions follow ================ */
02161 
02162 #include <cmath>
02163 #include <sstream>
02164 #include <numeric>
02165 //#include <algorithm>
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           // fabsl doesn't appear to be commonly available.
02220           // return std::fabsl(input);
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           // Strange.  fabsf shows up in the global namespace.  Is this
02256           // a bug in the g++ 3.3 stand library?
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           // Many compilers don't support C99.
02294           // return std::llabs(input);
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     // This function returns an array of the same size and element type
02322     // as its input argument, in which each element is set to the
02323     // absolute value of the corresponding element of the input array.
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     // This function returns an array of the same size and element type
02335     // as its input argument, in which each element is set to the
02336     // absolute value of the corresponding element of the input array.
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     // This function returns true if each element of its argument is
02349     // false, and returns false otherwise.
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     // This function returns true if each element of its argument is
02366     // true, and returns false otherwise.
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     // This function returns true if any element of its argument is
02383     // false, and returns false otherwise.
02384     template <class Type>
02385     inline bool
02386     anyFalse(const Array1D<Type>& array0)
02387     {
02388       return !allTrue(array0);
02389     }
02390 
02391   
02392     // This function returns true if any element of its argument is
02393     // true, and returns false otherwise.
02394     template <class Type>
02395     inline bool
02396     anyTrue(const Array1D<Type>& array0)
02397     {
02398       return !allFalse(array0);
02399     }
02400 
02401   
02402     // This function returns the index of the largest element of its input
02403     // array.
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     // This function returns the index of the largest element of its input
02413     // sequence.
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     // This function returns the index of the largest element of its
02423     // input array, where largeness is defined by the second argument.
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     // This function returns the index of the largest element of its
02433     // input sequence, where largeness is defined by the second argument.
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     // This function returns an Index2D instance indicating which is
02444     // the largest element of its input array.
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     // This function returns an Index2D instance indicating which is
02454     // the largest element of its input array, where largeness is
02455     // defined by the second argument.
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(); // Int division.
02462       size_t column = ravelIndex - row * array0.columns();
02463       return Index2D(row, column);
02464     }
02465   
02466 
02467     // This function returns the index of the smallest element of its input
02468     // array.
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     // This function returns the index of the smallest element of its
02478     // input array, where largeness is defined by the second argument.
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     // This function returns an Index2D instance indicating which is
02490     // the smallest element of its input array.
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     // This function returns an Index2D instance indicating which is
02500     // the smallest element of its input array, where smallness is
02501     // defined by the second argument.
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(); // Int division.
02508       size_t column = ravelIndex - row * array0.columns();
02509       return Index2D(row, column);
02510     }
02511 
02512 
02513     // This function returns an array of indices, result, so that the
02514     // sequence (array0[result[0]], array0[result[1]],
02515     // array0[result[2]], ...) is sorted from smallest to largest using
02516     // operator<().
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 //   // This function returns an array of indices, result, so that the
02535 //   // sequence (array0[result[0]], array0[result[1]],
02536 //   // array0[result[2]], ...) is sorted from smallest to largest using
02537 //   // the supplied comparison operator.
02538 //   template <class Type, class Functor>
02539 //   Array1D<size_t>
02540 //   argsort(const Array1D<Type>& array0, Functor comparator)
02541 //   {
02542 //     Array1D< std::pair<Type, size_t> > sortVector(array0.size());
02543 //     Array1D<size_t> resultVector(array0.size());
02544 
02545 //     for(size_t index0 = 0; index0 < array0.size(); ++index0) {
02546 //       sortVector[index0] = std::pair<Type, size_t>(array0[index0], index0);
02547 //     }
02548 //     std::sort(sortVector.begin(), sortVector.end(), comparator);
02549 //     std::transform(sortVector.begin(), sortVector.end(), resultVector.begin(),
02550 //                    ExtractSecondFunctor<Type, size_t>());
02551 //     return resultVector;    
02552 //   }
02553 
02554   
02555     // This function returns an Array1D in which each element has the
02556     // value of the largest element in one row or column of the input
02557     // Array2D.
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     // This function returns an Array1D in which each element has the
02567     // value of the largest element in one row or column of the input
02568     // Array2D, where largeness is defined by the third argument.
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     // This function returns an Array1D in which each element has the
02616     // value of the smallest element in one row or column of the input
02617     // Array2D.
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     // This function returns an Array1D in which each element has the
02627     // value of the smallest element in one row or column of the input
02628     // Array2D, where smallness is defined by the third argument.
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     // This function returns an Array1D in which each element has the
02676     // sum of one row or column of the input Array2D.
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     // This function returns an Array1D in which each element has the
02686     // sum of one row or column of the input Array2D.
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     // This function returns an Array1D in which each element has the
02696     // sum of one row or column of the input Array2D.
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     // columnIndices(rows, columns): Returns an Array2D in which each
02740     // element contains the index of its column.
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     // This function selects those elements of an input Array1D which
02756     // correspond to "true" values of a mask Array1D, and returns an
02757     // Array1D containing only those elements.
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     // This function behaves in exactly the same way as compress(const
02768     // Array1D&, const Array1D&), above, but it permits the user to
02769     // specify the number of true elements in the condition array.
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       // copy_mask(input.begin(), input.end(), condition.begin(),
02787       //           result.begin());
02788       //
02789       // Hmm... copy_mask was previously defined in dlrLib3 as follows:
02790       //
02791       // template <class In0, class In1, class Out>
02792       // Out copy_mask(In0 first, In0 last, In1 mask, Out res)
02793       // {
02794       //   while(first != last) {
02795       //     if(*mask) {
02796       //       *res++ = *first;
02797       //     }
02798       //     ++mask;
02799       //     ++first;
02800       //   }
02801       //   return res;
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     // This element counts the number of elements of the input array
02820     // which evaluate to true, and returns that number.
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     // This function computes the cross product of two Vector3D
02829     // instances.
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     // This function computes the inner product of two input arrays.
02838     // The computation is done using the ProductType specified by
02839     // the appropriate NumericTraits class.
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     // This function computes the inner product of two input arrays, and
02848     // allows the user to control which type is used to do the
02849     // calculation.
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     // This function computes the inner product of two Vector2D instances.
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     // This function computes the inner product of two Vector3D instances.
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     // This function computes the matrix X such that A * x = X * vec(A).
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     // This function returns a 2D matrix with the specified shape in
02896     // which the elements on the diagonal are set to 1, and all other
02897     // elements are set to 0.
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     // This function returns a 2D matrix with the specified shape in
02907     // which the elements on the diagonal are set to 1, and all other
02908     // elements are set to 0.
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     // This function returns an array in which each element is the
02923     // natural logarithm of the corresponding element of its input.
02924     template <class Type>
02925     Array1D<Type>
02926     ln(const Array1D<Type>& array0)
02927     {
02928       Array1D<Type> result(array0.size());
02929       // std::transform(array0.begin(), array0.end(), std::log);
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 //   // Strange troubles compiling these next two functions.
02938 //
02939 //   // This function returns an array in which each element is the
02940 //   // natural logarithm of the corresponding element of its input.
02941 //   template <>
02942 //   Array1D<float>
02943 //   ln(const Array1D<float>& array0)
02944 //   {
02945 //     Array1D<float> result(array0.size());
02946 //     std::transform(array0.begin(), array0.end(), std::logf);
02947 //     return result;
02948 //   }
02949 
02950 
02951 //   // This function returns an array in which each element is the
02952 //   // natural logarithm of the corresponding element of its input.
02953 //   template <>
02954 //   Array1D<long double>
02955 //   ln(const Array1D<long double>& array0)
02956 //   {
02957 //     Array1D<long double> result(array0.size());
02958 //     std::transform(array0.begin(), array0.end(), std::logl);
02959 //     return result;
02960 //   }
02961 
02962 
02963     // This function returns an array in which each element is the
02964     // natural logarithm of the corresponding element of its input.
02965     template <class Type>
02966     Array2D<Type>
02967     ln(const Array2D<Type>& array0)
02968     {
02969       Array2D<Type> result(array0.size());
02970       // std::transform(array0.begin(), array0.end(), std::log);
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     // This function returns an Array2D instance in which the value of
02979     // each element is the logical not of the corresponding element of
02980     // the input array.
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     // This function computes the magnitude of its input argument.
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     // This function computes the magnitude of its input argument.
03003     inline double
03004     magnitude(const Vector2D& vector0) {
03005       return std::sqrt(magnitudeSquared(vector0));
03006     }
03007 
03008   
03009     // This function computes the magnitude of its input argument.
03010     inline double
03011     magnitude(const Vector3D& vector0) {
03012       return std::sqrt(magnitudeSquared(vector0));
03013     }
03014 
03015 
03016     // This function computes the square of the magnitude of its input
03017     // argument.
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     // This function computes the square of the magnitude of its input
03029     // argument.
03030     inline double
03031     magnitudeSquared(const Vector2D& vector0) {
03032       return (vector0.x() * vector0.x() + vector0.y() * vector0.y());
03033     }
03034 
03035   
03036     // This function computes the square of the magnitude of its input
03037     // argument.
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     // This function computes a vector * matrix product.  Assuming the
03046     // first argument, vector0, represents a column vector and the
03047     // second argument, matrix0, represents a 2D array, then this function
03048     // computes the product:
03049     //   vector0' * matrix0
03050     // Where the single quote indicates transpose.
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     // This function computes a vector * matrix product just like
03060     // matrixMultiply(const Array1D<Type>&, const Array2D<Type>&), above.
03061     // This function differs in that the element type of the return
03062     // value is set explicitly using a third argument.
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     // This function computes a matrix * vector product.  Assuming the
03088     // first argument, matrix0, represents a matrix, and the and the
03089     // second argument, vector0, represents a column vector, then this
03090     // function computes the product:
03091     //  matrix0 * vector0
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     // This function computes a matrix * vector product just like
03102     // matrixMultiply(const Array2D<Type>&, const Array1D<Type>&), above.
03103     // This function differs in that the element type of the return
03104     // value is set explicitly using a third argument.
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     // This function computes a matrix * matrix product.  That is,
03125     // the elements of the resulting array are the dot products of the
03126     // rows of the first argument with the columns of the second argument.
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     // This function computes a matrix * matrix product, just like
03136     // matrixMultiply(const Array2D<Type>&, const Array2D<Type>&), above.
03137     // This function differs in that the element type of the return
03138     // value is set explicitly using a third argument.
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     // This function returns a copy of the largest element in the input
03169     // Array1D instance.
03170     template <class Type>
03171     inline Type
03172     maximum(const Array1D<Type>& array0)
03173     {
03174       return maximum(array0, std::less<Type>());
03175     }
03176 
03177     // This function returns a copy of the largest element in the input
03178     // Array1D instance, where largeness is defined by the return value
03179     // of the second argument.
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     // This function computes the average value, or geometric mean, of
03193     // the elements of input sequence.
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     // This function computes the average value, or geometric mean, of
03210     // the elements of its argument.
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     // This function computes the average value, or geometric mean, of
03227     // the elements of its argument, and allows the user to specify the
03228     // precision with which the computation is carried out.
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     // This function estimates the mean and covariance of a set of
03238     // vectors, which are represented by the rows (or columns) of the
03239     // input 2D array.
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       // Check input.
03248       if(sampleArray.size() == 0) {
03249         meanArray.reinit(0);
03250         covarianceArray.reinit(0, 0);
03251         return;
03252       }
03253 
03254       // Compute number of samples and sample dimensionality.
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       // Now compute the mean value of the sample points.
03266       // For efficiency, we simultaneously compute covariance.
03267       // We make use of the identity:
03268       // 
03269       //   covariance = E[(x - u)*(x - u)^T] = E[x * x^T] - u * u^T
03270       // 
03271       // where E[.] denotes expected value and u is the mean.
03272       //
03273       // Note that the sample mean is an unbiased estimator for u, but
03274       // substituting the sample mean for u in the covariance equation
03275       // gives a biased estimator of covariance.  We resolve this by
03276       // using the unbiased estimator:
03277       //
03278       //   covariance = (1/(N-1))sum(x * x^T) - (N/(N-1)) * uHat * uHat^T
03279       //
03280       // where uHat is the sample mean.
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             // Update accumulator for mean.
03294             meanArray[columnIndex] += *sampleIterator;
03295 
03296             // Update accumulator for covariance.  We only compute the top
03297             // half of the covariance matrix for now, since it's symmetric.
03298             typename Array2D<Type>::const_iterator sampleIterator2 =
03299               sampleIterator;
03300             for(size_t covarianceIndex = columnIndex;
03301                 covarianceIndex < dimensionality;
03302                 ++covarianceIndex) {
03303               // Note(xxx): Should fix this to run faster.
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             // Update accumulator for mean.
03316             meanArray[rowIndex] += *sampleIterator;
03317 
03318             // Update accumulator for covariance.  We only compute the top
03319             // half of the covariance matrix for now, since it's symmetric.
03320             typename Array2D<Type>::const_iterator sampleIterator2 =
03321               sampleIterator;
03322             for(size_t covarianceIndex = rowIndex;
03323                 covarianceIndex < dimensionality;
03324                 ++covarianceIndex) {
03325               // Note(xxx): Should fix this to run faster.
03326               covarianceArray(rowIndex, covarianceIndex) +=
03327                 *sampleIterator * *sampleIterator2;
03328               sampleIterator2 += numberOfSamples;
03329             }
03330             ++sampleIterator;
03331           }
03332         }
03333       }
03334     
03335       // Finish up the computation of the mean.
03336       meanArray /= static_cast<double>(numberOfSamples);
03337 
03338       // Finish up the computation of the covariance matrix.
03339       for(size_t index0 = 0; index0 < dimensionality; ++index0) {
03340         for(size_t index1 = index0; index1 < dimensionality; ++index1) {
03341           // Add correction to the top half of the covariance matrix.
03342           covarianceArray(index0, index1) -=
03343             numberOfSamples * meanArray[index0] * meanArray[index1];
03344           covarianceArray(index0, index1) /= (numberOfSamples - 1);
03345 
03346           // And copy to the bottom half.
03347           if(index0 != index1) {
03348             covarianceArray(index1, index0) = covarianceArray(index0, index1);
03349           }
03350         }
03351       }
03352 
03353     }
03354   
03355                     
03356     // This function returns a copy of the smallest element in the input
03357     // Array1D instance.
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     // This function returns a copy of the smallest element in the input
03370     // Array1D instance, where largeness is defined by the value of the
03371     // second argument.
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     // This function uses the iterative method of Newton and Raphson to
03381     // search for a zero crossing of the supplied functor.
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;  // keep compiler happy
03401     }
03402 
03403     // This function computes the normalized correlation of two Array1D
03404     // arguments.
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     // This function is equivalent to normalizedCorrelation(const
03414     // Array1D&, const Array1D&), except that the computation is carried
03415     // out using the type specified by the third argument.
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     // This function returns an Array1D of the specified size and type
03453     // in which the value of every element is initialized to 1.
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     // This function returns an Array2D of the specified size and type
03464     // in which the value of every element is initialized to one.
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     // This function computes the outer product of two input Array1D
03475     // instances. 
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     // This function computes the outer product of two input Array1D
03485     // instances and allows the user to control which type is used to do the
03486     // computation.
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     // This function returns an Array1D in which the first element has
03505     // value equal to argument "start," and each subsequent element has
03506     // value equal to the previous element plus argument "stride."
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     // This function computes the RMS (Root Mean Square) value of the
03535     // elements of its argument.
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     // This function computes the RMS (Root Mean Square) value of the
03547     // elements of its argument, and allows the user to specify the
03548     // precision with which the computation is carried out.
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     // rowIndices(rows, columns): Returns an Array2D in which each
03562     // element contains the index of its row.
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     // This function returns true if the two arrays have the same shape,
03576     // false otherwise.
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     // This function returns true if the two arrays have the same shape,
03585     // false otherwise.
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     // This function returns true if the two arrays have the same shape,
03595     // false otherwise.
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     // skewSymmetric(x): Returns a skew symmetric matrix X such
03606     // that matrixMultiply(X, y) = cross(x, y).
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     // This function computes the real roots of the quadratic polynomial
03634     // c0*x^2 + c1*x + c0 = 0.
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     // This function computes the standard deviation of a group of
03645     // scalar samples.
03646     template <class Type>
03647     inline double
03648     standardDeviation(const Array1D<Type>& array0)
03649     {
03650       return standardDeviation(array0, type_tag<double>());
03651     }
03652 
03653     // This function computes the standard deviation, of the elements of
03654     // its argument, and allows the user to specify the precision with
03655     // which the computation is carried out.
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     // This function computes the sum of all elements in the input
03664     // array.  The summation is done using the SumType specified by the
03665     // appropriate NumericTraits class.
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     // This function computes the sum of all elements in the input
03674     // array, and allows the user to control which type is used to do the
03675     // summation.
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     // This function computes the sum of those elements of its
03686     // argument which lie within a rectangular region of interest.
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     // This function computes the sum of those elements of its
03699     // argument which lie within a rectangular region of interest.
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     // This function returns an array made up of only those elements
03722     // of dataArray that correspond to indices in indexArray.
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     // This function returns an array made up of only those elements
03737     // of dataArray that correspond to indices in indexArray.
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     // This function works just like take(Array2D const&, Array1D
03748     // const&), with the exception that the input array is not
03749     // flattened, and entire rows (or columns) are selected.
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     // This function computes the standard deviation of a group of
03785     // scalar samples.
03786     template <class Type>
03787     inline double
03788     variance(const Array1D<Type>& array0)
03789     {
03790       return variance(array0, type_tag<double>());
03791     }
03792 
03793     // This function computes the standard deviation, of the elements of
03794     // its argument, and allows the user to specify the precision with
03795     // which the computation is carried out.
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     // This function returns an Array1D of the specified size and type
03810     // in which the value of every element is zero.
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     // This function returns an Array2D of the specified size and type
03821     // in which the value of every element is zero.
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     // This function returns an Array3D of the specified size and type
03833     // in which the value of every element is zero.
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   } // namespace numeric
03844 
03845 } // namespace dlr
03846 
03847 #endif /* #ifndef DLR_NUMERIC_UTILITIES_H */

Generated on Wed Nov 25 00:42:43 2009 for dlrUtilities Utility Library by  doxygen 1.5.8