Treat 3/4D SOC cone as dense

Cargo.toml

@@ -22,6 +22,7 @@ enum_dispatch  = "0.3.8"
 amd            = "0.2.2"
 thiserror      = "1.0"
 cfg-if         = "1.0"
+itertools      = "0.11"
 
 # -------------------------------
 # features

src/algebra/csc/core.rs

@@ -69,20 +69,18 @@ where
     J: IntoIterator<Item = &'a T>,
     T: FloatT,
 {
+    #[allow(clippy::needless_range_loop)]
     fn from(rows: I) -> CscMatrix<T> {
         let rows: Vec<Vec<T>> = rows
             .into_iter()
-            .map(|r| r.into_iter().map(|&v| v).collect())
+            .map(|r| r.into_iter().copied().collect())
             .collect();
 
         let m = rows.len();
         let n = rows.iter().map(|r| r.len()).next().unwrap_or(0);
+
         assert!(rows.iter().all(|r| r.len() == n));
-        let nnz = rows
-            .iter()
-            .flat_map(|r| r)
-            .filter(|&&v| v != T::zero())
-            .count();
+        let nnz = rows.iter().flatten().filter(|&&v| v != T::zero()).count();
 
         let mut colptr = Vec::with_capacity(n + 1);
         let mut rowval = Vec::with_capacity(nnz);

src/algebra/dense/core.rs

@@ -42,6 +42,7 @@ pub struct Matrix<T = f64> {
 ///        [0.0, 4.0]]);
 // ```
 //
+#[allow(clippy::needless_range_loop)]
 impl<'a, I, J, T> From<I> for Matrix<T>
 where
     I: IntoIterator<Item = J>,
@@ -51,7 +52,7 @@ where
     fn from(rows: I) -> Matrix<T> {
         let rows: Vec<Vec<T>> = rows
             .into_iter()
-            .map(|r| r.into_iter().map(|&v| v).collect())
+            .map(|r| r.into_iter().copied().collect())
             .collect();
 
         let m = rows.len();

src/algebra/floats.rs

@@ -1,3 +1,4 @@
+#![allow(non_snake_case)]
 use num_traits::{Float, FloatConst, FromPrimitive, NumAssign};
 use std::fmt::{Debug, Display, LowerExp};
 
@@ -79,8 +80,6 @@ cfg_if::cfg_if! {
 // NB: `AsFloatT` is a convenience trait for f32/64 and u32/64
 // so that we can do things like (2.0).as_T() everywhere on
 // constants, rather than the awful T::from_f32(2.0).unwrap()
-
-#[allow(non_snake_case)]
 pub trait AsFloatT<T>: 'static {
     fn as_T(&self) -> T;
 }

src/algebra/math_traits.rs

@@ -96,7 +96,10 @@ pub trait VectorMath {
     /// Standard Euclidian or 2-norm distance from `self` to `y`
     fn dist(&self, y: &Self) -> Self::T;
 
-    /// Sum of elements squared.
+    /// Sum of elements.
+    fn sum(&self) -> Self::T;
+
+    /// Sum of squares of the elements.
     fn sumsq(&self) -> Self::T;
 
     /// 2-norm

src/algebra/vecmath.rs

@@ -1,4 +1,5 @@
 use super::{FloatT, ScalarMath, VectorMath};
+use itertools::izip;
 use std::iter::zip;
 
 impl<T: FloatT> VectorMath for [T] {
@@ -88,7 +89,7 @@ impl<T: FloatT> VectorMath for [T] {
         assert_eq!(s.len(), ds.len());
 
         let mut out = T::zero();
-        for ((&s, &ds), (&z, &dz)) in zip(zip(s, ds), zip(z, dz)) {
+        for (&s, &ds, &z, &dz) in izip!(s, ds, z, dz) {
             let si = s + α * ds;
             let zi = z + α * dz;
             out += si * zi;
@@ -101,6 +102,10 @@ impl<T: FloatT> VectorMath for [T] {
         T::sqrt(dist2)
     }
 
+    fn sum(&self) -> T {
+        self.iter().fold(T::zero(), |acc, &x| acc + x)
+    }
+
     fn sumsq(&self) -> T {
         self.dot(self)
     }

src/solver/core/cones/compositecone.rs

@@ -52,7 +52,7 @@ where
         // create cones with the given dims
         for t in types.iter() {
             //make a new cone
-            let cone = make_cone(*t);
+            let cone = make_cone(t);
 
             //update global problem symmetry
             _is_symmetric = _is_symmetric && cone.is_symmetric();
@@ -157,7 +157,7 @@ where
     pub fn iter_mut(&mut self) -> std::slice::IterMut<'_, SupportedCone<T>> {
         self.cones.iter_mut()
     }
-    pub(crate) fn type_count(&self, tag: SupportedConeTag) -> usize {
+    pub(crate) fn get_type_count(&self, tag: SupportedConeTag) -> usize {
         if self.type_counts.contains_key(&tag) {
             self.type_counts[&tag]
         } else {
@@ -182,6 +182,18 @@ where
         self._is_symmetric
     }
 
+    fn is_sparse_expandable(&self) -> bool {
+        //This should probably never be called
+        //self.cones.iter().any(|cone| cone.is_sparse_expandable())
+        unreachable!();
+    }
+
+    fn allows_primal_dual_scaling(&self) -> bool {
+        self.cones
+            .iter()
+            .all(|cone| cone.allows_primal_dual_scaling())
+    }
+
     fn rectify_equilibration(&self, δ: &mut [T], e: &[T]) -> bool {
         let mut any_changed = false;
 
@@ -248,14 +260,7 @@ where
         //This function should probably never be called since
         //we only us it to interrogate the blocks, but we can
         //implement something reasonable anyway
-        let mut is_diag = true;
-        for cone in self.iter() {
-            is_diag &= cone.Hs_is_diagonal();
-            if !is_diag {
-                break;
-            }
-        }
-        is_diag
+        self.cones.iter().all(|cone| cone.Hs_is_diagonal())
     }
 
     #[allow(non_snake_case)]
@@ -348,9 +353,9 @@ where
         (α, α)
     }
 
-    fn compute_barrier(&self, z: &[T], s: &[T], dz: &[T], ds: &[T], α: T) -> T {
+    fn compute_barrier(&mut self, z: &[T], s: &[T], dz: &[T], ds: &[T], α: T) -> T {
         let mut barrier = T::zero();
-        for (cone, rng) in zip(&self.cones, &self.rng_cones) {
+        for (cone, rng) in zip(&mut self.cones, &self.rng_cones) {
             let zi = &z[rng.clone()];
             let si = &s[rng.clone()];
             let dzi = &dz[rng.clone()];