Merge branch 'main' of github.com:Anexen/pyxirr into feature/ffi

CHANGELOG.md

@@ -1,5 +1,9 @@
 # Changelog
 
+## not released
+
+- XIRR Performance improvements
+
 ## [0.10.2] - 2024-01-27
 
 - (X)IRR Performance improvements

README.md

@@ -32,9 +32,9 @@ and the [implementation from the Stack Overflow](https://stackoverflow.com/a/115
 
 ![bench](https://raw.githubusercontent.com/Anexen/pyxirr/main/docs/static/bench.png)
 
-PyXIRR is ~10-20x faster in XIRR calculation than the other implementations.
+PyXIRR is much faster than the other implementations.
 
-Powered by [github-action-benchmark](https://github.com/rhysd/github-action-benchmark) and [plotly.js](https://github.com/plotly/plotly.js).
+Powered by [github-action-benchmark](https://github.com/benchmark-action/github-action-benchmark) and [plotly.js](https://github.com/plotly/plotly.js).
 
 Live benchmarks are hosted on [Github Pages](https://anexen.github.io/pyxirr/bench).
 
@@ -61,16 +61,16 @@ xirr(['2020-01-01', '2021-01-01'], [-1000, 1200])
 
 # Multiple IRR problem
 
-The Multiple IRR problem occur when the signs of cash flows change more than
+The Multiple IRR problem occurs when the signs of cash flows change more than
 once. In this case, we say that the project has non-conventional cash flows.
 This leads to situation, where it can have more the one IRR or have no IRR at all.
 
-PyXIRR's approach to the Multiple IRR problem:
+PyXIRR addresses the Multiple IRR problem as follows:
 
 1. It looks for positive result around 0.1 (the same as Excel with the default guess=0.1).
 2. If it can't find a result, it uses several other attempts and selects the lowest IRR to be conservative.
 
-Here is an example of how to find multiple IRRs:
+Here is an example illustrating how to identify multiple IRRs:
 
 ```python
 import numpy as np

bench-requirements.txt

@@ -1,3 +1,2 @@
 scipy==1.*
 numpy-financial==1.*
-xirr==0.1.6

benches/comparison.rs

@@ -10,7 +10,7 @@ mod common;
 use common::PaymentsLoader;
 
 // https://stackoverflow.com/questions/8919718/financial-python-library-that-has-xirr-and-xnpv-function
-const TOP_STACK_OVERFLOW_ANSWER: &str = r#"
+const PURE_PYTHON_IMPL: &str = r#"
 def xirr(transactions):
     years = [(ta[0] - transactions[0][0]).days / 365.0 for ta in transactions]
     residual = 1
@@ -32,13 +32,29 @@ def xirr(transactions):
     return guess
 "#;
 
+const SCIPY_IMPL: &str = r#"
+import scipy.optimize
+
+def xnpv(rate, values, dates):
+    if rate <= -1.0:
+        return float('inf')
+    d0 = dates[0]    # or min(dates)
+    return sum([ vi / (1.0 + rate)**((di - d0).days / 365.0) for vi, di in zip(values, dates)])
+
+def xirr(values, dates):
+    try:
+        return scipy.optimize.newton(lambda r: xnpv(r, values, dates), 0.0)
+    except RuntimeError:    # Failed to converge?
+        return scipy.optimize.brentq(lambda r: xnpv(r, values, dates), -1.0, 1e10)
+"#;
+
 macro_rules! bench_rust {
     ($name:ident, $file:expr) => {
         #[bench]
         fn $name(b: &mut Bencher) {
             Python::with_gil(|py| {
                 let data = PaymentsLoader::from_csv(py, $file).to_records();
-                b.iter(|| pyxirr_call_impl!(py, "xirr", (data,)).unwrap());
+                b.iter(|| pyxirr_call_impl!(py, "xirr", black_box((data,))).unwrap());
             });
         }
     };
@@ -49,10 +65,12 @@ macro_rules! bench_scipy {
         #[bench]
         fn $name(b: &mut Bencher) {
             Python::with_gil(|py| {
-                let xirr =
-                    py.import("xirr").expect("xirr is not installed").getattr("xirr").unwrap();
-                let data = PaymentsLoader::from_csv(py, $file).to_dict();
-                b.iter(|| xirr.call1(black_box((data,))).unwrap())
+                let xirr = PyModule::from_code(py, SCIPY_IMPL, "xirr.py", "scipy_py_xirr")
+                    .unwrap()
+                    .getattr("xirr")
+                    .unwrap();
+                let data = PaymentsLoader::from_csv(py, $file).to_columns();
+                b.iter(|| xirr.call1(black_box((data.1, data.0))).unwrap())
             });
         }
     };
@@ -63,26 +81,28 @@ macro_rules! bench_python {
         #[bench]
         fn $name(b: &mut Bencher) {
             Python::with_gil(|py| {
-                let xirr =
-                    PyModule::from_code(py, TOP_STACK_OVERFLOW_ANSWER, "xirr.py", "pure_py_xirr")
-                        .unwrap()
-                        .getattr("xirr")
-                        .unwrap();
+                let xirr = PyModule::from_code(py, PURE_PYTHON_IMPL, "xirr.py", "pure_py_xirr")
+                    .unwrap()
+                    .getattr("xirr")
+                    .unwrap();
                 let data = PaymentsLoader::from_csv(py, $file).to_records();
                 b.iter(|| xirr.call1(black_box((data,))).unwrap())
             });
         }
     };
 }
 
-bench_rust!(bench_rust_100, "tests/samples/random_100.csv");
-bench_rust!(bench_rust_500, "tests/samples/random_500.csv");
-bench_rust!(bench_rust_1000, "tests/samples/random_1000.csv");
+bench_rust!(bench_rust_50, "tests/samples/rw-50.csv");
+bench_rust!(bench_rust_100, "tests/samples/rw-100.csv");
+bench_rust!(bench_rust_500, "tests/samples/rw-500.csv");
+bench_rust!(bench_rust_1000, "tests/samples/rw-1000.csv");
 
-bench_scipy!(bench_scipy_100, "tests/samples/random_100.csv");
-bench_scipy!(bench_scipy_500, "tests/samples/random_500.csv");
-bench_scipy!(bench_scipy_1000, "tests/samples/random_1000.csv");
+bench_scipy!(bench_scipy_50, "tests/samples/rw-50.csv");
+bench_scipy!(bench_scipy_100, "tests/samples/rw-100.csv");
+bench_scipy!(bench_scipy_500, "tests/samples/rw-500.csv");
+bench_scipy!(bench_scipy_1000, "tests/samples/rw-1000.csv");
 
-bench_python!(bench_python_100, "tests/samples/random_100.csv");
-bench_python!(bench_python_500, "tests/samples/random_500.csv");
-bench_python!(bench_python_1000, "tests/samples/random_1000.csv");
+bench_python!(bench_python_50, "tests/samples/rw-50.csv");
+bench_python!(bench_python_100, "tests/samples/rw-100.csv");
+bench_python!(bench_python_500, "tests/samples/rw-500.csv");
+bench_python!(bench_python_1000, "tests/samples/rw-1000.csv");

core/src/optimize.rs

@@ -1,5 +1,6 @@
 const MAX_ERROR: f64 = 1e-9;
 const MAX_ITERATIONS: u32 = 50;
+const MAX_FX_TOL: f64 = 1e-3;
 
 pub fn newton_raphson<Func, Deriv>(start: f64, f: &Func, d: &Deriv) -> f64
 where
@@ -11,13 +12,13 @@ where
     let mut x = start;
 
     for _ in 0..MAX_ITERATIONS {
-        let res = f(x);
+        let y = f(x);
 
-        if res.abs() < MAX_ERROR {
+        if y.abs() < MAX_ERROR {
             return x;
         }
 
-        let delta = res / d(x);
+        let delta = y / d(x);
 
         if delta.abs() < MAX_ERROR {
             return x - delta;
@@ -29,6 +30,35 @@ where
     f64::NAN
 }
 
+// a slightly modified version that accepts a callback function that
+// calculates the result and the derivative at once
+pub fn newton_raphson_2<Func>(start: f64, fd: &Func) -> f64
+where
+    Func: Fn(f64) -> (f64, f64),
+{
+    // x[n + 1] = x[n] - f(x[n])/f'(x[n])
+
+    let mut x = start;
+
+    for _ in 0..MAX_ITERATIONS {
+        let (y0, y1) = fd(x);
+
+        if y0.abs() < MAX_ERROR {
+            return x;
+        }
+
+        let delta = y0 / y1;
+
+        if delta.abs() < MAX_ERROR && y0.abs() < MAX_FX_TOL {
+            return x;
+        }
+
+        x -= delta;
+    }
+
+    f64::NAN
+}
+
 pub fn newton_raphson_with_default_deriv<Func>(start: f64, f: Func) -> f64
 where
     Func: Fn(f64) -> f64,
@@ -47,8 +77,8 @@ pub fn brentq<Func>(f: &Func, xa: f64, xb: f64, iter: usize) -> f64
 where
     Func: Fn(f64) -> f64,
 {
-    let xtol = 2e-14;
-    let rtol = 8.881784197001252e-16;
+    const XTOL: f64 = 2e-14;
+    const RTOL: f64 = 8.881784197001252e-16;
 
     let mut xpre = xa;
     let mut xcur = xb;
@@ -86,11 +116,15 @@ where
             fblk = fpre;
         }
 
-        let delta = (xtol + rtol * xcur.abs()) / 2.;
+        let delta = (XTOL + RTOL * xcur.abs()) / 2.;
         let sbis = (xblk - xcur) / 2.;
 
         if fcur == 0. || sbis.abs() < delta {
-            return xcur;
+            return if fcur.abs() < MAX_FX_TOL {
+                xcur
+            } else {
+                f64::NAN
+            };
         }
 
         if spre.abs() > delta && fcur.abs() < fpre.abs() {

core/src/periodic.rs

@@ -181,7 +181,6 @@ pub fn irr(values: &[f64], guess: Option<f64>) -> Result<f64, InvalidPaymentsErr
         self::npv(rate, values, Some(true))
     };
     let df = |rate| self::npv_deriv(rate, values);
-    let is_good_rate = |rate: f64| rate.is_finite() && f(rate).abs() < 1e-3;
 
     let guess = match guess {
         Some(g) => g,
@@ -193,13 +192,13 @@ pub fn irr(values: &[f64], guess: Option<f64>) -> Result<f64, InvalidPaymentsErr
 
     let rate = newton_raphson(guess, &f, &df);
 
-    if is_good_rate(rate) {
+    if utils::is_a_good_rate(rate, f) {
         return Ok(rate);
     }
 
     let rate = brentq(&f, -0.999999999999999, 100., 100);
 
-    if is_good_rate(rate) {
+    if utils::is_a_good_rate(rate, f) {
         return Ok(rate);
     }
 

core/src/scheduled/xirr.rs

@@ -1,8 +1,8 @@
 use super::{year_fraction, DayCount};
 use crate::{
     models::{validate, DateLike, InvalidPaymentsError},
-    optimize::{brentq, newton_raphson},
-    utils::{self, fast_pow},
+    optimize::{brentq, newton_raphson_2},
+    utils::fast_pow,
 };
 
 pub fn xirr(
@@ -15,32 +15,26 @@ pub fn xirr(
 
     let deltas = &day_count_factor(dates, day_count);
 
-    let f = |rate| {
-        if rate <= -1.0 {
-            // bound newton_raphson
-            return f64::INFINITY;
-        }
-        xnpv_result(amounts, deltas, rate)
-    };
-    let df = |rate| xnpv_result_deriv(amounts, deltas, rate);
+    let f = |rate| xnpv_result(amounts, deltas, rate);
+    let fd = |rate| xnpv_result_with_deriv(amounts, deltas, rate);
 
-    let rate = newton_raphson(guess.unwrap_or(0.1), &f, &df);
+    let rate = newton_raphson_2(guess.unwrap_or(0.1), &fd);
 
-    if utils::is_a_good_rate(rate, f) {
+    if rate.is_finite() {
         return Ok(rate);
     }
 
     let rate = brentq(&f, -0.999999999999999, 100., 100);
 
-    if utils::is_a_good_rate(rate, f) {
+    if rate.is_finite() {
         return Ok(rate);
     }
 
     let mut step = 0.01;
     let mut guess = -0.99999999999999;
     while guess < 1.0 {
-        let rate = newton_raphson(guess, &f, &df);
-        if utils::is_a_good_rate(rate, f) {
+        let rate = newton_raphson_2(guess, &fd);
+        if rate.is_finite() {
             return Ok(rate);
         }
         guess += step;
@@ -86,15 +80,32 @@ fn day_count_factor(dates: &[DateLike], day_count: Option<DayCount>) -> Vec<f64>
 
 // \sum_{i=1}^n \frac{P_i}{(1 + rate)^{(d_i - d_0)/365}}
 fn xnpv_result(payments: &[f64], deltas: &[f64], rate: f64) -> f64 {
+    if rate <= -1.0 {
+        // bound newton_raphson
+        return f64::INFINITY;
+    }
     payments.iter().zip(deltas).map(|(p, &e)| p * fast_pow(1.0 + rate, -e)).sum()
 }
 
 // XNPV first derivative
 // \sum_{i=1}^n P_i * (d_0 - d_i) / 365 * (1 + rate)^{((d_0 - d_i)/365 - 1)}}
 // simplify in order to reuse cached deltas (d_i - d_0)/365
 // \sum_{i=1}^n \frac{P_i * -(d_i - d_0) / 365}{(1 + rate)^{((d_i - d_0)/365 + 1)}}
-fn xnpv_result_deriv(payments: &[f64], deltas: &[f64], rate: f64) -> f64 {
-    payments.iter().zip(deltas).map(|(p, e)| p * -e * fast_pow(1.0 + rate, -e - 1.0)).sum()
+// fn xnpv_result_deriv(payments: &[f64], deltas: &[f64], rate: f64) -> f64 {
+//     payments.iter().zip(deltas).map(|(p, e)| p * -e * fast_pow(1.0 + rate, -e - 1.0)).sum()
+// }
+
+fn xnpv_result_with_deriv(payments: &[f64], deltas: &[f64], rate: f64) -> (f64, f64) {
+    if rate <= -1.0 {
+        return (f64::INFINITY, f64::INFINITY);
+    }
+    // pow is an expensive function.
+    // we can re-use the result of pow for derivative calculation
+    payments.iter().zip(deltas).fold((0.0, 0.0), |acc, (p, e)| {
+        let y0 = p * fast_pow(1.0 + rate, -e);
+        let y1 = y0 * -e / (1.0 + rate);
+        (acc.0 + y0, acc.1 + y1)
+    })
 }
 
 #[cfg(test)]