NumeRe v1.1.4
NumeRe: Framework für Numerische Rechnungen
command_implementations.cpp
Go to the documentation of this file.
1/*****************************************************************************
2 NumeRe: Framework fuer Numerische Rechnungen
3 Copyright (C) 2019 Erik Haenel et al.
4
5 This program is free software: you can redistribute it and/or modify
6 it under the terms of the GNU General Public License as published by
7 the Free Software Foundation, either version 3 of the License, or
8 (at your option) any later version.
9
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU General Public License for more details.
14
15 You should have received a copy of the GNU General Public License
16 along with this program. If not, see <http://www.gnu.org/licenses/>.
17******************************************************************************/
18
19#include <gsl/gsl_statistics.h>
20#include <gsl/gsl_sort.h>
21#include <algorithm>
22#include <memory>
23
24#include "command_implementations.hpp"
25#include "parser_functions.hpp"
26#include "matrixoperations.hpp"
27#include "spline.h"
28#include "wavelet.hpp"
29#include "filtering.hpp"
30#include "../AudioLib/audiofile.hpp"
31#include "../../kernel.hpp"
32#include "../../../common/http.h"
33
34#define TRAPEZOIDAL 1
35#define SIMPSON 2
36
37using namespace std;
38
39DefaultVariables _defVars;
40static mu::value_type localizeExtremum(string& sCmd, mu::value_type* dVarAdress, Parser& _parser, const Settings& _option, mu::value_type dLeft, mu::value_type dRight, double dEps = 1e-10, int nRecursion = 0);
41static mu::value_type localizeZero(string& sCmd, mu::value_type* dVarAdress, Parser& _parser, const Settings& _option, mu::value_type dLeft, mu::value_type dRight, double dEps = 1e-10, int nRecursion = 0);
42static std::vector<size_t> getSamplesForDatagrid(CommandLineParser& cmdParser, const std::vector<size_t>& nSamples);
43static void expandVectorToDatagrid(IntervalSet& ivl, std::vector<std::vector<mu::value_type>>& vZVals, size_t nSamples_x, size_t nSamples_y);
44
45
60static void integrationstep_trapezoidal(mu::value_type& x, mu::value_type x0, mu::value_type dx, vector<mu::value_type>& vResult, vector<mu::value_type>& vFunctionValues, bool bReturnFunctionPoints)
61{
62 x = x0;
63 int nResults;
64 mu::value_type* v = NumeReKernel::getInstance()->getParser().Eval(nResults);
65
66 // Evaluate the current integration step for each of the
67 // defined functions
68 for (int i = 0; i < nResults; i++)
69 {
70 if (isnan(v[i].real()) || isnan(v[i].imag()))
71 v[i] = 0.0;
72 }
73
74 // Now calculate the area below the curve
75 if (!bReturnFunctionPoints)
76 {
77 for (int i = 0; i < nResults; i++)
78 vResult[i] += dx * (vFunctionValues[i] + v[i]) * 0.5;
79 }
80 else
81 {
82 // Calculate the integral
83 if (vResult.size())
84 vResult.push_back(dx * (vFunctionValues[0] + v[0]) * 0.5 + vResult.back());
85 else
86 vResult.push_back(dx * (vFunctionValues[0] + v[0]) * 0.5);
87 }
88
89 // Set the last sample as the first one
90 vFunctionValues.assign(v, v+nResults);
91}
92
93
109static void integrationstep_simpson(mu::value_type& x, mu::value_type x0, mu::value_type x1, double dx, vector<mu::value_type>& vResult, vector<mu::value_type>& vFunctionValues, bool bReturnFunctionPoints)
110{
111 // Evaluate the intermediate function value
112 x = x0;
113 int nResults;
114 mu::value_type* v = NumeReKernel::getInstance()->getParser().Eval(nResults);
115
116 for (int i = 0; i < nResults; i++)
117 {
118 if (isnan(std::abs(v[i])))
119 v[i] = 0.0;
120 }
121
122 vector<mu::value_type> vInter(v, v+nResults);
123
124 // Evaluate the end function value
125 x = x1;
126 v = NumeReKernel::getInstance()->getParser().Eval(nResults);
127
128 for (int i = 0; i < nResults; i++)
129 {
130 if (isnan(std::abs(v[i])))
131 v[i] = 0.0;
132 }
133
134 // Now calculate the area below the curve
135 if (!bReturnFunctionPoints)
136 {
137 for (int i = 0; i < nResults; i++)
138 vResult[i] += dx / 6.0 * (vFunctionValues[i] + 4.0 * vInter[i] + v[i]); // b-a/6*(f(a)+4f(a+b/2)+f(b))
139 }
140 else
141 {
142 // Calculate the integral at the current x position
143 if (vResult.size())
144 vResult.push_back(dx / 6.0 * (vFunctionValues[0] + 4.0 * vInter[0] + v[0]) + vResult.back());
145 else
146 vResult.push_back(dx / 6.0 * (vFunctionValues[0] + 4.0 * vInter[0] + v[0]));
147 }
148
149 // Set the last sample as the first one
150 vFunctionValues.assign(v, v+nResults);
151}
152
153
162static vector<mu::value_type> integrateSingleDimensionData(CommandLineParser& cmdParser)
163{
164 bool bReturnFunctionPoints = cmdParser.hasParam("points");
165 bool bCalcXvals = cmdParser.hasParam("xvals");
166
167 vector<mu::value_type> vResult;
168
169 MemoryManager& _data = NumeReKernel::getInstance()->getMemoryManager();
170
171 // Extract the integration interval
172 IntervalSet ivl = cmdParser.parseIntervals();
173
174 // Get table name and the corresponding indices
175 DataAccessParser accessParser = cmdParser.getExprAsDataObject();
176 accessParser.evalIndices();
177 Indices& _idx = accessParser.getIndices();
178 std::string sDatatable = accessParser.getDataObject();
179
180 // The indices are vectors
181 //
182 // If it is a single column or row, then we simply
183 // summarize its contents, otherwise we calculate the
184 // integral with the trapezoidal method
185 if ((_idx.row.size() == 1 || _idx.col.size() == 1) && !bReturnFunctionPoints)
186 vResult.push_back(_data.sum(sDatatable, _idx.row, _idx.col));
187 else if (_idx.row.size() == 1 || _idx.col.size() == 1)
188 {
189 // cumulative sum
190 vResult.push_back(_data.getElement(_idx.row[0], _idx.col[0], sDatatable));
191
192 if (_idx.row.size() == 1)
193 {
194 for (size_t i = 1; i < _idx.col.size(); i++)
195 vResult.push_back(vResult.back() + _data.getElement(_idx.row[0], _idx.col[i], sDatatable));
196 }
197 else
198 {
199 for (size_t i = 1; i < _idx.row.size(); i++)
200 vResult.push_back(vResult.back() + _data.getElement(_idx.row[i], _idx.col[0], sDatatable));
201 }
202 }
203 else
204 {
205 MemoryManager _cache;
206
207 // Copy the data
208 for (size_t i = 0; i < _idx.row.size(); i++)
209 {
210 _cache.writeToTable(i, 0, "table", _data.getElement(_idx.row[i], _idx.col[0], sDatatable));
211 _cache.writeToTable(i, 1, "table", _data.getElement(_idx.row[i], _idx.col[1], sDatatable));
212 }
213
214 // Sort the data
215 _cache.sortElements("sort -table c=1[2]");
216 mu::value_type dResult = 0.0;
217 long long int j = 1;
218
219 // Calculate the integral by jumping over NaNs
220 for (long long int i = 0; i < _cache.getLines("table", false) - 1; i++) //nan-suche
221 {
222 j = 1;
223
224 if (!_cache.isValidElement(i, 1, "table"))
225 continue;
226
227 while (!_cache.isValidElement(i + j, 1, "table") && i + j < _cache.getLines("table", false) - 1)
228 j++;
229
230 if (!_cache.isValidElement(i + j, 0, "table") || !_cache.isValidElement(i + j, 1, "table"))
231 break;
232
233 if (ivl.intervals.size() >= 1 && ivl.intervals[0].front().real() > _cache.getElement(i, 0, "table").real())
234 continue;
235
236 if (ivl.intervals.size() >= 1 && ivl.intervals[0].back().real() < _cache.getElement(i + j, 0, "table").real())
237 break;
238
239 // Calculate either the integral, its samples or the corresponding x values
240 if (!bReturnFunctionPoints && !bCalcXvals)
241 dResult += (_cache.getElement(i, 1, "table") + _cache.getElement(i + j, 1, "table")) / 2.0 * (_cache.getElement(i + j, 0, "table") - _cache.getElement(i, 0, "table"));
242 else if (bReturnFunctionPoints && !bCalcXvals)
243 {
244 if (vResult.size())
245 vResult.push_back((_cache.getElement(i, 1, "table") + _cache.getElement(i + j, 1, "table")) / 2.0 * (_cache.getElement(i + j, 0, "table") - _cache.getElement(i, 0, "table")) + vResult.back());
246 else
247 vResult.push_back((_cache.getElement(i, 1, "table") + _cache.getElement(i + j, 1, "table")) / 2.0 * (_cache.getElement(i + j, 0, "table") - _cache.getElement(i, 0, "table")));
248 }
249 else
250 vResult.push_back(_cache.getElement(i + j, 0, "table"));
251 }
252
253 // If the integral was calculated, then there is a
254 // single result, which hasn't been stored yet
255 if (!bReturnFunctionPoints && !bCalcXvals)
256 vResult.push_back(dResult);
257 }
258
259 // Return the result of the integral
260 return vResult;
261}
262
263
272bool integrate(CommandLineParser& cmdParser)
273{
274 Parser& _parser = NumeReKernel::getInstance()->getParser();
275 MemoryManager& _data = NumeReKernel::getInstance()->getMemoryManager();
276 const Settings& _option = NumeReKernel::getInstance()->getSettings();
277 string sIntegrationExpression = cmdParser.getExprAsMathExpression();
278 mu::value_type* v = 0;
279 int nResults = 0;
280 vector<mu::value_type> vResult; // Ausgabe-Wert
281 vector<mu::value_type> vFunctionValues; // Werte an der Stelle n und n+1
282 bool bLargeInterval = false; // Boolean: TRUE, wenn ueber ein grosses Intervall integriert werden soll
283 bool bReturnFunctionPoints = cmdParser.hasParam("points");
284 bool bCalcXvals = cmdParser.hasParam("xvals");
285 unsigned int nMethod = TRAPEZOIDAL; // 1 = trapezoidal, 2 = simpson
286 size_t nSamples = 1e3;
287
288 mu::value_type& x = _defVars.vValue[0][0];
289 double range;
290
291 // It's not possible to calculate the integral of a string expression
292 if (NumeReKernel::getInstance()->getStringParser().isStringExpression(sIntegrationExpression))
293 throw SyntaxError(SyntaxError::STRINGS_MAY_NOT_BE_EVALUATED_WITH_CMD, cmdParser.getCommandLine(), SyntaxError::invalid_position, "integrate");
294
295 // Ensure that the function is available
296 if (!sIntegrationExpression.length())
297 throw SyntaxError(SyntaxError::NO_INTEGRATION_FUNCTION, cmdParser.getCommandLine(), SyntaxError::invalid_position);
298
299 // If the integration function contains a data object,
300 // the calculation is done way different from the usual
301 // integration
302 if (_data.isTable(sIntegrationExpression)
303 && getMatchingParenthesis(sIntegrationExpression) != string::npos
304 && sIntegrationExpression.find_first_not_of(' ', getMatchingParenthesis(sIntegrationExpression) + 1) == string::npos) // xvals
305 {
306 cmdParser.setReturnValue(integrateSingleDimensionData(cmdParser));
307 return true;
308 }
309
310 // Evaluate the parameters
311 IntervalSet ivl = cmdParser.parseIntervals();
312
313 if (ivl.size())
314 {
315 if (ivl[0].min() == ivl[0].max())
316 throw SyntaxError(SyntaxError::INVALID_INTEGRATION_RANGES, cmdParser.getCommandLine(), SyntaxError::invalid_position);
317
318 if (isinf(ivl[0].min()) || isnan(ivl[0].min())
319 || isinf(ivl[0].max()) || isnan(ivl[0].max()))
320 {
321 cmdParser.setReturnValue("nan");
322 return false;
323 }
324
325 range = ivl[0].range();
326 }
327 else
328 throw SyntaxError(SyntaxError::NO_INTEGRATION_RANGES, cmdParser.getCommandLine(), SyntaxError::invalid_position);
329
330 std::vector<mu::value_type> vParVal = cmdParser.getParameterValueAsNumericalValue("precision");
331
332 if (vParVal.size())
333 nSamples = std::rint(range / vParVal.front().real());
334 else
335 {
336 vParVal = cmdParser.getParameterValueAsNumericalValue("p");
337
338 if (vParVal.size())
339 nSamples = std::rint(range / vParVal.front().real());
340 else
341 {
342 vParVal = cmdParser.getParameterValueAsNumericalValue("eps");
343
344 if (vParVal.size())
345 nSamples = std::rint(range / vParVal.front().real());
346 }
347 }
348
349 vParVal = cmdParser.getParameterValueAsNumericalValue("steps");
350
351 if (vParVal.size())
352 nSamples = std::abs(intCast(vParVal.front()));
353 else
354 {
355 vParVal = cmdParser.getParameterValueAsNumericalValue("s");
356
357 if (vParVal.size())
358 nSamples = std::abs(intCast(vParVal.front()));
359 }
360
361 std::string sParVal = cmdParser.getParameterValue("method");
362
363 if (!sParVal.length())
364 sParVal = cmdParser.getParameterValue("m");
365
366 if (sParVal == "trapezoidal")
367 nMethod = TRAPEZOIDAL;
368 else if (sParVal == "simpson")
369 nMethod = SIMPSON;
370
371 // Check, whether the expression actual depends
372 // upon the integration variable
373 _parser.SetExpr(sIntegrationExpression);
374
375 _parser.Eval(nResults);
376 vResult.resize(nResults, 0.0);
377
378 // Set the calculation variables to their starting values
379 vFunctionValues.resize(nResults, 0.0);
380
381 // Calculate the numerical integration
382 // If the precision is invalid (e.g. due to a very
383 // small interval, simply guess a reasonable interval
384 // here
385 if (!nSamples)
386 nSamples = 100;
387
388 // Calculate the x values, if desired
389 if (bCalcXvals)
390 {
391 for (size_t i = 0; i < nSamples; i++)
392 {
393 vResult.push_back(ivl[0](i, nSamples));
394 }
395
396 cmdParser.setReturnValue(vResult);
397 return true;
398 }
399
400 // Set the expression in the parser
401 _parser.SetExpr(sIntegrationExpression);
402
403 // Is it a large interval (then it will need more time)
404 if ((nMethod == TRAPEZOIDAL && nSamples >= 9.9e6) || (nMethod == SIMPSON && nSamples >= 1e4))
405 bLargeInterval = true;
406
407 // Do not allow a very high number of integration steps
408 if (nSamples > 1e10)
409 throw SyntaxError(SyntaxError::INVALID_INTEGRATION_PRECISION, cmdParser.getCommandLine(), SyntaxError::invalid_position);
410
411 // Set the integration variable to the lower border
412 x = ivl[0](0);
413
414 // Calculate the first sample(s)
415 v = _parser.Eval(nResults);
416 vFunctionValues.assign(v, v+nResults);
417
418 double dx = range / (nSamples-1);
419
420 // Perform the actual numerical integration
421 for (size_t i = 1; i < nSamples; i++)
422 {
423 // Calculate the samples first
424 if (nMethod == TRAPEZOIDAL)
425 integrationstep_trapezoidal(x, ivl[0](i, nSamples), dx, vResult, vFunctionValues, bReturnFunctionPoints);
426 else if (nMethod == SIMPSON)
427 integrationstep_simpson(x, ivl[0](2*i-1, 2*nSamples-1), ivl[0](2*i, 2*nSamples-1), dx, vResult, vFunctionValues, bReturnFunctionPoints);
428
429 // Print a status value, if needed
430 if (_option.systemPrints() && bLargeInterval)
431 {
432 if ((int)(i / (double)nSamples * 100) > (int)((i-1) / (double)nSamples * 100))
433 NumeReKernel::printPreFmt("\r|INTEGRATE> " + _lang.get("COMMON_EVALUATING") + " ... " + toString((int)(i / (double)nSamples * 100)) + " %");
434
435 if (NumeReKernel::GetAsyncCancelState())
436 {
437 NumeReKernel::printPreFmt("\r|INTEGRATE> " + _lang.get("COMMON_EVALUATING") + " ... " + _lang.get("COMMON_CANCEL") + ".\n");
438 throw SyntaxError(SyntaxError::PROCESS_ABORTED_BY_USER, "", SyntaxError::invalid_position);
439 }
440 }
441 }
442
443 // Display a success message
444 if (_option.systemPrints() && bLargeInterval)
445 NumeReKernel::printPreFmt("\r|INTEGRATE> " + _lang.get("COMMON_EVALUATING") + " ... 100 %: " + _lang.get("COMMON_SUCCESS") + "!\n");
446
447 cmdParser.setReturnValue(vResult);
448 return true;
449}
450
451
464static void refreshBoundaries(IntervalSet& ivl, const string& sIntegrationExpression)
465{
466 // Refresh the y boundaries, if necessary
467 ivl[1].refresh();
468
469 NumeReKernel::getInstance()->getParser().SetExpr(sIntegrationExpression);
470}
471
472
481bool integrate2d(CommandLineParser& cmdParser)
482{
483 Parser& _parser = NumeReKernel::getInstance()->getParser();
484 const Settings& _option = NumeReKernel::getInstance()->getSettings();
485 string sIntegrationExpression = cmdParser.getExprAsMathExpression(); // string fuer die zu integrierende Funktion
486 mu::value_type* v = 0;
487 int nResults = 0;
488 vector<mu::value_type> vResult[3]; // value_type-Array, wobei vResult[0] das eigentliche Ergebnis speichert
489 // und vResult[1] fuer die Zwischenergebnisse des inneren Integrals ist
490 vector<mu::value_type> fx_n[2][3]; // value_type-Array fuer die jeweiligen Stuetzstellen im inneren und aeusseren Integral
491 bool bRenewBoundaries = false; // bool, der speichert, ob die Integralgrenzen von x oder y abhaengen
492 bool bLargeArray = false; // bool, der TRUE fuer viele Datenpunkte ist;
493 unsigned int nMethod = TRAPEZOIDAL; // trapezoidal = 1, simpson = 2
494 size_t nSamples = 1e3;
495
496 mu::value_type& x = _defVars.vValue[0][0];
497 mu::value_type& y = _defVars.vValue[1][0];
498
499 double range;
500
501 // Strings may not be integrated
502 if (NumeReKernel::getInstance()->getStringParser().isStringExpression(sIntegrationExpression))
503 throw SyntaxError(SyntaxError::STRINGS_MAY_NOT_BE_EVALUATED_WITH_CMD, cmdParser.getCommandLine(), SyntaxError::invalid_position, "integrate");
504
505 // Ensure that the integration function is available
506 if (!sIntegrationExpression.length())
507 throw SyntaxError(SyntaxError::NO_INTEGRATION_FUNCTION, cmdParser.getCommandLine(), SyntaxError::invalid_position);
508
509 // Evaluate the parameters
510 IntervalSet ivl = cmdParser.parseIntervals();
511
512 if (ivl.intervals.size() >= 2)
513 {
514 if (ivl[0].front() == ivl[0].back())
515 throw SyntaxError(SyntaxError::INVALID_INTEGRATION_RANGES, cmdParser.getCommandLine(), SyntaxError::invalid_position);
516
517 if (isinf(ivl[0].min()) || isnan(ivl[0].min())
518 || isinf(ivl[0].max()) || isnan(ivl[0].max())
519 || isinf(ivl[1].min()) || isnan(ivl[1].min())
520 || isinf(ivl[1].max()) || isnan(ivl[1].max()))
521 {
522 cmdParser.setReturnValue("nan");
523 return false;
524 }
525
526 range = std::min(ivl[0].range(), ivl[1].range());
527 }
528 else
529 throw SyntaxError(SyntaxError::NO_INTEGRATION_RANGES, cmdParser.getCommandLine(), SyntaxError::invalid_position);
530
531 std::vector<mu::value_type> vParVal = cmdParser.getParameterValueAsNumericalValue("precision");
532
533 if (vParVal.size())
534 nSamples = std::rint(range / vParVal.front().real());
535 else
536 {
537 vParVal = cmdParser.getParameterValueAsNumericalValue("p");
538
539 if (vParVal.size())
540 nSamples = std::rint(range / vParVal.front().real());
541 else
542 {
543 vParVal = cmdParser.getParameterValueAsNumericalValue("eps");
544
545 if (vParVal.size())
546 nSamples = std::rint(range / vParVal.front().real());
547 }
548 }
549
550 vParVal = cmdParser.getParameterValueAsNumericalValue("steps");
551
552 if (vParVal.size())
553 nSamples = std::abs(intCast(vParVal.front()));
554 else
555 {
556 vParVal = cmdParser.getParameterValueAsNumericalValue("s");
557
558 if (vParVal.size())
559 nSamples = std::abs(intCast(vParVal.front()));
560 }
561
562 std::string sParVal = cmdParser.getParameterValue("method");
563
564 if (!sParVal.length())
565 sParVal = cmdParser.getParameterValue("m");
566
567 if (sParVal == "trapezoidal")
568 nMethod = TRAPEZOIDAL;
569 else if (sParVal == "simpson")
570 nMethod = SIMPSON;
571
572 // Check, whether the expression depends upon one or both
573 // integration variables
574 _parser.SetExpr(sIntegrationExpression);
575
576 // Prepare the memory for integration
577 _parser.Eval(nResults);
578
579 for (int i = 0; i < 3; i++)
580 {
581 vResult[i].resize(nResults, 0.0);
582 fx_n[0][i].resize(nResults, 0.0);
583 fx_n[1][i].resize(nResults, 0.0);
584 }
585
586 if (ivl[1].contains(_defVars.sName[0]))
587 bRenewBoundaries = true; // Ja? Setzen wir den bool entsprechend
588
589 // Ensure that the precision is reasonble
590 // If the precision is invalid, we guess a reasonable value here
591 if (!nSamples)
592 {
593 // We use the smallest intervall and split it into
594 // 100 parts
595 nSamples = 100;
596 }
597
598 // Is it a very slow integration?
599 if ((nMethod == TRAPEZOIDAL && nSamples*nSamples >= 1e8) || (nMethod == SIMPSON && nSamples*nSamples >= 1e6))
600 bLargeArray = true;
601
602 // Avoid calculation with too many steps
603 if (nSamples*nSamples > 1e10)
604 throw SyntaxError(SyntaxError::INVALID_INTEGRATION_PRECISION, cmdParser.getCommandLine(), SyntaxError::invalid_position);
605
606 // --> Kleine Info an den Benutzer, dass der Code arbeitet <--
607 if (_option.systemPrints() && bLargeArray)
608 NumeReKernel::printPreFmt("\r|INTEGRATE> " + _lang.get("COMMON_EVALUATING") + " ... 0 %");
609
610 // --> Setzen wir "x" und "y" auf ihre Startwerte <--
611 x = ivl[0](0); // x = x_0
612
613 // y might depend on the starting value
614 if (bRenewBoundaries)
615 refreshBoundaries(ivl, sIntegrationExpression);
616
617 y = ivl[1](0); // y = y_0
618
619 double dx = ivl[0].range() / (nSamples-1);
620 double dy = ivl[1].range() / (nSamples-1);
621
622 // --> Werte mit den Startwerten die erste Stuetzstelle fuer die y-Integration aus <--
623 v = _parser.Eval(nResults);
624 fx_n[1][0].assign(v, v+nResults);
625
626 /* --> Berechne das erste y-Integral fuer die erste Stuetzstelle fuer x
627 * Die Schleife laeuft so lange wie y < y_1 <--
628 */
629 for (size_t j = 1; j < nSamples; j++)
630 {
631 if (nMethod == TRAPEZOIDAL)
632 integrationstep_trapezoidal(y, ivl[1](j, nSamples), dy, vResult[1], fx_n[1][0], false);
633 else if (nMethod == SIMPSON)
634 integrationstep_simpson(y, ivl[1](2*j-1, 2*nSamples-1), ivl[1](2*j, 2*nSamples-1), dy, vResult[1], fx_n[1][0], false);
635 }
636
637 fx_n[0][0] = vResult[1];
638
639 /* --> Das eigentliche, numerische Integral. Es handelt sich um nichts weiter als viele
640 * while()-Schleifendurchlaeufe.
641 * Die aeussere Schleife laeuft so lange x < x_1 ist. <--
642 */
643 for (size_t i = 1; i < nSamples; i++)
644 {
645 if (nMethod == TRAPEZOIDAL)
646 {
647 x = ivl[0](i, nSamples);
648
649 // Refresh the y boundaries, if necessary
650 if (bRenewBoundaries)
651 {
652 refreshBoundaries(ivl, sIntegrationExpression);
653 dy = ivl[1].range() / (nSamples-1);
654 }
655
656 // --> Setzen wir "y" auf den Wert, der von der unteren y-Grenze vorgegeben wird <--
657 y = ivl[1](0);
658 // --> Werten wir sofort die erste y-Stuetzstelle aus <--
659 v = _parser.Eval(nResults);
660 fx_n[1][0].assign(v, v+nResults);
661 vResult[1].assign(nResults, 0.0);
662
663 for (size_t j = 1; j < nSamples; j++)
664 integrationstep_trapezoidal(y, ivl[1](j, nSamples), dy, vResult[1], fx_n[1][0], false);
665
666 // --> Weise das Ergebnis der y-Integration an die zweite Stuetzstelle der x-Integration zu <--
667 for (int i = 0; i < nResults; i++)
668 {
669 if (isnan(std::abs(vResult[1][i])))
670 vResult[1][i] = 0.0;
671
672 vResult[0][i] += dx * (fx_n[0][0][i] + vResult[1][i]) * 0.5; // Berechne das Trapez zu x
673 fx_n[0][0][i] = vResult[1][i]; // Weise den Wert der zweiten Stuetzstelle an die erste Stuetzstelle zu
674 }
675 }
676 else if (nMethod == SIMPSON)
677 {
678 for (size_t n = 1; n <= 2; n++)
679 {
680 x = ivl[0](2*i+n-2, 2*nSamples-1);
681
682 // Refresh the y boundaries, if necessary
683 if (bRenewBoundaries)
684 {
685 refreshBoundaries(ivl, sIntegrationExpression);
686 dy = ivl[1].range() / (nSamples-1);
687 }
688
689 // Set y to the first position
690 y = ivl[1](0);
691
692 // Calculate the first position
693 if (n == 1)
694 {
695 v = _parser.Eval(nResults);
696 fx_n[1][0].assign(v, v+nResults);
697 }
698
699 vResult[n].assign(nResults, 0.0);
700
701 for (size_t j = 1; j < nSamples; j++)
702 integrationstep_simpson(y, ivl[1](2*j-1, 2*nSamples-1), ivl[1](2*j, 2*nSamples-1), dy, vResult[n], fx_n[1][0], false);
703
704 // --> Weise das Ergebnis der y-Integration an die zweite Stuetzstelle der x-Integration zu <--
705 for (int i = 0; i < nResults; i++)
706 {
707 if (isnan(std::abs(vResult[n][i])))
708 vResult[n][i] = 0.0;
709 }
710 }
711
712 for (int i = 0; i < nResults; i++)
713 {
714 vResult[0][i] += dx / 6.0 * (fx_n[0][0][i] + 4.0 * vResult[1][i] + vResult[2][i]); // Berechne das Trapez zu x
715 fx_n[0][0][i] = vResult[2][i]; // Weise den Wert der zweiten Stuetzstelle an die erste Stuetzstelle zu
716 }
717 }
718
719 // Show some progress
720 if (_option.systemPrints())
721 {
722 if (bLargeArray)
723 {
724 if ((int)(i / (double)nSamples * 100) > (int)((i-1) / (double)nSamples * 100))
725 NumeReKernel::printPreFmt("\r|INTEGRATE> " + _lang.get("COMMON_EVALUATING") + " ... " + toString((int)(i / (double)nSamples * 100)) + " %");
726 }
727
728 if (NumeReKernel::GetAsyncCancelState())//GetAsyncKeyState(VK_ESCAPE))
729 {
730 NumeReKernel::printPreFmt("\r|INTEGRATE> " + _lang.get("COMMON_EVALUATING") + " ... " + _lang.get("COMMON_CANCEL") + "!\n");
731 throw SyntaxError(SyntaxError::PROCESS_ABORTED_BY_USER, "", SyntaxError::invalid_position);
732 }
733 }
734 }
735
736 // Show a success message
737 if (_option.systemPrints() && bLargeArray)
738 NumeReKernel::printPreFmt("\r|INTEGRATE> " + _lang.get("COMMON_EVALUATING") + " ... 100 %: " + _lang.get("COMMON_SUCCESS") + "!\n");
739
740 // --> Fertig! Zurueck zur aufrufenden Funkton! <--
741 cmdParser.setReturnValue(vResult[0]);
742 return true;
743}
744
745
754bool differentiate(CommandLineParser& cmdParser)
755{
756 Parser& _parser = NumeReKernel::getInstance()->getParser();
757 MemoryManager& _data = NumeReKernel::getInstance()->getMemoryManager();
758 const Settings& _option = NumeReKernel::getInstance()->getSettings();
759 FunctionDefinitionManager& _functions = NumeReKernel::getInstance()->getDefinitions();
760 string sExpr = cmdParser.getExprAsMathExpression();
761 string sVar = "";
762 string sPos = "";
763 double dEps = 0.0;
764 mu::value_type dPos = 0.0;
765 mu::value_type* dVar = 0;
766 mu::value_type* v = 0;
767 int nResults = 0;
768 int nSamples = 100;
769 size_t order = 1;
770 std::vector<mu::value_type> vInterval;
771 std::vector<mu::value_type> vResult;
772 std::vector<mu::value_type> paramVal;
773
774 // Strings cannot be differentiated
775 if (NumeReKernel::getInstance()->getStringParser().isStringExpression(sExpr))
776 throw SyntaxError(SyntaxError::STRINGS_MAY_NOT_BE_EVALUATED_WITH_CMD, cmdParser.getCommandLine(), SyntaxError::invalid_position, "diff");
777
778 // Get the order of the differntiation
779 paramVal = cmdParser.getParameterValueAsNumericalValue("order");
780
781 if (paramVal.size())
782 {
783 order = intCast(paramVal.front());
784 order = std::min(order, 3u);
785 }
786
787 // Numerical expressions and data sets are handled differently
788 if (!_data.containsTablesOrClusters(sExpr) && cmdParser.getParameterList().length())
789 {
790 // Is the "eps" parameter available?
791 paramVal = cmdParser.getParameterValueAsNumericalValue("eps");
792
793 if (paramVal.size())
794 dEps = fabs(paramVal.front());
795
796 paramVal = cmdParser.getParameterValueAsNumericalValue("samples");
797
798 if (paramVal.size())
799 nSamples = intCast(fabs(paramVal.front()));
800
801 sVar = cmdParser.getParameterList();
802
803 if (!_functions.call(sVar))
804 throw SyntaxError(SyntaxError::FUNCTION_ERROR, cmdParser.getCommandLine(), sVar, sVar);
805
806 // Is a variable interval defined?
807 if (sVar.find('=') != string::npos ||
808 (sVar.find('[') != string::npos
809 && sVar.find(']', sVar.find('[')) != string::npos
810 && sVar.find(':', sVar.find('[')) != string::npos))
811 {
812 // Remove possible parameter list initializers
813 if (sVar.substr(0, 2) == "--")
814 sVar = sVar.substr(2);
815 else if (sVar.substr(0, 4) == "-set")
816 sVar = sVar.substr(4);
817
818 // Extract variable intervals or locations
819 if (sVar.find('[') != string::npos
820 && sVar.find(']', sVar.find('[')) != string::npos
821 && sVar.find(':', sVar.find('[')) != string::npos)
822 {
823 sPos = sVar.substr(sVar.find('[') + 1, getMatchingParenthesis(sVar.substr(sVar.find('['))) - 1);
824 sVar = "x";
825 StripSpaces(sPos);
826
827 if (sPos == ":")
828 sPos = "-10:10";
829 }
830 else
831 {
832 int nPos = sVar.find('=');
833 sPos = sVar.substr(nPos + 1, sVar.find(' ', nPos) - nPos - 1);
834 sVar = " " + sVar.substr(0, nPos);
835 sVar = sVar.substr(sVar.rfind(' '));
836 StripSpaces(sVar);
837 }
838
839 // Evaluate the position/range expression
840 if (isNotEmptyExpression(sPos))
841 {
842 if (_data.containsTablesOrClusters(sPos))
843 getDataElements(sPos, _parser, _data, _option);
844
845 if (sPos.find(':') != string::npos)
846 sPos.replace(sPos.find(':'), 1, ",");
847 _parser.SetExpr(sPos);
848 v = _parser.Eval(nResults);
849
850 if (isinf(std::abs(v[0])) || isnan(std::abs(v[0])))
851 {
852 cmdParser.setReturnValue("nan");
853 return true;
854 }
855
856 for (int i = 0; i < nResults; i++)
857 vInterval.push_back(v[i]);
858 }
859
860 // Set the expression for differentiation
861 // and evaluate it
862 _parser.SetExpr(sExpr);
863 _parser.Eval(nResults);
864
865 // Get the address of the variable
866 dVar = getPointerToVariable(sVar, _parser);
867 }
868
869 // Ensure that the address could be found
870 if (!dVar)
871 throw SyntaxError(SyntaxError::NO_DIFF_VAR, cmdParser.getCommandLine(), SyntaxError::invalid_position);
872
873 // Store the expression
874 string sCompl_Expr = sExpr;
875
876 // As long as the expression has a length
877 while (sCompl_Expr.length())
878 {
879 // Get the next subexpression and
880 // set it in the parser
881 sExpr = getNextArgument(sCompl_Expr, true);
882 _parser.SetExpr(sExpr);
883
884 // Evaluate the differential at the desired
885 // locations
886 if (vInterval.size() == 1 || vInterval.size() > 2)
887 {
888 // single point or a vector
889 for (unsigned int i = 0; i < vInterval.size(); i++)
890 {
891 dPos = vInterval[i];
892 vResult.push_back(_parser.Diff(dVar, dPos, dEps, order));
893 }
894 }
895 else
896 {
897 // a range -> use the samples
898 for (int i = 0; i < nSamples; i++)
899 {
900 dPos = vInterval[0] + (vInterval[1] - vInterval[0]) / (double)(nSamples - 1) * (double)i;
901 vResult.push_back(_parser.Diff(dVar, dPos, dEps, order));
902 }
903 }
904 }
905 }
906 else if (_data.containsTablesOrClusters(sExpr))
907 {
908 // This is a data set
909 //
910 // Get the indices first
911 DataAccessParser accessParser = cmdParser.getExprAsDataObject();
912 Indices& _idx = accessParser.getIndices();
913 std::string sTableName = accessParser.getDataObject();
914 size_t nFilterSize = 5;
915 paramVal = cmdParser.getParameterValueAsNumericalValue("points");
916
917 if (paramVal.size())
918 {
919 nFilterSize = intCast(paramVal.front());
920
921 if (!(nFilterSize % 2))
922 nFilterSize++;
923 }
924
925 NumeRe::SavitzkyGolayDiffFilter diff(nFilterSize, order);
926
927 // Validate the indices
928 if (!isValidIndexSet(_idx))
929 throw SyntaxError(SyntaxError::INVALID_INDEX, cmdParser.getCommandLine(), SyntaxError::invalid_position, _idx.row.to_string() + ", " + _idx.col.to_string());
930
931 if (_idx.row.isOpenEnd())
932 _idx.row.setRange(0, _data.getLines(sTableName, false)-1);
933
934 if (_idx.col.isOpenEnd())
935 _idx.col.setRange(0, _idx.col.front()+1);
936
937 // If shorter than filter's size return an invalid
938 // value
939 if (_idx.row.size() < nFilterSize)
940 {
941 vResult.push_back(NAN);
942 cmdParser.setReturnValue(vResult);
943 return true;
944 }
945
946 // Copy the data contents, sort the values
947 // and calculate the derivative
948
949 // Vectors as indices
950 //
951 // Depending on the number of selected columns, we either
952 // have to sort the data or we assume that the difference
953 // between two values is 1
954 if (_idx.col.size() == 1)
955 {
956 vResult.resize(_idx.row.size());
957
958 // No sorting, difference is 1
959 //
960 // Jump over NaNs and get the difference of the neighbouring
961 // values, which is identical to the derivative in this case
962 for (size_t i = nFilterSize/2; i < _idx.row.size() - nFilterSize/2; i++)
963 {
964 for (int j = 0; j < (int)nFilterSize; j++)
965 {
966 if (_data.isValidElement(_idx.row[i + j - nFilterSize/2], _idx.col.front(), sTableName))
967 vResult[i] += diff.apply(j, 0, _data.getElement(_idx.row[i + j - nFilterSize/2], _idx.col.front(), sTableName));
968 }
969 }
970
971 // Repeat the first and last values
972 for (size_t i = 0; i < nFilterSize/2; i++)
973 {
974 vResult[i] = vResult[nFilterSize/2];
975 vResult[vResult.size()-1-i] = vResult[vResult.size()-1-nFilterSize/2];
976 }
977 }
978 else
979 {
980 // We have to sort, because the difference is passed
981 // explicitly
982 //
983 // Copy the data first and sort afterwards
984 MemoryManager _cache;
985
986 for (size_t i = 0; i < _idx.row.size(); i++)
987 {
988 _cache.writeToTable(i, 0, "table", _data.getElement(_idx.row[i], _idx.col[0], sTableName));
989 _cache.writeToTable(i, 1, "table", _data.getElement(_idx.row[i], _idx.col[1], sTableName));
990 }
991
992 _cache.sortElements("sort -table c=1[2]");
993
994 // Shall the x values be calculated?
995 if (cmdParser.hasParam("xvals"))
996 {
997 NumeReKernel::getInstance()->issueWarning(_lang.get("COMMON_SYNTAX_DEPRECATED", cmdParser.getCommandLine()));
998
999 // The x values are approximated to be in the
1000 // middle of the two samplex
1001 for (int i = 0; i < _cache.getLines("table", false) - 1; i++)
1002 {
1003 if (_cache.isValidElement(i, 0, "table")
1004 && _cache.isValidElement(i + 1, 0, "table")
1005 && _cache.isValidElement(i, 1, "table")
1006 && _cache.isValidElement(i + 1, 1, "table"))
1007 vResult.push_back((_cache.getElement(i + 1, 0, "table") + _cache.getElement(i, 0, "table")) / 2.0);
1008 else
1009 vResult.push_back(NAN);
1010 }
1011 }
1012 else
1013 {
1014 vResult.resize(_cache.getLines("table", false));
1015
1016 // We calculate the derivative of the data
1017 // by approximating it linearily
1018 for (int i = nFilterSize/2; i < _cache.getLines("table", false) - (int)nFilterSize/2; i++)
1019 {
1020 std::pair<mu::value_type, size_t> avgDiff(0.0, 0);
1021
1022 for (int j = 0; j < (int)nFilterSize; j++)
1023 {
1024 if (_cache.isValidElement(i + j - nFilterSize/2, 0, "table")
1025 && _cache.isValidElement(i + j - nFilterSize/2, 1, "table"))
1026 {
1027 vResult[i] += diff.apply(j, 0, _cache.getElement(i + j - nFilterSize/2, 1, "table"));
1028
1029 // Calculate the average difference
1030 if (_cache.isValidElement(i + j - nFilterSize/2 - 1, 0, "table"))
1031 {
1032 avgDiff.first += _cache.getElement(i + j - nFilterSize/2, 0, "table")
1033 - _cache.getElement(i + j - nFilterSize/2 - 1, 0, "table");
1034 avgDiff.second++;
1035 }
1036 }
1037 }
1038
1039 if (!avgDiff.second)
1040 vResult[i] = NAN;
1041 else
1042 vResult[i] /= intPower(avgDiff.first/(double)avgDiff.second, order);
1043 }
1044
1045 // Repeat the first and last values
1046 for (size_t i = 0; i < nFilterSize/2; i++)
1047 {
1048 vResult[i] = vResult[nFilterSize/2];
1049 vResult[vResult.size()-1-i] = vResult[vResult.size()-1-nFilterSize/2];
1050 }
1051 }
1052 }
1053 }
1054 else
1055 {
1056 // Ensure that a parameter list is available
1057 throw SyntaxError(SyntaxError::NO_DIFF_OPTIONS, cmdParser.getCommandLine(), SyntaxError::invalid_position);
1058 }
1059
1060 cmdParser.setReturnValue(vResult);
1061 return true;
1062}
1063
1064
1075static string getIntervalForSearchFunctions(const string& sParams, string& sVar)
1076{
1077 string sInterval = "";
1078
1079 // Extract the interval definition from the
1080 // parameter string
1081 if (sParams.find('=') != string::npos)
1082 {
1083 int nPos = sParams.find('=');
1084 sInterval = getArgAtPos(sParams, nPos + 1);
1085
1086 if (sInterval.front() == '[' && sInterval.back() == ']')
1087 {
1088 sInterval.pop_back();
1089 sInterval.erase(0, 1);
1090 }
1091
1092 sVar = " " + sParams.substr(0, nPos);
1093 sVar = sVar.substr(sVar.rfind(' '));
1094 StripSpaces(sVar);
1095 }
1096 else
1097 {
1098 sVar = "x";
1099 sInterval = sParams.substr(sParams.find('[') + 1, getMatchingParenthesis(sParams.substr(sParams.find('['))) - 1);
1100 StripSpaces(sInterval);
1101
1102 if (sInterval == ":")
1103 sInterval = "-10:10";
1104 }
1105
1106 return sInterval;
1107}
1108
1109
1123static bool findExtremaInMultiResult(CommandLineParser& cmdParser, string& sExpr, string& sInterval, int nOrder, int nMode)
1124{
1125 Parser& _parser = NumeReKernel::getInstance()->getParser();
1126 _parser.SetExpr(sExpr);
1127 int nResults;
1128 mu::value_type* v = _parser.Eval(nResults);
1129 vector<mu::value_type> vResults;
1130 int nResults_x = 0;
1131 MemoryManager _cache;
1132
1133 // Store the results in the second column of a table
1134 for (int i = 0; i < nResults; i++)
1135 _cache.writeToTable(i, 1, "table", v[i]);
1136
1137 _parser.SetExpr(sInterval);
1138 v = _parser.Eval(nResults_x);
1139
1140 // Write the results for the x interval in the first column
1141 if (nResults_x > 1)
1142 {
1143 for (int i = 0; i < nResults; i++)
1144 {
1145 if (i >= nResults_x)
1146 _cache.writeToTable(i, 0, "table", 0.0);
1147 else
1148 _cache.writeToTable(i, 0, "table", v[i]);
1149 }
1150 }
1151 else
1152 return false;
1153
1154 std::string sSortingExpr = "sort -table cols=1[2]";
1155 _cache.sortElements(sSortingExpr);
1156
1157 double dMedian = 0.0, dExtremum = 0.0;
1158 double* data = new double[nOrder];
1159 double nDir = 0;
1160 int nanShift = 0;
1161
1162 if (nOrder >= nResults / 3)
1163 nOrder = nResults / 3;
1164
1165 // Ensure that the number of used points is reasonable
1166 if (nOrder < 3)
1167 {
1168 vResults.push_back(NAN);
1169 return false;
1170 }
1171
1172 // Find the first median and use it as starting point
1173 // for identifying the next extremum
1174 for (int i = 0; i + nanShift < _cache.getLines("table", true); i++)
1175 {
1176 if (i == nOrder)
1177 break;
1178
1179 while (isnan(std::abs(_cache.getElement(i + nanShift, 1, "table"))) && i + nanShift < _cache.getLines("table", true) - 1)
1180 nanShift++;
1181
1182 data[i] = _cache.getElement(i + nanShift, 1, "table").real();
1183 }
1184
1185 // Sort the data and find the median
1186 gsl_sort(data, 1, nOrder);
1187 dExtremum = gsl_stats_median_from_sorted_data(data, 1, nOrder);
1188
1189 // Go through the data points using sliding median to find the local
1190 // extrema in the data set
1191 for (int i = nOrder; i + nanShift < _cache.getLines("table", false) - nOrder; i++)
1192 {
1193 int currNanShift = 0;
1194 dMedian = 0.0;
1195
1196 for (int j = i; j < i + nOrder; j++)
1197 {
1198 while (isnan(std::abs(_cache.getElement(j + nanShift + currNanShift, 1, "table"))) && j + nanShift + currNanShift < _cache.getLines("table", true) - 1)
1199 currNanShift++;
1200
1201 data[j - i] = _cache.getElement(j + nanShift + currNanShift, 1, "table").real();
1202 }
1203
1204 gsl_sort(data, 1, nOrder);
1205 dMedian = gsl_stats_median_from_sorted_data(data, 1, nOrder);
1206
1207 if (!nDir)
1208 {
1209 if (dMedian > dExtremum)
1210 nDir = 1;
1211 else if (dMedian < dExtremum)
1212 nDir = -1;
1213
1214 dExtremum = dMedian;
1215 }
1216 else
1217 {
1218 if (nDir*dMedian < nDir*dExtremum)
1219 {
1220 if (!nMode || nMode == nDir)
1221 {
1222 int nExtremum = i + nanShift;
1223 double dExtremum = _cache.getElement(i + nanShift, 1, "table").real();
1224
1225 for (long long int k = i + nanShift; k >= 0; k--)
1226 {
1227 if (k == i - nOrder)
1228 break;
1229
1230 if (nDir*_cache.getElement(k, 1, "table").real() > nDir*dExtremum)
1231 {
1232 nExtremum = k;
1233 dExtremum = _cache.getElement(k, 1, "table").real();
1234 }
1235 }
1236
1237 vResults.push_back(_cache.getElement(nExtremum, 0, "table"));
1238 i = nExtremum + nOrder;
1239 }
1240
1241 nDir = 0;
1242 }
1243
1244 dExtremum = dMedian;
1245 }
1246
1247 nanShift += currNanShift;
1248 }
1249
1250 if (!vResults.size())
1251 vResults.push_back(NAN);
1252
1253 delete[] data;
1254 cmdParser.setReturnValue(vResults);
1255 return true;
1256}
1257
1258
1272static double calculateMedian(value_type* v, int nResults, int start, int nOrder, int nanShiftStart, int& nNewNanShift)
1273{
1274 vector<double> data;
1275
1276 for (int i = start; i < start + nOrder; i++)
1277 {
1278 while (isnan(v[i + nanShiftStart + nNewNanShift].real()) && i + nanShiftStart + nNewNanShift < nResults - 1)
1279 nNewNanShift++;
1280
1281 if (i + nanShiftStart + nNewNanShift >= nResults)
1282 break;
1283
1284 data.push_back(v[i + nanShiftStart + nNewNanShift].real());
1285 }
1286
1287 gsl_sort(&data[0], 1, data.size());
1288 return gsl_stats_median_from_sorted_data(&data[0], 1, data.size());
1289}
1290
1291
1303static bool findExtremaInData(CommandLineParser& cmdParser, string& sExpr, int nOrder, int nMode)
1304{
1305 mu::value_type* v;
1306 int nResults = 0;
1307 Parser& _parser = NumeReKernel::getInstance()->getParser();
1308 _parser.SetExpr(sExpr);
1309 v = _parser.Eval(nResults);
1310
1311 if (nResults > 1)
1312 {
1313 if (nOrder >= nResults / 3)
1314 nOrder = nResults / 3;
1315
1316 double dMedian = 0.0, dExtremum = 0.0;
1317 double nDir = 0;
1318 int nanShift = 0;
1319 vector<mu::value_type> vResults;
1320
1321 if (nOrder < 3)
1322 {
1323 vResults.push_back(NAN);
1324 return false;
1325 }
1326
1327 dExtremum = calculateMedian(v, nResults, 0, nOrder, 0, nanShift);
1328
1329 for (int i = nOrder; i + nanShift < nResults - nOrder; i++)
1330 {
1331 int currNanShift = 0;
1332 dMedian = calculateMedian(v, nResults, i, nOrder, nanShift, currNanShift);
1333
1334 if (!nDir)
1335 {
1336 if (dMedian > dExtremum)
1337 nDir = 1;
1338 else if (dMedian < dExtremum)
1339 nDir = -1;
1340
1341 dExtremum = dMedian;
1342 }
1343 else
1344 {
1345 if (nDir*dMedian < nDir*dExtremum)
1346 {
1347 if (!nMode || nMode == nDir)
1348 {
1349 int nExtremum = i + nanShift;
1350 double dLocalExtremum = v[i + nanShift].real();
1351
1352 for (long long int k = i + nanShift; k >= 0; k--)
1353 {
1354 if (k == i + nanShift - nOrder)
1355 break;
1356
1357 if (nDir*v[k].real() > nDir*dLocalExtremum)
1358 {
1359 nExtremum = k;
1360 dLocalExtremum = v[k].real();
1361 }
1362 }
1363
1364 vResults.push_back(nExtremum + 1);
1365 i = nExtremum + nOrder;
1366 nanShift = 0;
1367 }
1368
1369 nDir = 0;
1370 }
1371
1372 dExtremum = dMedian;
1373 }
1374
1375 nanShift += currNanShift;
1376 }
1377
1378 if (!vResults.size())
1379 vResults.push_back(NAN);
1380
1381 cmdParser.setReturnValue(vResults);
1382 return true;
1383 }
1384 else
1385 throw SyntaxError(SyntaxError::NO_EXTREMA_VAR, cmdParser.getCommandLine(), SyntaxError::invalid_position);
1386}
1387
1388
1398bool findExtrema(CommandLineParser& cmdParser)
1399{
1400 NumeReKernel* instance = NumeReKernel::getInstance();
1401 MemoryManager& _data = instance->getMemoryManager();
1402 Parser& _parser = instance->getParser();
1403
1404 unsigned int nSamples = 21;
1405 int nOrder = 5;
1406 mu::value_type dVal[2];
1407 mu::value_type dBoundaries[2] = {0.0, 0.0};
1408 int nMode = 0;
1409 mu::value_type* dVar = 0;
1410 string sExpr = "";
1411 string sParams = "";
1412 string sInterval = "";
1413 string sVar = "";
1414
1415 // We cannot search extrema in strings
1416 if (instance->getStringParser().isStringExpression(cmdParser.getExpr()))
1417 throw SyntaxError(SyntaxError::STRINGS_MAY_NOT_BE_EVALUATED_WITH_CMD, cmdParser.getCommandLine(), SyntaxError::invalid_position, "extrema");
1418
1419 if (!_data.containsTablesOrClusters(cmdParser.getExpr()) && !cmdParser.getParameterList().length())
1420 throw SyntaxError(SyntaxError::NO_EXTREMA_OPTIONS, cmdParser.getCommandLine(), SyntaxError::invalid_position);
1421
1422 // Isolate the expression
1423 StripSpaces(sExpr);
1424 sExpr = sExpr.substr(findCommand(sExpr).sString.length());
1425
1426 // Ensure that the expression is not empty
1427 // and that the custom functions don't throw
1428 // any errors
1429 sExpr = cmdParser.getExpr();
1430 sParams = cmdParser.getParameterList();
1431
1432 if (!isNotEmptyExpression(sExpr) || !instance->getDefinitions().call(sExpr))
1433 return false;
1434
1435 if (!instance->getDefinitions().call(sParams))
1436 return false;
1437
1438 StripSpaces(sParams);
1439
1440 // If the expression or the parameter list contains
1441 // data elements, get their values here
1442 if (_data.containsTablesOrClusters(sExpr))
1443 getDataElements(sExpr, _parser, _data, instance->getSettings(), false);
1444
1445 if (_data.containsTablesOrClusters(sParams))
1446 getDataElements(sParams, _parser, _data, instance->getSettings(), false);
1447
1448 // Evaluate the parameters
1449 if (findParameter(sParams, "min"))
1450 nMode = -1;
1451
1452 if (findParameter(sParams, "max"))
1453 nMode = 1;
1454
1455 if (findParameter(sParams, "samples", '='))
1456 {
1457 _parser.SetExpr(getArgAtPos(sParams, findParameter(sParams, "samples", '=') + 7));
1458 nSamples = intCast(_parser.Eval());
1459
1460 if (nSamples < 21)
1461 nSamples = 21;
1462
1463 sParams.erase(findParameter(sParams, "samples", '=') - 1, 8);
1464 }
1465
1466 if (findParameter(sParams, "points", '='))
1467 {
1468 _parser.SetExpr(getArgAtPos(sParams, findParameter(sParams, "points", '=') + 6));
1469 nOrder = intCast(_parser.Eval());
1470
1471 if (nOrder <= 3)
1472 nOrder = 3;
1473
1474 sParams.erase(findParameter(sParams, "points", '=') - 1, 7);
1475 }
1476
1477 // Extract the interval
1478 if (sParams.find('=') != string::npos
1479 || (sParams.find('[') != string::npos
1480 && sParams.find(']', sParams.find('['))
1481 && sParams.find(':', sParams.find('['))))
1482 {
1483 if (sParams.substr(0, 2) == "--")
1484 sParams = sParams.substr(2);
1485 else if (sParams.substr(0, 4) == "-set")
1486 sParams = sParams.substr(4);
1487
1488 int nResults = 0;
1489
1490 sInterval = getIntervalForSearchFunctions(sParams, sVar);
1491
1492 _parser.SetExpr(sExpr);
1493 _parser.Eval(nResults);
1494
1495 if (nResults > 1)
1496 return findExtremaInMultiResult(cmdParser, sExpr, sInterval, nOrder, nMode);
1497 else
1498 {
1499 if (findVariableInExpression(sExpr, sVar) == std::string::npos)
1500 {
1501 cmdParser.setReturnValue("nan");
1502 return true;
1503 }
1504
1505 dVar = getPointerToVariable(sVar, _parser);
1506
1507 if (!dVar)
1508 throw SyntaxError(SyntaxError::EXTREMA_VAR_NOT_FOUND, cmdParser.getCommandLine(), sVar, sVar);
1509
1510 if (sInterval.find(':') == string::npos || sInterval.length() < 3)
1511 return false;
1512
1513 auto indices = getAllIndices(sInterval);
1514
1515 for (size_t i = 0; i < 2; i++)
1516 {
1517 if (isNotEmptyExpression(indices[i]))
1518 {
1519 _parser.SetExpr(indices[i]);
1520 dBoundaries[i] = _parser.Eval();
1521
1522 if (isinf(std::abs(dBoundaries[i])) || isnan(std::abs(dBoundaries[i])))
1523 {
1524 cmdParser.setReturnValue("nan");
1525 return false;
1526 }
1527 }
1528 else
1529 return false;
1530 }
1531
1532 if (std::abs(dBoundaries[1]) < std::abs(dBoundaries[0]))
1533 {
1534 mu::value_type Temp = dBoundaries[1];
1535 dBoundaries[1] = dBoundaries[0];
1536 dBoundaries[0] = Temp;
1537 }
1538 }
1539 }
1540 else if (cmdParser.exprContainsDataObjects())
1541 return findExtremaInData(cmdParser, sExpr, nOrder, nMode);
1542 else
1543 throw SyntaxError(SyntaxError::NO_EXTREMA_VAR, cmdParser.getCommandLine(), SyntaxError::invalid_position);
1544
1545 // Calculate the number of samples depending on
1546 // the interval width
1547 if (intCast(std::abs(dBoundaries[1] - dBoundaries[0])))
1548 nSamples = (nSamples - 1) * intCast(std::abs(dBoundaries[1] - dBoundaries[0])) + 1;
1549
1550 // Ensure that we calculate a reasonable number of samples
1551 if (nSamples > 10001)
1552 nSamples = 10001;
1553
1554 // Set the expression and evaluate it once
1555 _parser.SetExpr(sExpr);
1556 _parser.Eval();
1557 vector<mu::value_type> vResults;
1558 dVal[0] = _parser.Diff(dVar, dBoundaries[0], 1e-7);
1559
1560 // Evaluate the extrema for all samples. We search for
1561 // a sign change in the derivative and examine these intervals
1562 // in more detail
1563 for (unsigned int i = 1; i < nSamples; i++)
1564 {
1565 // Evaluate the derivative at the current sample position
1566 dVal[1] = _parser.Diff(dVar, dBoundaries[0] + (double)i * (dBoundaries[1] - dBoundaries[0]) / (double)(nSamples - 1), 1e-7);
1567
1568 // Is it a sign change or a actual zero?
1569 if (dVal[0].real()*dVal[1].real() < 0)
1570 {
1571 if (!nMode
1572 || (nMode == 1 && (dVal[0].real() > 0 && dVal[1].real() < 0))
1573 || (nMode == -1 && (dVal[0].real() < 0 && dVal[1].real() > 0)))
1574 {
1575 // Examine the current interval in more detail
1576 vResults.push_back(localizeExtremum(sExpr, dVar, _parser, instance->getSettings(), dBoundaries[0] + (double)(i - 1) * (dBoundaries[1] - dBoundaries[0]) / (double)(nSamples - 1), dBoundaries[0] + (double)i * (dBoundaries[1] - dBoundaries[0]) / (double)(nSamples - 1)));
1577 }
1578 }
1579 else if (dVal[0].real()*dVal[1].real() == 0.0)
1580 {
1581 if (!nMode
1582 || (nMode == 1 && (dVal[0].real() > 0 || dVal[1].real() < 0))
1583 || (nMode == -1 && (dVal[0].real() < 0 || dVal[1].real() > 0)))
1584 {
1585 int nTemp = i - 1;
1586
1587 // Jump over multiple zeros due to constness
1588 if (dVal[0] != 0.0)
1589 {
1590 while (dVal[0].real()*dVal[1].real() == 0.0 && i + 1 < nSamples)
1591 {
1592 i++;
1593 dVal[1] = _parser.Diff(dVar, dBoundaries[0] + (double)i * (dBoundaries[1] - dBoundaries[0]) / (double)(nSamples - 1), 1e-7);
1594 }
1595 }
1596 else
1597 {
1598 while (dVal[1] == 0.0 && i + 1 < nSamples)
1599 {
1600 i++;
1601 dVal[1] = _parser.Diff(dVar, dBoundaries[0] + (double)i * (dBoundaries[1] - dBoundaries[0]) / (double)(nSamples - 1), 1e-7);
1602 }
1603 }
1604
1605 // Store the current location
1606 vResults.push_back(localizeExtremum(sExpr, dVar, _parser, instance->getSettings(), dBoundaries[0] + (double)nTemp * (dBoundaries[1] - dBoundaries[0]) / (double)(nSamples - 1), dBoundaries[0] + (double)i * (dBoundaries[1] - dBoundaries[0]) / (double)(nSamples - 1)));
1607 }
1608 }
1609 dVal[0] = dVal[1];
1610 }
1611
1612 // If we didn't find any results
1613 // examine the boundaries for possible extremas
1614 if (!vResults.size())
1615 {
1616 dVal[0] = _parser.Diff(dVar, dBoundaries[0]);
1617 dVal[1] = _parser.Diff(dVar, dBoundaries[1]);
1618 std::string sRetVal;
1619
1620 // Examine the left boundary
1621 if (std::abs(dVal[0])
1622 && (!nMode
1623 || (dVal[0].real() < 0 && nMode == 1)
1624 || (dVal[0].real() > 0 && nMode == -1)))
1625 sRetVal = toString(dBoundaries[0], instance->getSettings().getPrecision());
1626
1627 // Examine the right boundary
1628 if (std::abs(dVal[1])
1629 && (!nMode
1630 || (dVal[1].real() < 0 && nMode == -1)
1631 || (dVal[1].real() > 0 && nMode == 1)))
1632 {
1633 if (sRetVal.length())
1634 sRetVal += ", ";
1635
1636 sRetVal += toString(dBoundaries[1], instance->getSettings().getPrecision());
1637 }
1638
1639 // Still nothing found?
1640 if (!std::abs(dVal[0]) && ! std::abs(dVal[1]))
1641 sRetVal = "nan";
1642
1643 cmdParser.setReturnValue(sRetVal);
1644 }
1645 else
1646 {
1647 cmdParser.setReturnValue(vResults);
1648 }
1649
1650 return true;
1651}
1652
1653
1666static bool findZeroesInMultiResult(CommandLineParser& cmdParser, string& sExpr, string& sInterval, int nMode)
1667{
1668 Parser& _parser = NumeReKernel::getInstance()->getParser();
1669 _parser.SetExpr(sExpr);
1670 int nResults;
1671 mu::value_type* v = _parser.Eval(nResults);
1672 MemoryManager _cache;
1673
1674 vector<mu::value_type> vResults;
1675 int nResults_x = 0;
1676
1677 for (int i = 0; i < nResults; i++)
1678 _cache.writeToTable(i, 1, "table", v[i]);
1679
1680 _parser.SetExpr(sInterval);
1681 v = _parser.Eval(nResults_x);
1682
1683 if (nResults_x > 1)
1684 {
1685 for (int i = 0; i < nResults; i++)
1686 {
1687 if (i >= nResults_x)
1688 _cache.writeToTable(i, 0, "table", 0.0);
1689 else
1690 _cache.writeToTable(i, 0, "table", v[i]);
1691 }
1692 }
1693 else
1694 return false;
1695
1696 std::string sSortingExpr = "sort -table cols=1[2]";
1697 _cache.sortElements(sSortingExpr);
1698
1699 for (long long int i = 1; i < _cache.getLines("table", false); i++)
1700 {
1701 if (isnan(_cache.getElement(i - 1, 1, "table").real()))
1702 continue;
1703
1704 if (!nMode && _cache.getElement(i, 1, "table").real()*_cache.getElement(i - 1, 1, "table").real() <= 0.0)
1705 {
1706 if (_cache.getElement(i, 1, "table") == 0.0)
1707 {
1708 vResults.push_back(_cache.getElement(i, 0, "table"));
1709 i++;
1710 }
1711 else if (_cache.getElement(i - 1, 1, "table") == 0.0)
1712 vResults.push_back(_cache.getElement(i - 1, 0, "table"));
1713 else if (_cache.getElement(i, 1, "table").real()*_cache.getElement(i - 1, 1, "table").real() < 0.0)
1714 vResults.push_back(Linearize(_cache.getElement(i - 1, 0, "table").real(), _cache.getElement(i - 1, 1, "table").real(), _cache.getElement(i, 0, "table").real(), _cache.getElement(i, 1, "table").real()));
1715 }
1716 else if (nMode && _cache.getElement(i, 1, "table").real()*_cache.getElement(i - 1, 1, "table").real() <= 0.0)
1717 {
1718 if (_cache.getElement(i, 1, "table") == 0.0 && _cache.getElement(i - 1, 1, "table") == 0.0)
1719 {
1720 for (long long int j = i + 1; j < _cache.getLines("table", false); j++)
1721 {
1722 if (nMode * _cache.getElement(j, 1, "table").real() > 0.0)
1723 {
1724 for (long long int k = i - 1; k <= j; k++)
1725 vResults.push_back(_cache.getElement(k, 0, "table"));
1726
1727 break;
1728 }
1729 else if (nMode * _cache.getElement(j, 1, "table").real() < 0.0)
1730 break;
1731
1732 if (j + 1 == _cache.getLines("table", false) && i > 1 && nMode * _cache.getElement(i - 2, 1, "table").real() < 0.0)
1733 {
1734 for (long long int k = i - 1; k <= j; k++)
1735 vResults.push_back(_cache.getElement(k, 0, "table"));
1736
1737 break;
1738 }
1739 }
1740
1741 continue;
1742 }
1743 else if (_cache.getElement(i, 1, "table") == 0.0 && nMode * _cache.getElement(i - 1, 1, "table").real() < 0.0)
1744 vResults.push_back(_cache.getElement(i, 0, "table"));
1745 else if (_cache.getElement(i - 1, 1, "table") == 0.0 && nMode * _cache.getElement(i, 1, "table").real() > 0.0)
1746 vResults.push_back(_cache.getElement(i - 1, 0, "table"));
1747 else if (_cache.getElement(i, 1, "table").real()*_cache.getElement(i - 1, 1, "table").real() < 0.0 && nMode * _cache.getElement(i - 1, 1, "table").real() < 0.0)
1748 vResults.push_back(Linearize(_cache.getElement(i - 1, 0, "table").real(), _cache.getElement(i - 1, 1, "table").real(), _cache.getElement(i, 0, "table").real(), _cache.getElement(i, 1, "table").real()));
1749 }
1750 }
1751
1752 if (!vResults.size())
1753 vResults.push_back(NAN);
1754
1755 cmdParser.setReturnValue(vResults);
1756 return true;
1757}
1758
1759
1770static bool findZeroesInData(CommandLineParser& cmdParser, string& sExpr, int nMode)
1771{
1772 mu::value_type* v;
1773 int nResults = 0;
1774 Parser& _parser = NumeReKernel::getInstance()->getParser();
1775 _parser.SetExpr(sExpr);
1776 v = _parser.Eval(nResults);
1777
1778 if (nResults > 1)
1779 {
1780 vector<mu::value_type> vResults;
1781
1782 for (int i = 1; i < nResults; i++)
1783 {
1784 if (isnan(std::abs(v[i - 1])))
1785 continue;
1786
1787 if (!nMode && v[i].real()*v[i - 1].real() <= 0.0)
1788 {
1789 if (v[i] == 0.0)
1790 {
1791 vResults.push_back((double)i + 1);
1792 i++;
1793 }
1794 else if (v[i - 1] == 0.0)
1795 vResults.push_back((double)i);
1796 else if (fabs(v[i]) <= fabs(v[i - 1]))
1797 vResults.push_back((double)i + 1);
1798 else
1799 vResults.push_back((double)i);
1800 }
1801 else if (nMode && v[i].real()*v[i - 1].real() <= 0.0)
1802 {
1803 if (v[i] == 0.0 && v[i - 1] == 0.0)
1804 {
1805 for (int j = i + 1; j < nResults; j++)
1806 {
1807 if (nMode * v[j].real() > 0.0)
1808 {
1809 for (int k = i - 1; k <= j; k++)
1810 vResults.push_back((double)k);
1811
1812 break;
1813 }
1814 else if (nMode * v[j].real() < 0.0)
1815 break;
1816
1817 if (j + 1 == nResults && i > 2 && nMode * v[i - 2].real() < 0.0)
1818 {
1819 for (int k = i - 1; k <= j; k++)
1820 vResults.push_back((double)k);
1821
1822 break;
1823 }
1824 }
1825
1826 continue;
1827 }
1828 else if (v[i] == 0.0 && nMode * v[i - 1].real() < 0.0)
1829 vResults.push_back((double)i + 1);
1830 else if (v[i - 1] == 0.0 && nMode * v[i].real() > 0.0)
1831 vResults.push_back((double)i);
1832 else if (fabs(v[i]) <= fabs(v[i - 1]) && nMode * v[i - 1].real() < 0.0)
1833 vResults.push_back((double)i + 1);
1834 else if (nMode * v[i - 1].real() < 0.0)
1835 vResults.push_back((double)i);
1836 }
1837 }
1838
1839 if (!vResults.size())
1840 vResults.push_back(NAN);
1841
1842 cmdParser.setReturnValue(vResults);
1843 return true;
1844 }
1845 else
1846 throw SyntaxError(SyntaxError::NO_ZEROES_VAR, cmdParser.getCommandLine(), SyntaxError::invalid_position);
1847}
1848
1849
1859bool findZeroes(CommandLineParser& cmdParser)
1860{
1861 NumeReKernel* instance = NumeReKernel::getInstance();
1862 MemoryManager& _data = instance->getMemoryManager();
1863 Parser& _parser = instance->getParser();
1864
1865 unsigned int nSamples = 21;
1866 mu::value_type dVal[2];
1867 mu::value_type dBoundaries[2] = {0.0, 0.0};
1868 int nMode = 0;
1869 mu::value_type* dVar = 0;
1870 mu::value_type dTemp = 0.0;
1871 string sExpr = "";
1872 string sParams = "";
1873 string sInterval = "";
1874 string sVar = "";
1875
1876 // We cannot search zeroes in strings
1877 if (instance->getStringParser().isStringExpression(cmdParser.getExpr()))
1878 throw SyntaxError(SyntaxError::STRINGS_MAY_NOT_BE_EVALUATED_WITH_CMD, cmdParser.getCommandLine(), SyntaxError::invalid_position, "zeroes");
1879
1880 if (!_data.containsTablesOrClusters(cmdParser.getExpr()) && !cmdParser.getParameterList().length())
1881 throw SyntaxError(SyntaxError::NO_ZEROES_OPTIONS, cmdParser.getCommandLine(), SyntaxError::invalid_position);
1882
1883 // Ensure that the expression is not empty
1884 // and that the custom functions don't throw
1885 // any errors
1886 sExpr = cmdParser.getExpr();
1887 sParams = cmdParser.getParameterList();
1888
1889 if (!isNotEmptyExpression(sExpr) || !instance->getDefinitions().call(sExpr))
1890 return false;
1891
1892 if (!instance->getDefinitions().call(sParams))
1893 return false;
1894
1895 StripSpaces(sParams);
1896
1897 // If the expression or the parameter list contains
1898 // data elements, get their values here
1899 if (_data.containsTablesOrClusters(sExpr))
1900 getDataElements(sExpr, _parser, _data, NumeReKernel::getInstance()->getSettings(), false);
1901
1902 if (_data.containsTablesOrClusters(sParams))
1903 getDataElements(sParams, _parser, _data, NumeReKernel::getInstance()->getSettings(), false);
1904
1905 // Evaluate the parameter list
1906 if (findParameter(sParams, "min") || findParameter(sParams, "down"))
1907 nMode = -1;
1908 if (findParameter(sParams, "max") || findParameter(sParams, "up"))
1909 nMode = 1;
1910
1911 if (findParameter(sParams, "samples", '='))
1912 {
1913 _parser.SetExpr(getArgAtPos(sParams, findParameter(sParams, "samples", '=') + 7));
1914 nSamples = intCast(_parser.Eval());
1915
1916 if (nSamples < 21)
1917 nSamples = 21;
1918
1919 sParams.erase(findParameter(sParams, "samples", '=') - 1, 8);
1920 }
1921
1922 // Evaluate the interval
1923 if (sParams.find('=') != string::npos
1924 || (sParams.find('[') != string::npos
1925 && sParams.find(']', sParams.find('['))
1926 && sParams.find(':', sParams.find('['))))
1927 {
1928 if (sParams.substr(0, 2) == "--")
1929 sParams = sParams.substr(2);
1930 else if (sParams.substr(0, 4) == "-set")
1931 sParams = sParams.substr(4);
1932
1933 int nResults = 0;
1934
1935 sInterval = getIntervalForSearchFunctions(sParams, sVar);
1936
1937 _parser.SetExpr(sExpr);
1938 _parser.Eval(nResults);
1939
1940 if (nResults > 1)
1941 return findZeroesInMultiResult(cmdParser, sExpr, sInterval, nMode);
1942 else
1943 {
1944 if (findVariableInExpression(sExpr, sVar) == std::string::npos)
1945 {
1946 cmdParser.setReturnValue("nan");
1947 return true;
1948 }
1949
1950 dVar = getPointerToVariable(sVar, _parser);
1951
1952 if (!dVar)
1953 throw SyntaxError(SyntaxError::ZEROES_VAR_NOT_FOUND, cmdParser.getCommandLine(), sVar, sVar);
1954
1955 if (sInterval.find(':') == string::npos || sInterval.length() < 3)
1956 return false;
1957
1958 auto indices = getAllIndices(sInterval);
1959
1960 for (size_t i = 0; i < 2; i++)
1961 {
1962 if (isNotEmptyExpression(indices[i]))
1963 {
1964 _parser.SetExpr(indices[i]);
1965 dBoundaries[i] = _parser.Eval();
1966
1967 if (isinf(std::abs(dBoundaries[i])) || isnan(std::abs(dBoundaries[i])))
1968 {
1969 cmdParser.setReturnValue("nan");
1970 return false;
1971 }
1972 }
1973 else
1974 return false;
1975 }
1976
1977 if (std::abs(dBoundaries[1]) < std::abs(dBoundaries[0]))
1978 {
1979 mu::value_type Temp = dBoundaries[1];
1980 dBoundaries[1] = dBoundaries[0];
1981 dBoundaries[0] = Temp;
1982 }
1983 }
1984 }
1985 else if (cmdParser.exprContainsDataObjects())
1986 return findZeroesInData(cmdParser, sExpr, nMode);
1987 else
1988 throw SyntaxError(SyntaxError::NO_ZEROES_VAR, cmdParser.getCommandLine(), SyntaxError::invalid_position);
1989
1990 // Calculate the interval
1991 if (intCast(std::abs(dBoundaries[1] - dBoundaries[0])))
1992 nSamples = (nSamples - 1) * intCast(std::abs(dBoundaries[1] - dBoundaries[0])) + 1;
1993
1994 // Ensure that we calculate a reasonable
1995 // amount of samples
1996 if (nSamples > 10001)
1997 nSamples = 10001;
1998
1999 // Set the expression and evaluate it once
2000 _parser.SetExpr(sExpr);
2001 _parser.Eval();
2002
2003 dTemp = *dVar;
2004
2005 *dVar = dBoundaries[0];
2006 vector<mu::value_type> vResults;
2007 dVal[0] = _parser.Eval();
2008
2009 // Find near zeros to the left of the boundary
2010 // which are probably not located due toe rounding
2011 // errors
2012 if (dVal[0] != 0.0 && fabs(dVal[0]) < 1e-10)
2013 {
2014 *dVar = dBoundaries[0] - 1e-10;
2015 dVal[1] = _parser.Eval();
2016
2017 if (dVal[0].real()*dVal[1].real() < 0 && (nMode * dVal[0].real() <= 0.0))
2018 vResults.push_back(localizeExtremum(sExpr, dVar, _parser, instance->getSettings(), dBoundaries[0] - 1e-10, dBoundaries[0]));
2019 }
2020
2021 // Evaluate all samples. We try to find
2022 // sign changes and evaluate the intervals, which
2023 // contain the sign changes, further
2024 for (unsigned int i = 1; i < nSamples; i++)
2025 {
2026 // Evalute the current sample
2027 *dVar = dBoundaries[0] + (double)i * (dBoundaries[1] - dBoundaries[0]) / (double)(nSamples - 1);
2028 dVal[1] = _parser.Eval();
2029
2030 if (dVal[0].real()*dVal[1].real() < 0)
2031 {
2032 if (!nMode
2033 || (nMode == -1 && (dVal[0].real() > 0 && dVal[1].real() < 0))
2034 || (nMode == 1 && (dVal[0].real() < 0 && dVal[1].real() > 0)))
2035 {
2036 // Examine the current interval
2037 vResults.push_back((localizeZero(sExpr, dVar, _parser, instance->getSettings(), dBoundaries[0] + (double)(i - 1) * (dBoundaries[1] - dBoundaries[0]) / (double)(nSamples - 1), dBoundaries[0] + (double)i * (dBoundaries[1] - dBoundaries[0]) / (double)(nSamples - 1))));
2038 }
2039 }
2040 else if (dVal[0].real()*dVal[1].real() == 0.0)
2041 {
2042 if (!nMode
2043 || (nMode == -1 && (dVal[0].real() > 0 || dVal[1].real() < 0))
2044 || (nMode == 1 && (dVal[0].real() < 0 || dVal[1].real() > 0)))
2045 {
2046 int nTemp = i - 1;
2047
2048 // Ignore consecutive zeros due to
2049 // constness
2050 if (dVal[0] != 0.0)
2051 {
2052 while (dVal[0]*dVal[1] == 0.0 && i + 1 < nSamples)
2053 {
2054 i++;
2055 *dVar = dBoundaries[0] + (double)i * (dBoundaries[1] - dBoundaries[0]) / (double)(nSamples - 1);
2056 dVal[1] = _parser.Eval();
2057 }
2058 }
2059 else
2060 {
2061 while (dVal[1] == 0.0 && i + 1 < nSamples)
2062 {
2063 i++;
2064 *dVar = dBoundaries[0] + (double)i * (dBoundaries[1] - dBoundaries[0]) / (double)(nSamples - 1);
2065 dVal[1] = _parser.Eval();
2066 }
2067 }
2068
2069 // Store the result
2070 vResults.push_back(localizeZero(sExpr, dVar, _parser, instance->getSettings(), dBoundaries[0] + (double)nTemp * (dBoundaries[1] - dBoundaries[0]) / (double)(nSamples - 1), dBoundaries[0] + (double)i * (dBoundaries[1] - dBoundaries[0]) / (double)(nSamples - 1)));
2071 }
2072 }
2073
2074 dVal[0] = dVal[1];
2075 }
2076
2077 // Examine the right boundary, because there might be
2078 // a zero slightly right from the interval
2079 if (dVal[0] != 0.0 && fabs(dVal[0]) < 1e-10)
2080 {
2081 *dVar = dBoundaries[1] + 1e-10;
2082 dVal[1] = _parser.Eval();
2083
2084 if (dVal[0].real()*dVal[1].real() < 0 && nMode * dVal[0].real() <= 0.0)
2085 vResults.push_back(localizeZero(sExpr, dVar, _parser, instance->getSettings(), dBoundaries[1], dBoundaries[1] + 1e-10));
2086 }
2087
2088 *dVar = dTemp;
2089
2090 if (!vResults.size())
2091 cmdParser.setReturnValue("nan");
2092 else
2093 cmdParser.setReturnValue(vResults);
2094
2095 return true;
2096}
2097
2098
2118static mu::value_type localizeExtremum(string& sCmd, mu::value_type* dVarAdress, Parser& _parser, const Settings& _option, mu::value_type dLeft, mu::value_type dRight, double dEps, int nRecursion)
2119{
2120 const unsigned int nSamples = 101;
2121 mu::value_type dVal[2];
2122
2123 if (_parser.GetExpr() != sCmd)
2124 {
2125 _parser.SetExpr(sCmd);
2126 _parser.Eval();
2127 }
2128
2129 // Calculate the leftmost value
2130 dVal[0] = _parser.Diff(dVarAdress, dLeft, 1e-7);
2131
2132 // Separate the current interval in
2133 // nSamples steps and examine each step
2134 for (unsigned int i = 1; i < nSamples; i++)
2135 {
2136 // Calculate the next value
2137 dVal[1] = _parser.Diff(dVarAdress, dLeft + (double)i * (dRight - dLeft) / (double)(nSamples - 1), 1e-7);
2138
2139 // Multiply the values to find a sign change
2140 if (dVal[0].real()*dVal[1].real() < 0)
2141 {
2142 // Sign change
2143 // return, if precision is reached. Otherwise perform
2144 // a new recursion between the two values
2145 if (std::abs(dRight - dLeft) / (double)(nSamples - 1) <= dEps || fabs(log(dEps)) + 1 < nRecursion * 2)
2146 return dLeft + (double)(i - 1) * (dRight - dLeft) / (double)(nSamples - 1) + Linearize(0.0, dVal[0].real(), (dRight - dLeft).real() / (double)(nSamples - 1), dVal[1].real());
2147 else
2148 return localizeExtremum(sCmd, dVarAdress, _parser, _option, dLeft + (double)(i - 1) * (dRight - dLeft) / (double)(nSamples - 1), dLeft + (double)i * (dRight - dLeft) / (double)(nSamples - 1), dEps, nRecursion + 1);
2149 }
2150 else if (dVal[0]*dVal[1] == 0.0)
2151 {
2152 // One of the two vwlues is zero.
2153 // Jump over all following zeros due
2154 // to constness
2155 int nTemp = i - 1;
2156
2157 if (dVal[0] != 0.0)
2158 {
2159 while (dVal[0]*dVal[1] == 0.0 && i + 1 < nSamples)
2160 {
2161 i++;
2162 dVal[1] = _parser.Diff(dVarAdress, dLeft + (double)i * (dRight - dLeft) / (double)(nSamples - 1), 1e-7);
2163 }
2164 }
2165 else
2166 {
2167 while (dVal[1] == 0.0 && i + 1 < nSamples)
2168 {
2169 i++;
2170 dVal[1] = _parser.Diff(dVarAdress, dLeft + (double)i * (dRight - dLeft) / (double)(nSamples - 1), 1e-7);
2171 }
2172 }
2173
2174 // return, if precision is reached. Otherwise perform
2175 // a new recursion between the two values
2176 if ((i - nTemp) * std::abs(dRight - dLeft) / (double)(nSamples - 1) <= dEps || (!nTemp && i + 1 == nSamples) || fabs(log(dEps)) + 1 < nRecursion * 2)
2177 return dLeft + (double)nTemp * (dRight - dLeft) / (double)(nSamples - 1) + Linearize(0.0, dVal[0].real(), (i - nTemp) * (dRight - dLeft).real() / (double)(nSamples - 1), dVal[1].real());
2178 else
2179 return localizeExtremum(sCmd, dVarAdress, _parser, _option, dLeft + (double)nTemp * (dRight - dLeft) / (double)(nSamples - 1), dLeft + (double)i * (dRight - dLeft) / (double)(nSamples - 1), dEps, nRecursion + 1);
2180 }
2181
2182 dVal[0] = dVal[1];
2183 }
2184
2185 // If no explict sign change was found,
2186 // interpolate the position by linearisation
2187 *dVarAdress = dLeft;
2188 dVal[0] = _parser.Eval();
2189 *dVarAdress = dRight;
2190 dVal[1] = _parser.Eval();
2191 return Linearize(dLeft.real(), dVal[0].real(), dRight.real(), dVal[1].real());
2192}
2193
2194
2214static mu::value_type localizeZero(string& sCmd, mu::value_type* dVarAdress, Parser& _parser, const Settings& _option, mu::value_type dLeft, mu::value_type dRight, double dEps, int nRecursion)
2215{
2216 const unsigned int nSamples = 101;
2217 mu::value_type dVal[2];
2218
2219 if (_parser.GetExpr() != sCmd)
2220 {
2221 _parser.SetExpr(sCmd);
2222 _parser.Eval();
2223 }
2224
2225 // Calculate the leftmost value
2226 *dVarAdress = dLeft;
2227 dVal[0] = _parser.Eval();
2228
2229 // Separate the current interval in
2230 // nSamples steps and examine each step
2231 for (unsigned int i = 1; i < nSamples; i++)
2232 {
2233 // Calculate the next value
2234 *dVarAdress = dLeft + (double)i * (dRight - dLeft) / (double)(nSamples - 1);
2235 dVal[1] = _parser.Eval();
2236
2237 // Multiply the values to find a sign change
2238 if (dVal[0].real()*dVal[1].real() < 0)
2239 {
2240 // Sign change
2241 // return, if precision is reached. Otherwise perform
2242 // a new recursion between the two values
2243 if (std::abs(dRight - dLeft) / (double)(nSamples - 1) <= dEps || fabs(log(dEps)) + 1 < nRecursion * 2)
2244 return dLeft + (double)(i - 1) * (dRight - dLeft) / (double)(nSamples - 1) + Linearize(0.0, dVal[0].real(), (dRight - dLeft).real() / (double)(nSamples - 1), dVal[1].real());
2245 else
2246 return localizeZero(sCmd, dVarAdress, _parser, _option, dLeft + (double)(i - 1) * (dRight - dLeft) / (double)(nSamples - 1), dLeft + (double)i * (dRight - dLeft) / (double)(nSamples - 1), dEps, nRecursion + 1);
2247 }
2248 else if (dVal[0]*dVal[1] == 0.0)
2249 {
2250 // One of the two vwlues is zero.
2251 // Jump over all following zeros due
2252 // to constness
2253 int nTemp = i - 1;
2254
2255 if (dVal[0] != 0.0)
2256 {
2257 while (dVal[0]*dVal[1] == 0.0 && i + 1 < nSamples)
2258 {
2259 i++;
2260 *dVarAdress = dLeft + (double)i * (dRight - dLeft) / (double)(nSamples - 1);
2261 dVal[1] = _parser.Eval();
2262 }
2263 }
2264 else
2265 {
2266 while (dVal[1] == 0.0 && i + 1 < nSamples)
2267 {
2268 i++;
2269 *dVarAdress = dLeft + (double)i * (dRight - dLeft) / (double)(nSamples - 1);
2270 dVal[1] = _parser.Eval();
2271 }
2272 }
2273
2274 // return, if precision is reached. Otherwise perform
2275 // a new recursion between the two values
2276 if ((i - nTemp) * std::abs(dRight - dLeft) / (double)(nSamples - 1) <= dEps || (!nTemp && i + 1 == nSamples) || fabs(log(dEps)) + 1 < nRecursion * 2)
2277 return dLeft + (double)nTemp * (dRight - dLeft) / (double)(nSamples - 1) + Linearize(0.0, dVal[0].real(), (i - nTemp) * (dRight - dLeft).real() / (double)(nSamples - 1), dVal[1].real());
2278 else
2279 return localizeZero(sCmd, dVarAdress, _parser, _option, dLeft + (double)nTemp * (dRight - dLeft) / (double)(nSamples - 1), dLeft + (double)i * (dRight - dLeft) / (double)(nSamples - 1), dEps, nRecursion + 1);
2280 }
2281
2282 dVal[0] = dVal[1];
2283 }
2284
2285 // If no explict sign change was found,
2286 // interpolate the position by linearisation
2287 *dVarAdress = dLeft;
2288 dVal[0] = _parser.Eval();
2289 *dVarAdress = dRight;
2290 dVal[1] = _parser.Eval();
2291 return Linearize(dLeft.real(), dVal[0].real(), dRight.real(), dVal[1].real());
2292}
2293
2294
2305void taylor(CommandLineParser& cmdParser)
2306{
2307 Parser& _parser = NumeReKernel::getInstance()->getParser();
2308 const Settings& _option = NumeReKernel::getInstance()->getSettings();
2309
2310 const double dPRECISION = 1e-1;
2311 string sVarName = "";
2312 string sExpr = cmdParser.getExprAsMathExpression();
2313 string sExpr_cpy = "";
2314 string sArg = "";
2315 string sTaylor = "Taylor";
2316 string sPolynom = "";
2317 bool bUseUniqueName = cmdParser.hasParam("unique") || cmdParser.hasParam("u");
2318 size_t nth_taylor = 6;
2319 size_t nSamples = 0;
2320 mu::value_type* dVar = 0;
2321 mu::value_type dVarValue = 0.0;
2322 std::vector<mu::value_type> vCoeffs;
2323
2324 // We cannot approximate string expressions
2325 if (containsStrings(sExpr))
2326 throw SyntaxError(SyntaxError::STRINGS_MAY_NOT_BE_EVALUATED_WITH_CMD, cmdParser.getCommandLine(), SyntaxError::invalid_position, "taylor");
2327
2328 // Extract the parameter list
2329 if (!cmdParser.getParameterList().length())
2330 {
2331 NumeReKernel::print(LineBreak(_lang.get("PARSERFUNCS_TAYLOR_MISSINGPARAMS"), _option));
2332 return;
2333 }
2334
2335 // Evaluate the parameters
2336 auto vParVal = cmdParser.getParameterValueAsNumericalValue("n");
2337
2338 if (vParVal.size())
2339 nth_taylor = abs(intCast(vParVal.front()));
2340
2341 std::vector<std::string> vParams = cmdParser.getAllParametersWithValues();
2342
2343 for (const std::string& sPar : vParams)
2344 {
2345 if (sPar != "n")
2346 {
2347 sVarName = sPar;
2348 dVarValue = cmdParser.getParameterValueAsNumericalValue(sVarName).front();
2349
2350 // Ensure that the location was chosen reasonable
2351 if (isinf(std::abs(dVarValue)) || isnan(std::abs(dVarValue)))
2352 return;
2353
2354 // Create the string element, which is used
2355 // for the variable in the created funcction
2356 // string
2357 if (dVarValue == 0.0)
2358 sArg = "x";
2359 else if (dVarValue.real() < 0)
2360 sArg = "x+" + toString(-dVarValue, _option.getPrecision());
2361 else
2362 sArg = "x-" + toString(dVarValue, _option.getPrecision());
2363
2364 break;
2365 }
2366 }
2367
2368 // Extract the expression
2369 sExpr_cpy = sExpr;
2370
2371 // Create a unique function name, if it is desired
2372 if (bUseUniqueName)
2373 sTaylor += toString((int)nth_taylor) + "_" + cmdParser.getExpr();
2374
2375 StripSpaces(sExpr);
2376 _parser.SetExpr(sExpr);
2377
2378 // Ensure that the expression uses the selected variable
2379 if (findVariableInExpression(sExpr, sVarName) == std::string::npos)
2380 {
2381 NumeReKernel::print(LineBreak(_lang.get("PARSERFUNCS_TAYLOR_CONSTEXPR", sVarName), _option));
2382 return;
2383 }
2384
2385 // Get the address of the selected variable
2386 if (sVarName.length())
2387 dVar = getPointerToVariable(sVarName, _parser);
2388
2389 if (!dVar)
2390 return;
2391
2392 // If unique function names are desired,
2393 // generate them here by removing all operators
2394 // from the string
2395 if (bUseUniqueName)
2396 {
2397 string sOperators = " ,;-*/%^!<>&|?:=+[]{}()";
2398
2399 for (unsigned int i = 0; i < sTaylor.length(); i++)
2400 {
2401 if (sOperators.find(sTaylor[i]) != string::npos)
2402 {
2403 sTaylor.erase(i, 1);
2404 i--;
2405 }
2406 }
2407 }
2408
2409 sTaylor += "(x) := ";
2410
2411 // Generate the taylor polynomial
2412 if (!nth_taylor)
2413 {
2414 // zero order polynomial
2415 *dVar = dVarValue;
2416 vCoeffs.push_back(_parser.Eval());
2417 sTaylor += toString(vCoeffs.back(), _option.getPrecision());
2418 }
2419 else
2420 {
2421 // nth order polynomial
2422 *dVar = dVarValue;
2423
2424 // the constant term
2425 vCoeffs.push_back(_parser.Eval());
2426 sPolynom = toString(vCoeffs.back(), _option.getPrecision()) + ",";
2427
2428 nSamples = 6*nth_taylor + 1;
2429 const size_t nFILTERSIZE = 7;
2430 NumeRe::SavitzkyGolayDiffFilter filter(nFILTERSIZE, 1);
2431 std::vector<mu::value_type> vValues(nSamples, 0.0);
2432 std::vector<mu::value_type> vDiffValues;
2433 double dPrec = dPRECISION / nth_taylor;
2434
2435 // Prepare smoothing array
2436 for (size_t i = 0; i < vValues.size(); i++)
2437 {
2438 *dVar = dVarValue + ((int)i - (int)nSamples/2)*dPrec;
2439 vValues[i] = _parser.Eval();
2440 }
2441
2442 // Copy values for easier initialisation
2443 vDiffValues = vValues;
2444
2445 // Perform the derivation
2446 for (size_t n = 0; n < nth_taylor; n++)
2447 {
2448 for (size_t i = nFILTERSIZE/2; i < vValues.size()-nFILTERSIZE/2; i++)
2449 {
2450 vDiffValues[i] = 0.0;
2451
2452 for (size_t j = 0; j < nFILTERSIZE; j++)
2453 vDiffValues[i] += filter.apply(j, 0, vValues[i + j - nFILTERSIZE/2]);
2454
2455 vDiffValues[i] /= dPrec;
2456 }
2457
2458 vCoeffs.push_back(vDiffValues[vDiffValues.size()/2] / (double)integralFactorial(n+1));
2459 sPolynom += toString(vCoeffs.back(), _option.getPrecision()) + ",";
2460 vValues = vDiffValues;
2461 }
2462
2463 sTaylor += "polynomial(" + sArg + "," + sPolynom.substr(0, sPolynom.length()-1) + ")";
2464 }
2465
2466 //if (_option.systemPrints())
2467 // NumeReKernel::print(LineBreak(sTaylor, _option, true, 0, 8));
2468
2469 sTaylor += " " + _lang.get("PARSERFUNCS_TAYLOR_DEFINESTRING", sExpr_cpy, sVarName, toString(dVarValue, 4), toString((int)nth_taylor));
2470
2471 FunctionDefinitionManager& _functions = NumeReKernel::getInstance()->getDefinitions();
2472
2473 bool bDefinitionSuccess = false;
2474
2475 if (_functions.isDefined(sTaylor.substr(0, sTaylor.find(":="))))
2476 bDefinitionSuccess = _functions.defineFunc(sTaylor, true);
2477 else
2478 bDefinitionSuccess = _functions.defineFunc(sTaylor);
2479
2480 if (bDefinitionSuccess)
2481 NumeReKernel::print(_lang.get("DEFINE_SUCCESS"), _option.systemPrints());
2482 else
2483 NumeReKernel::issueWarning(_lang.get("DEFINE_FAILURE"));
2484
2485 cmdParser.setReturnValue(vCoeffs);
2486}
2487
2488
2499static bool detectPhaseOverflow(std::complex<double> cmplx[3])
2500{
2501 return (fabs(std::arg(cmplx[2]) - std::arg(cmplx[1])) >= M_PI) && ((std::arg(cmplx[2]) - std::arg(cmplx[1])) * (std::arg(cmplx[1]) - std::arg(cmplx[0])) < 0);
2502}
2503
2504
2515static std::vector<size_t> getShiftedAxis(size_t nElements, bool inverseTrafo)
2516{
2517 bool isOdd = nElements % 2;
2518 std::vector<size_t> vValues(nElements, 0u);
2519 size_t nOddVal;
2520
2521 // Prepare the axis values
2522 for (size_t i = 0; i < nElements; i++)
2523 {
2524 vValues[i] = i;
2525 }
2526
2527 // Extract the special odd value, if the axis
2528 // length is odd
2529 if (isOdd)
2530 {
2531 if (inverseTrafo)
2532 {
2533 nOddVal = vValues.back();
2534 vValues.pop_back();
2535 }
2536 else
2537 {
2538 nOddVal = vValues.front();
2539 vValues.erase(vValues.begin());
2540 }
2541 }
2542
2543 // Create the actual axis: first part
2544 std::vector<size_t> vAxis(vValues.begin() + nElements / 2+(isOdd && inverseTrafo), vValues.end());
2545
2546 // Insert the odd value, if necessay
2547 if (isOdd)
2548 vAxis.push_back(nOddVal);
2549
2550 // second part
2551 vAxis.insert(vAxis.end(), vValues.begin(), vValues.begin() + nElements / 2+(isOdd && inverseTrafo));
2552
2553 return vAxis;
2554}
2555
2556
2562struct FFTData
2563{
2564 int lines;
2565 int cols;
2566 bool bInverseTrafo;
2567 bool bShiftAxis;
2568 bool bComplex;
2569 double dFrequencyOffset;
2570 double dNyquistFrequency[2];
2571 double dTimeInterval[2];
2572};
2573
2574
2588static void calculate1dFFT(MemoryManager& _data, Indices& _idx, const std::string& sTargetTable, mglDataC& _fftData, std::vector<size_t>& vAxis, FFTData& _fft)
2589{
2590 _fftData.Save("D:/Software/NumeRe/save/fftdata.txt");
2591
2592 if (!_fft.bInverseTrafo)
2593 {
2594 _fftData.FFT("x");
2595
2596 double samples = _fft.lines/2.0;
2597
2598 _fftData.a[0] /= dual(2*samples, 0.0);
2599
2600 for (int i = 1; i < _fftData.GetNx(); i++)
2601 _fftData.a[i] /= dual(samples, 0.0);
2602 }
2603 else
2604 {
2605 double samples = _fftData.GetNx()/2.0;
2606
2607 _fftData.a[0] *= dual(2*samples, 0.0);
2608
2609 for (int i = 1; i < _fftData.GetNx(); i++)
2610 _fftData.a[i] *= dual(samples, 0.0);
2611
2612 _fftData.FFT("ix");
2613 }
2614
2615 if (_idx.col.isOpenEnd())
2616 _idx.col.setRange(0, _idx.col.front() + 3);
2617
2618 if (_fft.bShiftAxis && !_fft.bInverseTrafo)
2619 vAxis = getShiftedAxis(_fft.lines, _fft.bInverseTrafo);
2620
2621 // Store the results of the transformation in the target
2622 // table
2623 if (!_fft.bInverseTrafo)
2624 {
2625 size_t nElements = _fftData.GetNx();
2626 double dPhaseOffset = 0.0;
2627
2628 if (_idx.row.isOpenEnd())
2629 _idx.row.setRange(0, _idx.row.front() + nElements);
2630
2631 for (size_t i = 0; i < nElements; i++)
2632 {
2633 if (i > _idx.row.size())
2634 break;
2635
2636 if (_fft.bShiftAxis)
2637 _data.writeToTable(_idx.row[i], _idx.col.front(), sTargetTable,
2638 _fft.dFrequencyOffset + 2.0 * (double)(i)*_fft.dNyquistFrequency[0] / (double)(_fftData.GetNx()));
2639 else if (i <= nElements/2)
2640 _data.writeToTable(_idx.row[i], _idx.col.front(), sTargetTable,
2641 _fft.dFrequencyOffset + 2.0 * (double)(i)*_fft.dNyquistFrequency[0] / (double)(_fftData.GetNx()));
2642 else
2643 _data.writeToTable(_idx.row[i], _idx.col.front(), sTargetTable,
2644 _fft.dFrequencyOffset + 2.0 * (double)(-(int)nElements+(int)i)*_fft.dNyquistFrequency[0] / (double)(_fftData.GetNx()));
2645
2646 if (!_fft.bComplex)
2647 {
2648 _data.writeToTable(_idx.row[vAxis[i]], _idx.col[1], sTargetTable, std::abs(_fftData.a[i]));
2649
2650 // Stitch phase overflows into a continous array
2651 if (i > 2 && detectPhaseOverflow(&_fftData.a[i-2]))
2652 {
2653 if (std::arg(_fftData.a[i - 1]) - std::arg(_fftData.a[i - 2]) < 0.0)
2654 dPhaseOffset -= 2 * M_PI;
2655 else if (std::arg(_fftData.a[i - 1]) - std::arg(_fftData.a[i - 2]) > 0.0)
2656 dPhaseOffset += 2 * M_PI;
2657 }
2658
2659 _data.writeToTable(_idx.row[vAxis[i]], _idx.col[2], sTargetTable, std::arg(_fftData.a[i]) + dPhaseOffset);
2660 }
2661 else
2662 {
2663 // Only complex values
2664 _data.writeToTable(_idx.row[vAxis[i]], _idx.col[1], sTargetTable, _fftData.a[i]);
2665 }
2666 }
2667
2668 // Write headlines
2669 _data.setHeadLineElement(_idx.col.front(), sTargetTable, _lang.get("COMMON_FREQUENCY") + " [Hz]");
2670
2671 if (!_fft.bComplex)
2672 {
2673 _data.setHeadLineElement(_idx.col[1], sTargetTable, _lang.get("COMMON_AMPLITUDE"));
2674 _data.setHeadLineElement(_idx.col[2], sTargetTable, _lang.get("COMMON_PHASE") + " [rad]");
2675 }
2676 else
2677 _data.setHeadLineElement(_idx.col[1], sTargetTable, _lang.get("COMMON_AMPLITUDE"));
2678 }
2679 else
2680 {
2681 if (_idx.row.isOpenEnd())
2682 _idx.row.setRange(0, _idx.row.front() + _fftData.GetNx());
2683
2684 for (int i = 0; i < _fftData.GetNx(); i++)
2685 {
2686 if (i > (int)_idx.row.size())
2687 break;
2688
2689 _data.writeToTable(_idx.row[i], _idx.col[0], sTargetTable, (double)(i)*_fft.dTimeInterval[0] / (double)(_fftData.GetNx() - 1));
2690 _data.writeToTable(_idx.row[i], _idx.col[1], sTargetTable, _fftData.a[i]);
2691 }
2692
2693 // Write headlines
2694 _data.setHeadLineElement(_idx.col[0], sTargetTable, _lang.get("COMMON_TIME") + " [s]");
2695 _data.setHeadLineElement(_idx.col[1], sTargetTable, _lang.get("COMMON_SIGNAL"));
2696 }
2697}
2698
2699
2713static void calculate2dFFT(MemoryManager& _data, Indices& _idx, const std::string& sTargetTable, mglDataC& _fftData, std::vector<size_t>& vAxis, FFTData& _fft)
2714{
2715 if (!_fft.bInverseTrafo)
2716 {
2717 _fftData.FFT("xy");
2718
2719 double samples = _fft.lines*(_fft.cols-2)/2.0;
2720
2721 _fftData.a[0] /= dual(2*samples, 0.0);
2722
2723 for (long long int i = 1; i < _fftData.GetNN(); i++)
2724 _fftData.a[i] /= dual(samples, 0.0);
2725 }
2726 else
2727 {
2728 double samples = _fft.lines*(_fft.cols-2)/2.0;
2729
2730 _fftData.a[0] *= dual(2*samples, 0.0);
2731
2732 for (long long int i = 1; i < _fftData.GetNN(); i++)
2733 _fftData.a[i] *= dual(samples, 0.0);
2734
2735 _fftData.FFT("iy");
2736 _fftData.FFT("ix");
2737 }
2738
2739 if (_idx.col.isOpenEnd())
2740 _idx.col.setRange(0, _idx.col.front() + _fft.cols);
2741
2742 // Store the results of the transformation in the target
2743 // table
2744 if (!_fft.bInverseTrafo)
2745 {
2746 int nElemsLines = _fftData.GetNx();
2747 int nElemsCols = _fftData.GetNy();
2748
2749 if (_idx.row.isOpenEnd())
2750 _idx.row.setRange(0, _idx.row.front() + std::max(nElemsCols, nElemsLines));
2751
2752 for (int i = 0; i < nElemsLines; i++)
2753 {
2754 if ((size_t)i > _idx.row.size())
2755 break;
2756
2757 // Write x axis
2758 if (_fft.bShiftAxis)
2759 _data.writeToTable(_idx.row[i], _idx.col[0], sTargetTable, _fft.dNyquistFrequency[0]*(-1.0 + 2.0 * (double)(i) / nElemsLines));
2760 else
2761 _data.writeToTable(_idx.row[i], _idx.col[0], sTargetTable, 2.0 * (double)(i)*_fft.dNyquistFrequency[0] / (double)(nElemsLines));
2762
2763 for (int j = 0; j < nElemsCols; j++)
2764 {
2765 // Write the values
2766 if (_fft.bShiftAxis)
2767 _data.writeToTable(_idx.row[i + (i >= nElemsLines/2 ? -nElemsLines/2 : nElemsLines/2+nElemsLines % 2)],
2768 _idx.col[j+2 + (j >= nElemsCols/2 ? -nElemsCols/2 : nElemsCols/2+nElemsCols % 2)],
2769 sTargetTable,
2770 _fft.bComplex ? _fftData.a[i+j*nElemsLines] : std::abs(_fftData.a[i+j*nElemsLines]));
2771 else
2772 _data.writeToTable(_idx.row[i],
2773 _idx.col[j+2],
2774 sTargetTable,
2775 _fft.bComplex ? _fftData.a[i+j*nElemsLines] : std::abs(_fftData.a[i+j*nElemsLines]));
2776 }
2777 }
2778
2779 // Write y axis
2780 for (int i = 0; i < nElemsCols; i++)
2781 {
2782 if (i > (int)_idx.row.size())
2783 break;
2784
2785 if (_fft.bShiftAxis)
2786 _data.writeToTable(_idx.row[i], _idx.col[1], sTargetTable, _fft.dNyquistFrequency[1]*(-1.0 + 2.0 * (double)(i) / nElemsCols));
2787 else
2788 _data.writeToTable(_idx.row[i], _idx.col[1], sTargetTable, 2.0 * (double)(i)*_fft.dNyquistFrequency[1] / (double)(nElemsCols));
2789 }
2790
2791 // Write headlines
2792 _data.setHeadLineElement(_idx.col[0], sTargetTable, _lang.get("COMMON_FREQUENCY") + " [Hz]");
2793 _data.setHeadLineElement(_idx.col[1], sTargetTable, _lang.get("COMMON_FREQUENCY") + " [Hz]");
2794
2795 for (int j = 0; j < nElemsCols; j++)
2796 {
2797 _data.setHeadLineElement(_idx.col[j+2], sTargetTable, _lang.get("COMMON_AMPLITUDE") + "(:," + toString(j+1) + ")");
2798 }
2799 }
2800 else
2801 {
2802 int nElemsLines = _fftData.GetNx();
2803 int nElemsCols = _fftData.GetNy();
2804
2805 if (_idx.row.isOpenEnd())
2806 _idx.row.setRange(0, _idx.row.front() + std::max(nElemsCols, nElemsLines));
2807
2808 for (int i = 0; i < nElemsLines; i++)
2809 {
2810 if ((size_t)i > _idx.row.size())
2811 break;
2812
2813 // Write x axis
2814 _data.writeToTable(_idx.row[i], _idx.col[0], sTargetTable, (double)(i)*_fft.dTimeInterval[0] / (double)(nElemsLines - 1));
2815
2816 // Write the values
2817 for (int j = 0; j < nElemsCols; j++)
2818 {
2819 _data.writeToTable(_idx.row[i], _idx.col[j+2], sTargetTable, _fftData.a[i+j*nElemsLines]);
2820 }
2821 }
2822
2823 // Write y axis
2824 for (int i = 0; i < nElemsCols; i++)
2825 {
2826 if ((size_t)i > _idx.row.size())
2827 break;
2828
2829 _data.writeToTable(_idx.row[i], _idx.col[1], sTargetTable, (double)(i)*_fft.dTimeInterval[1] / (double)(nElemsCols - 1));
2830 }
2831
2832 // Write headlines
2833 _data.setHeadLineElement(_idx.col[0], sTargetTable, _lang.get("COMMON_TIME") + " [s]");
2834 _data.setHeadLineElement(_idx.col[1], sTargetTable, _lang.get("COMMON_TIME") + " [s]");
2835
2836 for (int j = 0; j < nElemsCols; j++)
2837 {
2838 _data.setHeadLineElement(_idx.col[j+2], sTargetTable, _lang.get("COMMON_SIGNAL") + "(:," + toString(j+1) + ")");
2839 }
2840 }
2841}
2842
2843
2855bool fastFourierTransform(CommandLineParser& cmdParser)
2856{
2857 MemoryManager& _data = NumeReKernel::getInstance()->getMemoryManager();
2858 const Settings& _option = NumeReKernel::getInstance()->getSettings();
2859
2860 mglDataC _fftData;
2861 Indices _idx;
2862 FFTData _fft;
2863
2864 _fft.dNyquistFrequency[0] = 1;
2865 _fft.dNyquistFrequency[1] = 1;
2866 _fft.dTimeInterval[0] = 0;
2867 _fft.dTimeInterval[1] = 0;
2868 _fft.dFrequencyOffset = 0.0;
2869 _fft.bInverseTrafo = cmdParser.hasParam("inverse");
2870 _fft.bComplex = cmdParser.hasParam("complex");
2871 _fft.bShiftAxis = cmdParser.hasParam("axisshift");
2872 bool bIs2DFFT = cmdParser.getCommand() == "fft2d";
2873 string sTargetTable = "fftdata";
2874
2875 // search for explicit "target" options and select the target cache
2876 sTargetTable = cmdParser.getTargetTable(_idx, sTargetTable);
2877
2878 DataAccessParser accessParser = cmdParser.getExprAsDataObject();
2879
2880 // get the data from the data object and sort only for the forward transformation
2881 std::unique_ptr<Memory> _mem(extractRange(cmdParser.getCommandLine(), accessParser, bIs2DFFT ? -1 : 2, !_fft.bInverseTrafo));
2882
2883 if (!_mem)
2884 throw SyntaxError(SyntaxError::TABLE_DOESNT_EXIST, cmdParser.getCommandLine(), accessParser.getDataObject() + "(", accessParser.getDataObject() + "()");
2885
2886 _mem->shrink();
2887
2888 _fft.lines = _mem->getElemsInColumn(0);
2889 _fft.cols = _mem->getCols();
2890
2891 if (_fft.lines % 2 && _fft.lines > 1e3)
2892 _fft.lines--;
2893
2894 _fft.dNyquistFrequency[0] = _fft.lines / (_mem->readMem(_fft.lines - 1, 0).real() - _mem->readMem(0, 0).real()) / 2.0;
2895 _fft.dTimeInterval[0] = (_fft.lines - 1) / (_mem->max(VectorIndex(0, VectorIndex::OPEN_END), VectorIndex(0)).real() - _mem->min(VectorIndex(0, VectorIndex::OPEN_END), VectorIndex(0)).real());
2896
2897 if (bIs2DFFT)
2898 {
2899 if (_fft.cols % 2 && _fft.cols > 1e3)
2900 _fft.cols--;
2901
2902 int collines = _mem->getElemsInColumn(1);
2903
2904 _fft.dNyquistFrequency[1] = collines / (_mem->readMem(collines - 1, 1).real() - _mem->readMem(0, 1).real()) / 2.0;
2905 _fft.dTimeInterval[1] = (collines - 1) / (_mem->readMem(collines - 1, 1).real());
2906 }
2907
2908 // Check the dimensions of the input data
2909 if (_fft.lines < 10 || _fft.cols < 2 || (bIs2DFFT && _fft.cols < _mem->getElemsInColumn(1)+2))
2910 throw SyntaxError(SyntaxError::WRONG_DATA_SIZE, cmdParser.getCommandLine(), cmdParser.getExpr());
2911
2912 // Adapt the values for the shifted axis
2913 if (_fft.bShiftAxis)
2914 {
2915 _fft.dFrequencyOffset = -_fft.dNyquistFrequency[0] * (1 + (_fft.lines % 2) * 1.0 / _fft.lines);
2916 _fft.dTimeInterval[0] = fabs((_fft.lines + (_fft.lines % 2)) / (_mem->readMem(0, 0).real())) * 0.5;
2917
2918 if (bIs2DFFT)
2919 _fft.dTimeInterval[1] = fabs((_fft.cols-2 + (_fft.cols % 2)) / (_mem->readMem(0, 1).real())) * 0.5;
2920 }
2921
2922 if (_option.systemPrints())
2923 {
2924 if (!_fft.bInverseTrafo)
2925 NumeReKernel::printPreFmt(LineBreak("|-> " + _lang.get("PARSERFUNCS_FFT_FOURIERTRANSFORMING", toString(_fft.cols)) + " ", _option, 0));
2926 else
2927 NumeReKernel::printPreFmt(LineBreak("|-> " + _lang.get("PARSERFUNCS_FFT_INVERSE_FOURIERTRANSFORMING", toString(_fft.cols)) + " ", _option, 0));
2928 }
2929
2930 if (bIs2DFFT)
2931 _fftData.Create(_fft.lines, _fft.cols-2);
2932 else
2933 _fftData.Create(_fft.lines);
2934
2935 std::vector<size_t> vAxis;
2936
2937 // Prepare the axis (shifted if necessary)
2938 if (_fft.bShiftAxis && _fft.bInverseTrafo)
2939 vAxis = getShiftedAxis(_fft.lines, _fft.bInverseTrafo);
2940 else
2941 {
2942 for (size_t i = 0; i < (size_t)_fft.lines; i++)
2943 vAxis.push_back(i);
2944 }
2945
2946 // Lambda expression for catching and converting NaNs into zeros
2947 auto nanguard = [](const mu::value_type& val) {return mu::isnan(val) ? mu::value_type(0.0) : val;};
2948
2949 // Copy the data
2950 for (int i = 0; i < _fft.lines; i++)
2951 {
2952 if (_fft.cols == 2)
2953 _fftData.a[i] = nanguard(_mem->readMem(vAxis[i], 1)); // Can be complex or not: does not matter
2954 else if (_fft.cols == 3 && _fft.bComplex)
2955 _fftData.a[i] = nanguard(dual(_mem->readMem(vAxis[i], 1).real(), _mem->readMem(vAxis[i], 2).real()));
2956 else if (_fft.cols == 3 && !_fft.bComplex)
2957 _fftData.a[i] = nanguard(dual(_mem->readMem(vAxis[i], 1).real() * cos(_mem->readMem(vAxis[i], 2).real()),
2958 _mem->readMem(vAxis[i], 1).real() * sin(_mem->readMem(vAxis[i], 2).real())));
2959 else if (bIs2DFFT)
2960 {
2961 int nLines = _fft.lines;
2962 int nCols = _fft.cols-2;
2963
2964 for (int j = 0; j < nCols; j++)
2965 {
2966 if (_fft.bShiftAxis && _fft.bInverseTrafo)
2967 _fftData.a[i+j*nLines] = nanguard(_mem->readMem(i + (i >= nLines/2 ? -nLines/2 : nLines/2 + nLines % 2),
2968 j+2 + (j >= nCols/2 ? -nCols/2 : nCols/2 + nCols % 2)));
2969 else
2970 _fftData.a[i+j*nLines] = nanguard(_mem->readMem(i, j+2));
2971 }
2972 }
2973
2974 }
2975
2976 // Calculate the actual transformation and apply some
2977 // normalisation
2978 if (bIs2DFFT)
2979 calculate2dFFT(_data, _idx, sTargetTable, _fftData, vAxis, _fft);
2980 else
2981 calculate1dFFT(_data, _idx, sTargetTable, _fftData, vAxis, _fft);
2982
2983
2984 if (_option.systemPrints())
2985 NumeReKernel::printPreFmt(toSystemCodePage(_lang.get("COMMON_DONE")) + ".\n");
2986
2987 return true;
2988}
2989
2990
3003bool fastWaveletTransform(CommandLineParser& cmdParser)
3004{
3005 MemoryManager& _data = NumeReKernel::getInstance()->getMemoryManager();
3006 const Settings& _option = NumeReKernel::getInstance()->getSettings();
3007
3008 vector<double> vWaveletData;
3009 vector<double> vAxisData;
3010 Indices _idx;
3011
3012 bool bInverseTrafo = cmdParser.hasParam("inverse");
3013 bool bTargetGrid = cmdParser.hasParam("grid");
3014 string sTargetTable = "fwtdata";
3015 string sType = "d"; // d = daubechies, cd = centered daubechies, h = haar, ch = centered haar, b = bspline, cb = centered bspline
3016 int k = 4;
3017
3018 std::string sParVal = cmdParser.getParameterValue("type");
3019
3020 if (sParVal.length())
3021 sType = sParVal;
3022
3023 sParVal = cmdParser.getParameterValue("method");
3024
3025 if (!sType.length() && sParVal.length())
3026 sType = sParVal;
3027
3028 std::vector<mu::value_type> vParVal = cmdParser.getParameterValueAsNumericalValue("k");
3029
3030 if (vParVal.size())
3031 k = intCast(vParVal.front());
3032
3033
3034 // search for explicit "target" options and select the target cache
3035 sTargetTable = cmdParser.getTargetTable(_idx, sTargetTable);
3036
3037 DataAccessParser accessParser = cmdParser.getExprAsDataObject();
3038
3039 // get the data from the data object
3040 std::unique_ptr<Memory> _mem(extractRange(cmdParser.getCommandLine(), accessParser, 2, true));
3041
3042 if (!_mem)
3043 throw SyntaxError(SyntaxError::TABLE_DOESNT_EXIST, cmdParser.getCommandLine(), accessParser.getDataObject() + "(", accessParser.getDataObject() + "()");
3044
3045 if (_option.systemPrints())
3046 {
3047 string sExplType = "";
3048
3049 if (sType.front() == 'c')
3050 sExplType = "Centered ";
3051
3052 if (sType.back() == 'd' || sType.find("daubechies") != string::npos)
3053 sExplType += "Daubechies";
3054 else if (sType.back() == 'h' || sType.find("haar") != string::npos)
3055 sExplType += "Haar";
3056 else if (sType.back() == 'b' || sType.find("bspline") != string::npos)
3057 sExplType += "BSpline";
3058
3059 if (!bInverseTrafo)
3060 NumeReKernel::printPreFmt(LineBreak("|-> " + _lang.get("PARSERFUNCS_WAVELET_TRANSFORMING", sExplType) + " ", _option, 0));
3061 else
3062 NumeReKernel::printPreFmt(LineBreak("|-> " + _lang.get("PARSERFUNCS_WAVELET_INVERSE_TRANSFORMING", sExplType) + " ", _option, 0));
3063 }
3064
3065 for (size_t i = 0; i < (size_t)_mem->getLines(); i++)
3066 {
3067 vWaveletData.push_back(_mem->readMem(i, 1).real());
3068
3069 if (bTargetGrid)
3070 vAxisData.push_back(_mem->readMem(i, 0).real());
3071 }
3072
3073 // calculate the wavelet:
3074 if (sType == "d" || sType == "daubechies")
3075 calculateWavelet(vWaveletData, Daubechies, k, !bInverseTrafo);
3076 else if (sType == "cd" || sType == "cdaubechies")
3077 calculateWavelet(vWaveletData, CenteredDaubechies, k, !bInverseTrafo);
3078 else if (sType == "h" || sType == "haar")
3079 calculateWavelet(vWaveletData, Haar, k, !bInverseTrafo);
3080 else if (sType == "ch" || sType == "chaar")
3081 calculateWavelet(vWaveletData, CenteredHaar, k, !bInverseTrafo);
3082 else if (sType == "b" || sType == "bspline")
3083 calculateWavelet(vWaveletData, BSpline, k, !bInverseTrafo);
3084 else if (sType == "cb" || sType == "cbspline")
3085 calculateWavelet(vWaveletData, CenteredBSpline, k, !bInverseTrafo);
3086
3087 // write the output as datagrid for plotting (only if not an inverse trafo)
3088 if (bTargetGrid && !bInverseTrafo)
3089 {
3090 NumeRe::Table tWaveletData = decodeWaveletData(vWaveletData, vAxisData);
3091
3092 if (_idx.col.isOpenEnd())
3093 _idx.col.setRange(0, _idx.col.front() + tWaveletData.getCols()-1);
3094
3095 if (_idx.row.isOpenEnd())
3096 _idx.row.setRange(0, _idx.row.front() + tWaveletData.getLines()-1);
3097
3098 for (size_t i = 0; i < tWaveletData.getLines(); i++)
3099 {
3100 if (_idx.row[i] == VectorIndex::INVALID)
3101 break;
3102
3103 for (size_t j = 0; j < tWaveletData.getCols(); j++)
3104 {
3105 // write the headlines
3106 if (!i)
3107 {
3108 string sHeadline = "";
3109
3110 if (!j)
3111 sHeadline = _lang.get("COMMON_TIME");
3112 else if (j == 1)
3113 sHeadline = _lang.get("COMMON_LEVEL");
3114 else
3115 sHeadline = _lang.get("COMMON_COEFFICIENT");
3116
3117 _data.setHeadLineElement(_idx.col[j], sTargetTable, sHeadline);
3118 }
3119
3120 if (_idx.col[j] == VectorIndex::INVALID)
3121 break;
3122
3123 _data.writeToTable(_idx.row[i], _idx.col[j], sTargetTable, tWaveletData.getValue(i, j).real());
3124 }
3125 }
3126
3127 if (_option.systemPrints())
3128 NumeReKernel::printPreFmt(toSystemCodePage(_lang.get("COMMON_DONE")) + ".\n");
3129
3130 return true;
3131 }
3132
3133 // write the output as usual data rows
3134 if (_idx.col.isOpenEnd())
3135 _idx.col.setRange(0, _idx.col.front() + 2);
3136
3137 if (_idx.row.isOpenEnd())
3138 _idx.row.setRange(0, _idx.row.front() + vWaveletData.size()-1);
3139
3140 for (long long int i = 0; i < vWaveletData.size(); i++)
3141 {
3142 if (_idx.row[i] == VectorIndex::INVALID)
3143 break;
3144
3145 _data.writeToTable(_idx.row[i], _idx.col[0], sTargetTable, (double)(i));
3146 _data.writeToTable(_idx.row[i], _idx.col[1], sTargetTable, vWaveletData[i]);
3147 }
3148
3149 if (!bInverseTrafo)
3150 {
3151 _data.setHeadLineElement(_idx.col[0], sTargetTable, _lang.get("COMMON_COEFFICIENT"));
3152 _data.setHeadLineElement(_idx.col[1], sTargetTable, _lang.get("COMMON_AMPLITUDE"));
3153 }
3154 else
3155 {
3156 _data.setHeadLineElement(_idx.col[0], sTargetTable, _lang.get("COMMON_TIME"));
3157 _data.setHeadLineElement(_idx.col[1], sTargetTable, _lang.get("COMMON_SIGNAL"));
3158 }
3159
3160 if (_option.systemPrints())
3161 NumeReKernel::printPreFmt(toSystemCodePage(_lang.get("COMMON_DONE")) + ".\n");
3162
3163 return true;
3164}
3165
3166
3175bool evalPoints(CommandLineParser& cmdParser)
3176{
3177 Parser& _parser = NumeReKernel::getInstance()->getParser();
3178 MemoryManager& _data = NumeReKernel::getInstance()->getMemoryManager();
3179 unsigned int nSamples = 100;
3180 mu::value_type* dVar = 0;
3181 mu::value_type dTemp = 0.0;
3182 string sExpr = cmdParser.getExprAsMathExpression();
3183 string sVar = "x";
3184 static string zero = "0.0";
3185 bool bLogarithmic = cmdParser.hasParam("logscale");
3186
3187 // Evaluate calls in the expression
3188 // to any table or cluster
3189 if (_data.containsTablesOrClusters(sExpr))
3190 {
3191 getDataElements(sExpr, _parser, _data, NumeReKernel::getInstance()->getSettings());
3192
3193 if (sExpr.find("{") != string::npos)
3194 convertVectorToExpression(sExpr, NumeReKernel::getInstance()->getSettings());
3195 }
3196
3197 IntervalSet ivl = cmdParser.parseIntervals();
3198
3199 // Extract the interval definition
3200 if (!ivl.size() && cmdParser.getParameterList().find('=') != string::npos)
3201 {
3202 std::vector<std::string> vParams = cmdParser.getAllParametersWithValues();
3203
3204 for (const std::string& sPar : vParams)
3205 {
3206 if (sPar != "samples")
3207 {
3208 sVar = sPar;
3209 std::string sInterval = cmdParser.getParameterValue(sPar);
3210
3211 if (sInterval.front() == '[' && sInterval.back() == ']')
3212 {
3213 sInterval.pop_back();
3214 sInterval.erase(0, 1);
3215 }
3216
3217 auto indices = getAllIndices(sInterval);
3218
3219 _parser.SetExpr(indices[0] + "," + indices[1]);
3220 int nIndices;
3221 mu::value_type* res = _parser.Eval(nIndices);
3222 ivl.intervals.push_back(Interval(res[0], res[1]));
3223
3224 break;
3225 }
3226 }
3227 }
3228
3229 if (!ivl.size())
3230 ivl.intervals.push_back(Interval(-10.0, 10.0));
3231
3232 std::vector<mu::value_type> vSamples = cmdParser.getParameterValueAsNumericalValue("samples");
3233
3234 if (vSamples.size())
3235 nSamples = intCast(vSamples.front());
3236 else if (ivl[0].getSamples())
3237 nSamples = ivl[0].getSamples();
3238
3239 if (isNotEmptyExpression(sExpr))
3240 _parser.SetExpr(sExpr);
3241 else
3242 _parser.SetExpr(sVar);
3243
3244 _parser.Eval();
3245 dVar = getPointerToVariable(sVar, _parser);
3246
3247 if (!dVar)
3248 throw SyntaxError(SyntaxError::EVAL_VAR_NOT_FOUND, cmdParser.getCommandLine(), sVar, sVar);
3249
3250 if (isnan(ivl[0].front().real()) && isnan(ivl[0].back().real()))
3251 ivl[0] = Interval(-10.0, 10.0);
3252 else if (isnan(ivl[0].front().real()) || isnan(ivl[0].back().real()) || isinf(ivl[0].front().real()) || isinf(ivl[0].back().real()))
3253 {
3254 cmdParser.setReturnValue("nan");
3255 return false;
3256 }
3257
3258 if (bLogarithmic && (ivl[0].min() <= 0.0))
3259 throw SyntaxError(SyntaxError::WRONG_PLOT_INTERVAL_FOR_LOGSCALE, cmdParser.getCommandLine(), SyntaxError::invalid_position);
3260
3261 // Set the corresponding expression
3262 if (isNotEmptyExpression(sExpr))
3263 _parser.SetExpr(sExpr);
3264 else if (dVar)
3265 _parser.SetExpr(sVar);
3266 else
3267 _parser.SetExpr(zero);
3268
3269 _parser.Eval();
3270 vector<mu::value_type> vResults;
3271
3272 // Evaluate the selected expression at
3273 // the selected samples
3274 if (dVar)
3275 {
3276 dTemp = *dVar;
3277 *dVar = ivl[0](0);
3278 vResults.push_back(_parser.Eval());
3279
3280 for (unsigned int i = 1; i < nSamples; i++)
3281 {
3282 // Is a logarithmic distribution needed?
3283 if (bLogarithmic)
3284 *dVar = ivl[0].log(i, nSamples);
3285 else
3286 *dVar = ivl[0](i, nSamples);
3287
3288 vResults.push_back(_parser.Eval());
3289 }
3290
3291 *dVar = dTemp;
3292 }
3293 else
3294 {
3295 for (unsigned int i = 0; i < nSamples; i++)
3296 vResults.push_back(_parser.Eval());
3297 }
3298
3299 cmdParser.setReturnValue(vResults);
3300 return true;
3301}
3302
3303
3312bool createDatagrid(CommandLineParser& cmdParser)
3313{
3314 std::vector<size_t> vSamples = {100, 100};
3315 bool bTranspose = cmdParser.hasParam("transpose");
3316
3317 Indices _iTargetIndex;
3318 vector<vector<mu::value_type> > vZVals;
3319
3320 MemoryManager& _data = NumeReKernel::getInstance()->getMemoryManager();
3321
3322 // search for explicit "target" options and select the target cache
3323 std::string sTargetCache = cmdParser.getTargetTable(_iTargetIndex, "grid");
3324
3325 cmdParser.clearReturnValue();
3326 cmdParser.setReturnValue(sTargetCache);
3327
3328 // Get the intervals
3329 IntervalSet ivl = cmdParser.parseIntervals();
3330
3331 // Add missing intervals
3332 while (ivl.size() < 2)
3333 {
3334 ivl.intervals.push_back(Interval(-10.0, 10.0));
3335 }
3336
3337 // Get the number of samples from the option list
3338 auto vParVal = cmdParser.getParameterValueAsNumericalValue("samples");
3339
3340 if (vParVal.size())
3341 {
3342 vSamples.front() = abs(intCast(vParVal.front()));
3343
3344 if (vParVal.size() >= 2)
3345 vSamples[1] = abs(intCast(vParVal[1]));
3346 else
3347 vSamples[1] = vSamples.front();
3348
3349 }
3350
3351 if (vSamples.front() < 2 || vSamples.back() < 2)
3352 throw SyntaxError(SyntaxError::TOO_FEW_DATAPOINTS, cmdParser.getCommandLine(), SyntaxError::invalid_position);
3353
3354 // Get the samples
3355 vSamples = getSamplesForDatagrid(cmdParser, vSamples);
3356
3357 // extract samples from the interval set
3358 if (ivl[0].getSamples())
3359 vSamples[bTranspose] = ivl[0].getSamples();
3360
3361 if (ivl[1].getSamples())
3362 vSamples[1-bTranspose] = ivl[1].getSamples();
3363
3364 //>> Z-Matrix
3365 if (cmdParser.exprContainsDataObjects())
3366 {
3367 // Get the datagrid from another table
3368 DataAccessParser _accessParser = cmdParser.getExprAsDataObject();
3369
3370 if (!_accessParser.getDataObject().length() || _accessParser.isCluster())
3371 throw SyntaxError(SyntaxError::TABLE_DOESNT_EXIST, cmdParser.getCommandLine(), SyntaxError::invalid_position);
3372
3373 Indices& _idx = _accessParser.getIndices();
3374
3375 // identify the table
3376 std::string& szDatatable = _accessParser.getDataObject();
3377
3378 // Check the indices
3379 if (!isValidIndexSet(_idx))
3380 throw SyntaxError(SyntaxError::INVALID_INDEX, cmdParser.getCommandLine(), SyntaxError::invalid_position, _idx.row.to_string() + ", " + _idx.col.to_string());
3381
3382 // the indices are vectors
3383 vector<mu::value_type> vVector;
3384
3385 if (_idx.col.isOpenEnd())
3386 _idx.col.setRange(0, _data.getCols(szDatatable)-1);
3387
3388 if (_idx.row.isOpenEnd())
3389 _idx.row.setRange(0, _data.getColElements(_idx.col.subidx(0), szDatatable)-1);
3390
3391 // Get the data. Choose the order of reading depending on the "transpose" command line option
3392 if (!bTranspose)
3393 {
3394 for (size_t i = 0; i < _idx.row.size(); i++)
3395 {
3396 vVector = _data.getElement(VectorIndex(_idx.row[i]), _idx.col, szDatatable);
3397 vZVals.push_back(vVector);
3398 vVector.clear();
3399 }
3400 }
3401 else
3402 {
3403 for (size_t j = 0; j < _idx.col.size(); j++)
3404 {
3405 vVector = _data.getElement(_idx.row, VectorIndex(_idx.col[j]), szDatatable);
3406 vZVals.push_back(vVector);
3407 vVector.clear();
3408 }
3409 }
3410
3411 // Check the content of the z matrix
3412 if (!vZVals.size() || (vZVals.size() == 1 && vZVals[0].size() == 1))
3413 throw SyntaxError(SyntaxError::TOO_FEW_DATAPOINTS, cmdParser.getCommandLine(), SyntaxError::invalid_position);
3414
3415 // Expand the z vector into a matrix for the datagrid if necessary
3416 expandVectorToDatagrid(ivl, vZVals, vSamples[bTranspose], vSamples[1 - bTranspose]);
3417 }
3418 else
3419 {
3420 Parser& _parser = NumeReKernel::getInstance()->getParser();
3421
3422 // Calculate the grid from formula
3423 _parser.SetExpr(cmdParser.getExprAsMathExpression());
3424
3425 vector<mu::value_type> vVector;
3426
3427 for (unsigned int x = 0; x < vSamples[bTranspose]; x++)
3428 {
3429 _defVars.vValue[0][0] = ivl[0](x, vSamples[bTranspose]);
3430
3431 for (unsigned int y = 0; y < vSamples[1-bTranspose]; y++)
3432 {
3433 _defVars.vValue[1][0] = ivl[1](y, vSamples[1-bTranspose]);
3434 vVector.push_back(_parser.Eval());
3435 }
3436
3437 vZVals.push_back(vVector);
3438 vVector.clear();
3439 }
3440 }
3441
3442 // Store the results in the target cache
3443 if (_iTargetIndex.row.isOpenEnd())
3444 _iTargetIndex.row.setRange(0, _iTargetIndex.row.front() + vSamples[bTranspose] - 1);
3445
3446 if (_iTargetIndex.col.isOpenEnd())
3447 _iTargetIndex.col.setRange(0, _iTargetIndex.col.front() + vSamples[1-bTranspose] + 1);
3448
3449 // Write the x axis
3450 for (size_t i = 0; i < vSamples[bTranspose]; i++)
3451 _data.writeToTable(i, _iTargetIndex.col[0], sTargetCache, ivl[0](i, vSamples[bTranspose]));
3452
3453 _data.setHeadLineElement(_iTargetIndex.col[0], sTargetCache, "x");
3454
3455 // Write the y axis
3456 for (size_t i = 0; i < vSamples[1-bTranspose]; i++)
3457 _data.writeToTable(i, _iTargetIndex.col[1], sTargetCache, ivl[1](i, vSamples[1-bTranspose]));
3458
3459 _data.setHeadLineElement(_iTargetIndex.col[1], sTargetCache, "y");
3460
3461 // Write the z matrix
3462 for (size_t i = 0; i < vZVals.size(); i++)
3463 {
3464 if (_iTargetIndex.row[i] == VectorIndex::INVALID)
3465 break;
3466
3467 for (size_t j = 0; j < vZVals[i].size(); j++)
3468 {
3469 if (_iTargetIndex.col[j+2] == VectorIndex::INVALID)
3470 break;
3471
3472 _data.writeToTable(_iTargetIndex.row[i], _iTargetIndex.col[j+2], sTargetCache, vZVals[i][j]);
3473
3474 if (!i)
3475 _data.setHeadLineElement(_iTargetIndex.col[j+2], sTargetCache, "z(x(:),y(" + toString((int)j + 1) + "))");
3476 }
3477 }
3478
3479 return true;
3480}
3481
3482
3492static vector<size_t> getSamplesForDatagrid(CommandLineParser& cmdParser, const std::vector<size_t>& nSamples)
3493{
3494 vector<size_t> vSamples = nSamples;
3495
3496 // If the z vals are inside of a table then obtain the correct number of samples here
3497 if (cmdParser.exprContainsDataObjects())
3498 {
3499 MemoryManager& _data = NumeReKernel::getInstance()->getMemoryManager();
3500
3501 // Get the indices and identify the table name
3502 DataAccessParser _accessParser = cmdParser.getExprAsDataObject();
3503
3504 if (!_accessParser.getDataObject().length() || _accessParser.isCluster())
3505 throw SyntaxError(SyntaxError::TABLE_DOESNT_EXIST, cmdParser.getCommandLine(), SyntaxError::invalid_position);
3506
3507 Indices& _idx = _accessParser.getIndices();
3508 std::string& sZDatatable = _accessParser.getDataObject();
3509
3510 // Check the indices
3511 if (!isValidIndexSet(_idx))
3512 throw SyntaxError(SyntaxError::INVALID_INDEX, cmdParser.getCommandLine(), SyntaxError::invalid_position, _idx.row.to_string() + ", " + _idx.col.to_string());
3513
3514 // The indices are vectors
3515 if (_idx.col.isOpenEnd())
3516 _idx.col.setRange(0, _data.getCols(sZDatatable)-1);
3517
3518 if (_idx.row.isOpenEnd())
3519 _idx.row.setRange(0, _data.getColElements(_idx.col, sZDatatable)-1);
3520
3521 vSamples.front() = _idx.row.size();
3522 vSamples.back() = _idx.col.size();
3523
3524 // Check for singletons
3525 if (vSamples[0] < 2 && vSamples[1] >= 2)
3526 vSamples[0] = vSamples[1];
3527 else if (vSamples[1] < 2 && vSamples[0] >= 2)
3528 vSamples[1] = vSamples[0];
3529 }
3530
3531 if (vSamples.size() < 2 || vSamples[0] < 2 || vSamples[1] < 2)
3532 throw SyntaxError(SyntaxError::TOO_FEW_DATAPOINTS, cmdParser.getCommandLine(), SyntaxError::invalid_position);
3533
3534 return vSamples;
3535}
3536
3537
3549static void expandVectorToDatagrid(IntervalSet& ivl, vector<vector<mu::value_type>>& vZVals, size_t nSamples_x, size_t nSamples_y)
3550{
3551 // Only if a dimension is a singleton
3552 if (vZVals.size() == 1 || vZVals[0].size() == 1)
3553 {
3554 vector<mu::value_type> vVector;
3555
3556 // construct the needed MGL objects
3557 mglData _mData[4];
3558 mglGraph _graph;
3559
3560 // Prepare the memory
3561 _mData[0].Create(nSamples_x, nSamples_y);
3562 _mData[1].Create(nSamples_x);
3563 _mData[2].Create(nSamples_y);
3564
3565 if (vZVals.size() != 1)
3566 _mData[3].Create(vZVals.size());
3567 else
3568 _mData[3].Create(vZVals[0].size());
3569
3570 // copy the x and y vectors
3571 for (unsigned int i = 0; i < nSamples_x; i++)
3572 _mData[1].a[i] = ivl[0](i, nSamples_x).real();
3573
3574 for (unsigned int i = 0; i < nSamples_y; i++)
3575 _mData[2].a[i] = ivl[1](i, nSamples_y).real();
3576
3577 // copy the z vector
3578 if (vZVals.size() != 1)
3579 {
3580 for (unsigned int i = 0; i < vZVals.size(); i++)
3581 _mData[3].a[i] = vZVals[i][0].real();
3582 }
3583 else
3584 {
3585 for (unsigned int i = 0; i < vZVals[0].size(); i++)
3586 _mData[3].a[i] = vZVals[0][i].real();
3587 }
3588
3589 // Set the ranges needed for the DataGrid function
3590 _graph.SetRanges(_mData[1], _mData[2], _mData[3]);
3591
3592 // Calculate the data grid using a triangulation
3593 _graph.DataGrid(_mData[0], _mData[1], _mData[2], _mData[3]);
3594
3595 vZVals.clear();
3596
3597 // Refill the x and y vectors
3598 ivl[0] = Interval(_mData[1].Minimal(), _mData[1].Maximal());
3599 ivl[1] = Interval(_mData[2].Minimal(), _mData[2].Maximal());
3600
3601 // Copy the z matrix
3602 for (unsigned int i = 0; i < nSamples_x; i++)
3603 {
3604 for (unsigned int j = 0; j < nSamples_y; j++)
3605 vVector.push_back(_mData[0].a[i + nSamples_x * j]);
3606
3607 vZVals.push_back(vVector);
3608 vVector.clear();
3609 }
3610 }
3611}
3612
3613
3622bool writeAudioFile(CommandLineParser& cmdParser)
3623{
3624 MemoryManager& _data = NumeReKernel::getInstance()->getMemoryManager();
3625
3626 string sAudioFileName = "<savepath>/audiofile.wav";
3627 int nSamples = 44100;
3628 int nChannels = 1;
3629 double dMax = 0.0;
3630 double dMin = 0.0;
3631
3632 // Samples lesen
3633 std::vector<mu::value_type> vVals = cmdParser.getParameterValueAsNumericalValue("samples");
3634
3635 if (vVals.size())
3636 nSamples = intCast(vVals.front());
3637
3638 // Dateiname lesen
3639 sAudioFileName = cmdParser.getFileParameterValueForSaving(".wav", "<savepath>", sAudioFileName);
3640
3641 // Indices lesen
3642 DataAccessParser _accessParser = cmdParser.getExprAsDataObject();
3643 Indices& _idx = _accessParser.getIndices();
3644
3645 if (_idx.col.isOpenEnd())
3646 _idx.col.setRange(0, _idx.col.front() + 1);
3647
3648 _accessParser.evalIndices();
3649
3650 if (_idx.col.size() > 2)
3651 return false;
3652
3653 // Find the absolute maximal value
3654 dMin = fabs(_data.min(_accessParser.getDataObject(), _idx.row, _idx.col));
3655 dMax = fabs(_data.max(_accessParser.getDataObject(), _idx.row, _idx.col));
3656
3657 dMax = std::max(dMin, dMax);
3658
3659 nChannels = _idx.col.size() > 1 ? 2 : 1;
3660
3661 std::unique_ptr<Audio::File> audiofile(Audio::getAudioFileByType(sAudioFileName));
3662
3663 if (!audiofile.get())
3664 return false;
3665
3666 audiofile.get()->setChannels(nChannels);
3667 audiofile.get()->setSampleRate(nSamples);
3668 audiofile.get()->newFile();
3669
3670 if (!audiofile.get()->isValid())
3671 return false;
3672
3673 for (size_t i = 0; i < _idx.row.size(); i++)
3674 {
3675 audiofile.get()->write(Audio::Sample(_data.getElement(_idx.row[i], _idx.col[0], _accessParser.getDataObject()).real() / dMax, nChannels > 1 ? _data.getElement(_idx.row[i], _idx.col[1], _accessParser.getDataObject()).real() / dMax : NAN));
3676 }
3677
3678 return true;
3679}
3680
3681
3690bool readAudioFile(CommandLineParser& cmdParser)
3691{
3692 MemoryManager& _data = NumeReKernel::getInstance()->getMemoryManager();
3693
3694 Indices _targetIdx;
3695 std::string sTarget = cmdParser.getTargetTable(_targetIdx, "");
3696 std::string sAudioFile = cmdParser.getExprAsFileName(".wav");
3697
3698 g_logger.info("Load audiofile '" + sAudioFile + "'.");
3699
3700 // Read the whole table or only the metadata
3701 if (sTarget.length())
3702 {
3703 std::unique_ptr<Audio::File> audiofile(Audio::getAudioFileByType(sAudioFile));
3704
3705 if (!audiofile.get() || !audiofile.get()->isValid())
3706 return false;
3707
3708 size_t nLen = audiofile.get()->getLength();
3709 size_t nChannels = audiofile.get()->getChannels();
3710
3711 // Try to read the entire file
3712 _data.resizeTable(nChannels > 1 && _targetIdx.col.size() > 1 ? _targetIdx.col.subidx(0, 2).max()+1 : _targetIdx.col.front()+1,
3713 sTarget);
3714
3715 Memory* _table = _data.getTable(sTarget);
3716 int rowmax = _targetIdx.row.subidx(0, nLen).max();
3717
3718 // Write the last row for preallocation
3719 _table->writeData(rowmax, _targetIdx.col.front(), 0.0);
3720
3721 if (nChannels > 1 && _targetIdx.col.size() > 1)
3722 _table->writeData(rowmax, _targetIdx.col[1], 0.0);
3723
3724 for (size_t i = 0; i < nLen; i++)
3725 {
3726 Audio::Sample sample = audiofile.get()->read();
3727
3728 if (_targetIdx.row.size() <= i)
3729 break;
3730
3731 _table->writeDataDirectUnsafe(_targetIdx.row[i], _targetIdx.col[0], sample.leftOrMono);
3732
3733 if (nChannels > 1 && _targetIdx.col.size() > 1)
3734 _table->writeDataDirectUnsafe(_targetIdx.row[i], _targetIdx.col[1], sample.right);
3735 }
3736
3737 _table->markModified();
3738
3739 // Create the storage indices
3740 std::vector<mu::value_type> vIndices = {_targetIdx.row.min()+1,
3741 _targetIdx.row.size() < nLen ? _targetIdx.row.max()+1 : _targetIdx.row.min()+nLen,
3742 (nChannels > 1 && _targetIdx.col.size() > 1 ? _targetIdx.col.subidx(0, 2).min()+1 : _targetIdx.col.front()+1),
3743 (nChannels > 1 && _targetIdx.col.size() > 1 ? _targetIdx.col.subidx(0, 2).max()+1 : _targetIdx.col.front()+1)};
3744
3745 cmdParser.setReturnValue(vIndices);
3746 g_logger.info("Audiofile read.");
3747 }
3748 else
3749 {
3750 std::unique_ptr<Audio::File> audiofile(Audio::getAudioFileByType(sAudioFile));
3751
3752 if (!audiofile.get() || !audiofile.get()->isValid())
3753 return false;
3754
3755 // Only read the metadata
3756 size_t nLen = audiofile.get()->getLength();
3757 size_t nChannels = audiofile.get()->getChannels();
3758 size_t nSampleRate = audiofile.get()->getSampleRate();
3759
3760 // Create the metadata
3761 std::vector<mu::value_type> vMetaData = {nLen, nChannels, nSampleRate, nLen / (double)nSampleRate};
3762
3763 cmdParser.setReturnValue(vMetaData);
3764 g_logger.info("Audiofile metadata read.");
3765 }
3766
3767 return true;
3768}
3769
3770
3779bool seekInAudioFile(CommandLineParser& cmdParser)
3780{
3781 MemoryManager& _data = NumeReKernel::getInstance()->getMemoryManager();
3782 Indices _targetIdx;
3783 std::string sTarget = cmdParser.getTargetTable(_targetIdx, "audiodata");
3784 std::vector<mu::value_type> vSeekIndices = cmdParser.parseExprAsNumericalValues();
3785
3786 if (vSeekIndices.size() < 2)
3787 return false;
3788
3789 std::string sAudioFile = cmdParser.getFileParameterValue(".wav");
3790
3791 g_logger.info("Load audiofile '" + sAudioFile + "'.");
3792 std::unique_ptr<Audio::File> audiofile(Audio::getAudioFileByType(sAudioFile));
3793
3794 if (!audiofile.get() || !audiofile.get()->isValid())
3795 return false;
3796
3797 size_t nLen = audiofile.get()->getLength();
3798 size_t nChannels = audiofile.get()->getChannels();
3799
3800 if (!audiofile.get()->isSeekable() || std::max(vSeekIndices[0].real()-1, 0.0) >= nLen)
3801 return false;
3802
3803 std::unique_ptr<Audio::SeekableFile> seekable(static_cast<Audio::SeekableFile*>(audiofile.release()));
3804
3805 seekable.get()->setPosition(std::max(vSeekIndices[0].real()-1, 0.0));
3806 nLen = std::min(nLen - seekable.get()->getPosition(), (size_t)(std::max(vSeekIndices[1].real(), 0.0)));
3807
3808 // Try to read the desired length from the file
3809 _data.resizeTable(nChannels > 1 && _targetIdx.col.size() > 1 ? _targetIdx.col.subidx(0, 2).max()+1 : _targetIdx.col.front()+1,
3810 sTarget);
3811
3812 Memory* _table = _data.getTable(sTarget);
3813 int rowmax = _targetIdx.row.subidx(0, nLen).max();
3814
3815 // Write the last row for pre-allocation
3816 _table->writeData(rowmax, _targetIdx.col.front(), 0.0);
3817
3818 if (nChannels > 1 && _targetIdx.col.size() > 1)
3819 _table->writeData(rowmax, _targetIdx.col[1], 0.0);
3820
3821 for (size_t i = 0; i < nLen; i++)
3822 {
3823 Audio::Sample sample = seekable.get()->read();
3824
3825 if (_targetIdx.row.size() <= i)
3826 break;
3827
3828 _table->writeDataDirectUnsafe(_targetIdx.row[i], _targetIdx.col[0], sample.leftOrMono);
3829
3830 if (nChannels > 1 && _targetIdx.col.size() > 1)
3831 _table->writeDataDirectUnsafe(_targetIdx.row[i], _targetIdx.col[1], sample.right);
3832 }
3833
3834 _table->markModified();
3835 cmdParser.setReturnValue(toString(nLen));
3836 g_logger.info("Seeked portion read.");
3837
3838 return true;
3839}
3840
3841
3851bool regularizeDataSet(CommandLineParser& cmdParser)
3852{
3853 MemoryManager& _data = NumeReKernel::getInstance()->getMemoryManager();
3854 int nSamples = 0;
3855 string sColHeaders[2] = {"", ""};
3856 double dXmin, dXmax;
3857
3858 // Samples lesen
3859 auto vParVal = cmdParser.getParameterValueAsNumericalValue("samples");
3860
3861 if (vParVal.size())
3862 nSamples = intCast(vParVal.front());
3863
3864 // Indices lesen
3865 DataAccessParser accessParser = cmdParser.getExprAsDataObject();
3866
3867 std::unique_ptr<Memory> _mem(extractRange(cmdParser.getCommandLine(), accessParser, 2, true));
3868
3869 if (!_mem)
3870 throw SyntaxError(SyntaxError::TABLE_DOESNT_EXIST, cmdParser.getCommandLine(), accessParser.getDataObject() + "(", accessParser.getDataObject() + "()");
3871
3872 sColHeaders[0] = _mem->getHeadLineElement(0) + "\n(regularized)";
3873 sColHeaders[1] = _mem->getHeadLineElement(1) + "\n(regularized)";
3874
3875 long long int nLines = _mem->getLines();
3876
3877 dXmin = _mem->min(VectorIndex(0, VectorIndex::OPEN_END), VectorIndex(0)).real();
3878 dXmax = _mem->max(VectorIndex(0, VectorIndex::OPEN_END), VectorIndex(0)).real();
3879
3880 // Create splines
3881 tk::spline _spline;
3882 _spline.set_points(mu::real(_mem->readMem(VectorIndex(0, nLines-1), VectorIndex(0))),
3883 mu::real(_mem->readMem(VectorIndex(0, nLines-1), VectorIndex(1))), false);
3884
3885 if (!nSamples)
3886 nSamples = nLines;
3887
3888 long long int nLastCol = _data.getCols(accessParser.getDataObject(), false);
3889
3890 // Interpolate the data points
3891 for (long long int i = 0; i < nSamples; i++)
3892 {
3893 _data.writeToTable(i, nLastCol, accessParser.getDataObject(), dXmin + i * (dXmax-dXmin) / (nSamples-1));
3894 _data.writeToTable(i, nLastCol + 1, accessParser.getDataObject(), _spline(dXmin + i*(dXmax-dXmin) / (nSamples-1)));
3895 }
3896
3897 _data.setHeadLineElement(nLastCol, accessParser.getDataObject(), sColHeaders[0]);
3898 _data.setHeadLineElement(nLastCol + 1, accessParser.getDataObject(), sColHeaders[1]);
3899 return true;
3900}
3901
3902
3915bool analyzePulse(CommandLineParser& cmdParser)
3916{
3917 mglData _v;
3918 vector<mu::value_type> vPulseProperties;
3919 double dXmin = NAN, dXmax = NAN;
3920 double dSampleSize = NAN;
3921
3922 // Indices lesen
3923 DataAccessParser accessParser = cmdParser.getExprAsDataObject();
3924
3925 std::unique_ptr<Memory> _mem(extractRange(cmdParser.getCommandLine(), accessParser, 2, true));
3926
3927 if (!_mem)
3928 throw SyntaxError(SyntaxError::TABLE_DOESNT_EXIST, cmdParser.getCommandLine(), accessParser.getDataObject() + "(", accessParser.getDataObject() + "()");
3929
3930 long long int nLines = _mem->getLines();
3931
3932 dXmin = _mem->min(VectorIndex(0, VectorIndex::OPEN_END), VectorIndex(0)).real();
3933 dXmax = _mem->max(VectorIndex(0, VectorIndex::OPEN_END), VectorIndex(0)).real();
3934
3935 _v.Create(nLines);
3936
3937 for (long long int i = 0; i < nLines; i++)
3938 _v.a[i] = _mem->readMem(i, 1).real();
3939
3940 dSampleSize = (dXmax - dXmin) / ((double)_v.GetNx() - 1.0);
3941 mglData _pulse(_v.Pulse('x'));
3942
3943 if (_pulse.nx >= 5)
3944 {
3945 vPulseProperties.push_back(_pulse[0]); // max Amp
3946 vPulseProperties.push_back(_pulse[1]*dSampleSize + dXmin); // pos max Amp
3947 vPulseProperties.push_back(2.0 * _pulse[2]*dSampleSize); // FWHM
3948 vPulseProperties.push_back(2.0 * _pulse[3]*dSampleSize); // Width near max
3949 vPulseProperties.push_back(_pulse[4]*dSampleSize); // Energy (Integral pulse^2)
3950 }
3951 else
3952 {
3953 cmdParser.setReturnValue("nan");
3954 return true;
3955 }
3956
3957 // Ausgabe
3958 if (NumeReKernel::getInstance()->getSettings().systemPrints())
3959 {
3960 NumeReKernel::toggleTableStatus();
3961 make_hline();
3962 NumeReKernel::print("NUMERE: " + toUpperCase(_lang.get("PARSERFUNCS_PULSE_HEADLINE")));
3963 make_hline();
3964
3965 for (unsigned int i = 0; i < vPulseProperties.size(); i++)
3966 NumeReKernel::printPreFmt("| " + _lang.get("PARSERFUNCS_PULSE_TABLE_" + toString((int)i + 1) + "_*", toString(vPulseProperties[i], NumeReKernel::getInstance()->getSettings().getPrecision())) + "\n");
3967
3968 NumeReKernel::toggleTableStatus();
3969 make_hline();
3970 }
3971
3972 cmdParser.setReturnValue(vPulseProperties);
3973 return true;
3974}
3975
3976
3985bool shortTimeFourierAnalysis(CommandLineParser& cmdParser)
3986{
3987 MemoryManager& _data = NumeReKernel::getInstance()->getMemoryManager();
3988
3989 Indices _target;
3990 mglData _real, _imag, _result;
3991 int nSamples = 0;
3992
3993 double dXmin = NAN, dXmax = NAN;
3994 double dFmin = 0.0, dFmax = 1.0;
3995 double dSampleSize = NAN;
3996
3997 auto vParVal = cmdParser.getParameterValueAsNumericalValue("samples");
3998
3999 if (vParVal.size())
4000 nSamples = std::max(intCast(vParVal.front()), 0LL);
4001
4002 std::string sTargetCache = cmdParser.getTargetTable(_target, "stfdat");
4003
4004 // Indices lesen
4005 DataAccessParser _accessParser = cmdParser.getExprAsDataObject();
4006 _accessParser.evalIndices();
4007 Indices& _idx = _accessParser.getIndices();
4008
4009 dXmin = _data.min(_accessParser.getDataObject(), _idx.row, _idx.col.subidx(0, 1)).real();
4010 dXmax = _data.max(_accessParser.getDataObject(), _idx.row, _idx.col.subidx(0, 1)).real();
4011
4012 _real.Create(_idx.row.size());
4013 _imag.Create(_idx.row.size());
4014
4015 for (size_t i = 0; i < _idx.row.size(); i++)
4016 {
4017 mu::value_type val = _data.getElement(_idx.row[i], _idx.col[1], _accessParser.getDataObject());
4018 _real.a[i] = val.real();
4019 _imag.a[i] = val.imag();
4020 }
4021
4022 if (!nSamples || nSamples > _real.GetNx())
4023 nSamples = _real.GetNx() / 32;
4024
4025 // Tatsaechliche STFA
4026 _result = mglSTFA(_real, _imag, nSamples);
4027
4028 dSampleSize = (dXmax - dXmin) / ((double)_result.GetNx() - 1.0);
4029
4030 // Nyquist: _real.GetNx()/(dXmax-dXmin)/2.0
4031 dFmax = _real.GetNx() / (dXmax - dXmin) / 2.0;
4032
4033 // Zielcache befuellen entsprechend der Fourier-Algorithmik
4034
4035 //if (_target.row.isOpenEnd())
4036 // _target.row.setRange(0, _target.row.front() + _result.GetNx() - 1);
4037
4038 if (_target.col.isOpenEnd())
4039 _target.col.setRange(0, _target.col.front() + _result.GetNy()/2 + 1);
4040
4041 _data.resizeTable(_target.col.max(), sTargetCache);
4042
4043 // Write the time axis
4044 for (int i = 0; i < _result.GetNx(); i++)
4045 _data.writeToTable(_target.row[i], _target.col.front(), sTargetCache, dXmin + i * dSampleSize);
4046
4047 // Define headline
4048 _data.setHeadLineElement(_target.col.front(), sTargetCache, _data.getHeadLineElement(_idx.col.front(), _accessParser.getDataObject()));
4049 dSampleSize = 2 * (dFmax - dFmin) / ((double)_result.GetNy() - 1.0);
4050
4051 // Write the frequency axis
4052 for (int i = 0; i < _result.GetNy() / 2; i++)
4053 _data.writeToTable(_target.row[i], _target.col[1], sTargetCache, dFmin + i * dSampleSize); // Fourier f
4054
4055 // Define headline
4056 _data.setHeadLineElement(_target.col[1], sTargetCache, "f [Hz]");
4057
4058 // Write the STFA map
4059 for (int i = 0; i < _result.GetNx(); i++)
4060 {
4061 if (_target.row[i] == VectorIndex::INVALID)
4062 break;
4063
4064 for (int j = 0; j < _result.GetNy() / 2; j++)
4065 {
4066 if (_target.col[j+2] == VectorIndex::INVALID)
4067 break;
4068
4069 _data.writeToTable(_target.row[i], _target.col[j+2], sTargetCache, _result[i + (j + _result.GetNy() / 2)*_result.GetNx()]);
4070
4071 // Update the headline
4072 if (!i)
4073 _data.setHeadLineElement(_target.col[j+2], sTargetCache, "A(" + toString((int)j + 1) + ")");
4074 }
4075 }
4076
4077 cmdParser.clearReturnValue();
4078 cmdParser.setReturnValue(sTargetCache);
4079 return true;
4080}
4081
4082
4094static double interpolateToGrid(const std::vector<double>& vAxis, double dInterpolVal, bool bExtent = false)
4095{
4096 if (isnan(dInterpolVal))
4097 return dInterpolVal;
4098
4099 // Find the base index
4100 int nBaseVal = intCast(dInterpolVal) + (dInterpolVal < 0 ? -1 : 0);
4101
4102 // Get the decimal part
4103 double x = dInterpolVal - nBaseVal;
4104
4105 if (nBaseVal >= 0 && nBaseVal+1 < (int)vAxis.size())
4106 return vAxis[nBaseVal] + (vAxis[nBaseVal+1] - vAxis[nBaseVal]) * x;
4107 else if (nBaseVal > 0 && nBaseVal < (int)vAxis.size()) // Repeat last distance
4108 return vAxis[nBaseVal] + (vAxis[nBaseVal] - vAxis[nBaseVal-1]) * x;
4109 else if (nBaseVal == -1 && nBaseVal+1 < (int)vAxis.size()) // Repeat first distance
4110 return vAxis.front() + (-1 + x) * (vAxis[1] - vAxis.front());
4111
4112 if (bExtent)
4113 {
4114 if (nBaseVal >= (int)vAxis.size())
4115 return vAxis.back() + (nBaseVal - vAxis.size() + 1 + x) * (vAxis.back() - vAxis[vAxis.size()-2]);
4116 else if (nBaseVal < -1)
4117 return vAxis.front() + (nBaseVal + x) * (vAxis[1] - vAxis.front());
4118 }
4119
4120 // nBaseVal below zero or not in the grid
4121 return NAN;
4122}
4123
4124
4133void boneDetection(CommandLineParser& cmdParser)
4134{
4135 MemoryManager& _data = NumeReKernel::getInstance()->getMemoryManager();
4136 mglData _mData;
4137 double dLevel = NAN, dAttrX = 0.0, dAttrY = 1.0, dMinLen = 0.0;
4138
4139 // detect TABLE() -set minval=LVL attract={x,y} minlen=MIN target=TARGET()
4140 auto vParVal = cmdParser.getParameterValueAsNumericalValue("minval");
4141
4142 if (vParVal.size())
4143 dLevel = vParVal.front().real();
4144
4145 vParVal = cmdParser.getParameterValueAsNumericalValue("attract");
4146
4147 if (vParVal.size() > 1)
4148 {
4149 dAttrX = fabs(vParVal[0]);
4150 dAttrY = fabs(vParVal[1]);
4151 }
4152 else if (vParVal.size())
4153 dAttrY = fabs(vParVal[0]);
4154
4155 vParVal = cmdParser.getParameterValueAsNumericalValue("minlen");
4156
4157 if (vParVal.size())
4158 dMinLen = fabs(vParVal.front());
4159
4160 Indices _target;
4161 std::string sTargetCache = cmdParser.getTargetTable(_target, "detectdat");
4162
4163 // Indices lesen
4164 DataAccessParser accessParser = cmdParser.getExprAsDataObject();
4165 Indices& _idx = accessParser.getIndices();
4166
4167 if (_idx.row.isOpenEnd())
4168 _idx.row.setRange(0, _data.getLines(accessParser.getDataObject())-1);
4169
4170 if (_idx.col.isOpenEnd())
4171 _idx.col.setRange(0, _idx.col.front() + _data.getLines(accessParser.getDataObject(), true) - _data.getAppendedZeroes(_idx.col[1], accessParser.getDataObject()) + 1);
4172
4173 // Get x and y axis for the final scaling
4174 std::vector<mu::value_type> vX = _data.getElement(_idx.row, _idx.col.subidx(0, 1), accessParser.getDataObject());
4175 std::vector<mu::value_type> vY = _data.getElement(_idx.row, _idx.col.subidx(1, 1), accessParser.getDataObject());
4176
4177 _idx.col = _idx.col.subidx(2);
4178
4179 // Get the data
4180 std::unique_ptr<Memory> _mem(extractRange(cmdParser.getCommandLine(), accessParser, 100));
4181
4182 if (!_mem)
4183 throw SyntaxError(SyntaxError::TABLE_DOESNT_EXIST, cmdParser.getCommandLine(), accessParser.getDataObject() + "(", accessParser.getDataObject() + "()");
4184
4185 long long int nLines = _mem->getLines();
4186 long long int nCols = _mem->getCols();
4187
4188 // Restrict the attraction to not exceed the axis range
4189 dAttrX = min(dAttrX, (double)nLines);
4190 dAttrY = min(dAttrY, (double)nCols);
4191
4192 // Use minimal data value if level is NaN
4193 if (isnan(dLevel))
4194 dLevel = _mem->min(VectorIndex(0, VectorIndex::OPEN_END), VectorIndex(0, VectorIndex::OPEN_END)).real();
4195
4196 _mData.Create(nLines, nCols);
4197
4198 // Copy the data to the mglData object
4199 for (long long int i = 0; i < nLines; i++)
4200 {
4201 for (long long int j = 0; j < nCols; j++)
4202 {
4203 _mData.a[i+j*nLines] = _mem->readMem(i, j).real();
4204 }
4205 }
4206
4207 // Perform the actual bone detection
4208 mglData _res = _mData.Detect(dLevel, dAttrY, dAttrX, dMinLen);
4209
4210 if (_target.row.isOpenEnd())
4211 _target.row.setOpenEndIndex(_res.GetNy() + _target.row.front());
4212
4213 if (_target.col.isOpenEnd())
4214 _target.col.setOpenEndIndex(_res.GetNx() + _target.col.front());
4215
4216 // Copy the results to the target table
4217 for (int i = 0; i < _res.GetNy(); i++)
4218 {
4219 if (_target.row.size() <= (size_t)i)
4220 break;
4221
4222 for (int j = 0; j < _res.GetNx(); j++)
4223 {
4224 if (_target.col.size() <= (size_t)j)
4225 break;
4226
4227 if (!j)
4228 _data.writeToTable(_target.row[i], _target.col[j], sTargetCache, interpolateToGrid(mu::real(vX), _res.a[j+i*_res.GetNx()]));
4229 else
4230 _data.writeToTable(_target.row[i], _target.col[j], sTargetCache, interpolateToGrid(mu::real(vY), _res.a[j+i*_res.GetNx()]));
4231 }
4232 }
4233
4234 // Write axis labels
4235 _data.setHeadLineElement(_target.col[0], sTargetCache, "Structure_x");
4236 _data.setHeadLineElement(_target.col[1], sTargetCache, "Structure_y");
4237
4238 if (NumeReKernel::getInstance()->getSettings().systemPrints())
4239 NumeReKernel::print(_lang.get("BUILTIN_CHECKKEYWORD_DETECT_SUCCESS", sTargetCache));
4240}
4241
4242
4253bool calculateSplines(CommandLineParser& cmdParser)
4254{
4255 FunctionDefinitionManager& _functions = NumeReKernel::getInstance()->getDefinitions();
4256
4257 tk::spline _spline;
4258 vector<double> xVect, yVect;
4259
4260 DataAccessParser accessParser = cmdParser.getExprAsDataObject();
4261
4262 std::unique_ptr<Memory> _mem(extractRange(cmdParser.getCommandLine(), accessParser, 2, true));
4263
4264 if (!_mem)
4265 throw SyntaxError(SyntaxError::TABLE_DOESNT_EXIST, cmdParser.getCommandLine(), accessParser.getDataObject() + "(", accessParser.getDataObject() + "()");
4266
4267 int nLines = _mem->getLines();
4268
4269 if (nLines < 2)
4270 throw SyntaxError(SyntaxError::TOO_FEW_DATAPOINTS, cmdParser.getCommandLine(), accessParser.getDataObject());
4271
4272 for (int i = 0; i < nLines; i++)
4273 {
4274 xVect.push_back(_mem->readMem(i, 0).real());
4275 yVect.push_back(_mem->readMem(i, 1).real());
4276 }
4277
4278 // Set the points for the spline to calculate
4279 _spline.set_points(xVect, yVect);
4280
4281 string sDefinition = "Spline(x) := ";
4282
4283 // Create the polynomial, which will be defined for future use
4284 for (size_t i = 0; i < xVect.size() - 1; i++)
4285 {
4286 string sRange = "polynomial(";
4287
4288 if (xVect[i] == 0)
4289 sRange += "x,";
4290 else if (xVect[i] < 0)
4291 sRange += "x+" + toString(fabs(xVect[i]), 4) + ",";
4292 else
4293 sRange += "x-" + toString(xVect[i], 4) + ",";
4294
4295 vector<double> vCoeffs = _spline[i];
4296 sRange += toString(vCoeffs[0], 4) + "," + toString(vCoeffs[1], 4) + "," + toString(vCoeffs[2], 4) + "," + toString(vCoeffs[3], 4) + ")";
4297
4298 if (i == xVect.size() - 2)
4299 sRange += "*ivl(x," + toString(xVect[i], 4) + "," + toString(xVect[i + 1], 4) + ",1,1)";
4300 else
4301 sRange += "*ivl(x," + toString(xVect[i], 4) + "," + toString(xVect[i + 1], 4) + ",1,2)";
4302
4303 sDefinition += sRange;
4304
4305 if (i < xVect.size() - 2)
4306 sDefinition += " + ";
4307 }
4308
4309 if (NumeReKernel::getInstance()->getSettings().systemPrints() && !NumeReKernel::bSupressAnswer)
4310 NumeReKernel::print(sDefinition);
4311
4312 bool bDefinitionSuccess = false;
4313
4314 if (_functions.isDefined(sDefinition.substr(0, sDefinition.find(":="))))
4315 bDefinitionSuccess = _functions.defineFunc(sDefinition, true);
4316 else
4317 bDefinitionSuccess = _functions.defineFunc(sDefinition);
4318
4319 if (bDefinitionSuccess)
4320 NumeReKernel::print(_lang.get("DEFINE_SUCCESS"), NumeReKernel::getInstance()->getSettings().systemPrints() && !NumeReKernel::bSupressAnswer);
4321 else
4322 NumeReKernel::issueWarning(_lang.get("DEFINE_FAILURE"));
4323
4324 return true;
4325}
4326
4327
4336void rotateTable(CommandLineParser& cmdParser)
4337{
4338 MemoryManager& _data = NumeReKernel::getInstance()->getMemoryManager();
4339 DataAccessParser _accessParser = cmdParser.getExprAsDataObject();
4340
4341 Indices _idx;
4342 double dAlpha = 0.0; // radians
4343
4344 std::vector<double> source_x;
4345 std::vector<double> source_y;
4346
4347 // Find the target of this operation
4348 std::string sTargetTable = cmdParser.getTargetTable(_idx, "rotdata");
4349
4350 auto vParVal = cmdParser.getParameterValueAsNumericalValue("alpha");
4351
4352 if (vParVal.size())
4353 dAlpha = -vParVal.front().real() / 180.0 * M_PI; // deg2rad and change orientation for mathematical positive rotation
4354
4355 _accessParser.getIndices().row.setOpenEndIndex(_data.getLines(_accessParser.getDataObject())-1);
4356 _accessParser.getIndices().col.setOpenEndIndex(_data.getCols(_accessParser.getDataObject())-1);
4357
4358 // Handle the image and datagrid cases
4359 if (cmdParser.getCommand() == "imrot" || cmdParser.getCommand() == "gridrot")
4360 {
4361 if (cmdParser.getCommand() == "gridrot")
4362 {
4363 source_x = mu::real(_data.getElement(_accessParser.getIndices().row, _accessParser.getIndices().col.subidx(0, 1), _accessParser.getDataObject()));
4364 source_y = mu::real(_data.getElement(_accessParser.getIndices().row, _accessParser.getIndices().col.subidx(1, 1), _accessParser.getDataObject()));
4365
4366 // Remove trailing NANs
4367 while (isnan(source_x.back()))
4368 source_x.pop_back();
4369
4370 while (isnan(source_y.back()))
4371 source_y.pop_back();
4372
4373 // Determine the interval range of both axes
4374 // and calculate the relations between size and
4375 // interval ranges
4376 auto source_x_minmax = std::minmax_element(source_x.begin(), source_x.end());
4377 auto source_y_minmax = std::minmax_element(source_y.begin(), source_y.end());
4378 double dSizeRelation = source_x.size() / (double)source_y.size();
4379 double dRangeRelation = fabs(*source_x_minmax.second - *source_x_minmax.first) / fabs(*source_y_minmax.second - *source_y_minmax.first);
4380
4381 // Warn, if the relations are too different
4382 if (fabs(dSizeRelation - dRangeRelation) > 1e-3)
4383 NumeReKernel::issueWarning(_lang.get("BUILTIN_CHECKKEYWORD_ROTATETABLE_WARN_AXES_NOT_PRESERVED"));
4384 }
4385
4386 _accessParser.getIndices().col = _accessParser.getIndices().col.subidx(2);
4387 }
4388
4389 // Extract the range (will pass the ownership)
4390 Memory* _source = _data.getTable(_accessParser.getDataObject())->extractRange(_accessParser.getIndices().row, _accessParser.getIndices().col);
4391
4392 // Remove obsolete entries (needed since new memory model)
4393 _source->shrink();
4394
4395 // Get the edges
4396 Point topLeft(0, 0);
4397 Point topRightS(_source->getCols(false), 0);
4398 Point topRightP(_source->getCols(false)-1, 0);
4399 Point bottomLeftS(0, _source->getLines(false));
4400 Point bottomLeftP(0, _source->getLines(false)-1);
4401 Point bottomRightS(_source->getCols(false), _source->getLines(false));
4402 Point bottomRightP(_source->getCols(false)-1, _source->getLines(false)-1);
4403
4404 // get the rotation origin
4405 Point origin = (bottomRightS + topLeft) / 2.0;
4406
4407 // Calculate their final positions
4408 topLeft.rotate(dAlpha, origin);
4409 topRightS.rotate(dAlpha, origin);
4410 topRightP.rotate(dAlpha, origin);
4411 bottomLeftS.rotate(dAlpha, origin);
4412 bottomLeftP.rotate(dAlpha, origin);
4413 bottomRightS.rotate(dAlpha, origin);
4414 bottomRightP.rotate(dAlpha, origin);
4415
4416 // Calculate the final image extent
4417 double topS = std::min(topLeft.y, std::min(topRightS.y, std::min(bottomLeftS.y, bottomRightS.y)));
4418 double topP = std::min(topLeft.y, std::min(topRightP.y, std::min(bottomLeftP.y, bottomRightP.y)));
4419 double bot = std::max(topLeft.y, std::max(topRightS.y, std::max(bottomLeftS.y, bottomRightS.y)));
4420 double leftS = std::min(topLeft.x, std::min(topRightS.x, std::min(bottomLeftS.x, bottomRightS.x)));
4421 double leftP = std::min(topLeft.x, std::min(topRightP.x, std::min(bottomLeftP.x, bottomRightP.x)));
4422 double right = std::max(topLeft.x, std::max(topRightS.x, std::max(bottomLeftS.x, bottomRightS.x)));
4423
4424 int rows = ceil(bot - topS);
4425 int cols = ceil(right - leftS);
4426
4427 // Insert the axes, if necessary
4428 if (cmdParser.getCommand() == "imrot")
4429 {
4430 // Write the x axis
4431 for (int i = 0; i < rows; i++)
4432 {
4433 if (_idx.row.size() <= (size_t)i)
4434 break;
4435
4436 _data.writeToTable(_idx.row[i], _idx.col[0], sTargetTable, i+1);
4437 }
4438
4439 // Write the y axis
4440 for (int i = 0; i < cols; i++)
4441 {
4442 if (_idx.row.size() <= (size_t)i)
4443 break;
4444
4445 _data.writeToTable(_idx.row[i], _idx.col[1], sTargetTable, i+1);
4446 }
4447
4448 _idx.col = _idx.col.subidx(2);
4449 }
4450 else if (cmdParser.getCommand() == "gridrot")
4451 {
4452 Point origin_axis(interpolateToGrid(source_x, origin.x, true), interpolateToGrid(source_y, origin.y, true));
4453
4454 // Rotate "cell" coordinates to old coord sys,
4455 // interpolate values and rotate back to obtain
4456 // the values of the new coordinates
4457 for (int i = 0; i < rows; i++)
4458 {
4459 if (_idx.row.size() <= (size_t)i)
4460 break;
4461
4462 Point p(i+topP, origin.y);
4463 p.rotate(-dAlpha, origin);
4464
4465 p.x = interpolateToGrid(source_x, p.x, true);
4466 p.y = interpolateToGrid(source_y, p.y, true);
4467
4468 p.rotate(dAlpha, origin_axis);
4469
4470 _data.writeToTable(_idx.row[i], _idx.col[0], sTargetTable, p.x);
4471 }
4472
4473 for (int j = 0; j < cols; j++)
4474 {
4475 if (_idx.row.size() <= (size_t)j)
4476 break;
4477
4478 Point p(origin.x, j+leftP);
4479 p.rotate(-dAlpha, origin);
4480
4481 p.x = interpolateToGrid(source_x, p.x, true);
4482 p.y = interpolateToGrid(source_y, p.y, true);
4483
4484 p.rotate(dAlpha, origin_axis);
4485
4486 _data.writeToTable(_idx.row[j], _idx.col[1], sTargetTable, p.y);
4487 }
4488
4489 _idx.col = _idx.col.subidx(2);
4490 }
4491
4492 Memory* _table = _data.getTable(sTargetTable);
4493
4494 // Prepare the needed number of columns
4495 _table->resizeMemory(-1, _idx.col.subidx(0, cols).max()+1);
4496
4497 // Calculate the rotated grid
4498 #pragma omp parallel for
4499 for (int j = 0; j < cols; j++)
4500 {
4501 if (_idx.col.size() <= (size_t)j)
4502 continue;
4503
4504 for (int i = 0; i < rows; i++)
4505 {
4506 if (_idx.row.size() <= (size_t)i)
4507 break;
4508
4509 // Create a point in rotated source coordinates
4510 // and rotate it backwards
4511 Point p(j + leftP, i + topP);
4512 p.rotate(-dAlpha, origin);
4513
4514 // Store the interpolated value in target coordinates
4515 _table->writeDataDirect(_idx.row[i], _idx.col[j], _source->readMemInterpolated(p.y, p.x));
4516 }
4517 }
4518
4519 // clear the memory instance
4520 delete _source;
4521 _table->markModified();
4522 _data.shrink(sTargetTable);
4523 g_logger.debug("Dataset rotated.");
4524
4525 if (NumeReKernel::getInstance()->getSettings().systemPrints())
4526 NumeReKernel::print(_lang.get("BUILTIN_CHECKKEYWORD_ROTATETABLE_SUCCESS", sTargetTable));
4527}
4528
4529
4530// Forward declaration to make this function usable by the
4531// particle swarm optimizer (needs randomness)
4532value_type parser_Random(const value_type& vRandMin, const value_type& vRandMax);
4533
4534
4547void particleSwarmOptimizer(CommandLineParser& cmdParser)
4548{
4549 Parser& _parser = NumeReKernel::getInstance()->getParser();
4550
4551 size_t nGlobalBest = 0;
4552 size_t nNumParticles = 100;
4553 size_t nMaxIterations = 100;
4554 size_t nDims = 1;
4555
4556 // Extract the interval information
4557 IntervalSet ivl = cmdParser.parseIntervals();
4558
4559 // Handle parameters
4560 std::vector<mu::value_type> vParVal = cmdParser.getParameterValueAsNumericalValue("particles");
4561
4562 if (vParVal.size())
4563 nNumParticles = intCast(vParVal.front());
4564
4565 vParVal = cmdParser.getParameterValueAsNumericalValue("iter");
4566
4567 if (vParVal.size())
4568 nMaxIterations = intCast(vParVal.front());
4569
4570 // Set the expression
4571 _parser.SetExpr(cmdParser.getExprAsMathExpression(true));
4572
4573 // Determine intervals and dimensionality
4574 if (!ivl.size())
4575 ivl.intervals.push_back(Interval(-10.0, 10.0));
4576
4577 nDims = ivl.size();
4578
4579 // Restrict to 4 dimensions, because there are
4580 // only 4 default variables
4581 nDims = std::min(4u, nDims);
4582
4583 // Determine the random range for the velocity vector
4584 double minRange = fabs(ivl[0].max() - ivl[0].min());
4585
4586 for (size_t i = 1; i < nDims; i++)
4587 {
4588 if (fabs(ivl[i].max() - ivl[i].min()) < minRange)
4589 minRange = fabs(ivl[i].max() - ivl[i].min());
4590 }
4591
4592 // The random range is a 10th of the smallest interval
4593 // range to avoid too large jumps in the particle velocity
4594 // compared to the interval range.
4595 double fRandRange = minRange / 10.0;
4596
4597 std::vector<std::vector<double>> vPos;
4598 std::vector<std::vector<double>> vBest;
4599 std::vector<std::vector<double>> vVel;
4600 std::vector<double> vFunctionValues;
4601 std::vector<double> vBestValues;
4602
4603 vPos.resize(nDims);
4604 vBest.resize(nDims);
4605 vVel.resize(nDims);
4606
4607 // Prepare initial vectors
4608 for (size_t i = 0; i < nNumParticles; i++)
4609 {
4610 for (size_t j = 0; j < nDims; j++)
4611 {
4612 vPos[j].push_back(parser_Random(ivl[j].min(), ivl[j].max()).real());
4613 vVel[j].push_back(parser_Random(ivl[j].min(), ivl[j].max()).real()/5.0);
4614
4615 _defVars.vValue[j][0] = vPos[j].back();
4616 }
4617
4618 vFunctionValues.push_back(_parser.Eval().real());
4619 }
4620
4621 for (size_t j = 0; j < nDims; j++)
4622 {
4623 vBest[j] = vPos[j];
4624 }
4625
4626 vBestValues = vFunctionValues;
4627
4628 // Find global best
4629 nGlobalBest = std::min_element(vBestValues.begin(), vBestValues.end()) - vBestValues.begin();
4630
4631 // Iterate
4632 for (size_t i = 0; i < nMaxIterations; i++)
4633 {
4634 // Create an adaptive factor to reduce
4635 // particle position variation over time
4636 double fAdaptiveVelFactor = (nMaxIterations - i) / (double)nMaxIterations;
4637
4638 for (size_t j = 0; j < nNumParticles; j++)
4639 {
4640 for (size_t n = 0; n < nDims; n++)
4641 {
4642 // Update velocities
4643 vVel[n][j] += parser_Random(0, fRandRange).real() * (vBest[n][j] - vPos[n][j]) + parser_Random(0, fRandRange).real() * (vBest[n][nGlobalBest] - vPos[n][j]);
4644
4645 // Update positions
4646 vPos[n][j] += fAdaptiveVelFactor * vVel[n][j];
4647
4648 // Restrict to interval boundaries
4649 vPos[n][j] = std::max(ivl[n].min(), std::min(vPos[n][j], ivl[n].max()));
4650
4651 // Update the corresponding default variable
4652 _defVars.vValue[n][0] = vPos[n][j];
4653 }
4654
4655 // Recalculate the function value
4656 vFunctionValues[j] = _parser.Eval().real();
4657
4658 // Update the best positions
4659 if (vFunctionValues[j] < vBestValues[j])
4660 {
4661 vBestValues[j] = vFunctionValues[j];
4662
4663 for (size_t n = 0; n < nDims; n++)
4664 {
4665 vBest[n][j] = vPos[n][j];
4666 }
4667 }
4668 }
4669
4670 // Update best global position
4671 nGlobalBest = std::min_element(vBestValues.begin(), vBestValues.end()) - vBestValues.begin();
4672 }
4673
4674 // Create return value
4675 std::vector<mu::value_type> vRes;
4676
4677 for (size_t j = 0; j < nDims; j++)
4678 {
4679 vRes.push_back(vBest[j][nGlobalBest]);
4680 }
4681
4682 // Create a temporary vector variable
4683 cmdParser.setReturnValue(vRes);
4684}
4685
4686
4687// Forward declaration for urlExecute
4688std::string removeQuotationMarks(const std::string&);
4689
4690
4700void urlExecute(CommandLineParser& cmdParser)
4701{
4702 // Try to get the contents of the desired URL
4703 try
4704 {
4705 // Get URL, username and password
4706 std::string sUrl = cmdParser.parseExprAsString();
4707 std::string sUserName = cmdParser.getParameterValueAsString("usr", "", true);
4708 std::string sPassword = cmdParser.getParameterValueAsString("pwd", "", true);
4709
4710 // Push the response into a file, if necessary
4711 if (cmdParser.hasParam("file"))
4712 {
4713 std::string sFileName = sUrl.substr(sUrl.rfind("/"));
4714
4715 if (sFileName.find('.') == std::string::npos)
4716 sFileName = "index.html";
4717
4718 if (cmdParser.hasParam("up"))
4719 {
4720 // Get the file parameter value
4721 sFileName = cmdParser.getFileParameterValue(sFileName.substr(sFileName.rfind('.')), "<savepath>", sFileName);
4722
4723 // Upload the file
4724 size_t bytes = url::put(sUrl, sFileName, sUserName, sPassword);
4725 cmdParser.setReturnValue(toString(bytes));
4726 }
4727 else
4728 {
4729 // Get the response from the server
4730 std::string sUrlResponse = url::get(sUrl, sUserName, sPassword);
4731
4732 // Get the file parameter value
4733 sFileName = cmdParser.getFileParameterValueForSaving(sFileName.substr(sFileName.rfind('.')), "<savepath>", sFileName);
4734
4735 // Open the file binary and clean it
4736 std::ofstream file(sFileName, std::ios_base::trunc | std::ios_base::binary);
4737
4738 // Stream the response to the file
4739 if (file.good())
4740 {
4741 file << sUrlResponse;
4742 cmdParser.setReturnValue(toString(sUrlResponse.length()));
4743 }
4744 else
4745 throw SyntaxError(SyntaxError::CANNOT_OPEN_TARGET, cmdParser.getCommandLine(), sFileName, sFileName);
4746 }
4747 }
4748 else
4749 {
4750 // Get the response from the server
4751 std::string sUrlResponse = url::get(sUrl, sUserName, sPassword);
4752
4753 // Replace all masked characters in the return value
4754 replaceAll(sUrlResponse, "\r\n", "\n");
4755// replaceAll(sUrlResponse, "\\", "\\ ");
4756// replaceAll(sUrlResponse, "\"", "\\\"");
4757 cmdParser.setReturnValue(std::vector<std::string>({"\"" + sUrlResponse + "\""}));
4758 }
4759 }
4760 catch (url::Error& e)
4761 {
4762 throw SyntaxError(SyntaxError::URL_ERROR, cmdParser.getCommandLine(), cmdParser.parseExprAsString(), e.what());
4763 }
4764}
4765
4766
4767
Audio::SeekableFile
This class extends the generic audio file with seeking functionalities as well as possibilities for r...
Definition: audiofile.hpp:85
CommandLineParser
This class provides the functionality to extract the different components of a command line into the ...
Definition: commandlineparser.hpp:39
CommandLineParser::parseExprAsNumericalValues
std::vector< mu::value_type > parseExprAsNumericalValues() const
Parses the expression into numerical values, returned as a vector of doubles.
Definition: commandlineparser.cpp:426
CommandLineParser::getExprAsDataObject
DataAccessParser getExprAsDataObject() const
Parses the expression to a DataAccessParser, which will extract the needed information for the curren...
Definition: commandlineparser.cpp:334
CommandLineParser::getFileParameterValue
std::string getFileParameterValue(std::string sFileExt, const std::string &sBaseFolder="", const std::string &sDefaultName="") const
Parses the value of the common "file" command line parameter and returns a valid filename.
Definition: commandlineparser.cpp:558
CommandLineParser::getCommandLine
const std::string & getCommandLine() const
Returns the command line used for constructing this instance (e.g. for errors).
Definition: commandlineparser.hpp:142
CommandLineParser::getExprAsMathExpression
std::string getExprAsMathExpression(bool parseDataObjects=false) const
Prepares the expression by calling custom function definitions and resolving vector braces so that th...
Definition: commandlineparser.cpp:359
CommandLineParser::getParameterValueAsNumericalValue
std::vector< mu::value_type > getParameterValueAsNumericalValue(const std::string &sParameter) const
Parses the selected parameter as (one or more) numerical value(s) and returns them as a vector of dou...
Definition: commandlineparser.cpp:668
CommandLineParser::getParameterValueAsString
std::string getParameterValueAsString(const std::string &sParameter, const std::string &sDefaultValue, bool stripAlways=false, bool onlyStringEvaluation=false) const
Parses the selected parameter value as a string and returns it. If the parameter is not found,...
Definition: commandlineparser.cpp:637
CommandLineParser::getExpr
const std::string & getExpr() const
Returns the expression as plain value.
Definition: commandlineparser.hpp:79
CommandLineParser::getParameterList
const std::string & getParameterList() const
Returns the parameter list.
Definition: commandlineparser.hpp:90
CommandLineParser::getFileParameterValueForSaving
std::string getFileParameterValueForSaving(std::string sFileExt, const std::string &sBaseFolder="", const std::string &sDefaultName="") const
Parses the value of the common "file" command line parameter and returns a valid filename....
Definition: commandlineparser.cpp:616
CommandLineParser::getTargetTable
std::string getTargetTable(Indices &_targetIndices, const std::string &sDefaultTableName)
Evaluates any target=TABLE() statements in the parameter list and returns the needed information....
Definition: commandlineparser.cpp:486
CommandLineParser::setReturnValue
void setReturnValue(const std::string &sRetVal)
Sets the return value of the current command by simply appending it to the return value statement.
Definition: commandlineparser.hpp:116
CommandLineParser::clearReturnValue
void clearReturnValue()
Removes the return value statement.
Definition: commandlineparser.hpp:130
CommandLineParser::getCommand
const std::string & getCommand() const
Returns the command.
Definition: commandlineparser.hpp:68
CommandLineParser::parseExprAsString
std::string parseExprAsString() const
Prepares the expression by handling all string operations and removing the surrounding quotation mark...
Definition: commandlineparser.cpp:392
CommandLineParser::getParameterValue
std::string getParameterValue(const std::string &sParameter) const
Simply returns the parameter value or an empty string. Does not do any parsing steps.
Definition: commandlineparser.cpp:536
CommandLineParser::getAllParametersWithValues
std::vector< std::string > getAllParametersWithValues() const
Returns a vector containing all parameters with values in the current parameter list....
Definition: commandlineparser.cpp:502
CommandLineParser::getExprAsFileName
std::string getExprAsFileName(std::string sFileExt, const std::string &sBasePath="") const
Converts the expression to a file name and removes the surrounding quotation marks,...
Definition: commandlineparser.cpp:290
CommandLineParser::exprContainsDataObjects
bool exprContainsDataObjects() const
Simply returns, whether the expression contains any data objects.
Definition: commandlineparser.cpp:273
CommandLineParser::parseIntervals
IntervalSet parseIntervals(bool bErase=false)
Parses intervals in the parameter list and returns them as a vector. The interval may be deleted,...
Definition: commandlineparser.cpp:469
CommandLineParser::hasParam
bool hasParam(const std::string &sParameter) const
Simple wrapper around findParameter(), if used as a boolean flag.
Definition: commandlineparser.cpp:695
DataAccessParser
This class is defined to abstrahize the determination of the correct data object and the calculation ...
Definition: dataaccess.hpp:38
DataAccessParser::getIndices
Indices & getIndices()
Returns a reference to the stored indices.
Definition: dataaccess.cpp:196
DataAccessParser::evalIndices
void evalIndices()
Evaluates open end indices using the identified data object size.
Definition: dataaccess.cpp:141
DataAccessParser::isCluster
bool isCluster() const
Determines, whether the data access references a cluster.
Definition: dataaccess.cpp:209
DataAccessParser::getDataObject
std::string & getDataObject()
Returns a reference to the data object identifier.
Definition: dataaccess.cpp:170
DetachedLogger::info
void info(const std::string &sMessage)
Convenience member function.
Definition: logger.hpp:106
DetachedLogger::debug
void debug(const std::string &sMessage)
Convenience member function.
Definition: logger.hpp:94
FunctionDefinitionManager
This class implements the function definition managing instance.
Definition: define.hpp:74
FunctionDefinitionManager::defineFunc
bool defineFunc(const std::string &sExpr, bool bRedefine=false, bool bFallback=false)
This function defines a custom function, by passing it to a new FunctionDefinition class instance.
Definition: define.cpp:655
FunctionDefinitionManager::call
bool call(std::string &sExpr, int nRecursion=0)
This function searches for known custom definitions in the passed expression and replaces them with t...
Definition: define.cpp:801
FunctionDefinitionManager::isDefined
bool isDefined(const std::string &sFunc)
This method checks, whether the passed function is already defined.
Definition: define.cpp:625
Interval
This class represents a single interval in code providing reading access functionality.
Definition: interval.hpp:34
Language::get
std::string get(const std::string &sMessage, const std::vector< std::string > &vTokens) const
This member function returns the language string for the passed language identifier and replaces all ...
Definition: language.cpp:292
Memory
This class represents a single table in memory, or a - so to say - single memory page to be handled b...
Definition: memory.hpp:68
Memory::writeDataDirect
void writeDataDirect(int _nLine, int _nCol, const mu::value_type &_dData)
This member function provides an unsafe but direct way of writing data to the table....
Definition: memory.cpp:1249
Memory::readMemInterpolated
mu::value_type readMemInterpolated(double _dLine, double _dCol) const
This member function returns a (bilinearily) interpolated element at the selected double positions.
Definition: memory.cpp:398
Memory::getLines
int getLines(bool _bFull=false) const
This member function will return the number of lines, which are currently available in this table.
Definition: memory.cpp:258
Memory::markModified
void markModified()
Mark this table as modified.
Definition: memory.cpp:1121
Memory::getCols
int getCols(bool _bFull=false) const
This member function will return the number of columns, which are currently available in this table.
Definition: memory.cpp:241
Memory::writeData
void writeData(int _nLine, int _nCol, const mu::value_type &_dData)
This member function writes the passed value to the selected position. The table is automatically enl...
Definition: memory.cpp:1219
Memory::writeDataDirectUnsafe
void writeDataDirectUnsafe(int _nLine, int _nCol, const mu::value_type &_dData)
This member function provides an even more unsafe but direct way of writing data to the table....
Definition: memory.cpp:1271
Memory::resizeMemory
bool resizeMemory(size_t _nLines, size_t _nCols)
This member function will handle all memory grow operations by doubling the base size,...
Definition: memory.cpp:221
Memory::shrink
bool shrink()
This member function shrinks the table memory to the smallest possible dimensions reachable in powers...
Definition: memory.cpp:847
Memory::extractRange
Memory * extractRange(const VectorIndex &_vLine, const VectorIndex &_vCol) const
This member function extracts a range of this table and returns it as a new Memory instance.
Definition: memory.cpp:723
MemoryManager
This class represents the central memory managing instance. It will handle all tables and clusters,...
Definition: memorymanager.hpp:45
MemoryManager::getElement
mu::value_type getElement(int _nLine, int _nCol, const std::string &_sTable) const
Definition: memorymanager.hpp:478
MemoryManager::getAppendedZeroes
int getAppendedZeroes(int _i, const std::string &_sTable) const
Definition: memorymanager.hpp:560
MemoryManager::sortElements
std::vector< int > sortElements(const std::string &sLine)
This member function wraps the sorting functionality and evaluates the passed parameter string before...
Definition: memorymanager.cpp:204
MemoryManager::isValidElement
bool isValidElement(long long int _nLine, long long int _nCol, const std::string &_sTable) const
Definition: memorymanager.hpp:151
MemoryManager::min
std::vector< mu::value_type > min(const std::string &sTable, std::string sDir) const
Definition: memorymanager.hpp:633
MemoryManager::getLines
int getLines(StringView sTable, bool _bFull=false) const
Definition: memorymanager.hpp:422
MemoryManager::shrink
void shrink(const std::string &_sCache)
Definition: memorymanager.hpp:205
MemoryManager::isTable
bool isTable(const std::string &sTable) const
This member function returns, whether the passed table name corresponds to a known table.
Definition: memorymanager.cpp:975
MemoryManager::getHeadLineElement
std::string getHeadLineElement(int _i, const std::string &_sTable) const
Definition: memorymanager.hpp:539
MemoryManager::max
std::vector< mu::value_type > max(const std::string &sTable, std::string sDir) const
Definition: memorymanager.hpp:628
MemoryManager::sum
std::vector< mu::value_type > sum(const std::string &sTable, std::string sDir) const
Definition: memorymanager.hpp:643
MemoryManager::resizeTable
bool resizeTable(int _nCols, const std::string &_sTable)
Definition: memorymanager.hpp:185
MemoryManager::getTable
Memory * getTable(const std::string &sTable)
This member function returns a pointer to an existing Memory instance representing the selected table...
Definition: memorymanager.cpp:826
MemoryManager::containsTablesOrClusters
bool containsTablesOrClusters(const std::string &sCmdLine)
This member function evaluates, whether the passed command line contains tables or clusters.
Definition: memorymanager.cpp:960
MemoryManager::writeToTable
void writeToTable(int _nLine, int _nCol, const std::string &_sCache, const mu::value_type &_dData)
Definition: memorymanager.hpp:578
MemoryManager::setHeadLineElement
bool setHeadLineElement(int _i, const std::string &_sTable, std::string _sHead)
Definition: memorymanager.hpp:598
MemoryManager::getCols
int getCols(StringView sTable, bool _bFull=false) const
Definition: memorymanager.hpp:442
MemoryManager::getColElements
int getColElements(const VectorIndex &cols, const std::string &_sTable) const
Returns the maximal number of elements in the selected column range.
Definition: memorymanager.cpp:618
NumeRe::SavitzkyGolayDiffFilter
This class implements a Savitzky-Golay filter for differentiation providing a derivative up to degree...
Definition: filtering.hpp:692
NumeRe::SavitzkyGolayDiffFilter::apply
virtual mu::value_type apply(size_t i, size_t j, const mu::value_type &val) const override
Override for the abstract apply method of the base class. Applies the filter to the value at the sele...
Definition: filtering.hpp:788
NumeRe::StringParser::isStringExpression
virtual bool isStringExpression(const std::string &sExpression) override
Returns true, if the passed expression is an expression containing strings, string variables or strin...
Definition: stringparser.cpp:2297
NumeRe::Table
This data container is a copy- efficient table to interchange data between Kernel and GUI.
Definition: table.hpp:87
NumeRe::Table::getCols
size_t getCols() const
Get the number of columns.
Definition: table.cpp:656
NumeRe::Table::getLines
size_t getLines() const
Get the number of lines.
Definition: table.cpp:636
NumeRe::Table::getValue
mu::value_type getValue(size_t i, size_t j) const
Getter function for the value of the selected cell.
Definition: table.cpp:561
NumeReKernel
This class provides the interface to the core of NumeRe. It provides all functionalities,...
Definition: kernel.hpp:97
NumeReKernel::getInstance
static NumeReKernel * getInstance()
This static member function returns a a pointer to the singleton instance of the kernel.
Definition: kernel.hpp:221
NumeReKernel::getStringParser
NumeRe::StringParser & getStringParser()
Definition: kernel.hpp:286
NumeReKernel::bSupressAnswer
static bool bSupressAnswer
Definition: kernel.hpp:199
NumeReKernel::printPreFmt
static void printPreFmt(const std::string &__sLine, bool printingEnabled=true)
This member function appends the pre- formatted string to the buffer and informs the terminal that we...
Definition: kernel.cpp:2683
NumeReKernel::getParser
mu::Parser & getParser()
Definition: kernel.hpp:281
NumeReKernel::getMemoryManager
MemoryManager & getMemoryManager()
Definition: kernel.hpp:263
NumeReKernel::print
static void print(const std::string &__sLine, bool printingEnabled=true)
This member function appends the passed string as a new output line to the buffer and informs the ter...
Definition: kernel.cpp:2636
NumeReKernel::getDefinitions
FunctionDefinitionManager & getDefinitions()
Definition: kernel.hpp:291
NumeReKernel::issueWarning
static void issueWarning(std::string sWarningMessage)
This static function may be used to issue a warning to the user. The warning will be printed by the t...
Definition: kernel.cpp:2833
NumeReKernel::GetAsyncCancelState
static bool GetAsyncCancelState()
This function is used by the kernel to get informed, when the user pressed ESC or used other means of...
Definition: kernel.cpp:3703
NumeReKernel::toggleTableStatus
static void toggleTableStatus()
Toggles the table writing status, which will reduce the number or events send to the terminal.
Definition: kernel.cpp:3671
NumeReKernel::getSettings
Settings & getSettings()
Definition: kernel.hpp:296
Settings
This class manages the setting values of the internal (kernel) settings of this application.
Definition: settings.hpp:663
Settings::systemPrints
bool systemPrints() const
Returns, whether system messages shall be printed to the terminal.
Definition: settings.hpp:1140
Settings::getPrecision
size_t getPrecision() const
Returns the precision for displaying floating numbers in the terminal. This value determines the numb...
Definition: settings.hpp:1000
SyntaxError
Common exception class for all exceptions thrown in NumeRe.
Definition: error.hpp:32
SyntaxError::NO_DIFF_OPTIONS
@ NO_DIFF_OPTIONS
Definition: error.hpp:162
SyntaxError::INVALID_INTEGRATION_RANGES
@ INVALID_INTEGRATION_RANGES
Definition: error.hpp:131
SyntaxError::TABLE_DOESNT_EXIST
@ TABLE_DOESNT_EXIST
INSERT HERE.
Definition: error.hpp:211
SyntaxError::NO_EXTREMA_VAR
@ NO_EXTREMA_VAR
Definition: error.hpp:168
SyntaxError::EVAL_VAR_NOT_FOUND
@ EVAL_VAR_NOT_FOUND
Definition: error.hpp:98
SyntaxError::WRONG_DATA_SIZE
@ WRONG_DATA_SIZE
Definition: error.hpp:229
SyntaxError::NO_INTEGRATION_RANGES
@ NO_INTEGRATION_RANGES
Definition: error.hpp:172
SyntaxError::NO_ZEROES_OPTIONS
@ NO_ZEROES_OPTIONS
Definition: error.hpp:182
SyntaxError::NO_ZEROES_VAR
@ NO_ZEROES_VAR
Definition: error.hpp:183
SyntaxError::NO_INTEGRATION_FUNCTION
@ NO_INTEGRATION_FUNCTION
Definition: error.hpp:171
SyntaxError::CANNOT_OPEN_TARGET
@ CANNOT_OPEN_TARGET
Definition: error.hpp:73
SyntaxError::NO_DIFF_VAR
@ NO_DIFF_VAR
Definition: error.hpp:163
SyntaxError::TOO_FEW_DATAPOINTS
@ TOO_FEW_DATAPOINTS
Definition: error.hpp:214
SyntaxError::FUNCTION_ERROR
@ FUNCTION_ERROR
Definition: error.hpp:110
SyntaxError::EXTREMA_VAR_NOT_FOUND
@ EXTREMA_VAR_NOT_FOUND
Definition: error.hpp:100
SyntaxError::INVALID_INTEGRATION_PRECISION
@ INVALID_INTEGRATION_PRECISION
Definition: error.hpp:130
SyntaxError::ZEROES_VAR_NOT_FOUND
@ ZEROES_VAR_NOT_FOUND
INSERT HERE.
Definition: error.hpp:231
SyntaxError::STRINGS_MAY_NOT_BE_EVALUATED_WITH_CMD
@ STRINGS_MAY_NOT_BE_EVALUATED_WITH_CMD
Definition: error.hpp:206
SyntaxError::INVALID_INDEX
@ INVALID_INDEX
Definition: error.hpp:129
SyntaxError::WRONG_PLOT_INTERVAL_FOR_LOGSCALE
@ WRONG_PLOT_INTERVAL_FOR_LOGSCALE
Definition: error.hpp:228
SyntaxError::NO_EXTREMA_OPTIONS
@ NO_EXTREMA_OPTIONS
Definition: error.hpp:167
SyntaxError::PROCESS_ABORTED_BY_USER
@ PROCESS_ABORTED_BY_USER
Definition: error.hpp:202
SyntaxError::invalid_position
static size_t invalid_position
Definition: error.hpp:235
VectorIndex
This class abstracts all the index logics, i.e. the logical differences between single indices and in...
Definition: structures.hpp:42
VectorIndex::subidx
VectorIndex subidx(size_t pos, size_t nLen=std::string::npos) const
This member function returns a subset of the internal stored index just like the std::string::substr(...
Definition: structures.hpp:238
VectorIndex::setOpenEndIndex
void setOpenEndIndex(int nLast) const
This member function can be used to replace the open end state with a defined index value although th...
Definition: structures.hpp:756
VectorIndex::max
int max() const
This function calculates the maximal index value obtained from the values stored internally.
Definition: structures.hpp:558
VectorIndex::size
size_t size() const
This member function returns the size of the indices stored in this class.
Definition: structures.hpp:314
VectorIndex::isOpenEnd
bool isOpenEnd() const
This member function determines, whether the internal index set has an open end.
Definition: structures.hpp:614
VectorIndex::setRange
void setRange(int nMin, int nMax)
This member function can be used to force the indices stored internally to be in a defined interval....
Definition: structures.hpp:712
VectorIndex::to_string
std::string to_string() const
This member function converts the vector indexes contents into a human-readable string representation...
Definition: structures.hpp:770
VectorIndex::min
int min() const
This member function calculates the minimal index value obtained from the values stored internally.
Definition: structures.hpp:576
VectorIndex::front
int & front()
This member function returns a reference to the first index value stored internally.
Definition: structures.hpp:640
mu::ParserBase::GetExpr
const string_type & GetExpr() const
Retrieve the formula.
Definition: muParserBase.cpp:1366
mu::ParserBase::SetExpr
void SetExpr(StringView a_sExpr)
Set the expression. Triggers first time calculation thus the creation of the bytecode and scanning of...
Definition: muParserBase.cpp:616
mu::ParserBase::Eval
value_type Eval()
Single-value wrapper around the vectorized overload of this member function.
Definition: muParserBase.cpp:3073
mu::Parser
Mathematical expressions parser.
Definition: muParser.h:51
mu::Parser::Diff
value_type Diff(value_type *a_Var, value_type a_fPos, value_type a_fEpsilon=0, size_t order=1)
Numerically differentiate with regard to a variable.
Definition: muParser.cpp:383
url::Error
A class for URL exceptions.
Definition: http.h:32
url::Error::what
virtual const char * what() const noexcept
Definition: http.h:38
findZeroes
bool findZeroes(CommandLineParser &cmdParser)
This function is a wrapper to the actual zeros localisation function localizeZero() further below.
Definition: command_implementations.cpp:1859
urlExecute
void urlExecute(CommandLineParser &cmdParser)
This function implements the url command providing an interface to http(s) and (s)ftp URLs.
Definition: command_implementations.cpp:4700
calculateMedian
static double calculateMedian(value_type *v, int nResults, int start, int nOrder, int nanShiftStart, int &nNewNanShift)
Static helper function for findExtremaInData.
Definition: command_implementations.cpp:1272
writeAudioFile
bool writeAudioFile(CommandLineParser &cmdParser)
This function creates a WAVE file from the selected data set.
Definition: command_implementations.cpp:3622
shortTimeFourierAnalysis
bool shortTimeFourierAnalysis(CommandLineParser &cmdParser)
This function performs the short-time fourier analysis on the passed data set.
Definition: command_implementations.cpp:3985
integrateSingleDimensionData
static vector< mu::value_type > integrateSingleDimensionData(CommandLineParser &cmdParser)
This static function integrates single dimension data.
Definition: command_implementations.cpp:162
findExtremaInMultiResult
static bool findExtremaInMultiResult(CommandLineParser &cmdParser, string &sExpr, string &sInterval, int nOrder, int nMode)
This static function finds extrema in a multi-result expression, i.e. an expression containing a tabl...
Definition: command_implementations.cpp:1123
findZeroesInMultiResult
static bool findZeroesInMultiResult(CommandLineParser &cmdParser, string &sExpr, string &sInterval, int nMode)
This static function finds zeroes in a multi-result expression, i.e. an expression containing a table...
Definition: command_implementations.cpp:1666
integrate2d
bool integrate2d(CommandLineParser &cmdParser)
Calculate the integral of a function in two dimensions.
Definition: command_implementations.cpp:481
rotateTable
void rotateTable(CommandLineParser &cmdParser)
This function rotates a table, an image or a datagrid around a specified angle.
Definition: command_implementations.cpp:4336
localizeZero
static mu::value_type localizeZero(string &sCmd, mu::value_type *dVarAdress, Parser &_parser, const Settings &_option, mu::value_type dLeft, mu::value_type dRight, double dEps=1e-10, int nRecursion=0)
This function searches for the positions of all zeroes (roots), which are located in the selected int...
Definition: command_implementations.cpp:2214
analyzePulse
bool analyzePulse(CommandLineParser &cmdParser)
This function performs a pulse analysis on the selected data set.
Definition: command_implementations.cpp:3915
expandVectorToDatagrid
static void expandVectorToDatagrid(IntervalSet &ivl, std::vector< std::vector< mu::value_type > > &vZVals, size_t nSamples_x, size_t nSamples_y)
This function will expand the z vector into a z matrix using triangulation.
Definition: command_implementations.cpp:3549
createDatagrid
bool createDatagrid(CommandLineParser &cmdParser)
This function calculates a datagrid from passed functions or (x-y-z) data values.
Definition: command_implementations.cpp:3312
integrationstep_simpson
static void integrationstep_simpson(mu::value_type &x, mu::value_type x0, mu::value_type x1, double dx, vector< mu::value_type > &vResult, vector< mu::value_type > &vFunctionValues, bool bReturnFunctionPoints)
This static function performs an integration step using the Simpson approximation algorithm.
Definition: command_implementations.cpp:109
readAudioFile
bool readAudioFile(CommandLineParser &cmdParser)
Reads either the audio file meta information or the whole audio file to memory.
Definition: command_implementations.cpp:3690
fastWaveletTransform
bool fastWaveletTransform(CommandLineParser &cmdParser)
This function calculates the fast wavelet transform of the passed data set.
Definition: command_implementations.cpp:3003
removeQuotationMarks
std::string removeQuotationMarks(const std::string &)
This function simply removes the surrounding quotation marks.
Definition: stringfunctions.hpp:79
_defVars
DefaultVariables _defVars
Definition: command_implementations.cpp:39
localizeExtremum
static mu::value_type localizeExtremum(string &sCmd, mu::value_type *dVarAdress, Parser &_parser, const Settings &_option, mu::value_type dLeft, mu::value_type dRight, double dEps=1e-10, int nRecursion=0)
This function searches for the positions of all extrema, which are located in the selected interval.
Definition: command_implementations.cpp:2118
getIntervalForSearchFunctions
static string getIntervalForSearchFunctions(const string &sParams, string &sVar)
This static function extracts the interval definition for the extrema and root search functions.
Definition: command_implementations.cpp:1075
findExtremaInData
static bool findExtremaInData(CommandLineParser &cmdParser, string &sExpr, int nOrder, int nMode)
This static function finds extrema in the selected data sets.
Definition: command_implementations.cpp:1303
interpolateToGrid
static double interpolateToGrid(const std::vector< double > &vAxis, double dInterpolVal, bool bExtent=false)
This static function is a helper for boneDetection to map/interpolate the calculated values to the fi...
Definition: command_implementations.cpp:4094
boneDetection
void boneDetection(CommandLineParser &cmdParser)
This function is the implementation of the detect command.
Definition: command_implementations.cpp:4133
differentiate
bool differentiate(CommandLineParser &cmdParser)
Calculate the numerical differential of the passed expression or data set.
Definition: command_implementations.cpp:754
SIMPSON
#define SIMPSON
Definition: command_implementations.cpp:35
getShiftedAxis
static std::vector< size_t > getShiftedAxis(size_t nElements, bool inverseTrafo)
Calculates an axis index, which performs the necessary data flips used for the shifted fft axis.
Definition: command_implementations.cpp:2515
calculate2dFFT
static void calculate2dFFT(MemoryManager &_data, Indices &_idx, const std::string &sTargetTable, mglDataC &_fftData, std::vector< size_t > &vAxis, FFTData &_fft)
This static function calculates a 2D FFT and stores the result in the target table.
Definition: command_implementations.cpp:2713
integrationstep_trapezoidal
static void integrationstep_trapezoidal(mu::value_type &x, mu::value_type x0, mu::value_type dx, vector< mu::value_type > &vResult, vector< mu::value_type > &vFunctionValues, bool bReturnFunctionPoints)
This static function performs an integration step using a trapezoidal approximation algorithm.
Definition: command_implementations.cpp:60
findExtrema
bool findExtrema(CommandLineParser &cmdParser)
This function is a wrapper to the actual extrema localisation function localizeExtremum() further bel...
Definition: command_implementations.cpp:1398
seekInAudioFile
bool seekInAudioFile(CommandLineParser &cmdParser)
Seek a position in an audiofile and extract a length of samples from it.
Definition: command_implementations.cpp:3779
taylor
void taylor(CommandLineParser &cmdParser)
This function approximates the passed expression using Taylor's method.
Definition: command_implementations.cpp:2305
calculateSplines
bool calculateSplines(CommandLineParser &cmdParser)
This function approximates the passed data set using cubic splines.
Definition: command_implementations.cpp:4253
parser_Random
value_type parser_Random(const value_type &vRandMin, const value_type &vRandMax)
This function returns a uniformly distributed random number between both boundaries.
Definition: functionimplementation.cpp:1834
regularizeDataSet
bool regularizeDataSet(CommandLineParser &cmdParser)
This function regularizes the samples of a defined x-y-data array such that DeltaX is equal for every...
Definition: command_implementations.cpp:3851
calculate1dFFT
static void calculate1dFFT(MemoryManager &_data, Indices &_idx, const std::string &sTargetTable, mglDataC &_fftData, std::vector< size_t > &vAxis, FFTData &_fft)
This static function calculates a 1D FFT and stores the result in the target table.
Definition: command_implementations.cpp:2588
TRAPEZOIDAL
#define TRAPEZOIDAL
Definition: command_implementations.cpp:34
integrate
bool integrate(CommandLineParser &cmdParser)
Calculate the integral of a function or a data set in a single dimension.
Definition: command_implementations.cpp:272
findZeroesInData
static bool findZeroesInData(CommandLineParser &cmdParser, string &sExpr, int nMode)
This static function finds zeroes in the selected data set.
Definition: command_implementations.cpp:1770
particleSwarmOptimizer
void particleSwarmOptimizer(CommandLineParser &cmdParser)
This function implements a particle swarm optimizer in up to four dimensions (depending on the number...
Definition: command_implementations.cpp:4547
evalPoints
bool evalPoints(CommandLineParser &cmdParser)
This function samples a defined expression in an array of discrete values.
Definition: command_implementations.cpp:3175
refreshBoundaries
static void refreshBoundaries(IntervalSet &ivl, const string &sIntegrationExpression)
This static function re-evaluates the boundary expression and updates the internal variables correspo...
Definition: command_implementations.cpp:464
detectPhaseOverflow
static bool detectPhaseOverflow(std::complex< double > cmplx[3])
This static function is a helper function for fastFourierTransform() and detects phase overflows,...
Definition: command_implementations.cpp:2499
fastFourierTransform
bool fastFourierTransform(CommandLineParser &cmdParser)
This function calculates the fast fourier transform of the passed data set.
Definition: command_implementations.cpp:2855
getSamplesForDatagrid
static std::vector< size_t > getSamplesForDatagrid(CommandLineParser &cmdParser, const std::vector< size_t > &nSamples)
This function will obtain the samples of the datagrid for each spatial direction.
Definition: command_implementations.cpp:3492
command_implementations.hpp
getDataElements
string getDataElements(string &sLine, Parser &_parser, MemoryManager &_data, const Settings &_option, int options)
Searches the passed string for calls to any table or cluster and replaces them with internal vectors ...
Definition: dataaccess.cpp:276
isNotEmptyExpression
bool isNotEmptyExpression(const string &sExpr)
This function checks, whether the passed expression is non-empty (i.e. it contains more than white sp...
Definition: dataaccess.cpp:2275
extractRange
Memory * extractRange(const std::string &sCmd, DataAccessParser &_accessParser, int nDesiredCols, bool bSort)
This function extracts a portion of a table and returns it to the calling function....
Definition: dataaccess.cpp:2161
isValidIndexSet
bool isValidIndexSet(const Indices &_idx)
Definition: dataaccess.hpp:81
_lang
Language _lang
Definition: kernel.cpp:39
toSystemCodePage
std::string toSystemCodePage(std::string)
Converts an internal to an external string. Does nothing currently.
Definition: stringtools.cpp:1253
findVariableInExpression
size_t findVariableInExpression(const std::string &sExpr, const std::string &sVarName, size_t nPosStart)
This function searches for the selected variable in the passed string and returns the position of the...
Definition: parser_functions.cpp:42
g_logger
DetachedLogger g_logger
Definition: logger.cpp:23
getMatchingParenthesis
unsigned int getMatchingParenthesis(const StringView &)
Returns the position of the closing parenthesis.
Definition: tools.cpp:414
intPower
std::complex< double > intPower(const std::complex< double > &, int)
This function calculates the power of a value with the specialization that the exponent is an integer...
Definition: tools.cpp:3640
Audio::getAudioFileByType
File * getAudioFileByType(const std::string &sFileName)
Return a audio file type depending on the file extension or a nullptr if the file type is not support...
Definition: audiofile.cpp:34
date::detail::trunc
CONSTCD11 std::enable_if<!std::chrono::treat_as_floating_point< T >::value, T >::type trunc(T t) NOEXCEPT
Definition: date.h:1113
date::abs
CONSTCD11 std::chrono::duration< Rep, Period > abs(std::chrono::duration< Rep, Period > d)
Definition: date.h:1317
date::ceil
CONSTCD14 To ceil(const std::chrono::duration< Rep, Period > &d)
Definition: date.h:1302
mu::value_type
MUP_BASETYPE value_type
The numeric datatype used by the parser.
Definition: muParserDef.h:251
mu::isnan
bool isnan(const value_type &v)
Definition: muParserDef.h:379
mu::real
std::vector< double > real(const std::vector< value_type > &vVec)
Definition: muParserBase.cpp:72
mu::imag
std::vector< double > imag(const std::vector< value_type > &vVec)
Definition: muParserBase.cpp:83
mu::isinf
bool isinf(const value_type &v)
Definition: muParserDef.h:374
mu::rint
value_type rint(value_type v)
Definition: muParserBase.cpp:94
url::put
size_t put(const std::string &sUrl, const std::string &sFileName, const std::string &sUserName, const std::string &sPassWord)
Upload a file to a destination and return the transmitted bytes.
Definition: http.cpp:286
url::get
std::string get(const std::string &sUrl, const std::string &sUserName, const std::string &sPassWord)
Get the contents of a URL.
Definition: http.cpp:251
convertVectorToExpression
void convertVectorToExpression(string &sLine, const Settings &_option)
This function replaces vector expressions with their corresponding multi- expression equation.
Definition: parser_functions.cpp:82
integralFactorial
int integralFactorial(int nNumber)
This function returns the factorial of the passed integral value.
Definition: parser_functions.cpp:535
getPointerToVariable
mu::value_type * getPointerToVariable(const string &sVarName, Parser &_parser)
This function returns the pointer to the passed variable.
Definition: parser_functions.cpp:511
parser_functions.hpp
min
#define min(a, b)
Definition: resampler.cpp:34
M_PI
#define M_PI
Definition: resampler.cpp:47
max
#define max(a, b)
Definition: resampler.cpp:30
StripSpaces
void StripSpaces(std::string &)
Removes leading and trailing white spaces and tabulator characters.
Definition: stringtools.cpp:463
findParameter
int findParameter(const std::string &sCmd, const std::string &sParam, const char cFollowing)
This function searches the passed parameter in the passed command string. If something is found,...
Definition: tools.cpp:113
toUpperCase
std::string toUpperCase(const std::string &sLowerCase)
Converts lowercase letters to uppercase ones.
Definition: stringtools.cpp:932
getNextArgument
std::string getNextArgument(std::string &sArgList, bool bCut)
Definition: tools.cpp:2294
replaceAll
void replaceAll(std::string &sToModify, const char *sToRep, const char *sNewValue, size_t nStart, size_t nEnd)
This function replaces all occurences of the string sToRep in the string sToModify with the new value...
Definition: stringtools.cpp:1443
Audio::Sample
Defines a sample in the audio stream, i.e. a single value per channel.
Definition: audiofile.hpp:33
Audio::Sample::right
double right
Definition: audiofile.hpp:35
Audio::Sample::leftOrMono
double leftOrMono
Definition: audiofile.hpp:34
DefaultVariables
Structure for the four standard variables.
Definition: structures.hpp:2502
DefaultVariables::sName
std::string sName[4]
Definition: structures.hpp:2503
DefaultVariables::vValue
mu::value_type vValue[4][4]
Definition: structures.hpp:2504
FFTData
This structure gathers all information needed for calculating a FFT in one or two dimensions.
Definition: command_implementations.cpp:2563
FFTData::dNyquistFrequency
double dNyquistFrequency[2]
Definition: command_implementations.cpp:2570
Indices
This structure is central for managing the indices of a table or cluster read or write data access....
Definition: structures.hpp:2341
Indices::col
VectorIndex col
Definition: structures.hpp:2343
Indices::row
VectorIndex row
Definition: structures.hpp:2342
IntervalSet
This class represents a set of intervals used together for calculations and simulations.
Definition: interval.hpp:84
IntervalSet::size
size_t size() const
Return the number of intervals.
Definition: interval.cpp:643
IntervalSet::intervals
std::vector< Interval > intervals
Definition: interval.hpp:85
Point
This represents a point in 2D space.
Definition: structures.hpp:2381
Point::rotate
void rotate(double dAlpha, const Point &origin=Point(0, 0))
Definition: structures.hpp:2424
Point::x
double x
Definition: structures.hpp:2382
Point::y
double y
Definition: structures.hpp:2383
intCast
long long int intCast(const std::complex< double > &)
Casts the real part of the complex number to an integer and avoids rounding errors.
Definition: tools.cpp:1824
toString
std::string toString(int)
Converts an integer to a string without the Settings bloat.
Definition: stringtools.cpp:121
findCommand
Match findCommand(StringView sCmd, const std::string &sCommand)
This function is very important for the command handler.
Definition: tools.cpp:1275
containsStrings
bool containsStrings(const string &sLine)
This function checks, whether the passed expression contains strings or valtostring parser.
Definition: tools.cpp:2470
getArgAtPos
string getArgAtPos(const string &sCmd, unsigned int nPos, int extraction)
Extracts a options value at the selected position and applies automatic parsing, if necessary.
Definition: tools.cpp:1598
Linearize
double Linearize(double x_0, double y_0, double x_1, double y_1)
Determines the root of the line connecting the two passed points.
Definition: tools.cpp:2222
LineBreak
std::string LineBreak(std::string sOutput, const Settings &_option, bool bAllowDashBreaks, int nFirstIndent, int nIndent)
This function takes a string, splits it into multiple lines if it is too long and returns the result.
Definition: tools.cpp:2205
getAllIndices
EndlessVector< string > getAllIndices(string sArgList)
Splits up the complete index list and returns them as an EndlessVector.
Definition: tools.cpp:2384
make_hline
void make_hline(int nLength=-1)
This function prints a horizontal line to the terminal using either minus or equal signs.
Definition: kernel.cpp:3720
calculateWavelet
void calculateWavelet(std::vector< double > &data, WaveletType _type, int k, int dir)
Definition: wavelet.cpp:26
decodeWaveletData
NumeRe::Table decodeWaveletData(const std::vector< double > &vWaveletData, const std::vector< double > &vAxisData)
Definition: wavelet.cpp:99
Daubechies
@ Daubechies
Definition: wavelet.hpp:28
CenteredDaubechies
@ CenteredDaubechies
Definition: wavelet.hpp:29
CenteredBSpline
@ CenteredBSpline
Definition: wavelet.hpp:33
BSpline
@ BSpline
Definition: wavelet.hpp:32
CenteredHaar
@ CenteredHaar
Definition: wavelet.hpp:31
Haar
@ Haar
Definition: wavelet.hpp:30