Skip to content

Commit

Permalink
under development [no ci]
Browse files Browse the repository at this point in the history
  • Loading branch information
liborty committed Mar 7, 2024
1 parent 0bc9cc3 commit 754b05a
Show file tree
Hide file tree
Showing 3 changed files with 73 additions and 65 deletions.
104 changes: 61 additions & 43 deletions src/algos.rs
Original file line number Diff line number Diff line change
@@ -1,8 +1,39 @@
use crate::{merror, Me};
use core::cmp::{Ordering, Ordering::*};
use indxvec::Mutops;
use std::ops::Range;

/// Partitions `s: &mut [&T]` within range `rng`, using comparator `c`.
/// Returns the boundaries of the rearranged partitions, (eqstart,gtstart), where
/// `rng.start..eqstart` (may be empty) contains references to items lesser than the pivot,
/// `eqstart..gtstart`are items equal to the pivot,
/// `gtstart..rng.end` (may be empty) contains references to items greater than the pivot.
pub fn partit<'a,T>(
s: &mut [&'a T],
pivot: &'a T,
rng: &Range<usize>,
c: &mut impl FnMut(&T, &T) -> Ordering,
) -> (usize, usize) {
let mut eqstart = rng.start;
let mut gtstart = eqstart;
for t in rng.start..rng.end {
match c(s[t], pivot) {
Less => {
s[eqstart] = s[t];
eqstart += 1;
s[t] = s[gtstart];
gtstart += 1;
}
Equal => {
s[t] = s[gtstart];
gtstart += 1;
}
Greater => (),
}
}
for e in eqstart..gtstart { s[e] = pivot; }
(eqstart, gtstart)
}

/// Scan a slice of f64s for NANs
pub fn nans(v: &[f64]) -> bool {
for &f in v {
Expand Down Expand Up @@ -186,31 +217,26 @@ pub fn medianu8(s: &[u8]) -> Result<f64, Me> {
}

/// Median of odd sized generic data with Odering comparisons by custom closure
pub(super) fn oddmed_by<'a, T>(s: &mut [&'a T], c: &mut impl FnMut(&T, &T) -> Ordering) -> &'a T {
let mut rng = 0..s.len();
let need = s.len() / 2; // median target position in fully partitioned set
pub(super) fn oddmed_by<'a, T>(
s: &mut [&'a T],
mut rng: Range<usize>,
c: &mut impl FnMut(&T, &T) -> Ordering,
) -> &'a T {
let need = rng.len() / 2; // median target position in fully partitioned set
loop {
let len = rng.len();
if len > 99 {
let mut pivotref = oddmed_by(
&mut s[rng.start..rng.start+len/20],
c,
);
std::mem::swap(&mut s[rng.start],&mut pivotref);
let n = rng.len();
let pivotref = if n > 99 {
oddmed_by(s, rng.start..rng.start + n / 20, c)
} else {
let pivotsub = middling(
rng.start,
s[middling(
rng.start,
rng.start + 1,
rng.end - 2,
rng.end - 1,
&mut |a, b| c(s[a], s[b]),
);
if pivotsub != rng.start {
s.swap(rng.start, pivotsub);
};
)]
};

let (eqsub, gtsub) = <&mut [T]>::part(s, &rng, c);
let (eqsub, gtsub) = partit(s, pivotref, &rng, c);
// well inside lt partition, iterate on it
if need + 2 < eqsub {
rng.end = eqsub;
Expand All @@ -228,7 +254,7 @@ pub(super) fn oddmed_by<'a, T>(s: &mut [&'a T], c: &mut impl FnMut(&T, &T) -> Or
};
if need < gtsub {
// within equals partition, return the pivot
return s[rng.start];
return pivotref;
};
// first place in gt partition, the solution is its minimum
if need == gtsub {
Expand All @@ -245,31 +271,24 @@ pub(super) fn oddmed_by<'a, T>(s: &mut [&'a T], c: &mut impl FnMut(&T, &T) -> Or
/// Median of even sized generic data with Odering comparisons by custom closure
pub(super) fn evenmed_by<'a, T>(
s: &mut [&'a T],
mut rng: Range<usize>,
c: &mut impl FnMut(&T, &T) -> Ordering,
) -> (&'a T, &'a T) {
let mut rng = 0..s.len();
let need = s.len() / 2 - 1; // median target position in fully partitioned set
loop {
let len = rng.len();
if len > 99 {
let mut pivotref = oddmed_by(
&mut s[rng.start..rng.start+len/20],
c,
);
std::mem::swap(&mut s[rng.start],&mut pivotref);
} else {
let pivotsub = middling(
rng.start,
rng.start + 1,
rng.end - 2,
rng.end - 1,
&mut |a, b| c(s[a], s[b]),
);
if pivotsub != rng.start {
s.swap(rng.start, pivotsub);
loop {
let n = rng.len();
let pivotref = if n > 99 {
oddmed_by(s, rng.start..rng.start + n / 20, c)
} else {
s[middling(
rng.start,
rng.start + 1,
rng.end - 2,
rng.end - 1,
&mut |a, b| c(s[a], s[b]),
)]
};
};
let (eqsub, gtsub) = <&mut [T]>::part(s, &rng, c);
let (eqsub, gtsub) = partit(s, pivotref, &rng, c);
// well inside lt partition, iterate on it narrowing the range
if need + 2 < eqsub {
rng.end = eqsub;
Expand All @@ -284,8 +303,7 @@ pub(super) fn evenmed_by<'a, T>(
// first place in gt partition, the solution are its two minima
if need == gtsub {
return min2(s, gtsub..rng.end, c);
};
let pivotref = s[rng.start];
};
// last place in the lt partition:
if need + 1 == eqsub {
// swapped comparison arguments to get maximum
Expand Down
30 changes: 10 additions & 20 deletions src/implementations.rs
Original file line number Diff line number Diff line change
Expand Up @@ -58,10 +58,10 @@ impl Medianf64 for &[f64] {
})
.collect::<Result<Vec<&f64>, Me>>()?;
if (n & 1) == 1 {
let oddm = oddmed_by(&mut s, &mut <f64>::total_cmp);
let oddm = oddmed_by(&mut s, 0..n, &mut <f64>::total_cmp);
Ok(*oddm)
} else {
let (&med1, &med2) = evenmed_by(&mut s, &mut <f64>::total_cmp);
let (&med1, &med2) = evenmed_by(&mut s, 0..n, &mut <f64>::total_cmp);
Ok((med1+med2) / 2.0)
}
}
Expand All @@ -76,12 +76,12 @@ impl Medianf64 for &[f64] {
2 => return (self[0] + self[1]) / 2.0,
_ => (),
};
let mut s = self.ref_vec(0..self.len());
let mut s = self.ref_vec(0..n);
if (n & 1) == 1 {
let oddm = oddmed_by(&mut s, &mut <f64>::total_cmp);
let oddm = oddmed_by(&mut s, 0..n,&mut <f64>::total_cmp);
*oddm
} else {
let (&med1, &med2) = evenmed_by(&mut s, &mut <f64>::total_cmp);
let (&med1, &med2) = evenmed_by(&mut s, 0..n,&mut <f64>::total_cmp);
(med1 + med2) / 2.0
}
}
Expand Down Expand Up @@ -175,9 +175,9 @@ impl<'a, T> Median<'a, T> for &'a [T] {
};
let mut s = self.ref_vec(0..self.len());
if (n & 1) == 1 {
Ok(q(oddmed_by(&mut s, c)))
Ok(q(oddmed_by(&mut s, 0..n,c)))
} else {
let (med1, med2) = evenmed_by(&mut s, c);
let (med1, med2) = evenmed_by(&mut s, 0..n, c);
Ok((q(med1) + q(med2)) / 2.0)
}
}
Expand All @@ -193,9 +193,9 @@ impl<'a, T> Median<'a, T> for &'a [T] {
};
let mut s = self.ref_vec(0..self.len());
if (n & 1) == 1 {
Ok(Medians::Odd(oddmed_by(&mut s, c)))
Ok(Medians::Odd(oddmed_by(&mut s, 0..n, c)))
} else {
Ok(Medians::Even(evenmed_by(&mut s, c)))
Ok(Medians::Even(evenmed_by(&mut s, 0..n, c)))
}
}

Expand All @@ -204,17 +204,7 @@ impl<'a, T> Median<'a, T> for &'a [T] {
Ok(self.iter().map(|s| q(s) - centre).collect())
}
/// We define median based correlation as cosine of an angle between two
/// zero median vectors (analogously to Pearson's zero mean vectors)
/// # Example
/// ```
/// use medians::{Medianf64,Median};
/// use core::convert::identity;
/// use core::cmp::Ordering::*;
/// let v1 = vec![1_f64,2.,3.,4.,5.,6.,7.,8.,9.,10.,11.,12.,13.,14.];
/// let v2 = vec![14_f64,1.,13.,2.,12.,3.,11.,4.,10.,5.,9.,6.,8.,7.];
/// assert_eq!(v1.medf_correlation(&v2).unwrap(),-0.08131842079168768);
/// assert_eq!(v1.med_correlation(&v2,&mut |a,b| a.total_cmp(b),|&a| identity(a)).unwrap(),-0.08131842079168768);
/// ```
/// zero median vectors (analogously to Pearson's zero mean vectors)
fn med_correlation(
self,
v: Self,
Expand Down
4 changes: 2 additions & 2 deletions tests/tests.rs
Original file line number Diff line number Diff line change
Expand Up @@ -18,9 +18,9 @@ fn parting() -> Result<(), Me> {
// println!("To f64s: {}",to_f64s(&to_u64s(&data)).gr());
// println!("Scrubbed: {}", scrub_nans(&to_f64s(&to_u64s(&data))).gr());
let len = data.len();
println!("Pivot {}: {}", data[0].yl(), data.gr());
println!("Pivot {}: {}", 7.0.yl(), data.gr());
let mut refdata = data.ref_vec(0..len);
let (eqsub, gtsub) = <&mut [f64]>::part(&mut refdata, &(0..len), &mut <f64>::total_cmp);
let (eqsub, gtsub) = partit(&mut refdata, &7.0, &(0..len), &mut <f64>::total_cmp);
println!(
"Result: {}\nCommas show the subranges:\n\
{GR}[{}, {}, {}]{UN}\n{} items equal to pivot {}",
Expand Down

0 comments on commit 754b05a

Please sign in to comment.