Changeset 2730


Ignore:
Timestamp:
01/16/13 16:00:28 (11 years ago)
Author:
Takeshi Nakazato
Message:

New Development: No

JIRA Issue: Yes CAS-4770

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No?

Module(s): Module Names change impacts.

Description: Describe your changes here...

Rewrite implementations for locator and interpolator.
Documentation (doxygen format) is added to header files.


Location:
trunk/src
Files:
18 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/BisectionLocator.cpp

    r2727 r2730  
    1515
    1616namespace asap {
     17BisectionLocator::BisectionLocator()
     18  : Locator()
     19{
     20}
    1721
    18 BisectionLocator::BisectionLocator(double *v, unsigned int n)
    19   : Locator(v, n)
     22
     23BisectionLocator::BisectionLocator(double *v, unsigned int n, bool copystorage)
     24  : Locator(v, n, copystorage)
    2025{}
    2126
    … …  
    2732  if (n_ == 1)
    2833    return 0;
    29   bool ascending = (x_[n_-1] >= x_[0]);
    30   if (ascending) {
    31     if (x <= x_[0])
    32       return 0;
    33     else if (x > x_[n_-1])
    34       return n_;
    3534
    36     unsigned int jl = 0;
    37     unsigned int ju = n_;
    38     unsigned int jm;
    39     while (ju - jl > 1) {
    40       jm = (ju + jl) >> 1;
    41       if (x > x_[jm])
    42         jl = jm;
    43       else
    44         ju = jm;
    45     }
    46     return ju;
    47   }
    48   else {
    49     if (x >= x_[0])
    50       return 0;
    51     else if (x < x_[n_-1])
    52       return n_;
    53 
    54     unsigned int jl = 0;
    55     unsigned int ju = n_;
    56     unsigned int jm;
    57     while (ju - jl > 1) {
    58       jm = (ju + jl) >> 1;
    59       if (x < x_[jm])
    60         jl = jm;
    61       else
    62         ju = jm;
    63     }
    64     return ju;
    65   }
     35  return bisection(x, 0, n_);
    6636}
    6737

  • trunk/src/BisectionLocator.h

    r2727 r2730  
    2323class BisectionLocator : public Locator {
    2424public:
    25   BisectionLocator() {;}
    26   BisectionLocator(double *v, unsigned int n);
     25  // Default constructor.
     26  BisectionLocator();
    2727
     28  // Construct with data
     29  // @param[in] v pointer to input data.
     30  // @param[in] n length of the data.
     31  // @param[in] copystorage whether allocate internal memory or not.
     32  // @see Locator::set()
     33  BisectionLocator(double *v, unsigned int n, bool copystorage=true);
     34
     35  // Destructor.
    2836  virtual ~BisectionLocator();
    2937
     38  // Return right hand side index of location using bisection search.
     39  // @param[in] x input value to be located.
     40  // @return location as an index j.
     41  // @see Locator::locate()
    3042  unsigned int locate(double x);
    3143};

  • trunk/src/BufferedBisectionLocator.cpp

    r2727 r2730  
    1515
    1616namespace asap {
     17BufferedBisectionLocator::BufferedBisectionLocator()
     18  : Locator(),
     19    prev_(0)
     20{}
    1721
    18 BufferedBisectionLocator::BufferedBisectionLocator(double *v, unsigned int n)
    19   : Locator(v, n),
     22BufferedBisectionLocator::BufferedBisectionLocator(double *v, unsigned int n,
     23                                                   bool copystorage)
     24  : Locator(v, n, copystorage),
    2025    prev_(0)
    2126{}
    … …  
    2833  if (n_ == 1)
    2934    return 0;
    30   bool ascending = (x_[n_-1] >= x_[0]);
    31   if (ascending) {
     35
     36  unsigned int jl = 0;
     37  unsigned int ju = n_;
     38  if (ascending_) {
     39    // ascending order
    3240    if (x <= x_[0])
    3341      return 0;
    … …  
    3543      return n_;
    3644
    37     unsigned int jl = 0;
    38     unsigned int ju = n_;
    39     unsigned int jm;
     45    if (x < x_[prev_]) {
     46      ju = bisection(x, jl, prev_);
     47    }
     48    else {
     49      ju = bisection(x, prev_, ju);
     50    }
    4051
    41     if (x < x_[prev_]) {
    42       ju = prev_;
    43       prev_ = 0;
    44     }
    45     else
    46       jl = prev_;
    47 
    48     while (ju - jl > 1) {
    49       jm = (ju + jl) >> 1;
    50       if (x > x_[jm])
    51         jl = jm;
    52       else
    53         ju = jm;
    54     }
    55     prev_ = jl;
    56     return ju;
    5752  }
    5853  else {
     54    // descending order
    5955    if (x >= x_[0])
    6056      return 0;
    … …  
    6258      return n_;
    6359
    64     unsigned int jl = 0;
    65     unsigned int ju = n_;
    66     unsigned int jm;
     60    if (x > x_[prev_]) {
     61      ju = bisection(x, jl, prev_);
     62    }
     63    else {
     64      ju = bisection(x, prev_, ju);
     65    }
     66  }
    6767
    68     if (x > x_[prev_]) {
    69       ju = prev_;
    70       prev_ = 0;
    71     }
    72     else
    73       jl = prev_;
    74 
    75     while (ju - jl > 1) {
    76       jm = (ju + jl) >> 1;
    77       if (x < x_[jm])
    78         jl = jm;
    79       else
    80         ju = jm;
    81     }
    82     prev_ = jl;
    83     return ju;
    84   }
     68  prev_ = (ju > 0) ? ju - 1 : 0;
     69  return ju;   
    8570}
    8671

  • trunk/src/BufferedBisectionLocator.h

    r2727 r2730  
    1818
    1919/**
    20  * Implementation of locate operation by bisection search
     20 * Implementation of locate operation by bisection search with
     21 * some buffer.
    2122 * @author TakeshiNakazato
    2223 */
    2324class BufferedBisectionLocator : public Locator {
    2425public:
    25   BufferedBisectionLocator() {;}
    26   BufferedBisectionLocator(double *v, unsigned int n);
     26  // Default constructor.
     27  BufferedBisectionLocator();
    2728
     29  // Construct with data
     30  // @param[in] v pointer to input data.
     31  // @param[in] n length of the data.
     32  // @param[in] copystorage whether allocate internal memory or not.
     33  // @see Locator::set()
     34  BufferedBisectionLocator(double *v, unsigned int n, bool copystorage=true);
     35
     36  // Destructor.
    2837  virtual ~BufferedBisectionLocator();
    2938
     39  // Return right hand side index of location using bisection search.
     40  // @param[in] x input value to be located.
     41  // @return location as an index j.
     42  // @see Locator::locate()
    3043  unsigned int locate(double x);
    3144private:
     45
     46  // Previous location index.
    3247  unsigned int prev_;
    3348};

  • trunk/src/BufferedLinearInterpolator1D.cpp

    r2727 r2730  
    2323BufferedLinearInterpolator1D::~BufferedLinearInterpolator1D()
    2424{}
     25
     26void BufferedLinearInterpolator1D::setData(double *x, float *y, unsigned int n)
     27{
     28  Interpolator1D::setData(x, y, n);
     29  reusable_ = false;
     30}
    2531
    2632void BufferedLinearInterpolator1D::setX(double *x, unsigned int n)

  • trunk/src/BufferedLinearInterpolator1D.h

    r2727 r2730  
    1818
    1919/**
    20  * Linear interpolation with some buffers
     20 * Linear interpolation with some buffers for acceleration.
    2121 * @author TakeshiNakazato
    2222 */
    2323class BufferedLinearInterpolator1D : public Interpolator1D {
    2424public:
     25  // Default constructor.
    2526  BufferedLinearInterpolator1D();
    2627
     28  // Destructor.
    2729  virtual ~BufferedLinearInterpolator1D();
    2830
     31  // Set horizontal (x) and vertical (y) data.
     32  // @param[in] x pointer to horizontal data.
     33  // @param[in] y pointer to vertical data.
     34  // @param[in] n number of data.
     35  // @see Interpolator1D::setData()
     36  void setData(double *x, float *y, unsigned int n);
     37
     38  // Set horizontal data (x).
     39  // @param[in] x pointer to horizontal data.
     40  // @param[in] n number of data.
     41  // @see Interpolator1D::setX()
    2942  void setX(double *x, unsigned int n);
     43
     44  // Perform interpolation.
     45  // @param[in] x horizontal location where the value is evaluated
     46  //              by interpolation.
     47  // @return interpolated value at x.
     48  // @see Interpolator1D::interpolate()
    3049  float interpolate(double x);
    3150
    3251private:
     52  // Numerical factor for linear interpolation.
    3353  double factor_;
     54
     55  // Previous location.
    3456  double xold_;
     57
     58  // Previous location as an index
    3559  unsigned int prev_;
     60
     61  // Boolean parameter whether buffered values are effective or not.
    3662  bool reusable_;
    3763};

  • trunk/src/CubicSplineInterpolator1D.cpp

    r2727 r2730  
    3232}
    3333
     34void CubicSplineInterpolator1D::setData(double *x, float *y, unsigned int n)
     35{
     36  Interpolator1D::setData(x, y, n);
     37  reusable_ = false;
     38}
     39
    3440void CubicSplineInterpolator1D::setY(float *y, unsigned int n)
    3541{
    … …  
    5662  // determine second derivative of each point
    5763  if (!reusable_) {
    58     spline();
     64    evaly2();
    5965    reusable_ = true;
    6066  }
    6167
    6268  // cubic spline interpolation
    63   float y = splint(x, i);
     69  float y = dospline(x, i);
    6470  return y;
    6571}
    6672
    67 void CubicSplineInterpolator1D::spline()
     73void CubicSplineInterpolator1D::evaly2()
    6874{
    6975  if (n_ > ny2_) {
    … …  
    7581
    7682  float *u = new float[ny2_-1];
     83
     84  // Natural cubic spline.
    7785  y2_[0] = 0.0;
     86  y2_[ny2_-1] = 0.0;
    7887  u[0] = 0.0;
    7988
     89  // Solve tridiagonal system.
     90  // Here, tridiagonal matrix is decomposed to triangular matrix
     91  // u stores upper triangular components while y2_ stores
     92  // right-hand side vector.
     93  double a1 = x_[1] - x_[0];
     94  double a2, bi;
    8095  for (unsigned int i = 1; i < ny2_ - 1; i++) {
    81     double sig = (x_[i] - x_[i-1]) / (x_[i+1] - x_[i-1]);
    82     double p = sig * y2_[i-1] + 2.0;
    83     y2_[i] = (sig - 1.0) / p;
    84     u[i] = (y_[i+1] - y_[i]) / (x_[i+1] - x_[i])
    85       - (y_[i] - y_[i-1]) / (x_[i] - x_[i-1]);
    86     u[i] = (6.0 * u[i] / (x_[i+1] - x_[i-1]) - sig * u[i-1]) / p;
     96    a2 = x_[i+1] - x_[i];
     97    bi = 1.0 / (x_[i+1] - x_[i-1]);
     98    y2_[i] = 3.0 * bi * ((y_[i+1] - y_[i]) / a2 - (y_[i] - y_[i-1]) / a1
     99                         - y2_[i-1] * 0.5 * a1);
     100    a1 = 1.0 / (1.0 - u[i-1] * 0.5 * a1 * bi);
     101    y2_[i] *= a1;
     102    u[i] = 0.5 * a2 * bi * a1;
     103    a1 = a2;
    87104  }
    88 
    89   y2_[ny2_-1] = 0.0;
    90 
     105 
     106  // Then, solve the system by backsubstitution and store solution
     107  // vector to y2_.
    91108  for (int k = ny2_ - 2; k >= 0; k--)
    92     y2_[k] = y2_[k] * y2_[k+1] + u[k];
     109    y2_[k] -= u[k] * y2_[k+1];
    93110
    94111  delete[] u;
    95112}
    96113
    97 float CubicSplineInterpolator1D::splint(double x, unsigned int i)
     114float CubicSplineInterpolator1D::dospline(double x, unsigned int i)
    98115{
    99116  unsigned int j = i - 1;

  • trunk/src/CubicSplineInterpolator1D.h

    r2727 r2730  
    1818
    1919/**
    20  * Cubic spline interpolation
     20 * Implementation of (natural) cubic spline interpolation.
    2121 * @author TakeshiNakazato
    2222 */
    2323class CubicSplineInterpolator1D : public Interpolator1D {
    2424public:
     25  // Default constructor.
    2526  CubicSplineInterpolator1D();
    2627
     28  // Destructor.
    2729  virtual ~CubicSplineInterpolator1D();
    2830
     31  // Override Interpolator1D::setData.
     32  // @see Interpolator1D::setData
     33  void setData(double *x, float *y, unsigned int n);
     34
     35  // Override Interpolator1D::setY.
     36  // @see Interpolator1D::setY()
    2937  void setY(float *y, unsigned int n);
     38
     39  // Perform interpolation.
     40  // @param[in] x horizontal location where the value is evaluated
     41  //              by interpolation.
     42  // @return interpolated value at x.
    3043  float interpolate(double x);
    3144private:
    32   // determine second derivatives of each point based on
    33   // natural cubic spline condition
    34   void spline();
     45  // Determine second derivatives of each point based on
     46  // natural cubic spline condition (second derivative at each
     47  // end is zero).
     48  void evaly2();
    3549
    36   // do interpolation using second derivatives from spline()
    37   float splint(double x, unsigned int i);
     50  // Do interpolation using second derivatives determined by evaly2().
     51  // @param[in] x horizontal location where the value is evaluated
     52  //              by interpolation.
     53  // @param[in] i location index for x.
     54  // @return interpolated value at x.
     55  float dospline(double x, unsigned int i);
    3856 
     57  // Array to store second derivatives on the data points.
    3958  float *y2_;
     59
     60  // number of data points for second derivatives
    4061  unsigned int ny2_;
     62
     63  // Boolean parameter whether buffered values are effective or not.
    4164  bool reusable_;
    4265};

  • trunk/src/HuntLocator.cpp

    r2727 r2730  
    1616namespace asap {
    1717
    18 HuntLocator::HuntLocator(double *v, unsigned int n)
    19   : Locator(v, n),
     18HuntLocator::HuntLocator()
     19  : Locator(),
     20    prev_(0)
     21{}
     22
     23HuntLocator::HuntLocator(double *v, unsigned int n, bool copystorage)
     24  : Locator(v, n, copystorage),
    2025    prev_(0)
    2126{}
    … …  
    2833  if (n_ == 1)
    2934    return 0;
    30   bool ascending = (x_[n_-1] >= x_[0]);
    31   if (ascending) {
     35
     36  if (ascending_) {
    3237    if (x <= x_[0])
    3338      return 0;
    … …  
    4449  unsigned int jl = 0;
    4550  unsigned int ju = n_;
     51
     52  // hunt phase
    4653  if (prev_ > 0 && prev_ < n_) {
    47     // do hunt
    4854    jl = prev_;
    49     unsigned int inc = 1;
    50     if ((x >= x_[jl]) == ascending) {
    51       // hunt up
    52       if (jl >= n_ - 1)
    53         return jl;
    54       ju = jl + inc;
    55       while ((x >= x_[ju]) == ascending) {
    56         jl = ju;
    57         inc <<= 1;
    58         ju = jl + inc;
    59         if (ju > n_ - 1) {
    60           ju = n_;
    61           break;
    62         }
    63       }
    64     }
    65     else {
    66       // hunt down
    67       if (jl == 0)
    68         return jl;
    69       ju = jl;
    70       jl -= inc;
    71       while ((x < x_[jl]) == ascending) {
    72         ju = jl;
    73         inc <<= 1;
    74         if (inc >= ju) {
    75           jl = 0;
    76           break;
    77         }
    78         else
    79           jl = ju - inc;
    80       }
    81     }
     55    hunt(x, jl, ju);
    8256  }
    8357
    8458  // final bisection phase
    85   unsigned int jm;
    86   if (ascending) {
    87     while (ju - jl > 1) {
    88       jm = (ju + jl) >> 1;
    89       if (x > x_[jm])
    90         jl = jm;
    91       else
    92         ju = jm;
    93     }
    94   }
    95   else {
    96     while (ju - jl > 1) {
    97       jm = (ju + jl) >> 1;
    98       if (x < x_[jm])
    99         jl = jm;
    100       else
    101         ju = jm;
    102     }
    103   }
    104   prev_ = jl;
    105   return ju;
     59  unsigned int j = bisection(x, jl, ju);
     60  prev_ = (j > 0) ? j - 1 : 0;
     61  return j;
    10662}
    10763

  • trunk/src/HuntLocator.h

    r2727 r2730  
    2323class HuntLocator : public Locator {
    2424public:
    25   HuntLocator() {;}
    26   HuntLocator(double *v, unsigned int n);
     25  // Default constructor.
     26  HuntLocator();
    2727
     28  // Construct with data
     29  // @param[in] v pointer to input data.
     30  // @param[in] n length of the data.
     31  // @param[in] copystorage whether allocate internal memory or not.
     32  // @see Locator::set()
     33  HuntLocator(double *v, unsigned int n, bool copystorage=true);
     34
     35  // Destructor.
    2836  virtual ~HuntLocator();
    2937
     38  // Return right hand side index of location using bisection search
     39  // plus hunt algorithm.
     40  // @param[in] x input value to be located.
     41  // @return location as an index j.
     42  // @see Locator::locate()
    3043  unsigned int locate(double x);
    3144private:
     45  // Storage for previous result.
    3246  unsigned int prev_;
    3347};

  • trunk/src/Interpolator1D.cpp

    r2727 r2730  
    7979}
    8080
    81 int Interpolator1D::locate(double x)
     81unsigned int Interpolator1D::locate(double x)
    8282{
    8383  return locator_->locate(x);

  • trunk/src/Interpolator1D.h

    r2727 r2730  
    2323class Interpolator1D {
    2424public:
     25  // Default constructor.
    2526  Interpolator1D();
    2627
     28  // Destructor.
    2729  virtual ~Interpolator1D();
    2830
     31  // Set horizontal (x) and vertical (y) data.
     32  // @param[in] x pointer to horizontal data.
     33  // @param[in] y pointer to vertical data.
     34  // @param[in] n number of data.
    2935  void setData(double *x, float *y, unsigned int n);
     36
     37  // Set horizontal data (x).
     38  // @param[in] x pointer to horizontal data.
     39  // @param[in] n number of data.
    3040  void setX(double *x, unsigned int n);
     41
     42  // Set vertical data (y).
     43  // @param[in] y pointer to vertical data.
     44  // @param[in] n number of data.
    3145  void setY(float *y, unsigned int n);
     46
     47  // Reset object.
    3248  void reset();
    3349
    34   // currently only effective for polynomial interpolation
     50  // Set order for polynomial interpolation.
     51  // @param order order of the polynomial.
     52  //
     53  // This method is effective only for polynomial interpolation.
    3554  void setOrder(unsigned int order) {order_ = order;}
    3655
     56  // Perform interpolation.
     57  // @param[in] x horizontal location where the value is evaluated
     58  //              by interpolation.
     59  // @return interpolated value at x.
    3760  virtual float interpolate(double x) = 0;
    3861
    3962protected:
    40   int locate(double x);
     63  // Locate x.
     64  // @param[in] x horizontal location.
     65  // @return location as an index.
     66  // @see Locator::locate()
     67  unsigned int locate(double x);
     68
     69  // Query function whether the object is ready to interpolate.
     70  // @return true if object is ready else false.
    4171  bool isready();
     72
     73  // Fuctory method to create appropriate Locator object.
    4274  void createLocator();
    4375
     76  // Order of the polynomial (only effective for polynomial interpolation).
    4477  unsigned int order_;
     78
     79  // Number of data points.
    4580  unsigned int n_;
     81
     82  // Horizontal data.
    4683  double *x_;
     84
     85  // Vertical data.
    4786  float *y_;
     87
     88  // Pointer to the Locator object.
    4889  Locator *locator_;
    4990};

  • trunk/src/LinearInterpolator1D.h

    r2727 r2730  
    1818
    1919/**
    20  * Linear interpolation
     20 * Implementation of linear interpolation
    2121 * @author TakeshiNakazato
    2222 */
    2323class LinearInterpolator1D : public Interpolator1D {
    2424public:
     25  // Default constructor.
    2526  LinearInterpolator1D();
    2627
     28  // Destructor.
    2729  virtual ~LinearInterpolator1D();
    2830
     31  // Perform interpolation.
     32  // @param[in] x horizontal location where the value is evaluated
     33  //              by interpolation.
     34  // @return interpolated value at x.
     35  // @see Interpolator1D::interpolate()
    2936  float interpolate(double x);
    3037};

  • trunk/src/Locator.cpp

    r2727 r2730  
    1515
    1616namespace asap {
     17Locator::Locator()
     18  : x_(0),
     19    n_(0),
     20    ascending_(true),
     21    copy_(false)
     22{
     23}
    1724
    18 Locator::Locator(double *v, unsigned int n)
     25Locator::Locator(double *v, unsigned int n, bool copystorage)
     26  : x_(0),
     27    n_(0),
     28    ascending_(true),
     29    copy_(false)
    1930{
    20   set(v, n);
     31  set(v, n, copystorage);
    2132}
    2233
    2334Locator::~Locator()
    24 {}
    25 
    26 void Locator::set(double *v, unsigned int n)
    2735{
    28   x_ = v;
    29   n_ = n;
     36  if (copy_ && x_)
     37    delete[] x_;
    3038}
    3139
     40void Locator::set(double *v, unsigned int n, bool copystorage)
     41{
     42  if (copy_) {
     43    if (!copystorage || n > n_) {
     44      delete[] x_;
     45      x_ = 0;
     46    }
     47  }
     48  copy_ = copystorage;
     49  n_ = n;
     50  if (copy_) {
     51    if (!x_)
     52      x_ = new double[n_];
     53    for (unsigned int i = 0; i < n_; i++)
     54      x_[i] = v[i];
     55  }
     56  else {
     57    x_ = v;
     58  }
     59  ascending_ = (x_[0] <= x_[n_-1]);
    3260}
     61
     62unsigned int Locator::bisection(double x, unsigned int left, unsigned int right)
     63{
     64  unsigned int jl = left;
     65  unsigned int ju = right;
     66
     67  if (ascending_) {
     68    // ascending order
     69    if (x <= x_[0])
     70      return 0;
     71    else if (x > x_[n_-1])
     72      return n_;
     73
     74    unsigned int jm;
     75    while (ju - jl > 1) {
     76      jm = (ju + jl) / 2;
     77      if (x > x_[jm])
     78        jl = jm;
     79      else
     80        ju = jm;
     81    }
     82  }
     83  else {
     84    // descending order
     85    if (x >= x_[0])
     86      return 0;
     87    else if (x < x_[n_-1])
     88      return n_;
     89
     90    unsigned int jm;
     91    while (ju - jl > 1) {
     92      jm = (ju + jl) / 2;
     93      if (x < x_[jm])
     94        jl = jm;
     95      else
     96        ju = jm;
     97    }
     98  }
     99
     100  return ju;
     101}
     102
     103void Locator::hunt(double x, unsigned int &left, unsigned int &right)
     104{
     105  unsigned int inc = 1;
     106  if (ascending_) {
     107    // ascending order
     108    if (x >= x_[left]) {
     109      // forward hunt
     110      if (left >= n_ - 1) {
     111        right = n_;
     112        return;
     113      }
     114      right = left + inc;
     115      while (x >= x_[right]) {
     116        left = right;
     117        inc *= 2;
     118        right = left + inc;
     119        if (right > n_ - 1) {
     120          right = n_;
     121          break;
     122        }
     123      }
     124    }
     125    else {
     126      // backward hunt
     127      if (left == 0) {
     128        right = 0;
     129        return;
     130      }
     131      right = left;
     132      left -= inc;
     133      while (x < x_[left]) {
     134        right = left;
     135        inc *= 2;
     136        if (inc >= right) {
     137          left = 0;
     138          break;
     139        }
     140        left = right - inc;
     141      }
     142    }
     143  }
     144  else {
     145    // descending order
     146    if (x < x_[left]) {
     147      // forward hunt
     148      if (left >= n_ - 1) {
     149        right = n_;
     150        return;
     151      }
     152      right = left + inc;
     153      while (x < x_[right]) {
     154        left = right;
     155        inc *= 2;
     156        right = left + inc;
     157        if (right > n_ - 1) {
     158          right = n_;
     159          break;
     160        }
     161      }
     162    }
     163    else {
     164      // backward hunt
     165      if (left == 0)
     166        return;
     167      right = left;
     168      left -= inc;
     169      while (x >= x_[left]) {
     170        right = left;
     171        inc *= 2;
     172        if (inc >= right) {
     173          left = 0;
     174          break;
     175        }
     176        left = right - inc;
     177      }
     178    }
     179  }
     180}
     181}

  • trunk/src/Locator.h

    r2727 r2730  
    2121class Locator {
    2222public:
    23   Locator() {;}
    24   Locator(double *v, unsigned int n);
    25   void set(double *v, unsigned int n);
     23  // Default constructor.
     24  Locator();
     25 
     26  // Construct with data.
     27  // @param[in] v pointer to input data.
     28  // @param[in] n length of the data.
     29  // @param[in] copystorage whether allocate internal memory or not.
     30  // @see set()
     31  Locator(double *v, unsigned int n, bool copystorage=true);
    2632
     33  // Set data. The data must be sorted in either ascending or descending
     34  // order, and must not have any duplicate elements.
     35  // @param[in] v pointer to input data.
     36  // @param[in] n length of the data.
     37  // @param[in] copystorage whether allocate internal memory or not.
     38  //
     39  // Setting copystorage=false is bit faster since it directly points
     40  // to the input array instead to allocate memory and copy input array.
     41  // However, you have to be careful to set copystorage to false since
     42  // it will take a risk to allow to edit the data to be searched from
     43  // outside the class.
     44  void set(double *v, unsigned int n, bool copystorage=true);
     45
     46  // Destructor.
    2747  virtual ~Locator();
    2848
    29   // return right hand side index of location
    30   // (return j+1 if x[j] < x <= x[j+1])
    31   // return value 0 or x.nelements() indicates out of range
     49  // Return right hand side index of location.
     50  // @param[in] x input value to be located.
     51  // @return location as an index j.
     52  //
     53  // Returned index j satisfies x_[j-1] < x <= x_[j] for ascending
     54  // case while x_[j-1] > x >= x_[j] for descending case.
     55  // Returned value 0 or x.nelements() indicates out of range.
    3256  virtual unsigned int locate(double x) = 0;
    3357
    3458protected:
     59  // Bisection search.
     60  // @param[in] x input value to be located.
     61  // @param[in] left the leftmost index to search.
     62  // @param[in] right the rightmost index to search.
     63  unsigned int bisection(double x, unsigned int left, unsigned int right);
     64 
     65  // Hunt algorithm
     66  // @param[in] x input value to be located.
     67  // @param[in,out] left input: the starting point for hunt.
     68  //                     output: the left index of hunted region.
     69  // @param[out] right the right index of hunted region.
     70  void hunt(double x, unsigned int &left, unsigned int &right);
     71
     72  // Pointer to the data.
    3573  double *x_;
     74
     75  // Length of the data.
    3676  unsigned int n_;
     77
     78  // Is data ascending or descending?
     79  bool ascending_;
     80
     81  // Is internal storage allocated?
     82  bool copy_;
    3783};
    3884

  • trunk/src/NearestInterpolator1D.h

    r2720 r2730  
    2828
    2929/**
    30  * Base class for interpolation operation
     30 * Implementation of nearest interpolation.
    3131 * @author TakeshiNakazato
    3232 */
    3333class NearestInterpolator1D : public Interpolator1D {
    3434public:
     35  // Default constructor.
    3536  NearestInterpolator1D();
    3637
     38  // Destructor.
    3739  virtual ~NearestInterpolator1D();
    3840
     41  // Perform interpolation.
     42  // @param[in] x horizontal location where the value is evaluated
     43  //              by interpolation.
     44  // @return interpolated value at x.
     45  // @see Interpolator1D::interpolate()
    3946  float interpolate(double x);
    4047};

  • trunk/src/PolynomialInterpolator1D.cpp

    r2727 r2730  
    1212#include <assert.h>
    1313#include <math.h>
     14#include <iostream>
     15using namespace std;
     16
     17#include <casa/Exceptions/Error.h>
    1418
    1519#include "PolynomialInterpolator1D.h"
    … …  
    4347  float y;
    4448  if (order_ >= n_ - 1) {
    45     y = polint(x, i, 0, n_);
     49    // use full region
     50    y = dopoly(x, 0, n_);
    4651  }
    4752  else {
     53    // use partial region
    4854    int j = i - 1 - order_ / 2;
    4955    unsigned int m = n_ - 1 - order_;
    5056    unsigned int k = (unsigned int)((j > 0) ? j : 0);
    5157    k = ((k > m) ? m : k);
    52     y = polint(x, i, k, order_ + 1);
     58    y = dopoly(x, k, order_ + 1);
    5359  }
    5460
    … …  
    5662}
    5763
    58 float PolynomialInterpolator1D::polint(double x, unsigned int loc,
    59                                        unsigned int left, unsigned int n)
     64float PolynomialInterpolator1D::dopoly(double x, unsigned int left,
     65                                       unsigned int n)
    6066{
    61   int ns = loc - left;
    62   if (fabs(x - x_[loc]) >= fabs(x - x_[loc-1])) {
    63     ns--;
    64   }
    65 
    6667  double *xa = &x_[left];
    6768  float *ya = &y_[left];
    6869
    69   float y = ya[ns];
    70  
    71   // c stores C11, C21, C31, ..., C[n-1]1
    72   // d stores D11, D21, D31, ..., D[n-1]1
     70  // storage for C and D in Neville's algorithm
    7371  float *c = new float[n];
    7472  float *d = new float[n];
    … …  
    7876  }
    7977
     78  // Neville's algorithm
     79  float y = c[0];
    8080  for (unsigned int m = 1; m < n; m++) {
     81    // Evaluate Cm1, Cm2, Cm3, ... Cm[n-m] and Dm1, Dm2, Dm3, ... Dm[n-m].
     82    // Those are stored to c[0], c[1], ..., c[n-m-1] and d[0], d[1], ...,
     83    // d[n-m-1].
    8184    for (unsigned int i = 0; i < n - m; i++) {
    82       double ho = xa[i] - x;
    83       double hp = xa[i+m] - x;
    84       float w = c[i+1] - d[i];
    85       double denom = ho - hp;
    86       if (denom == 0.0) {
     85      float cd = c[i+1] - d[i];
     86      double dx = xa[i] - xa[i+m];
     87      try {
     88        cd /= dx;
     89      }
     90      catch (...) {
    8791        delete[] c;
    8892        delete[] d;
    89         assert(denom != 0.0);
     93        throw casa::AipsError("x_ has duplicate elements");
    9094      }
    91       denom = w / denom;
    92      
    93       d[i] = hp * denom;
    94       c[i] = ho * denom;
     95      c[i] = (xa[i] - x) * cd;
     96      d[i] = (xa[i+m] - x) * cd;
    9597    }
    9698
    97     float dy;
    98     if (2 * ns < (int)(n - m)) {
    99       dy = c[ns];
    100     }
    101     else {
    102       dy = d[ns-1];
    103       ns--;
    104     }
    105     y += dy;
     99    // In each step, c[0] holds Cm1 which is a correction between
     100    // P12...m and P12...[m+1]. Thus, the following repeated update
     101    // corresponds to the route P1 -> P12 -> P123 ->...-> P123...n.
     102    y += c[0];
    106103  }
    107  
     104
    108105  delete[] c;
    109106  delete[] d;

  • trunk/src/PolynomialInterpolator1D.h

    r2727 r2730  
    1818
    1919/**
    20  * Polynomial interpolation
     20 * Implementation of polynomial interpolation.
    2121 * @author TakeshiNakazato
    2222 */
    2323class PolynomialInterpolator1D : public Interpolator1D {
    2424public:
     25  // Default constructor.
    2526  PolynomialInterpolator1D();
    2627
     28  // Destructor.
    2729  virtual ~PolynomialInterpolator1D();
    2830
     31  // Perform interpolation.
     32  // @param[in] x horizontal location where the value is evaluated
     33  //              by interpolation.
     34  // @return interpolated value at x.
    2935  float interpolate(double x);
    3036private:
    31   float polint(double x, unsigned int loc, unsigned int left, unsigned int n);
     37  // Perform polynomial interpolation.
     38  // If (number of data points) >  (polynomial order + 1), polynomial
     39  // interpolation must be done in the sub-region that contains x.
     40  // This method takes arguments that specifies sub-region to be used.
     41  // @param[in] x horizontal location where the value is evaluated
     42  //              by interpolation.
     43  // @param[in] left the leftmost index of sub-region.
     44  // @param[in] n number of data points of sub-region.
     45  // @return interpolated value at x.
     46  float dopoly(double x, unsigned int left, unsigned int n);
    3247};
    3348
Note: See TracChangeset for help on using the changeset viewer.