OpenShot Audio Library | OpenShotAudio  0.3.0
juce_Matrix.cpp
1 /*
2  ==============================================================================
3 
4  This file is part of the JUCE library.
5  Copyright (c) 2017 - ROLI Ltd.
6 
7  JUCE is an open source library subject to commercial or open-source
8  licensing.
9 
10  By using JUCE, you agree to the terms of both the JUCE 5 End-User License
11  Agreement and JUCE 5 Privacy Policy (both updated and effective as of the
12  27th April 2017).
13 
14  End User License Agreement: www.juce.com/juce-5-licence
15  Privacy Policy: www.juce.com/juce-5-privacy-policy
16 
17  Or: You may also use this code under the terms of the GPL v3 (see
18  www.gnu.org/licenses).
19 
20  JUCE IS PROVIDED "AS IS" WITHOUT ANY WARRANTY, AND ALL WARRANTIES, WHETHER
21  EXPRESSED OR IMPLIED, INCLUDING MERCHANTABILITY AND FITNESS FOR PURPOSE, ARE
22  DISCLAIMED.
23 
24  ==============================================================================
25 */
26 
27 namespace juce
28 {
29 namespace dsp
30 {
31 
32 template <typename ElementType>
34 {
35  Matrix result (size, size);
36 
37  for (size_t i = 0; i < size; ++i)
38  result(i, i) = 1;
39 
40  return result;
41 }
42 
43 template <typename ElementType>
45 {
46  jassert (vector.isOneColumnVector());
47  jassert (size <= vector.rows);
48 
49  Matrix result (size, size);
50 
51  for (size_t i = 0; i < size; ++i)
52  result (i, i) = vector (0, 0);
53 
54  for (size_t i = 1; i < size; ++i)
55  for (size_t j = i; j < size; ++j)
56  result (j, j - i) = result (j - i, j) = vector (i, 0);
57 
58  return result;
59 }
60 
61 template <typename ElementType>
62 Matrix<ElementType> Matrix<ElementType>::hankel (const Matrix& vector, size_t size, size_t offset)
63 {
64  jassert(vector.isOneColumnVector());
65  jassert(vector.rows >= (2 * (size - 1) + 1));
66 
67  Matrix result (size, size);
68 
69  for (size_t i = 0; i < size; ++i)
70  result (i, i) = vector ((2 * i) + offset, 0);
71 
72  for (size_t i = 1; i < size; ++i)
73  for (size_t j = i; j < size; ++j)
74  result (j, j - i) = result (j - i, j) = vector (i + 2 * (j - i) + offset, 0);
75 
76  return result;
77 }
78 
79 //==============================================================================
80 template <typename ElementType>
81 Matrix<ElementType>& Matrix<ElementType>::swapColumns (size_t columnOne, size_t columnTwo) noexcept
82 {
83  jassert (columnOne < columns && columnTwo < columns);
84 
85  auto* p = data.getRawDataPointer();
86 
87  for (size_t i = 0; i < rows; ++i)
88  {
89  auto offset = dataAcceleration.getUnchecked (static_cast<int> (i));
90  std::swap (p[offset + columnOne], p[offset + columnTwo]);
91  }
92 
93  return *this;
94 }
95 
96 template <typename ElementType>
97 Matrix<ElementType>& Matrix<ElementType>::swapRows (size_t rowOne, size_t rowTwo) noexcept
98 {
99  jassert (rowOne < rows && rowTwo < rows);
100 
101  auto offset1 = rowOne * columns;
102  auto offset2 = rowTwo * columns;
103 
104  auto* p = data.getRawDataPointer();
105 
106  for (size_t i = 0; i < columns; ++i)
107  std::swap (p[offset1 + i], p[offset2 + i]);
108 
109  return *this;
110 }
111 
112 //==============================================================================
113 template <typename ElementType>
115 {
116  auto n = getNumRows(), m = other.getNumColumns(), p = getNumColumns();
117  Matrix result (n, m);
118 
119  jassert (p == other.getNumRows());
120 
121  size_t offsetMat = 0, offsetlhs = 0;
122 
123  auto* dst = result.getRawDataPointer();
124  auto* a = getRawDataPointer();
125  auto* b = other.getRawDataPointer();
126 
127  for (size_t i = 0; i < n; ++i)
128  {
129  size_t offsetrhs = 0;
130 
131  for (size_t k = 0; k < p; ++k)
132  {
133  auto ak = a[offsetlhs++];
134 
135  for (size_t j = 0; j < m; ++j)
136  dst[offsetMat + j] += ak * b[offsetrhs + j];
137 
138  offsetrhs += m;
139  }
140 
141  offsetMat += m;
142  }
143 
144  return result;
145 }
146 
147 //==============================================================================
148 template <typename ElementType>
149 bool Matrix<ElementType>::compare (const Matrix& a, const Matrix& b, ElementType tolerance) noexcept
150 {
151  if (a.rows != b.rows || a.columns != b.columns)
152  return false;
153 
154  tolerance = std::abs (tolerance);
155 
156  auto* bPtr = b.begin();
157  for (auto aValue : a)
158  if (std::abs (aValue - *bPtr++) > tolerance)
159  return false;
160 
161  return true;
162 }
163 
164 //==============================================================================
165 template <typename ElementType>
166 bool Matrix<ElementType>::solve (Matrix& b) const noexcept
167 {
168  auto n = columns;
169  jassert (n == n && n == b.rows && b.isOneColumnVector());
170 
171  auto* x = b.getRawDataPointer();
172  const auto& A = *this;
173 
174  switch (n)
175  {
176  case 1:
177  {
178  auto denominator = A (0,0);
179 
180  if (denominator == 0)
181  return false;
182 
183  b (0, 0) /= denominator;
184  }
185  break;
186 
187  case 2:
188  {
189  auto denominator = A (0, 0) * A (1, 1) - A (0, 1) * A (1, 0);
190 
191  if (denominator == 0)
192  return false;
193 
194  auto factor = (1 / denominator);
195  auto b0 = x[0], b1 = x[1];
196 
197  x[0] = factor * (A (1, 1) * b0 - A (0, 1) * b1);
198  x[1] = factor * (A (0, 0) * b1 - A (1, 0) * b0);
199  }
200  break;
201 
202  case 3:
203  {
204  auto denominator = A (0, 0) * (A (1, 1) * A (2, 2) - A (1, 2) * A (2, 1))
205  + A (0, 1) * (A (1, 2) * A (2, 0) - A (1, 0) * A (2, 2))
206  + A (0, 2) * (A (1, 0) * A (2, 1) - A (1, 1) * A (2, 0));
207 
208  if (denominator == 0)
209  return false;
210 
211  auto factor = 1 / denominator;
212  auto b0 = x[0], b1 = x[1], b2 = x[2];
213 
214  x[0] = ( ( A (0, 1) * A (1, 2) - A (0, 2) * A (1, 1)) * b2
215  + (-A (0, 1) * A (2, 2) + A (0, 2) * A (2, 1)) * b1
216  + ( A (1, 1) * A (2, 2) - A (1, 2) * A (2, 1)) * b0) * factor;
217 
218  x[1] = -( ( A (0, 0) * A (1, 2) - A (0, 2) * A (1, 0)) * b2
219  + (-A (0, 0) * A (2, 2) + A (0, 2) * A (2, 0)) * b1
220  + ( A (1, 0) * A (2, 2) - A (1, 2) * A (2, 0)) * b0) * factor;
221 
222  x[2] = ( ( A (0, 0) * A (1, 1) - A (0, 1) * A (1, 0)) * b2
223  + (-A (0, 0) * A (2, 1) + A (0, 1) * A (2, 0)) * b1
224  + ( A (1, 0) * A (2, 1) - A (1, 1) * A (2, 0)) * b0) * factor;
225  }
226  break;
227 
228 
229  default:
230  {
231  Matrix<ElementType> M (A);
232 
233  for (size_t j = 0; j < n; ++j)
234  {
235  if (M (j, j) == 0)
236  {
237  auto i = j;
238  while (i < n && M (i, j) == 0)
239  ++i;
240 
241  if (i == n)
242  return false;
243 
244  for (size_t k = 0; k < n; ++k)
245  M (j, k) += M (i, k);
246 
247  x[j] += x[i];
248  }
249 
250  auto t = 1 / M (j, j);
251 
252  for (size_t k = 0; k < n; ++k)
253  M (j, k) *= t;
254 
255  x[j] *= t;
256 
257  for (size_t k = j + 1; k < n; ++k)
258  {
259  auto u = -M (k, j);
260 
261  for (size_t l = 0; l < n; ++l)
262  M (k, l) += u * M (j, l);
263 
264  x[k] += u * x[j];
265  }
266  }
267 
268  for (int k = static_cast<int> (n) - 2; k >= 0; --k)
269  for (size_t i = static_cast<size_t> (k) + 1; i < n; ++i)
270  x[k] -= M (static_cast<size_t> (k), i) * x[i];
271  }
272  }
273 
274  return true;
275 }
276 
277 //==============================================================================
278 template <typename ElementType>
280 {
281  StringArray entries;
282  int sizeMax = 0;
283 
284  auto* p = data.begin();
285 
286  for (size_t i = 0; i < rows; ++i)
287  {
288  for (size_t j = 0; j < columns; ++j)
289  {
290  String entry (*p++, 4);
291  sizeMax = jmax (sizeMax, entry.length());
292 
293  entries.add (entry);
294  }
295  }
296 
297  sizeMax = ((sizeMax + 1) / 4 + 1) * 4;
298 
299  MemoryOutputStream result;
300 
301  auto n = static_cast<size_t> (entries.size());
302 
303  for (size_t i = 0; i < n; ++i)
304  {
305  result << entries[(int) i].paddedRight (' ', sizeMax);
306 
307  if (i % columns == (columns - 1))
308  result << newLine;
309  }
310 
311  return result.toString();
312 }
313 
314 template class Matrix<float>;
315 template class Matrix<double>;
316 
317 } // namespace dsp
318 } // namespace juce
String * begin() noexcept
int size() const noexcept
void add(String stringToAdd)
int length() const noexcept
Matrix & swapRows(size_t rowOne, size_t rowTwo) noexcept
Definition: juce_Matrix.cpp:97
size_t getNumRows() const noexcept
Definition: juce_Matrix.h:93
static bool compare(const Matrix &a, const Matrix &b, ElementType tolerance=0) noexcept
ElementType * getRawDataPointer() noexcept
Definition: juce_Matrix.h:131
bool isOneColumnVector() const noexcept
Definition: juce_Matrix.h:185
Matrix & swapColumns(size_t columnOne, size_t columnTwo) noexcept
Definition: juce_Matrix.cpp:81
static Matrix hankel(const Matrix &vector, size_t size, size_t offset=0)
Definition: juce_Matrix.cpp:62
Matrix operator*(ElementType scalar) const
Definition: juce_Matrix.h:159
static Matrix identity(size_t size)
Definition: juce_Matrix.cpp:33
bool solve(Matrix &b) const noexcept
static Matrix toeplitz(const Matrix &vector, size_t size)
Definition: juce_Matrix.cpp:44
size_t getNumColumns() const noexcept
Definition: juce_Matrix.h:96
String toString() const