Skip to content

Commit

Permalink
Avoid allocs in Lagrange + cleanup CSPICE comments
Browse files Browse the repository at this point in the history
  • Loading branch information
@ChristopherRabotin
ChristopherRabotin committed Jun 28, 2024
1 parent 0da9ac2 commit d7a2ce3
Show file tree
Hide file tree
Showing 3 changed files with 67 additions and 67 deletions.
122 changes: 60 additions & 62 deletions anise/src/math/interpolation/hermite.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,37 +24,37 @@

/* $ Restrictions */

/* None. */
/* None. */

/* $ Literature_References */

/* [1] "Numerical Recipes---The Art of Scientific Computing" by */
/* William H. Press, Brian P. Flannery, Saul A. Teukolsky, */
/* William T. Vetterling (see sections 3.0 and 3.1). */
/* [1] "Numerical Recipes---The Art of Scientific Computing" by */
/* William H. Press, Brian P. Flannery, Saul A. Teukolsky, */
/* William T. Vetterling (see sections 3.0 and 3.1). */

/* [2] "Elementary Numerical Analysis---An Algorithmic Approach" */
/* by S. D. Conte and Carl de Boor. See p. 64. */
/* [2] "Elementary Numerical Analysis---An Algorithmic Approach" */
/* by S. D. Conte and Carl de Boor. See p. 64. */

/* $ Author_and_Institution */

/* N.J. Bachman (JPL) */
/* N.J. Bachman (JPL) */

/* $ Version */

/* - SPICELIB Version 1.2.1, 28-JAN-2014 (NJB) */

/* Fixed a few comment typos. */
/* Fixed a few comment typos. */

/* - SPICELIB Version 1.2.0, 01-FEB-2002 (NJB) (EDW) */

/* Bug fix: declarations of local variables XI and XIJ */
/* were changed from DOUBLE PRECISION to INTEGER. */
/* Note: bug had no effect on behavior of this routine. */
/* Bug fix: declarations of local variables XI and XIJ */
/* were changed from DOUBLE PRECISION to INTEGER. */
/* Note: bug had no effect on behavior of this routine. */

/* - SPICELIB Version 1.1.0, 28-DEC-2001 (NJB) */

/* Blanks following final newline were truncated to */
/* suppress compilation warnings on the SGI-N32 platform. */
/* Blanks following final newline were truncated to */
/* suppress compilation warnings on the SGI-N32 platform. */

/* - SPICELIB Version 1.0.0, 01-MAR-2000 (NJB) */

Expand Down Expand Up @@ -95,24 +95,24 @@ pub fn hermite_eval(
let work: &mut [f64] = &mut [0.0; 8 * MAX_SAMPLES];
let n: usize = xs.len();

/* Copy the input array into WORK. After this, the first column */
/* of WORK represents the first column of our triangular */
/* interpolation table. */
/* Copy the input array into WORK. After this, the first column */
/* of WORK represents the first column of our triangular */
/* interpolation table. */

for i in 0..n {
work[2 * i] = ys[i];
work[2 * i + 1] = ydots[i];
}

/* Compute the second column of the interpolation table: this */
/* consists of the N-1 values obtained by evaluating the */
/* first-degree interpolants at X. We'll also evaluate the */
/* derivatives of these interpolants at X and save the results in */
/* the second column of WORK. Because the derivative computations */
/* depend on the function computations from the previous column in */
/* the interpolation table, and because the function interpolation */
/* overwrites the previous column of interpolated function values, */
/* we must evaluate the derivatives first. */
/* Compute the second column of the interpolation table: this */
/* consists of the N-1 values obtained by evaluating the */
/* first-degree interpolants at X. We'll also evaluate the */
/* derivatives of these interpolants at X and save the results in */
/* the second column of WORK. Because the derivative computations */
/* depend on the function computations from the previous column in */
/* the interpolation table, and because the function interpolation */
/* overwrites the previous column of interpolated function values, */
/* we must evaluate the derivatives first. */

for i in 1..=n - 1 {
let c1 = xs[i] - x_eval;
Expand All @@ -127,51 +127,50 @@ pub fn hermite_eval(
});
}

/* The second column of WORK contains interpolated derivative */
/* values. */
/* The second column of WORK contains interpolated derivative */
/* values. */

/* The odd-indexed interpolated derivatives are simply the input */
/* derivatives. */
/* The odd-indexed interpolated derivatives are simply the input */
/* derivatives. */

let prev = 2 * i - 1;
let curr = 2 * i;
work[prev + 2 * n - 1] = work[prev];

/* The even-indexed interpolated derivatives are the slopes of */
/* the linear interpolating polynomials for adjacent input */
/* abscissa/ordinate pairs. */
/* The even-indexed interpolated derivatives are the slopes of */
/* the linear interpolating polynomials for adjacent input */
/* abscissa/ordinate pairs. */

work[prev + 2 * n] = (work[curr] - work[prev - 1]) / denom;

/* The first column of WORK contains interpolated function values. */
/* The odd-indexed entries are the linear Taylor polynomials, */
/* for each input abscissa value, evaluated at X. */
/* The first column of WORK contains interpolated function values. */
/* The odd-indexed entries are the linear Taylor polynomials, */
/* for each input abscissa value, evaluated at X. */

let temp = work[prev] * (x_eval - xs[i - 1]) + work[prev - 1];
work[prev] = (c1 * work[prev - 1] + c2 * work[curr]) / denom;
work[prev - 1] = temp;
}

/* The last column entries were not computed by the preceding loop; */
/* compute them now. */
/* The last column entries were not computed by the preceding loop; */
/* compute them now. */

work[4 * n - 2] = work[(2 * n) - 1];
work[2 * (n - 1)] += work[(2 * n) - 1] * (x_eval - xs[n - 1]);

/* Compute columns 3 through 2*N of the table. */
/* Compute columns 3 through 2*N of the table. */

for j in 2..=(2 * n) - 1 {
for i in 1..=(2 * n) - j {
/* In the theoretical construction of the interpolation table,
*/
/* there are 2*N abscissa values, since each input abcissa */
/* value occurs with multiplicity two. In this theoretical */
/* construction, the Jth column of the interpolation table */
/* contains results of evaluating interpolants that span J+1 */
/* consecutive abscissa values. The indices XI and XIJ below */
/* are used to pick the correct abscissa values out of the */
/* physical XVALS array, in which the abscissa values are not */
/* repeated. */
/* In the theoretical construction of the interpolation table, */
/* there are 2*N abscissa values, since each input abcissa */
/* value occurs with multiplicity two. In this theoretical */
/* construction, the Jth column of the interpolation table */
/* contains results of evaluating interpolants that span J+1 */
/* consecutive abscissa values. The indices XI and XIJ below */
/* are used to pick the correct abscissa values out of the */
/* physical XVALS array, in which the abscissa values are not */
/* repeated. */

let xi = (i + 1) / 2;
let xij = (i + j + 1) / 2;
Expand All @@ -186,30 +185,29 @@ pub fn hermite_eval(
});
}

/* Compute the interpolated derivative at X for the Ith */
/* interpolant. This is the derivative with respect to X of */
/* the expression for the interpolated function value, which */
/* is the second expression below. This derivative computation
*/
/* is done first because it relies on the interpolated */
/* function values from the previous column of the */
/* interpolation table. */
/* Compute the interpolated derivative at X for the Ith */
/* interpolant. This is the derivative with respect to X of */
/* the expression for the interpolated function value, which */
/* is the second expression below. This derivative computation */
/* is done first because it relies on the interpolated */
/* function values from the previous column of the */
/* interpolation table. */

/* The derivative expression here corresponds to equation */
/* 2.35 on page 64 in reference [2]. */
/* The derivative expression here corresponds to equation */
/* 2.35 on page 64 in reference [2]. */

work[i + 2 * n - 1] =
(c1 * work[i + 2 * n - 1] + c2 * work[i + 2 * n] + (work[i] - work[i - 1])) / denom;

/* Compute the interpolated function value at X for the Ith */
/* interpolant. */
/* Compute the interpolated function value at X for the Ith */
/* interpolant. */

work[i - 1] = (c1 * work[i - 1] + c2 * work[i]) / denom;
}
}

/* Our interpolated function value is sitting in WORK(1,1) at this */
/* point. The interpolated derivative is located in WORK(1,2). */
/* Our interpolated function value is sitting in WORK(1,1) at this */
/* point. The interpolated derivative is located in WORK(1,2). */

let f = work[0];
let df = work[2 * n];
Expand Down
11 changes: 6 additions & 5 deletions anise/src/math/interpolation/lagrange.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@

use crate::errors::MathError;

use super::InterpolationError;
use super::{InterpolationError, MAX_SAMPLES};

pub fn lagrange_eval(
xs: &[f64],
Expand All @@ -29,8 +29,11 @@ pub fn lagrange_eval(

// At this point, we know that the lengths of items is correct, so we can directly address them without worry for overflowing the array.

let mut work = ys.to_vec();
let mut dwork = vec![0.0; ys.len()];
let work: &mut [f64] = &mut [0.0; MAX_SAMPLES];
let dwork: &mut [f64] = &mut [0.0; MAX_SAMPLES];
for (ii, y) in ys.iter().enumerate() {
work[ii] = *y;
}

let n = xs.len();

Expand All @@ -41,8 +44,6 @@ pub fn lagrange_eval(

let denom = xi - xij;
if denom.abs() < f64::EPSILON {
dbg!(&xs);
dbg!(&ys);
return Err(InterpolationError::InterpMath {
source: MathError::DivisionByZero {
action: "lagrange data contains duplicate states",
Expand Down
1 change: 1 addition & 0 deletions anise/tests/ephemerides/validation/type13_hermite.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,7 @@ fn validate_hermite_type13_from_gmat() {
#[ignore = "Requires Rust SPICE -- must be executed serially"]
#[test]
fn validate_hermite_type13_with_varying_segment_sizes() {
// ISSUE: This file is corrupt, cf. https://github.com/nyx-space/anise/issues/262
let file_name = "spk-type13-validation-variable-seg-size".to_string();
let comparator = CompareEphem::new(
vec!["../data/variable-seg-size-hermite.bsp".to_string()],
Expand Down

0 comments on commit d7a2ce3

Please sign in to comment.