Entropy Piano Tuner  1.1.3 (documentation not yet complete) An open-source experimental software for piano tuning by entropy minimization
mathtools.cpp
Go to the documentation of this file.
1 /*****************************************************************************
2  * Copyright 2015 Haye Hinrichsen, Christoph Wick
3  *
4  * This file is part of Entropy Piano Tuner.
5  *
6  * Entropy Piano Tuner is free software: you can redistribute it and/or modify
8  * Free Software Foundation, either version 3 of the License, or (at your
9  * option) any later version.
10  *
11  * Entropy Piano Tuner is distributed in the hope that it will be useful, but
12  * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
13  * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
14  * more details.
15  *
16  * You should have received a copy of the GNU General Public License along with
17  * Entropy Piano Tuner. If not, see http://www.gnu.org/licenses/.
18  *****************************************************************************/
19
20 //=============================================================================
21 // Various mathematical tools
22 //=============================================================================
23
24 #include "mathtools.h"
25
26 #include <algorithm>
27 #include <numeric>
28
29 #include "../system/eptexception.h"
30
31 //-----------------------------------------------------------------------------
32 // Determine the raw moments of a distribution
33 //-----------------------------------------------------------------------------
34
37
38 double MathTools::computeMoment (const std::vector<double> &v, const int n)
39 {
40  EptAssert (n>=0 and v.size()>0, "Parameters inconsistent");
41  double norm=0,sum=0;
42  for (size_t i=0; i<v.size(); ++i)
43  {
44  norm += v[i];
45  sum += v[i] * pow(1.0*i,1.0*n);
46  }
47  EptAssert (norm>0, "The norm should be positive");
48  return sum / norm;
49 }
50
51
52 //-----------------------------------------------------------------------------
53 // Compute Shannon entropy
54 //-----------------------------------------------------------------------------
55
56 double MathTools::computeEntropy (const std::vector<double> &v)
57 {
58  assert(v.size()>0);
59  double sum=0;
60  for (auto x : v) sum -= (x>0 ? x*log(x) : 0);
61  return sum;
62 }
63
64
65 //-----------------------------------------------------------------------------
66 // Compute Renyi entropy
67 //-----------------------------------------------------------------------------
68
69 double MathTools::computeRenyiEntropy (const std::vector<double> &v, const double q)
70 {
71  assert(v.size()>0 and q != 1);
72  double sum=0;
73  for (auto x : v) sum += pow(x,q);
74  return log(sum) / (1-q);
75 }
76
77
78 //-----------------------------------------------------------------------------
79 // Compute the norm
80 //-----------------------------------------------------------------------------
81
82 double MathTools::computeNorm (std::vector<double> &vec)
83 {
84  return std::accumulate(vec.begin(), vec.end(), 0.0);
85 }
86
87
88 //-----------------------------------------------------------------------------
89 // Normalize a distribution
90 //-----------------------------------------------------------------------------
91
92 void MathTools::normalize (std::vector<double> &vec)
93 {
94  double norm = computeNorm(vec);
95  EptAssert (norm!=0,"Vectors with norm zero cannot be normalized");
96  for (auto &x : vec) x /= norm;
97 }
98
99
100 //-----------------------------------------------------------------------------
101 // Map a probability distribution from one vector to another
102 //-----------------------------------------------------------------------------
103
118
119
120 void MathTools::coarseGrainSpectrum (const std::vector<double> &X,
121  std::vector<double> &Y,
122  std::function<double(double y)> f,
123  double exponent)
124 {
125  assert(X.size()>0 and Y.size()>0);
126  double xs1 = f(-0.5);
127  int x1 = std::max<int>(0, roundToInteger(xs1));
128  double leftarea = (x1-xs1+0.5)*X[x1];
129  for (int y=0; y<static_cast<int>(Y.size()); ++y)
130  {
131  double xs2 = f(y+0.5);
132  int x2 = std::min<int>(roundToInteger(xs2), X.size() - 1);
133  double sum=0;
134  for (int x=x1+1; x<=x2; ++x) sum += X[x];
135  double rightarea = (x2-xs2+0.5)*X[x2];
136  Y[y] = (sum + leftarea - rightarea) * pow(xs1*xs2,exponent);
137  x1=x2; xs1=xs2; leftarea=rightarea;
138  }
139 }
140
141 //-----------------------------------------------------------------------------
142 // Find the maximum in a vector
143 //-----------------------------------------------------------------------------
144
150
151 int MathTools::findMaximum (const std::vector<double> &X, int i, int j)
152 {
153
154  int N = X.size();
155  assert (i>=0 and i<N and j>i and j<=N);
156  return std::distance(X.begin(), std::max_element(X.begin()+i,X.begin()+j));
157 }
158
159 int MathTools::findMaximum (const std::vector<double> &X)
160 {
161  return std::distance(X.begin(), std::max_element(X.begin(),X.end()));
162 }
163
164 double MathTools::findSmoothedMaximum (const std::vector<double> &x)
165 {
166  auto maxElem = std::max_element(std::next(x.begin()), std::prev(x.end()));
167  int maxIndex = std::distance(x.begin(), maxElem);
168  // interpolate with neighbours using a parabolas maximum
169  double y1 = *std::prev(maxElem);
170  double y2 = *maxElem;
171  double y3 = *std::next(maxElem);
172  double x1 = maxIndex - 1, x2 = maxIndex, x3 = maxIndex + 1;
173
174  // compute max of parabola
175  return (x3 * x3 * (y2 - y1) + x2 * x2 * (y1 - y3) + x1 * x1 * (y3 - y2)) / 2.0 / (x3 * (y2 - y1) + x2 * (y1 - y3) + x1 * (y3 - y2));
176 }
177
178 double MathTools::weightedArithmetricMean(const std::vector<double> &Y, size_t start, size_t end)
179 {
180  double norm = 0;
181  double mean = 0;
182  end = std::min(end, Y.size());
183  for (size_t i = start; i < end; ++i)
184  {
185  norm += Y[i] * Y[i];
186  mean += Y[i] * Y[i] * i;
187  }
188
189  return mean / norm;
190 }
191
192 //-----------------------------------------------------------------------------
193 // Restrict floating point value to an interval
194 //-----------------------------------------------------------------------------
195
196 double MathTools::restrictToInterval(double x, double xmin, double xmax)
197 {
198  EptAssert(xmax > xmin, "xmax should be larger than xmin");
199  return std::max(xmin, std::min(xmax, x));
200 }
double restrictToInterval(double x, double xmin, double xmax)
Restrict floating point value to an interval.
Definition: mathtools.cpp:196
double computeMoment(const std::vector< double > &v, const int n)
Determine the n-th moment of a distribution stored in a vector.
Definition: mathtools.cpp:38
double computeRenyiEntropy(const std::vector< double > &v, const double q)
Compute the Renyi entropy of a normalized probability distribution.
Definition: mathtools.cpp:69
double weightedArithmetricMean(const std::vector< double > &Y, size_t start=0, size_t end=std::numeric_limits< size_t >::max())
Computes the weighted arithmetric mean index of the given Y data.
Definition: mathtools.cpp:178
int roundToInteger(T x)
Round a floating point number to an integer.
Definition: mathtools.h:43
void coarseGrainSpectrum(const std::vector< double > &X, std::vector< double > &Y, std::function< double(double y)> f, double exponent=0)
Definition: mathtools.cpp:120
int findMaximum(const std::vector< double > &X, int i, int j)
Find the component where the vector has its maximum.
Definition: mathtools.cpp:151
double findSmoothedMaximum(const std::vector< double > &x)
Use a parabola to fit the maximum.
Definition: mathtools.cpp:164
double computeEntropy(const std::vector< double > &v)
Compute the Shannon entropy of a normalized probability distribution.
Definition: mathtools.cpp:56
void normalize(std::vector< double > &vec)
Normalize a distribution stored in a vector.
Definition: mathtools.cpp:92
#define EptAssert(a, b)
Definition: eptexception.h:47
double computeNorm(std::vector< double > &vec)
Compute the norm of a vector.
Definition: mathtools.cpp:82