UpsilonNumworks/Upsilon
1
0
Fork 0
mirror of https://github.com/UpsilonNumworks/Upsilon.git synced 2026-01-18 16:27:34 +01:00
Code Issues Packages Projects Releases Wiki Activity

[poincare/zoom] Move automatic zoom to Poincare

Browse Source 
Change-Id: I7a70ab54b33b9150682cc4f999ff4999a8288b7f
This commit is contained in:
Gabriel Ozouf
2020-09-28 15:18:17 +02:00
committed by Émilie Feral
parent 7499a9e0a4
commit d2b36be846
9 changed files with 299 additions and 256 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
apps/shared/continuous_function.cpp
View File

@@ -263,7 +263,7 @@ void ContinuousFunction::setTMax(float tMax) {
void ContinuousFunction::rangeForDisplay(float * xMin, float * xMax, float * yMin, float * yMax, Poincare::Context * context, bool tuneXRange) const {
if (plotType() == PlotType::Cartesian) {
Function::rangeForDisplay(xMin, xMax, yMin, yMax, context, tuneXRange);
protectedRangeForDisplay(xMin, xMax, yMin, yMax, context, tuneXRange, true);
} else {
fullXYRange(xMin, xMax, yMin, yMax, context);
}

1
apps/shared/continuous_function.h
View File

@@ -92,7 +92,6 @@ private:
template <typename T> Poincare::Coordinate2D<T> privateEvaluateXYAtParameter(T t, Poincare::Context * context) const;
void fullXYRange(float * xMin, float * xMax, float * yMin, float * yMax, Poincare::Context * context) const;
void refinedYRangeForDisplay(float xMin, float xMax, float * yMin, float * yMax, Poincare::Context * context) const override { protectedRefinedYRangeForDisplay(xMin, xMax, yMin, yMax, context, true); }
/* RecordDataBuffer is the layout of the data buffer of Record
* representing a ContinuousFunction. See comment on

261
apps/shared/function.cpp
View File

@@ -2,6 +2,7 @@
#include "range_1D.h"
#include "poincare_helpers.h"
#include "poincare/src/parsing/parser.h"
#include <poincare/zoom.h>
#include <ion/display.h>
#include <ion/unicode/utf8_decoder.h>
#include <algorithm>
@@ -78,257 +79,21 @@ Function::RecordDataBuffer * Function::recordData() const {
return reinterpret_cast<RecordDataBuffer *>(const_cast<void *>(d.buffer));
}
// Ranges
static float evaluateAndRound(const Function * f, float x, Context * context, float precision = 1e-5) {
void Function::protectedRangeForDisplay(float * xMin, float * xMax, float * yMin, float * yMax, Poincare::Context * context, bool tuneXRange, bool boundByMagnitude) const {
Zoom::ValueAtAbscissa evaluation = [](float x, Context * context, const void * auxiliary) {
/* When evaluating sin(x)/x close to zero using the standard sine function,
* one can detect small varitions, while the cardinal sine is supposed to be
* one can detect small variations, while the cardinal sine is supposed to be
* locally monotonous. To smooth our such variations, we round the result of
* the evaluations. As we are not interested in precise results but only in
* ordering, this approximation is sufficient. */
return precision * std::round(f->evaluateXYAtParameter(x, context).x2() / precision);
constexpr float precision = 1e-5;
return precision * std::round(static_cast<const Function *>(auxiliary)->evaluateXYAtParameter(x, context).x2() / precision);
};
Zoom::InterestingRangesForDisplay(evaluation, xMin, xMax, yMin, yMax, tMin(), tMax(), context, this, tuneXRange);
evaluation = [](float x, Context * context, const void * auxiliary) {
return static_cast<const Function *>(auxiliary)->evaluateXYAtParameter(x, context).x2();
};
Zoom::RefinedYRangeForDisplay(evaluation, *xMin, *xMax, yMin, yMax, context, this, boundByMagnitude);
}
/* TODO : These three methods perform checks that will also be relevant for the
* equation solver. Remember to factorize this code when integrating the new
* solver. */
static bool boundOfIntervalOfDefinitionIsReached(float y1, float y2) {
return std::isfinite(y1) && !std::isinf(y2) && std::isnan(y2);
}
static bool rootExistsOnInterval(float y1, float y2) {
return ((y1 < 0.f && y2 > 0.f) || (y1 > 0.f && y2 < 0.f));
}
static bool extremumExistsOnInterval(float y1, float y2, float y3) {
return (y1 < y2 && y2 > y3) || (y1 > y2 && y2 < y3);
}
/* This function checks whether an interval contains an extremum or an
* asymptote, by recursively computing the slopes. In case of an extremum, the
* slope should taper off toward the center. */
static bool isExtremum(const Function * f, float x1, float x2, float x3, float y1, float y2, float y3, Context * context, int iterations = 3) {
if (iterations <= 0) {
return false;
}
float x[2] = {x1, x3}, y[2] = {y1, y3};
float xm, ym;
for (int i = 0; i < 2; i++) {
xm = (x[i] + x2) / 2.f;
ym = evaluateAndRound(f, xm, context);
bool res = ((y[i] < ym) != (ym < y2)) ? isExtremum(f, x[i], xm, x2, y[i], ym, y2, context, iterations - 1) : std::fabs(ym - y[i]) >= std::fabs(y2 - ym);
if (!res) {
return false;
}
}
return true;
}
enum class PointOfInterest : uint8_t {
None,
Bound,
Extremum,
Root
};
void Function::rangeForDisplay(float * xMin, float * xMax, float * yMin, float * yMax, Context * context, bool tuneXRange) const {
assert(xMin && xMax && yMin && yMax);
/* Constants of the algorithm. */
constexpr float defaultMaxInterval = 2e5f;
constexpr float minDistance = 1e-2f;
constexpr float asymptoteThreshold = 2e-1f;
constexpr float stepFactor = 1.1f;
constexpr int maxNumberOfPoints = 3;
constexpr float breathingRoom = 0.3f;
constexpr float maxRatioBetweenPoints = 100.f;
const bool hasIntervalOfDefinition = std::isfinite(tMin()) && std::isfinite(tMax());
float center, maxDistance;
if (!tuneXRange) {
center = (*xMax + *xMin) / 2.f;
maxDistance = (*xMax - *xMin) / 2.f;
} else if (hasIntervalOfDefinition) {
center = (tMax() + tMin()) / 2.f;
maxDistance = (tMax() - tMin()) / 2.f;
} else {
center = 0.f;
maxDistance = defaultMaxInterval / 2.f;
}
float resultX[2] = {FLT_MAX, - FLT_MAX};
float resultYMin = FLT_MAX, resultYMax = - FLT_MAX;
float asymptote[2] = {FLT_MAX, - FLT_MAX};
int numberOfPoints;
float xFallback, yFallback[2] = {NAN, NAN};
float firstResult;
float dXOld, dXPrev, dXNext, yOld, yPrev, yNext;
/* Look for a point of interest at the center. */
const float a = center - minDistance - FLT_EPSILON, b = center + FLT_EPSILON, c = center + minDistance + FLT_EPSILON;
const float ya = evaluateAndRound(this, a, context), yb = evaluateAndRound(this, b, context), yc = evaluateAndRound(this, c, context);
if (boundOfIntervalOfDefinitionIsReached(ya, yc) ||
boundOfIntervalOfDefinitionIsReached(yc, ya) ||
rootExistsOnInterval(ya, yc) ||
extremumExistsOnInterval(ya, yb, yc) || ya == yc)
{
resultX[0] = resultX[1] = center;
if (extremumExistsOnInterval(ya, yb, yc) && isExtremum(this, a, b, c, ya, yb, yc, context)) {
resultYMin = resultYMax = yb;
}
}
/* We search for points of interest by exploring the function leftward from
* the center and then rightward, hence the two iterations. */
for (int i = 0; i < 2; i++) {
/* Initialize the search parameters. */
numberOfPoints = 0;
firstResult = NAN;
xFallback = NAN;
dXPrev = i == 0 ? - minDistance : minDistance;
dXNext = dXPrev * stepFactor;
yPrev = evaluateAndRound(this, center + dXPrev, context);
yNext = evaluateAndRound(this, center + dXNext, context);
while(std::fabs(dXPrev) < maxDistance) {
/* Update the slider. */
dXOld = dXPrev;
dXPrev = dXNext;
dXNext *= stepFactor;
yOld = yPrev;
yPrev = yNext;
yNext = evaluateAndRound(this, center + dXNext, context);
if (std::isinf(yNext)) {
continue;
}
/* Check for a change in the profile. */
const PointOfInterest variation = boundOfIntervalOfDefinitionIsReached(yPrev, yNext) ? PointOfInterest::Bound :
rootExistsOnInterval(yPrev, yNext) ? PointOfInterest::Root :
extremumExistsOnInterval(yOld, yPrev, yNext) ? PointOfInterest::Extremum :
PointOfInterest::None;
switch (static_cast<uint8_t>(variation)) {
/* The fallthrough is intentional, as we only want to update the Y
* range when an extremum is detected, but need to update the X range
* in all cases. */
case static_cast<uint8_t>(PointOfInterest::Extremum):
if (isExtremum(this, center + dXOld, center + dXPrev, center + dXNext, yOld, yPrev, yNext, context)) {
resultYMin = std::min(resultYMin, yPrev);
resultYMax = std::max(resultYMax, yPrev);
}
case static_cast<uint8_t>(PointOfInterest::Bound):
/* We only count extrema / discontinuities for limiting the number
* of points. This prevents cos(x) and cos(x)+2 from having different
* profiles. */
if (++numberOfPoints == maxNumberOfPoints) {
/* When too many points are encountered, we prepare their erasure by
* setting a restore point. */
xFallback = dXNext + center;
yFallback[0] = resultYMin;
yFallback[1] = resultYMax;
}
case static_cast<uint8_t>(PointOfInterest::Root):
asymptote[i] = i == 0 ? FLT_MAX : - FLT_MAX;
resultX[0] = std::min(resultX[0], center + (i == 0 ? dXNext : dXPrev));
resultX[1] = std::max(resultX[1], center + (i == 1 ? dXNext : dXPrev));
if (std::isnan(firstResult)) {
firstResult = dXNext;
}
break;
default:
const float slopeNext = (yNext - yPrev) / (dXNext - dXPrev), slopePrev = (yPrev - yOld) / (dXPrev - dXOld);
if ((std::fabs(slopeNext) < asymptoteThreshold) && (std::fabs(slopePrev) > asymptoteThreshold)) {
// Horizontal asymptote begins
asymptote[i] = (i == 0) ? std::min(asymptote[i], center + dXNext) : std::max(asymptote[i], center + dXNext);
} else if ((std::fabs(slopeNext) < asymptoteThreshold) && (std::fabs(slopePrev) > asymptoteThreshold)) {
// Horizontal asymptote invalidates : it might be an asymptote when
// going the other way.
asymptote[1 - i] = (i == 1) ? std::min(asymptote[1 - i], center + dXPrev) : std::max(asymptote[1 - i], center + dXPrev);
}
}
}
if (std::fabs(resultX[i]) > std::fabs(firstResult) * maxRatioBetweenPoints && !std::isnan(xFallback)) {
/* When there are too many points, cut them if their orders are too
* different. */
resultX[i] = xFallback;
resultYMin = yFallback[0];
resultYMax = yFallback[1];
}
}
if (tuneXRange) {
/* Cut after horizontal asymptotes. */
resultX[0] = std::min(resultX[0], asymptote[0]);
resultX[1] = std::max(resultX[1], asymptote[1]);
if (resultX[0] >= resultX[1]) {
/* Fallback to default range. */
resultX[0] = - Range1D::k_default;
resultX[1] = Range1D::k_default;
} else {
/* Add breathing room around points of interest. */
float xRange = resultX[1] - resultX[0];
resultX[0] -= breathingRoom * xRange;
resultX[1] += breathingRoom * xRange;
/* Round to the next integer. */
resultX[0] = std::floor(resultX[0]);
resultX[1] = std::ceil(resultX[1]);
}
*xMin = std::min(resultX[0], *xMin);
*xMax = std::max(resultX[1], *xMax);
}
*yMin = std::min(resultYMin, *yMin);
*yMax = std::max(resultYMax, *yMax);
refinedYRangeForDisplay(*xMin, *xMax, yMin, yMax, context);
}
void Function::protectedRefinedYRangeForDisplay(float xMin, float xMax, float * yMin, float * yMax, Context * context, bool boundByMagnitude) const {
/* This methods computes the Y range that will be displayed for cartesian
* functions and sequences, given an X range (xMin, xMax) and bounds yMin and
* yMax that must be inside the Y range.*/
assert(yMin && yMax);
constexpr int sampleSize = Ion::Display::Width / 4;
constexpr float boundNegligigbleThreshold = 0.2f;
float sampleYMin = FLT_MAX, sampleYMax = -FLT_MAX;
const float step = (xMax - xMin) / (sampleSize - 1);
float x, y;
float sum = 0.f;
int pop = 0;
for (int i = 1; i < sampleSize; i++) {
x = xMin + i * step;
y = evaluateXYAtParameter(x, context).x2();
if (!std::isfinite(y)) {
continue;
}
sampleYMin = std::min(sampleYMin, y);
sampleYMax = std::max(sampleYMax, y);
if (std::fabs(y) > FLT_EPSILON) {
sum += std::log(std::fabs(y));
pop++;
}
}
/* sum/pop is the log mean value of the function, which can be interpreted as
* its average order of magnitude. Then, bound is the value for the next
* order of magnitude and is used to cut the Y range. */
if (boundByMagnitude) {
float bound = (pop > 0) ? std::exp(sum / pop + 1.f) : FLT_MAX;
sampleYMin = std::max(sampleYMin, - bound);
sampleYMax = std::min(sampleYMax, bound);
}
*yMin = std::min(*yMin, sampleYMin);
*yMax = std::max(*yMax, sampleYMax);
if (*yMin == *yMax) {
float d = (*yMin == 0.f) ? 1.f : *yMin * 0.2f;
*yMin -= d;
*yMax += d;
}
/* Round out the smallest bound to 0 if it is negligible compare to the
* other one. This way, we can display the X axis for positive functions such
* as sqrt(x) even if we do not sample close to 0. */
if (*yMin > 0.f && *yMin / *yMax < boundNegligigbleThreshold) {
*yMin = 0.f;
} else if (*yMax < 0.f && *yMax / *yMin < boundNegligigbleThreshold) {
*yMax = 0.f;
}
}
}

6
apps/shared/function.h
View File

@@ -55,7 +55,7 @@ public:
virtual Poincare::Expression sumBetweenBounds(double start, double end, Poincare::Context * context) const = 0;
// Range
virtual void rangeForDisplay(float * xMin, float * xMax, float * yMin, float * yMax, Poincare::Context * context, bool tuneXRange = true) const;
virtual void rangeForDisplay(float * xMin, float * xMax, float * yMin, float * yMax, Poincare::Context * context, bool tuneXRange = true) const = 0;
protected:
/* RecordDataBuffer is the layout of the data buffer of Record
@@ -93,11 +93,9 @@ protected:
bool m_active;
};
void protectedRefinedYRangeForDisplay(float xMin, float xMax, float * yMin, float * yMax, Poincare::Context * context, bool boundByMagnitude) const;
void protectedRangeForDisplay(float * xMin, float * xMax, float * yMin, float * yMax, Poincare::Context * context, bool tuneXRange, bool boundByMagnitude) const;
private:
virtual void refinedYRangeForDisplay(float xMin, float xMax, float * yMin, float * yMax, Poincare::Context * context) const = 0;
RecordDataBuffer * recordData() const;
};

3
apps/shared/range_1D.h
View File

@@ -3,6 +3,7 @@
#include <cmath>
#include <float.h>
#include <poincare/zoom.h>
#if __EMSCRIPTEN__
#include <emscripten.h>
@@ -19,7 +20,7 @@ public:
/* If m_min and m_max are too close, we cannot divide properly the range by
* the number of pixels, which creates a drawing problem. */
constexpr static float k_minFloat = 1E-4f;
constexpr static float k_default = 10.0f;
constexpr static float k_default = Poincare::Zoom::k_defaultHalfRange;
Range1D(float min = -k_default, float max = k_default) :
m_min(min),
m_max(max)

5
apps/shared/sequence.h
View File

@@ -74,6 +74,10 @@ public:
Poincare::Expression sumBetweenBounds(double start, double end, Poincare::Context * context) const override;
constexpr static int k_initialRankNumberOfDigits = 3; // m_initialRank is capped by 999
//Range
void rangeForDisplay(float * xMin, float * xMax, float * yMin, float * yMax, Poincare::Context * context, bool tuneXRange = true) const override { protectedRangeForDisplay(xMin, xMax, yMin, yMax, context, tuneXRange, false); };
private:
constexpr static const KDFont * k_layoutFont = KDFont::LargeFont;
@@ -153,7 +157,6 @@ private:
};
template<typename T> T templatedApproximateAtAbscissa(T x, SequenceContext * sqctx) const;
void refinedYRangeForDisplay(float xMin, float xMax, float * yMin, float * yMax, Poincare::Context * context) const override { protectedRefinedYRangeForDisplay(xMin, xMax, yMin, yMax, context, false); }
size_t metaDataSize() const override { return sizeof(RecordDataBuffer); }
const Shared::ExpressionModel * model() const override { return &m_definition; }
RecordDataBuffer * recordData() const;

1
poincare/Makefile
View File

@@ -156,6 +156,7 @@ poincare_src += $(addprefix poincare/src/,\
vector_cross.cpp \
vector_dot.cpp \
vector_norm.cpp \
zoom.cpp \
)
poincare_src += $(addprefix poincare/src/parsing/,\

50
poincare/include/poincare/zoom.h Normal file
View File

@@ -0,0 +1,50 @@
#ifndef POINCARE_ZOOM_H
#define POINCARE_ZOOM_H
#include <poincare/context.h>
#include <ion/display.h>
namespace Poincare {
class Zoom {
public:
static constexpr float k_defaultHalfRange = 10.f;
typedef float (*ValueAtAbscissa)(float abscissa, Context * context, const void * auxiliary);
static void InterestingRangesForDisplay(ValueAtAbscissa evaluation, float * xMin, float * xMax, float * yMin, float * yMax, float tMin, float tMax, Context * context, const void * auxiliary, bool tuneXRange = false);
static void RefinedYRangeForDisplay(ValueAtAbscissa evaluation, float xMin, float xMax, float * yMin, float * yMax, Context * context, const void * auxiliary, bool boundByMagnitude = false);
private:
static constexpr int k_peakNumberOfPointsOfInterest = 3;
static constexpr int k_sampleSize = Ion::Display::Width / 4;
static constexpr float k_defaultMaxInterval = 2e5f;
static constexpr float k_minimalDistance = 1e-2f;
static constexpr float k_asymptoteThreshold = 2e-1f;
static constexpr float k_stepFactor = 1.1f;
static constexpr float k_breathingRoom = 0.3f;
static constexpr float k_forceXAxisThreshold = 0.2f;
static constexpr float k_maxRatioBetweenPointsOfInterest = 100.f;
enum class PointOfInterest : uint8_t {
None,
Bound,
Extremum,
Root
};
/* TODO : These methods perform checks that will also be relevant for the
* equation solver. Remember to factorize this code when integrating the new
* solver. */
static bool BoundOfIntervalOfDefinitionIsReached(float y1, float y2) { return std::isfinite(y1) && !std::isinf(y2) && std::isnan(y2); }
static bool RootExistsOnInterval(float y1, float y2) { return ((y1 < 0.f && y2 > 0.f) || (y1 > 0.f && y2 < 0.f)); }
static bool ExtremumExistsOnInterval(float y1, float y2, float y3) { return (y1 < y2 && y2 > y3) || (y1 > y2 && y2 < y3); }
/* IsConvexAroundExtremum checks whether an interval contains an extremum or
* an asymptote, by recursively computing the slopes. In case of an extremum,
* the slope should taper off toward the center. */
static bool IsConvexAroundExtremum(ValueAtAbscissa evaluation, float x1, float x2, float x3, float y1, float y2, float y3, Context * context, const void * auxiliary, int iterations = 3);
};
}
#endif

226
poincare/src/zoom.cpp Normal file
View File

@@ -0,0 +1,226 @@
#include <poincare/zoom.h>
#include <float.h>
#include <algorithm>
namespace Poincare {
constexpr int
Zoom::k_peakNumberOfPointsOfInterest,
Zoom::k_sampleSize;
constexpr float
Zoom::k_defaultMaxInterval,
Zoom::k_minimalDistance,
Zoom::k_asymptoteThreshold,
Zoom::k_stepFactor,
Zoom::k_breathingRoom,
Zoom::k_forceXAxisThreshold,
Zoom::k_defaultHalfRange,
Zoom::k_maxRatioBetweenPointsOfInterest;
void Zoom::InterestingRangesForDisplay(ValueAtAbscissa evaluation, float * xMin, float * xMax, float * yMin, float * yMax, float tMin, float tMax, Context * context, const void * auxiliary, bool tuneXRange) {
assert(xMin && xMax && yMin && yMax);
const bool hasIntervalOfDefinition = std::isfinite(tMin) && std::isfinite(tMax);
float center, maxDistance;
if (!tuneXRange) {
center = (*xMax + *xMin) / 2.f;
maxDistance = (*xMax - *xMin) / 2.f;
} else if (hasIntervalOfDefinition) {
center = (tMax + tMin) / 2.f;
maxDistance = (tMax - tMin) / 2.f;
} else {
center = 0.f;
maxDistance = k_defaultMaxInterval / 2.f;
}
float resultX[2] = {FLT_MAX, - FLT_MAX};
float resultYMin = FLT_MAX, resultYMax = - FLT_MAX;
float asymptote[2] = {FLT_MAX, - FLT_MAX};
int numberOfPoints;
float xFallback, yFallback[2] = {NAN, NAN};
float firstResult;
float dXOld, dXPrev, dXNext, yOld, yPrev, yNext;
/* Look for a point of interest at the center. */
const float a = center - k_minimalDistance - FLT_EPSILON, b = center + FLT_EPSILON, c = center + k_minimalDistance + FLT_EPSILON;
const float ya = evaluation(a, context, auxiliary), yb = evaluation(b, context, auxiliary), yc = evaluation(c, context, auxiliary);
if (BoundOfIntervalOfDefinitionIsReached(ya, yc) ||
BoundOfIntervalOfDefinitionIsReached(yc, ya) ||
RootExistsOnInterval(ya, yc) ||
ExtremumExistsOnInterval(ya, yb, yc) || ya == yc)
{
resultX[0] = resultX[1] = center;
if (ExtremumExistsOnInterval(ya, yb, yc) && IsConvexAroundExtremum(evaluation, a, b, c, ya, yb, yc, context, auxiliary)) {
resultYMin = resultYMax = yb;
}
}
/* We search for points of interest by exploring the function leftward from
* the center and then rightward, hence the two iterations. */
for (int i = 0; i < 2; i++) {
/* Initialize the search parameters. */
numberOfPoints = 0;
firstResult = NAN;
xFallback = NAN;
dXPrev = i == 0 ? - k_minimalDistance : k_minimalDistance;
dXNext = dXPrev * k_stepFactor;
yPrev = evaluation(center + dXPrev, context, auxiliary);
yNext = evaluation(center + dXNext, context, auxiliary);
while(std::fabs(dXPrev) < maxDistance) {
/* Update the slider. */
dXOld = dXPrev;
dXPrev = dXNext;
dXNext *= k_stepFactor;
yOld = yPrev;
yPrev = yNext;
yNext = evaluation(center + dXNext, context, auxiliary);
if (std::isinf(yNext)) {
continue;
}
/* Check for a change in the profile. */
const PointOfInterest variation = BoundOfIntervalOfDefinitionIsReached(yPrev, yNext) ? PointOfInterest::Bound :
RootExistsOnInterval(yPrev, yNext) ? PointOfInterest::Root :
ExtremumExistsOnInterval(yOld, yPrev, yNext) ? PointOfInterest::Extremum :
PointOfInterest::None;
switch (static_cast<uint8_t>(variation)) {
/* The fallthrough is intentional, as we only want to update the Y
* range when an extremum is detected, but need to update the X range
* in all cases. */
case static_cast<uint8_t>(PointOfInterest::Extremum):
if (IsConvexAroundExtremum(evaluation, center + dXOld, center + dXPrev, center + dXNext, yOld, yPrev, yNext, context, auxiliary)) {
resultYMin = std::min(resultYMin, yPrev);
resultYMax = std::max(resultYMax, yPrev);
}
case static_cast<uint8_t>(PointOfInterest::Bound):
/* We only count extrema / discontinuities for limiting the number
* of points. This prevents cos(x) and cos(x)+2 from having different
* profiles. */
if (++numberOfPoints == k_peakNumberOfPointsOfInterest) {
/* When too many points are encountered, we prepare their erasure by
* setting a restore point. */
xFallback = dXNext + center;
yFallback[0] = resultYMin;
yFallback[1] = resultYMax;
}
case static_cast<uint8_t>(PointOfInterest::Root):
asymptote[i] = i == 0 ? FLT_MAX : - FLT_MAX;
resultX[0] = std::min(resultX[0], center + (i == 0 ? dXNext : dXPrev));
resultX[1] = std::max(resultX[1], center + (i == 1 ? dXNext : dXPrev));
if (std::isnan(firstResult)) {
firstResult = dXNext;
}
break;
default:
const float slopeNext = (yNext - yPrev) / (dXNext - dXPrev), slopePrev = (yPrev - yOld) / (dXPrev - dXOld);
if ((std::fabs(slopeNext) < k_asymptoteThreshold) && (std::fabs(slopePrev) > k_asymptoteThreshold)) {
// Horizontal asymptote begins
asymptote[i] = (i == 0) ? std::min(asymptote[i], center + dXNext) : std::max(asymptote[i], center + dXNext);
} else if ((std::fabs(slopeNext) < k_asymptoteThreshold) && (std::fabs(slopePrev) > k_asymptoteThreshold)) {
// Horizontal asymptote invalidates : it might be an asymptote when
// going the other way.
asymptote[1 - i] = (i == 1) ? std::min(asymptote[1 - i], center + dXPrev) : std::max(asymptote[1 - i], center + dXPrev);
}
}
}
if (std::fabs(resultX[i]) > std::fabs(firstResult) * k_maxRatioBetweenPointsOfInterest && !std::isnan(xFallback)) {
/* When there are too many points, cut them if their orders are too
* different. */
resultX[i] = xFallback;
resultYMin = yFallback[0];
resultYMax = yFallback[1];
}
}
if (tuneXRange) {
/* Cut after horizontal asymptotes. */
resultX[0] = std::min(resultX[0], asymptote[0]);
resultX[1] = std::max(resultX[1], asymptote[1]);
if (resultX[0] >= resultX[1]) {
/* Fallback to default range. */
resultX[0] = - k_defaultHalfRange;
resultX[1] = k_defaultHalfRange;
} else {
/* Add breathing room around points of interest. */
float xRange = resultX[1] - resultX[0];
resultX[0] -= k_breathingRoom * xRange;
resultX[1] += k_breathingRoom * xRange;
/* Round to the next integer. */
resultX[0] = std::floor(resultX[0]);
resultX[1] = std::ceil(resultX[1]);
}
*xMin = std::min(resultX[0], *xMin);
*xMax = std::max(resultX[1], *xMax);
}
*yMin = std::min(resultYMin, *yMin);
*yMax = std::max(resultYMax, *yMax);
}
void Zoom::RefinedYRangeForDisplay(ValueAtAbscissa evaluation, float xMin, float xMax, float * yMin, float * yMax, Context * context, const void * auxiliary, bool boundByMagnitude) {
/* This methods computes the Y range that will be displayed for cartesian
* functions and sequences, given an X range (xMin, xMax) and bounds yMin and
* yMax that must be inside the Y range.*/
assert(yMin && yMax);
float sampleYMin = FLT_MAX, sampleYMax = -FLT_MAX;
const float step = (xMax - xMin) / (k_sampleSize - 1);
float x, y;
float sum = 0.f;
int pop = 0;
for (int i = 1; i < k_sampleSize; i++) {
x = xMin + i * step;
y = evaluation(x, context, auxiliary);
if (!std::isfinite(y)) {
continue;
}
sampleYMin = std::min(sampleYMin, y);
sampleYMax = std::max(sampleYMax, y);
if (std::fabs(y) > FLT_EPSILON) {
sum += std::log(std::fabs(y));
pop++;
}
}
/* sum/pop is the log mean value of the function, which can be interpreted as
* its average order of magnitude. Then, bound is the value for the next
* order of magnitude and is used to cut the Y range. */
if (boundByMagnitude) {
float bound = (pop > 0) ? std::exp(sum / pop + 1.f) : FLT_MAX;
sampleYMin = std::max(sampleYMin, - bound);
sampleYMax = std::min(sampleYMax, bound);
}
*yMin = std::min(*yMin, sampleYMin);
*yMax = std::max(*yMax, sampleYMax);
if (*yMin == *yMax) {
float d = (*yMin == 0.f) ? 1.f : *yMin * 0.2f;
*yMin -= d;
*yMax += d;
}
/* Round out the smallest bound to 0 if it is negligible compare to the
* other one. This way, we can display the X axis for positive functions such
* as sqrt(x) even if we do not sample close to 0. */
if (*yMin > 0.f && *yMin / *yMax < k_forceXAxisThreshold) {
*yMin = 0.f;
} else if (*yMax < 0.f && *yMax / *yMin < k_forceXAxisThreshold) {
*yMax = 0.f;
}
}
bool Zoom::IsConvexAroundExtremum(ValueAtAbscissa evaluation, float x1, float x2, float x3, float y1, float y2, float y3, Context * context, const void * auxiliary, int iterations) {
if (iterations <= 0) {
return false;
}
float x[2] = {x1, x3}, y[2] = {y1, y3};
float xm, ym;
for (int i = 0; i < 2; i++) {
xm = (x[i] + x2) / 2.f;
ym = evaluation(xm, context, auxiliary);
bool res = ((y[i] < ym) != (ym < y2)) ? IsConvexAroundExtremum(evaluation, x[i], xm, x2, y[i], ym, y2, context, auxiliary, iterations - 1) : std::fabs(ym - y[i]) >= std::fabs(y2 - ym);
if (!res) {
return false;
}
}
return true;
}
}