1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
#![allow(non_snake_case)]
use crate::algebra::FloatT;
use crate::solver::{core::ScalingStrategy, CoreSettings};
use enum_dispatch::*;
// the supported cone wrapper type for primitives
// and the composite cone
mod compositecone;
mod supportedcone;
// primitive cone types
mod expcone;
mod genpowcone;
mod nonnegativecone;
mod powcone;
mod socone;
mod zerocone;
// partially specialized traits and blanket implementataions
mod nonsymmetric_common;
mod symmetric_common;
//re-export everything to appear as one module
use nonsymmetric_common::*;
pub use {
compositecone::*, expcone::*, genpowcone::*, nonnegativecone::*, powcone::*, socone::*,
supportedcone::*, symmetric_common::*, zerocone::*,
};
// only use PSD cones with SDP/Blas enabled
#[cfg(feature = "sdp")]
mod psdtrianglecone;
#[cfg(feature = "sdp")]
pub use psdtrianglecone::*;
// marker for primal / dual distinctions
#[derive(Eq, PartialEq, Clone, Debug, Copy)]
pub enum PrimalOrDualCone {
PrimalCone,
DualCone,
}
#[enum_dispatch]
pub trait Cone<T>
where
T: FloatT,
{
// functions relating to basic sizing
fn degree(&self) -> usize;
fn numel(&self) -> usize;
//Can the cone provide a sparse expanded representation?
fn is_sparse_expandable(&self) -> bool;
// is the cone symmetric? NB: zero cone still reports true
fn is_symmetric(&self) -> bool;
// report false here if only dual scaling is implemented (e.g. GenPowerCone)
fn allows_primal_dual_scaling(&self) -> bool;
// converts an elementwise scaling into
// a scaling that preserves cone memership
fn rectify_equilibration(&self, δ: &mut [T], e: &[T]) -> bool;
// returns (α,β) such that:
// z - α⋅e is just on the cone boundary, with value
// α >=0 indicates z \in cone, i.e. negative margin ===
// outside of the cone.
//
// β is the sum of the margins that are positive. For most
// cones this will just be β = max(0.,α), but for cones that
// are composites (e.g. the R_n^+), it is the sum of all of
// the positive margin terms.
fn margins(&mut self, z: &mut [T], pd: PrimalOrDualCone) -> (T, T);
// functions relating to unit vectors and cone initialization
fn scaled_unit_shift(&self, z: &mut [T], α: T, pd: PrimalOrDualCone);
fn unit_initialization(&self, z: &mut [T], s: &mut [T]);
// Compute scaling points
fn set_identity_scaling(&mut self);
fn update_scaling(
&mut self, s: &[T], z: &[T], μ: T, scaling_strategy: ScalingStrategy
) -> bool;
// operations on the Hessian of the centrality condition
// : W^TW for symmmetric cones
// : μH(s) for nonsymmetric cones
fn Hs_is_diagonal(&self) -> bool;
fn get_Hs(&self, Hsblock: &mut [T]);
fn mul_Hs(&mut self, y: &mut [T], x: &[T], work: &mut [T]);
// ---------------------------------------------------------
// Linearized centrality condition functions
//
// For nonsymmetric cones:
// -----------------------
//
// The centrality condition is : s = -μg(z)
//
// The linearized version is :
// Δs + μH(z)Δz = -ds = -(affine_ds + combined_ds_shift)
//
// The affine term (computed in affine_ds!) is s
// The shift term is μg(z) plus any higher order corrections
//
// # To recover Δs from Δz, we can write
// Δs = - (ds + μHΔz)
// The "offset" in Δs_from_Δz_offset is then just ds
//
// For symmetric cones:
// --------------------
//
// The centrality condition is : (W(z + Δz) ∘ W⁻ᵀ(s + Δs) = μe
//
// The linearized version is :
// λ ∘ (WΔz + WᵀΔs) = -ds = - (affine_ds + combined_ds_shift)
//
// The affine term (computed in affine_ds!) is λ ∘ λ
// The shift term is W⁻¹Δs_aff ∘ WΔz_aff - σμe, where the terms
// Δs_aff an Δz_aff are from the affine KKT solve, i.e. they
// are the Mehrotra correction terms.
//
// To recover Δs from Δz, we can write
// Δs = - ( Wᵀ(λ \ ds) + WᵀW Δz)
// The "offset" in Δs_from_Δz_offset is then Wᵀ(λ \ ds)
//
// Note that the Δs_from_Δz_offset function is only needed in the
// general combined step direction. In the affine step direction,
// we have the identity Wᵀ(λ \ (λ ∘ λ )) = s. The symmetric and
// nonsymmetric cases coincide and offset is taken directly as s.
//
// The affine step directions terms steps_z and step_s are
// passed to combined_ds_shift as mutable. Once they have been
// used to compute the combined ds shift they are no longer needed,
// so may be modified in place as workspace.
// ---------------------------------------------------------
fn affine_ds(&self, ds: &mut [T], s: &[T]);
fn combined_ds_shift(&mut self, shift: &mut [T], step_z: &mut [T], step_s: &mut [T], σμ: T);
fn Δs_from_Δz_offset(&mut self, out: &mut [T], ds: &[T], work: &mut [T], z: &[T]);
// Find the maximum step length in some search direction
fn step_length(
&mut self,
dz: &[T],
ds: &[T],
z: &[T],
s: &[T],
settings: &CoreSettings<T>,
αmax: T,
) -> (T, T);
// return the barrier function at (z+αdz,s+αds)
fn compute_barrier(&mut self, z: &[T], s: &[T], dz: &[T], ds: &[T], α: T) -> T;
}
#[test]
fn numel_degree() {
use crate::solver::core::cones::*;
let zcone = ZeroCone::<f64>::new(5);
let nncone = NonnegativeCone::<f64>::new(5);
let scone = SecondOrderCone::<f64>::new(5);
let expcone = ExponentialCone::<f64>::new();
let powcone = PowerCone::<f64>::new(0.5);
assert_eq!(zcone.numel(), 5);
assert_eq!(zcone.degree(), 0);
assert_eq!(nncone.numel(), 5);
assert_eq!(nncone.degree(), 5);
assert_eq!(scone.numel(), 5);
assert_eq!(scone.degree(), 1);
assert_eq!(expcone.numel(), 3);
assert_eq!(expcone.degree(), 3);
assert_eq!(powcone.numel(), 3);
assert_eq!(powcone.degree(), 3);
#[cfg(feature = "sdp")]
{
let sdpcone = PSDTriangleCone::<f64>::new(5);
assert_eq!(sdpcone.numel(), 15);
assert_eq!(sdpcone.degree(), 5);
}
}