quantwave_core/portfolio/
projection.rs

1//! Constraint projection: capped-simplex projection for box + budget=1.
2
3use super::{Bounds, PortfolioError, PortfolioResult};
4use nalgebra::DVector;
5use serde::{Deserialize, Serialize};
6
7/// Directional constraint applied to portfolio weights during projection.
8#[derive(Debug, Clone, Copy, PartialEq, Eq, Default, Serialize, Deserialize)]
9pub enum PortfolioConstraint {
10    /// Long-only weights (non-negative lower bound of 0.0).
11    #[default]
12    LongOnly,
13    /// Long-short weights (unbounded negative lower bound of -inf).
14    LongShort,
15}
16
17impl PortfolioConstraint {
18    pub const fn is_long_only(self) -> bool {
19        matches!(self, Self::LongOnly)
20    }
21
22    pub const fn is_long_short(self) -> bool {
23        matches!(self, Self::LongShort)
24    }
25}
26
27impl From<bool> for PortfolioConstraint {
28    fn from(long_only: bool) -> Self {
29        if long_only {
30            Self::LongOnly
31        } else {
32            Self::LongShort
33        }
34    }
35}
36
37impl From<PortfolioConstraint> for bool {
38    fn from(constraint: PortfolioConstraint) -> Self {
39        constraint.is_long_only()
40    }
41}
42
43/// Project raw weights onto the box-constrained, budget=1 feasible set.
44///
45/// Computes the Euclidean projection of `w` onto
46/// `{x : lo <= x <= hi, sum(x) = 1}` via bisection on a scalar shift
47/// `lambda`, using the standard capped-simplex projection identity
48/// `x_i = clip(w_i - lambda, lo_i, hi_i)`. `sum(clip(w - lambda, lo, hi))`
49/// is monotonically non-increasing in `lambda`, so bisection converges to
50/// the unique `lambda` with `sum(x) = 1` whenever
51/// `sum(lo) <= 1 <= sum(hi)` (validated below). This guarantees bounds are
52/// respected *exactly* (unlike clip-then-renormalize, which can push
53/// values back out of bounds).
54///
55/// Direct port of `_apply_bounds_and_budget` in the Python v1
56/// (`quantwave/portfolio.py`).
57pub fn apply_bounds_and_budget(
58    w: &DVector<f64>,
59    bounds: Option<&[Bounds]>,
60    constraint: PortfolioConstraint,
61) -> PortfolioResult<DVector<f64>> {
62    let n = w.len();
63
64    let (lo, hi): (DVector<f64>, DVector<f64>) = match bounds {
65        None => {
66            let lo = if constraint.is_long_only() {
67                DVector::zeros(n)
68            } else {
69                DVector::from_element(n, f64::NEG_INFINITY)
70            };
71            let hi = DVector::from_element(n, 1.0);
72            (lo, hi)
73        }
74        Some(b) => {
75            if b.len() != n {
76                return Err(PortfolioError::BoundsLengthMismatch(b.len(), n));
77            }
78            let lo = DVector::from_iterator(n, b.iter().map(|x| x.0));
79            let hi = DVector::from_iterator(n, b.iter().map(|x| x.1));
80            if lo.iter().zip(hi.iter()).any(|(l, h)| l > h) {
81                return Err(PortfolioError::BoundsLowGreaterThanHigh);
82            }
83            (lo, hi)
84        }
85    };
86
87    if hi.sum() < 1.0 - 1e-9 {
88        return Err(PortfolioError::BoundsUpperSumTooLow);
89    }
90    let finite_lo_sum: f64 = lo.iter().filter(|v| v.is_finite()).sum();
91    if finite_lo_sum > 1.0 + 1e-9 {
92        return Err(PortfolioError::BoundsLowerSumTooHigh);
93    }
94
95    let lo_b: DVector<f64> =
96        DVector::from_iterator(n, lo.iter().map(|v| if v.is_finite() { *v } else { -1e9 }));
97    let hi_b: DVector<f64> =
98        DVector::from_iterator(n, hi.iter().map(|v| if v.is_finite() { *v } else { 1e9 }));
99
100    let excess = |lam: f64| -> f64 {
101        let sum: f64 = (0..n).map(|i| (w[i] - lam).clamp(lo_b[i], hi_b[i])).sum();
102        sum - 1.0
103    };
104
105    let mut lam_lo = -1e9_f64;
106    let mut lam_hi = 1e9_f64;
107    for _ in 0..200 {
108        let mid = (lam_lo + lam_hi) / 2.0;
109        if excess(mid) > 0.0 {
110            lam_lo = mid;
111        } else {
112            lam_hi = mid;
113        }
114        if lam_hi - lam_lo < 1e-14 {
115            break;
116        }
117    }
118    let lam = (lam_lo + lam_hi) / 2.0;
119
120    let mut w_final =
121        DVector::from_iterator(n, (0..n).map(|i| (w[i] - lam).clamp(lo_b[i], hi_b[i])));
122    let total = w_final.sum();
123    if !w_final.iter().all(|v| v.is_finite()) || (total - 1.0).abs() > 1e-6 {
124        return Err(PortfolioError::ProjectionFailed);
125    }
126    // Nudge for the last bit of float error so sum is exactly 1.
127    w_final /= total;
128    Ok(w_final)
129}
130
131#[cfg(test)]
132mod tests {
133    use super::*;
134
135    #[test]
136    fn unconstrained_long_only_sums_to_one() {
137        let w = DVector::from_vec(vec![0.5, 0.3, 0.2]);
138        let out = apply_bounds_and_budget(&w, None, PortfolioConstraint::LongOnly).unwrap();
139        assert!((out.sum() - 1.0).abs() < 1e-9);
140        assert!(out.iter().all(|v| *v >= 0.0));
141    }
142
143    #[test]
144    fn negative_weight_clipped_to_zero_and_redistributed() {
145        let w = DVector::from_vec(vec![-0.1, 0.6, 0.5]);
146        let out = apply_bounds_and_budget(&w, None, PortfolioConstraint::LongOnly).unwrap();
147        assert!((out.sum() - 1.0).abs() < 1e-9);
148        assert!(out.iter().all(|v| *v >= -1e-9));
149    }
150
151    #[test]
152    fn infeasible_upper_bounds_error() {
153        let w = DVector::from_vec(vec![0.5, 0.5]);
154        let bounds = vec![(0.0, 0.3), (0.0, 0.3)];
155        let err =
156            apply_bounds_and_budget(&w, Some(&bounds), PortfolioConstraint::LongOnly).unwrap_err();
157        assert_eq!(err, PortfolioError::BoundsUpperSumTooLow);
158    }
159}
quantwave_core/portfolio/projection.rs

quantwave_core/portfolio/
projection.rs
