From eada2cac7d42a9262381e3f54c2d4aaf4d7fa6c5 Mon Sep 17 00:00:00 2001
From: RunasSudo <runassudo@yingtongli.me>
Date: Tue, 26 Dec 2023 19:28:01 +1100
Subject: [PATCH] =?UTF-8?q?turnbull:=20Use=20Anderson=E2=80=93Bj=C3=B6rck?=
 =?UTF-8?q?=20method=20for=20likelihood-ratio=20confidence=20intervals?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

7.1% speedup compared with Illinois method
36% speedup compared with interval bisection
---
 src/root_finding.rs | 64 +++++++++++++++++++++++++++++----------------
 src/turnbull.rs     |  6 ++---
 2 files changed, 45 insertions(+), 25 deletions(-)

diff --git a/src/root_finding.rs b/src/root_finding.rs
index d712b14..9062eaf 100644
--- a/src/root_finding.rs
+++ b/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 {
-	pub fn new(bound_lower: f64, bound_upper: f64, value_lower: f64, value_upper: f64) -> IllinoisRootFinder {
-		return IllinoisRootFinder {
+impl AndersonBjorckRootFinder {
+	pub fn new(bound_lower: f64, bound_upper: f64, value_lower: f64, value_upper: f64) -> AndersonBjorckRootFinder {
+		return AndersonBjorckRootFinder {
 			bound_lower: bound_lower,
 			bound_upper: bound_upper,
 			value_lower: value_lower,
@@ -33,40 +33,60 @@ 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 {
+				if self.last_sign == 1.0 {
+					// 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;
+					}
+				}
+				
 				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
-					self.value_upper *= 0.5;
-				}
 			} 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;
-				
-				if self.last_sign == 1.0 {
-					self.value_lower *= 0.5;
-				}
 			}
 			self.last_sign = 1.0;
 		} else {
 			if self.value_lower < 0.0 || self.value_upper > 0.0 {
+				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;
+					}
+				}
+				
 				self.bound_lower = guess;
 				self.value_lower = value;
-				
-				if self.last_sign == -1.0 {
-					self.value_upper *= 0.5;
-				}
 			} 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;
-				
-				if self.last_sign == -1.0 {
-					self.value_lower *= 0.5;
-				}
 			}
 			self.last_sign = -1.0;
 		}
diff --git a/src/turnbull.rs b/src/turnbull.rs
index 3c32b61..ace0568 100644
--- a/src/turnbull.rs
+++ b/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
 	);