// ArrayMechanisms.cpp // // 2003-8-1 Initial code DD // 2003-80-6 DD major code review and changes // 2004-4-3 DD testing etc. // 2004-4-9 DD Code review and final fiat (for time being) // 2005-6-19 DD print function for Tensor // 2005-1-8 DD big fix in l2Norm (wrong name) // 2007-6-25 DD major new update, code // // (C) Datasim Education BV 2003-2007 // #ifndef ArrayMechanisms_cpp #define ArrayMechanisms_cpp #include "ArrayMechanisms.hpp" #include #include #include #include using namespace std; /////////////////////////////////////////////////////////////////////////////////////////////// // Sums and averages template V sum(const Vector& x) { // Sum of elements V ans = x[x.MinIndex()]; for (I j = x.MinIndex()+1; j <= x.MaxIndex(); ++j) { ans += x[j]; } return ans; } template V product(const Vector& x) { // Product of elements V ans = x[x.MinIndex()]; for (I j = x.MinIndex() + 1; j <= x.MaxIndex(); ++j) { ans *= x[j]; } return ans; } template V sumReciprocals(const Vector& x) { // Sum of reciprocals // Precondition (PREC): x is (strictly) positive V ans = x[x.MinIndex()]; for (I j = x.MinIndex() + 1; j <= x.MaxIndex(); ++j) { ans += 1.0/x[j]; } return ans; } template V sumAbsoluteValues(const Vector& x) { // Sum of the absolute values of a vector V ans = fabs(x[x.MinIndex()]); for (I j = x.MinIndex()+1; j <= x.MaxIndex(); ++j) { ans += fabs(x[j]); } return ans; } // Mean value == sum() / N template V arithmeticMean(const Vector& x) { return sum(x) / V(x.Size()); } // Weighted arithmetic mean template V weightedArithmeticMean(const Vector& x, const Vector& w) { // PREC: x and w have the same size; start indexes not necessarily the same // PREC: sum(w) is not zero V ans = w[w.MinIndex()]* x[x.MinIndex()]; for (I j = x.MinIndex() +1; j <= x.MaxIndex(); ++j) { ans += w[j] * x[j]; } return ans / sum(w); } // Geometric mean or geometric average template V geometricMean(const Vector& x) { V ans = product(x); return pow(ans, 1.0 / x.Size() ); } // Harmonic mean template V harmonicMean(const Vector& x) { // PREC: sumReciprocals(x) not zero return V(x.Size()) / sumReciprocals(x); } // Root mean square (RMS) template V quadraticMean(const Vector& x) { return sqrt( sumSquares(x) / V(x.Size()) ); } // Sum of squares template V sumSquares(const Vector& x) { V ans = x[x.MinIndex()]; for (I j = x.MinIndex()+1; j <= x.MaxIndex(); ++j) { ans += (x[j] * x[j]); } return ans; } // A function returning all of the above values in one foul swoop (performance) template SimplePropertySet allAverages(const Vector& x) { SimplePropertySet result; // Empty list result.add(Property ("SUM", sum(x))); result.add(Property ("PRODUCT", product(x))); result.add(Property ("SUMREC", sumReciprocals(x))); result.add(Property ("SUMABSVAL", sumAbsoluteValues(x))); result.add(Property ("MEAN", arithmeticMean(x))); result.add(Property ("GMEAN", geometricMean(x))); result.add(Property ("HMEAN", harmonicMean(x))); result.add(Property ("RMS", quadraticMean(x))); result.add(Property ("SUMSQ", sumSquares(x))); return result; } ///////////////////////////////////////////////////////////////////////////// // Measures of Dispersion template V deviationFromMean(const Vector& x) { V avg = mean(x); V ans = x[x.MinIndex()]; for (I j = x.MinIndex() + 1; j <= x.MaxIndex(); ++j) { ans += fabs(x[j] - avg); } return ans / x.Size(); } template V standardDeviation(const Vector& x) { V avg = mean(x); V tmp; V ans = x[x.MinIndex()]; for (I j = x.MinIndex() + 1; j <= x.MaxIndex(); ++j) { tmp = x[j] - avg; ans += tmp*tmp; } return sqrt(ans / V(x.Size())); } template V variance(const Vector& x) { V s = standardDeviation(x); return s*s; } // A function returning all of the above values in one foul swoop (performance) template SimplePropertySet allDispersions(const Vector& x) { SimplePropertySet result; // Empty list result.add(Property ("MDEV", deviationFromMean(x))); result.add(Property ("STD", standardDeviation())); result.add(Property ("VARIANCE", variance(x))); return result; } //////////////////////////////////////////////////////////////////////////////////////////////// // Moments, Skewness and Kurtosis // The rth moment about the value 0.0 template V rthMoment(const Vector& x, const I& r) { return rthMoment(x, r, 0.0); } // The rth moment about the Mean m(r) as a special origin template V rthMomentMean(const Vector& x, const I& r) { return rthMoment(x, r, mean(x)); } template V rthMomentMean(const Vector& x, const Vector& freq, const I& r) { return rthMoment(x, freq, r, mean(x)); } // The rth moment about an origin A template V rthMoment(const Vector& x, const I& r, const V& A) { V ans = pow(x[x.MinIndex()] - A, pr); V pr = V(r); for (I j = x.MinIndex() + 1; j <= x.MaxIndex(); ++j) { ans += pow(x[j] - A, pr); } return ans / x.Size(); } template V rthMoment(const Vector& x, const Vector& freq, const I& r, const V& A) { V ans = freq[freq.MinIndex()] * pow(x[x.MinIndex()] - A, pr); V vr = V(r); for (I j = x.MinIndex() + 1; j <= x.MaxIndex(); ++j) { ans += freq[j] * pow(x[j] - A, vr); } return ans / sum(freq); } // Mode: either the middel element or the mean of the two middle elements template V median(const Vector& x) { // Work around; ideally Vector should have a sort() vector v = createSTLvector(x); stable_sort(v.begin(), v.end()); // Now v is sorted int N = v.size(); if ((N/2)*2 == N) // Even list { return (v[N/2] + v[(N/2) - 1]) * V(0.5); } else { return v[N/2]; } } // Number of occurrences of value d in vector x template I occurs(const Vector& x, const V& d) { I result = I(0); for (I j = x.MinIndex(); j <= x.MaxIndex(); ++j) { if (d == x[j]) result++; } return result; } // The element in x that occurs with the greatest frequency template V mode(const Vector& x) { // The value that occurs with the greatest frequency // Exx. Code can be optimised map elements; // Create unique list of elements for (I j = x.MinIndex(); j <= x.MaxIndex(); ++j) { elements[x[j]] = occurs(x, x[j]); } map::const_iterator i = elements.begin(); I ans = (*i).second; V result = (*i).first; while (i != elements.end()) { if (ans < (*i).second) result = (*i).second; i++; } return result; } template V skewness(const Vector& x) { //return (mean(x) - mode(x)) / standardDeviation(x); // Good example for parallel processing return 3.0 * (arithmeticMean(x) - median(x))/ standardDeviation(x); } //////////////////////////////////////////////////////////////////////////////////////////////////// // Extremum operations on vectors template V maxValue(const Vector& x) { V ans = x[x.MinIndex()]; for (I j = x.MinIndex() + 1; j <= x.MaxIndex(); ++j) { if (ans < x[j]) ans = x[j]; } return ans; } template V minValue(const Vector& x) { V ans = x[x.MinIndex()]; for (I j = x.MinIndex() + 1; j <= x.MaxIndex(); ++j) { if (ans > x[j]) ans = x[j]; } return ans; } // Max and min of the absolute values template V maxAbsValue(const Vector& x) { V ans = fabs(x[x.MinIndex()]); for (I j = x.MinIndex() + 1; j <= x.MaxIndex(); ++j) { if (ans < fabs(x[j])) ans = fabs(x[j]); } return ans; } template V minAbsValue(const Vector& x) { V ans = fabs(x[x.MinIndex()]); for (I j = x.MinIndex() + 1; j <= x.MaxIndex(); ++j) { if (ans > fabs(x[j])) ans = fabs(x[j]); } return ans; } // Index of max and min values template I indexMaxValue(const Vector& x) { I index = x.MinIndex(); V ans = x[x.MinIndex()]; for (I j = x.MinIndex() + 1; j <= x.MaxIndex(); ++j) { if (ans < x[j]) { index = j; ans = x[j]; } } return index; } template V indexMinValue(const Vector& x) { I index = x.MinIndex(); V ans = x[x.MinIndex()]; for (I j = x.MinIndex() + 1; j <= x.MaxIndex(); ++j) { if (ans > x[j]) { index = j; ans = x[j]; } } return index; } template I indexMaxAbsValue(const Vector& x) { I index = x.MinIndex(); V ans = fabs(x[x.MinIndex()]); for (I j = x.MinIndex() + 1; j <= x.MaxIndex(); ++j) { if (ans < fabs(x[j])) { index = j; ans = fabs(x[j]); } } return index; } template V indexMinAbsValue(const Vector& x) { I index = x.MinIndex(); V ans = fabs(x[x.MinIndex()]); for (I j = x.MinIndex() + 1; j <= x.MaxIndex(); ++j) { if (ans > fabs(x[j])) { index = j; ans = fabs(x[j]); } } return index; } // Vector-vector extremum (difference of two vectors) template V maxValue(const Vector& vectorA, const Vector& vectorB) { // PREC: A and B have same size (holds for all the following functions too) Vector vecDiff = vectorA - vectorB; return maxValue(vecDiff); } template V minValue(const Vector& vectorA, const Vector& vectorB) { Vector vecDiff = vectorA - vectorB; return minValue(vecDiff); } template V maxAbsValue(const Vector& vectorA, const Vector& vectorB) { Vector vecDiff = vectorA - vectorB; return maxAbsValue(vecDiff); } template V minAbsValue(const Vector& vectorA, const Vector& vectorB) { Vector vecDiff = vectorA - vectorB; return minAbsValue(vecDiff); } /////////////////////////////////////////////////////////////////////////////////////////////////////// // Vector and matrix norms template V innerProduct(const Vector& x, const Vector& y) { // PREC: x and y have same size V ans = x[x.MinIndex()] * y[y.MinIndex()]; for (I j = x.MinIndex() + 1; j <= x.MaxIndex(); ++j) { ans += x[j] * y[j]; } return ans; } template V l1Norm(const Vector& x) { return sumAbsoluteValues(x); } // template V ManyToManyRelationorm(const Vector& x) template V l2Norm(const Vector& x) { return sqrt(sumSquares(x)); } template V lpNorm(const Vector& x, const I& p) { V myPower = V(p); V ans = pow(x[x.MinIndex()], myPower); for (I j = x.MinIndex()+ 1; j <= x.MaxIndex(); ++j) { ans += pow(x[j], myPower);; } return pow(ans, V(1.0)/ myPower); } template V lInfinityNorm(const Vector& x) { return maxAbsValue(x); } template SimplePropertySet allNorms(const Vector& x) { SimplePropertySet result; // Empty list result.add(Property ("l1", l1Norm(x))); result.add(Property ("l2", ManyToManyRelationorm(x))); result.add(Property ("linf", lInfinityNorm(x))); return result; } template SimplePropertySet allNorms(const Vector& vectorA, const Vector& vectorB) { Vector vecDiff = vectorA - vectorB; return allNorms(vecDiff); } // Same vector morms as above except for the difference of two vectors template V l1Norm(const Vector& vectorA, const Vector& vectorB) { Vector vecDiff = vectorA - vectorB; return l1Norm(vecDiff); } template V l2Norm(const Vector& vectorA, const Vector& vectorB) { Vector vecDiff = vectorA - vectorB; return l2Norm(vecDiff); } template V lpNorm(const Vector& vectorA, const Vector& vectorB, const I& p) { Vector vecDiff = vectorA - vectorB; return lpNorm(vecDiff, p); } template V lInfinityNorm(const Vector& vectorA, const Vector& vectorB) { Vector vecDiff = vectorA - vectorB; return lInfinityNorm(vecDiff); } //////////////////////////////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////////////////////////// // Functions for Operations Research // Comparing vectors with each other // Are all elements of a vector positive? template bool positive(const Vector& x) { V zero = V(0.0); for (I j = x.MinIndex(); j <= x.MaxIndex(); ++j) { if (x[j] <= zero) return false; } return true; } template bool negative(const Vector& x) { V zero = V(0.0); for (I j = x.MinIndex(); j <= x.MaxIndex(); ++j) { if (x[j] >= zero) return false; } return true; } // Is v1 < v2? etc. template bool operator < (const Vector& v1, const Vector& v2) { // Iterate in the matrix; when the condition is NOT true // (the inverse of the inequality) then exit and return false for (I j = v1.MinIndex(); j <= v1.MaxIndex(); ++j) { if (v1[j] >= v2[j]) return false; } return true; } template bool operator <= (const Vector& v1, const Vector& v2) { // Iterate in the matrix; when the condition is NOT true // (the inverse of the inequality) then exit and return false for (I j = v1.MinIndex(); j <= v1.MaxIndex(); ++j) { if (v1[j] > v2[j]) return false; } return true; } template bool operator > (const Vector& v1, const Vector& v2) { // Iterate in the matrix; when the condition is NOT true // (the inverse of the inequality) then exit and return false for (I j = v1.MinIndex(); j <= v1.MaxIndex(); ++j) { if (v1[j] <= v2[j]) return false; } return true; } template bool operator >= (const Vector& v1, const Vector& v2) { // Iterate in the matrix; when the condition is NOT true // (the inverse of the inequality) then exit and return false for (I j = v1.MinIndex(); j <= v1.MaxIndex(); ++j) { if (v1[j] < v2[j]) return false; } return true; } template bool operator == (const Vector& v1, const Vector& v2) { // Iterate in the matrix; when the condition is NOT true // (the inverse of the inequality) then exit and return false for (I j = v1.MinIndex(); j <= v1.MaxIndex(); ++j) { if (v1[j] != v2[j]) return false; } return true; } template bool operator != (const Vector& v1, const Vector& v2) { if (v1 == v2) return false; return true; } // Utility functions template vector createSTLvector (const Vector& myVector) { // Create an STL vector from a general Vector vector result(myVector.Size()); for (int i = 0; i < result.size(); i++) { result[i] = myVector[i+myVector.MinIndex()]; } return result; } template Vector createDatasimVector (const vector& mySTLvector) { // Create a general Vector from an STL vector Vector result(mySTLvector.size()); for (int i = 0, j = result.MinIndex(); i < mySTLvector.size(); i++, ++j) { result[j] = mySTLvector[i]; } return result; } template Vector cumulativeVector (const Vector& x) { // Cumulative vector c[j] = c[j-1] + x[j] Vector result(x.Size(), x.MaxIndex()); result[x.MinIndex()] = x[x.MinIndex()]; for (int i = x.MinIndex() + 1; i <= x.MaxIndex(); i++) { result[i] = result[i-1] + x[i]; } return result; } template Vector reverse (const Vector& x) { Vector result(x); I min = x.MInIndex(); for (I k = result.MinIndex(); k <= result.MaxIndex(); k++) { result[k] = x[x.MaxIndex() - k + min]; } return result; } ////////////////////// Print Functions ///////////////////////////////////////// template void print(const Array& v) { cout << "\n\nMinIndex: " << v.MinIndex() << " , MaxIndex: " << v.MaxIndex() << endl; cout << "\nARR:["; for (I j = v.MinIndex(); j <= v.MaxIndex(); ++j) { cout << v[j] << ", "; } cout << "]"; } template void print(const Vector& v) { cout << "\n\nMinIndex: " << v.MinIndex() << " , MaxIndex: " << v.MaxIndex() << endl; cout << "\nARR:["; for (I j = v.MinIndex(); j <= v.MaxIndex(); ++j) { cout << v[j] << ", "; } cout << "]"; } template void print (Tensor& tensor) { cout << "Tensor, Rows " << tensor.Rows() << ", Columns " << tensor.Columns() << ", Third Dimension " << tensor.sizeThird() << endl; for (long k = tensor.MinThirdIndex(); k<= tensor.MaxThirdIndex(); k++) { print (tensor[k]); } } #endif