turnbull: Use Anderson–Björck method for likelihood-ratio confidence intervals

7.1% speedup compared with Illinois method
36% speedup compared with interval bisection

src/root_finding.rs

@@ -14,7 +14,7 @@
 // You should have received a copy of the GNU Affero General Public License
 // along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
-pub struct IllinoisRootFinder {
+pub struct AndersonBjorckRootFinder {
 	bound_lower: f64,
 	bound_upper: f64,
 	value_lower: f64,
@@ -22,9 +22,9 @@ pub struct IllinoisRootFinder {
 	last_sign: f64  // Sign of the function at last evaluation (1.0 or -1.0) or 0.0 if first iteration
 }
 
-impl IllinoisRootFinder {
+impl AndersonBjorckRootFinder {
-	pub fn new(bound_lower: f64, bound_upper: f64, value_lower: f64, value_upper: f64) -> IllinoisRootFinder {
+	pub fn new(bound_lower: f64, bound_upper: f64, value_lower: f64, value_upper: f64) -> AndersonBjorckRootFinder {
-		return IllinoisRootFinder {
+		return AndersonBjorckRootFinder {
 			bound_lower: bound_lower,
 			bound_upper: bound_upper,
 			value_lower: value_lower,
@@ -33,41 +33,61 @@ impl IllinoisRootFinder {
 		}
 	}
 	
-	pub fn update(&mut self, guess: f64, mut value: f64) {
+	pub fn update(&mut self, guess: f64, value: f64) {
 		if value > 0.0 {
 			if self.value_lower > 0.0 || self.value_upper < 0.0 {
-				self.bound_lower = guess;
-				self.value_lower = value;
-				
 				if self.last_sign == 1.0 {
-					// Illinois adjustment: Halve the y-value of the retained end point (the other end point) when the new y-value has the same sign as the previous one
+					// Anderson–Björck adjustment
+					let m = 1.0 - value / self.value_lower;
+					if m > 0.0 {
+						self.value_upper *= m;
+					} else {
 						self.value_upper *= 0.5;
 					}
-			} else {
+				}
-				self.bound_upper = guess;
-				self.value_upper = value;
 				
+				self.bound_lower = guess;
+				self.value_lower = value;
+			} else {
 				if self.last_sign == 1.0 {
+					let m = 1.0 - value / self.value_upper;
+					if m > 0.0 {
+						self.value_lower *= m;
+					} else {
 						self.value_lower *= 0.5;
 					}
 				}
+				
+				self.bound_upper = guess;
+				self.value_upper = value;
+			}
 			self.last_sign = 1.0;
 		} else {
 			if self.value_lower < 0.0 || self.value_upper > 0.0 {
-				self.bound_lower = guess;
-				self.value_lower = value;
-				
 				if self.last_sign == -1.0 {
+					let m = 1.0 - value / self.value_lower;
+					if m > 0.0 {
+						self.value_upper *= m;
+					} else {
 						self.value_upper *= 0.5;
 					}
-			} else {
+				}
-				self.bound_upper = guess;
-				self.value_upper = value;
 				
+				self.bound_lower = guess;
+				self.value_lower = value;
+			} else {
 				if self.last_sign == -1.0 {
+					let m = 1.0 - value / self.value_upper;
+					if m > 0.0 {
+						self.value_lower *= m;
+					} else {
 						self.value_lower *= 0.5;
 					}
 				}
+				
+				self.bound_upper = guess;
+				self.value_upper = value;
+			}
 			self.last_sign = -1.0;
 		}
 	}

src/turnbull.rs

@@ -29,7 +29,7 @@ use serde::{Serialize, Deserialize};
 
 use crate::csv::read_csv;
 use crate::pava::monotonic_regression_pava;
-use crate::root_finding::IllinoisRootFinder;
+use crate::root_finding::AndersonBjorckRootFinder;
 use crate::term::UnconditionalTermLike;
 
 #[derive(Args)]
@@ -622,7 +622,7 @@ fn compute_hessian(data: &TurnbullData, p: &Vec<f64>) -> DMatrix<f64> {
 
 fn survival_prob_likelihood_ratio_ci(data: &TurnbullData, progress_bar: ProgressBar, max_iterations: u32, ll_tolerance: f64, ci_precision: f64, p: &Vec<f64>, ll_model: f64, s: &Vec<f64>, oim_se: &Vec<f64>, time_index: usize) -> (f64, f64) {
 	// Compute lower confidence limit
-	let mut root_finder = IllinoisRootFinder::new(
+	let mut root_finder = AndersonBjorckRootFinder::new(
 		0.0, s[time_index],
 		f64::NAN, -CHI2_1DF_95  // Value of (lr_statistic - CHI2_1DF_95), which we are seeking the roots of
 	);
@@ -660,7 +660,7 @@ fn survival_prob_likelihood_ratio_ci(data: &TurnbullData, progress_bar: Progress
 	let ci_lower = ci_estimate;
 	
 	// Compute upper confidence limit
-	root_finder = IllinoisRootFinder::new(
+	root_finder = AndersonBjorckRootFinder::new(
 		s[time_index], 1.0,
 		-CHI2_1DF_95, f64::NAN  // Value of (lr_statistic - CHI2_1DF_95), which we are seeking the roots of
 	);