SinEqLine Class Reference

Helper class to calculate roots for the function f(x) = sin x - slope * x - intercept. More...

#include <SinEqLine.h>

Public Member Functions
	SinEqLine ()
	Default constructor initializing slope and intercept to zero.
	SinEqLine (const Line2D &line2D)
	Constructor taking the line that shall be superimposed with the sin curve.
	SinEqLine (const double slope, const double intercept)
	Constructor taking the slope and intercept of the line that shall be superimposed with the sin curve.
double	map (const double x) const
	Interpreting as the function f this method carries out the translation from x to y coordinates.
double	gradient (const double x) const
	Interpreting as the function f this method calculates the gradient as need in Newtons algorithms.
int	getIHalfPeriod (const double x) const
	Returns the half period index in which the x position is located.
double	computeSmallestPositiveRoot (int maxIHalfPeriod=5) const
double	computeRootLargerThanExtemumInHalfPeriod (int iHalfPeriod) const
	Computes the solution that is addressed by the given half period index.
double	computeRootForLargeSlope () const
	Compute single solution in the case that fabs(slope) >= 1.
double	computeRootInInterval (double lowerX, double upperX) const
	Computes the solution in between the given x values. The x values are generally chosen consecutive local extremes.
double	computeExtremumXInHalfPeriod (int iHalfPeriod) const
	Get the local extremum that is located in the half period indicated by the given index.
bool	hasLargeSlope () const
	Indicates that the slope is so large such that the function has no local extrema.
double	getSlope () const
	Getter for the slope.
double	getIntercept () const
	Getter for the intercept.

Static Public Member Functions
static bool	changesSign (const Vector2D &lower, const Vector2D &upper)
	Checks if the function changes sign in the interval.
static EIncDec	getEIncDec (const Vector2D &lower, const Vector2D &upper)
	Determines if the function is increasing or decreasing in the interval.
static int	getIPeriodFromIHalfPeriod (int iHalfPeriod)
	Helper function to translate the index of the half period to index of the containing period.

Private Member Functions
double	newtonX (const Vector2D &pos) const
	Shrinking method of the newton algorithm return the next candidate root.

Static Private Member Functions
static double	secantX (const Vector2D &lower, const Vector2D &upper)
	Fall back shrinking method to the secant algorithm.
static double	middleX (const Vector2D &lower, const Vector2D &upper)
	Simple fall back shrinking method using trivial division of the interval.
static bool	updateBounds (Vector2D &lower, Vector2D &upper, const Vector2D &next)
	Replaces the lower or upper bound inplace if the next candidate position is valid and within the interval. Returns true on success.
static bool	isBetween (const Vector2D &lower, const Vector2D &next, const Vector2D &upper)
	Check is next position is within the interval given by lower and upper.
static bool	isConverged (const Vector2D &lower, const Vector2D &upper)
	Check if the interval has shrunk close enough to the solution.
static double	getConvergedBound (const Vector2D &lower, const Vector2D &upper)
	Returns the better solution x from the bounds of the interval.

Private Attributes
double	m_slope
	Memory for the slope.
double	m_intercept
	Memory for the intercept.

Detailed Description

Helper class to calculate roots for the function f(x) = sin x - slope * x - intercept.

Solves the equation sin x = slope * x + intercept by Newton's method on the function f(x) = sin x - slope * x - intercept. There can be zero, infinitely many solutions and anything in-between depending on the value of slope and intercept. To find solutions we exploit that the local extrema of the function f(x) can be found easily and that maximally one solutions can be found between consecutive extrema. There is exactly one local extremum in each half period of sin(x) and local extrema are therefore addressed by there half period index. Each possible solution of the equation is than addressed by the index of the closest smaller local extrema. However only a range of indices corresponds to realized solutions. For non existent solutions NAN is returned. Note that for fabs(slope) >= 1 there are no more local maxima. In this case only a single solution exists, which is always returned for any index.

Definition at line 36 of file SinEqLine.h.

Constructor & Destructor Documentation

SinEqLine() [1/3]

SinEqLine ( )

inline

Default constructor initializing slope and intercept to zero.

Definition at line 40 of file SinEqLine.h.

                  :
        m_slope(0.0),
        m_intercept(0.0)
      {}

SinEqLine() [2/3]

SinEqLine ( const Line2D &  line2D )

inlineexplicit

Constructor taking the line that shall be superimposed with the sin curve.

Definition at line 47 of file SinEqLine.h.

                                               :
        m_slope(line2D.slope()),
        m_intercept(line2D.intercept())
      {}

SinEqLine() [3/3]

SinEqLine ( const double  slope, const double  intercept  )

inline

Constructor taking the slope and intercept of the line that shall be superimposed with the sin curve.

Definition at line 53 of file SinEqLine.h.

                                                            :
        m_slope(slope),
        m_intercept(intercept)
      {}

Member Function Documentation

changesSign()

static bool changesSign ( const Vector2D &  lower, const Vector2D &  upper  )

inlinestatic

Checks if the function changes sign in the interval.

Definition at line 129 of file SinEqLine.h.

      { return (lower.y() > 0 and upper.y() < 0) or (lower.y() < 0 and upper.y() > 0); }

computeExtremumXInHalfPeriod()

double computeExtremumXInHalfPeriod

(

int

iHalfPeriod

)

const

Get the local extremum that is located in the half period indicated by the given index.

Definition at line 185 of file SinEqLine.cc.

{
  const double slope = getSlope();
  double extremumInFirstHalfPeriod = acos(slope);
 
  double extremumInFirstPeriod =
    isEven(iHalfPeriod) ? extremumInFirstHalfPeriod : 2 * M_PI - extremumInFirstHalfPeriod;
 
  int iPeriod = getIPeriodFromIHalfPeriod(iHalfPeriod);
  return extremumInFirstPeriod + 2 * M_PI * iPeriod;
}

computeRootForLargeSlope()

double computeRootForLargeSlope

(

)

const

Compute single solution in the case that fabs(slope) >= 1.

Definition at line 50 of file SinEqLine.cc.

{
  if (not hasLargeSlope()) return NAN;
 
  double xLineIsPlusOne = (1.0 - getIntercept()) / getSlope();
  double xLineIsMinusOne = (-1.0 - getIntercept()) / getSlope();
 
  if (xLineIsPlusOne > xLineIsMinusOne) {
    return computeRootInInterval(xLineIsMinusOne, xLineIsPlusOne);
 
  } else if (xLineIsPlusOne < xLineIsMinusOne) {
    return computeRootInInterval(xLineIsPlusOne, xLineIsMinusOne);
  } else {
    return NAN;
  }
}

computeRootInInterval()

double computeRootInInterval	(	double	lowerX,
		double	upperX
	)		const

Computes the solution in between the given x values. The x values are generally chosen consecutive local extremes.

Checks if convergence criterium has been met. For instance if one bound is already exactly at the root.

Checks if the function changes sign in the interval

Definition at line 67 of file SinEqLine.cc.

{
  if (not(lowerX < upperX)) return NAN;
 
  Vector2D lower(lowerX, map(lowerX));
  Vector2D upper(upperX, map(upperX));
 
  if (isConverged(lower, upper)) {
    B2INFO("Coverage early");
    return getConvergedBound(lower, upper);
  }
 
  if (not changesSign(lower, upper)) {
    return NAN;
  }
 
  Vector2D last(lower);
  Vector2D current(upper);
 
  Vector2D next;
  next.setX(secantX(last, current));
  next.setY(map(next.x()));
 
  // Should  always succeed since we checked everything before
  bool updatedBound = updateBounds(lower, upper, next);
  if (not updatedBound) return NAN;
 
  while (not isConverged(lower, upper)) {
    // swap accepted values
    last.set(current);
    current.set(next);
 
    next.setX(newtonX(current));
    next.setY(map(next.x()));
 
    updatedBound = updateBounds(lower, upper, next);
 
    if (not updatedBound) {
 
      // fall back to secant.
      next.setX(secantX(last, current));
      next.setY(map(next.x()));
 
      updatedBound = updateBounds(lower, upper, next);
 
      if (not updatedBound) {
        // fallback to interval division.
        next.setX(middleX(lower, upper));
        next.setY(map(next.x()));
        updatedBound = updateBounds(lower, upper, next);
 
        if (not updatedBound) break;
      }
    }
  }
 
  if (isConverged(lower, upper)) {
    return getConvergedBound(lower, upper);
 
  } else {
    return NAN;
  }
}

computeRootLargerThanExtemumInHalfPeriod()

double computeRootLargerThanExtemumInHalfPeriod

(

int

iHalfPeriod

)

const

Computes the solution that is addressed by the given half period index.

Definition at line 38 of file SinEqLine.cc.

{
  if (hasLargeSlope()) {
    return computeRootForLargeSlope();
 
  } else {
    double lowerX = computeExtremumXInHalfPeriod(iHalfPeriod);
    double upperX = computeExtremumXInHalfPeriod(iHalfPeriod + 1);
    return computeRootInInterval(lowerX, upperX);
  }
}

computeSmallestPositiveRoot()

double computeSmallestPositiveRoot

(

int

maxIHalfPeriod = 5

)

const

Smallest positive root might before the first positive extremum

Most of the time to should be sufficient to look for the root in the first two half periods. So we spell out this case explicitly.

Check if the solution returned is positive, that implies that it is not NAN.

Definition at line 17 of file SinEqLine.cc.

{
  double root = computeRootLargerThanExtemumInHalfPeriod(-1);
  if (root > 0) return root;
 
  root = computeRootLargerThanExtemumInHalfPeriod(0);
  if (root > 0) return root;
 
  for (int iHalfPeriod = 1; iHalfPeriod <= maxIHalfPeriod; ++iHalfPeriod) {
    root = computeRootLargerThanExtemumInHalfPeriod(iHalfPeriod);
 
    if (root > 0) return root;
  }
 
  return NAN;
}

getConvergedBound()

static double getConvergedBound ( const Vector2D &  lower, const Vector2D &  upper  )

inlinestaticprivate

Returns the better solution x from the bounds of the interval.

Definition at line 110 of file SinEqLine.h.

      {
        if (not std::isfinite(lower.y()) or not std::isfinite(upper.y())) {
          return NAN;
        }
 
        if (fabs(lower.y()) <= fabs(upper.y())) {
          return lower.x();
        }
 
        if (fabs(lower.y()) > fabs(upper.y())) {
          return upper.x();
        }
 
        return NAN;
      }

getEIncDec()

static EIncDec getEIncDec ( const Vector2D &  lower, const Vector2D &  upper  )

inlinestatic

Determines if the function is increasing or decreasing in the interval.

Definition at line 133 of file SinEqLine.h.

      {
        if (lower.y() < upper.y()) {
          return EIncDec::c_Increasing;
        } else if (lower.y() > upper.y()) {
          return EIncDec::c_Decreasing;
        } else if (lower.y() == upper.y()) {
          return EIncDec::c_Constant;
        } else {
          return EIncDec::c_Invalid;
        }
      }

getIHalfPeriod()

int getIHalfPeriod ( const double  x ) const

inline

Returns the half period index in which the x position is located.

Definition at line 68 of file SinEqLine.h.

      { return floor(x / M_PI); }

getIntercept()

double getIntercept ( ) const

inline

Getter for the intercept.

Definition at line 164 of file SinEqLine.h.

      { return m_intercept; }

getIPeriodFromIHalfPeriod()

static int getIPeriodFromIHalfPeriod ( int  iHalfPeriod )

inlinestatic

Helper function to translate the index of the half period to index of the containing period.

Definition at line 151 of file SinEqLine.h.

      { return isEven(iHalfPeriod) ? iHalfPeriod / 2 : (iHalfPeriod - 1) / 2; }

getSlope()

double getSlope ( ) const

inline

Getter for the slope.

Definition at line 160 of file SinEqLine.h.

      { return m_slope; }

gradient()

double gradient ( const double  x ) const

inline

Interpreting as the function f this method calculates the gradient as need in Newtons algorithms.

Definition at line 64 of file SinEqLine.h.

      { return cos(x) - getSlope(); }

hasLargeSlope()

bool hasLargeSlope ( ) const

inline

Indicates that the slope is so large such that the function has no local extrema.

Definition at line 156 of file SinEqLine.h.

      { return fabs(getSlope()) >= 1; }

isBetween()

static bool isBetween ( const Vector2D &  lower, const Vector2D &  next, const Vector2D &  upper  )

inlinestaticprivate

Check is next position is within the interval given by lower and upper.

Definition at line 100 of file SinEqLine.h.

      { return lower.x() < next.x() and next.x() < upper.x(); }

isConverged()

static bool isConverged ( const Vector2D &  lower, const Vector2D &  upper  )

inlinestaticprivate

Check if the interval has shrunk close enough to the solution.

Definition at line 104 of file SinEqLine.h.

      {
        return fabs(lower.y()) < 10e-7 or fabs(upper.y()) < 10e-7;
      }

map()

double map ( const double  x ) const

inline

Interpreting as the function f this method carries out the translation from x to y coordinates.

Definition at line 60 of file SinEqLine.h.

      { return sin(x) - getSlope() * x - getIntercept(); }

middleX()

double middleX ( const Vector2D &  lower, const Vector2D &  upper  )

staticprivate

Simple fall back shrinking method using trivial division of the interval.

Definition at line 134 of file SinEqLine.cc.

{
  return (lower.x() + upper.x()) / 2.0;
}

newtonX()

double newtonX ( const Vector2D &  pos ) const

private

Shrinking method of the newton algorithm return the next candidate root.

Definition at line 144 of file SinEqLine.cc.

{
  return pos.x() - pos.y() / gradient(pos.x());
}

secantX()

double secantX ( const Vector2D &  lower, const Vector2D &  upper  )

staticprivate

Fall back shrinking method to the secant algorithm.

Definition at line 139 of file SinEqLine.cc.

{
  return (lower.x() * upper.y() - upper.x() * lower.y()) / (upper.y() - lower.y());
}

updateBounds()

bool updateBounds ( Vector2D &  lower, Vector2D &  upper, const Vector2D &  next  )

staticprivate

Replaces the lower or upper bound inplace if the next candidate position is valid and within the interval. Returns true on success.

Only update if the next point is in-between the lower and upper bound

Definition at line 149 of file SinEqLine.cc.

{
  if (not isBetween(lower, next, upper)) return false;
 
  EIncDec incDecInfo = getEIncDec(lower, upper);
  if (incDecInfo == EIncDec::c_Increasing) {
    if (next.y() > 0.0 and upper.y() > 0.0) {
      upper.set(next);
      return true;
 
    } else if (next.y() <= 0.0 and lower.y() <= 0.0) {
      lower.set(next);
      return true;
 
    } else {
      return false;
    }
 
  } else if (incDecInfo == EIncDec::c_Decreasing) {
    if (next.y() >= 0 and lower.y() >= 0.0) {
      lower.set(next);
      return true;
 
    } else if (next.y() < 0 and upper.y() < 0) {
      upper.set(next);
      return true;
 
    } else {
      return false;
    }
  } else {
    return false;
  }
}

Member Data Documentation

m_intercept

double m_intercept

private

Memory for the intercept.

Definition at line 172 of file SinEqLine.h.

m_slope

double m_slope

private

Memory for the slope.

Definition at line 169 of file SinEqLine.h.

---

The documentation for this class was generated from the following files:

• tracking/trackFindingCDC/numerics/include/SinEqLine.h
• tracking/trackFindingCDC/numerics/src/SinEqLine.cc