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 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321
use std::iter::FromIterator;
use crate::{ Med, MStats, MinMax, MutVecg, MutVecf64, Stats, Vecg, Vecf64, VecVec};
pub use indxvec::{here,tof64,Printing,merge::*,Indices};
impl<T> VecVec<T> for &[Vec<T>] where T: Copy+PartialOrd+std::fmt::Display,
f64: From<T> {
/// Transpose vec of vecs like a matrix
fn transpose(self) -> Vec<Vec<T>> {
let n = self.len();
let d = self[0].len();
let mut transp:Vec<Vec<T>> = Vec::with_capacity(d);
for i in 0..d {
let mut column = Vec::with_capacity(n);
for v in self {
column.push(v[i]);
}
transp.push(column); // column becomes row
}
transp
}
/// Joint probability density function of n matched slices of the same length
fn jointpdfn(self) -> Vec<f64> {
let d = self[0].len(); // their common dimensionality (length)
for v in self.iter().skip(1) {
if v.len() != d { panic!("{} all vectors must be of equal length!",here!()) };
}
let mut res:Vec<f64> = Vec::with_capacity(d);
let mut tuples = self.transpose();
let df = tuples.len() as f64; // for turning counts to probabilities
println!("{}",df);
// lexical sort to group together occurrences of identical tuples
tuples.sort_unstable_by(|a, b| a.partial_cmp(b).unwrap());
let mut count = 1_usize; // running count
let mut lastindex = 0; // initial index of the last unique tuple
tuples.iter().enumerate().skip(1).for_each(|(i,ti)|
if ti > &tuples[lastindex] { // new tuple ti (Vec<T>) encountered
res.push((count as f64)/df); // save frequency count as probability
lastindex = i; // current index becomes the new one
count = 1_usize; // reset counter
}
else { count += 1; } );
res.push((count as f64)/df); // flush the rest!
res
}
/// Joint entropy of vectors of the same length
fn jointentropyn(self) -> f64 {
let jpdf = self.jointpdfn();
jpdf.iter().map(|&x| -x*(x.ln()) ).sum()
}
/// Dependence (component wise) of a set of vectors.
/// i.e. `dependencen` returns 0 iff they are statistically independent
/// bigger values when they are dependent
fn dependencen(self) -> f64 {
self.iter().map(|v| v.entropy()).sum::<f64>()/self.jointentropyn() - 1.0
}
/// Flattened lower triangular part of a symmetric matrix for column vectors in self.
/// The upper triangular part can be trivially generated for all j>i by: c(j,i) = c(i,j).
/// Applies closure f which computes a scalar relationship between two vectors,
/// that is different features stored in columns of self.
/// The closure typically invokes one of the methods from Vecg trait (in vecg.rs),
/// such as dependencies or correlations.
/// Example call: `pts.transpose().crossfeatures(|v1,v2| v1.mediancorr(v2))`
/// computes correlations between all column vectors (features) in pts.
fn crossfeatures<F>(self,f:F) -> Vec<f64> where F: Fn(&[T],&[T]) -> f64 {
let n = self.len(); // number of the vector(s)
let mut codp:Vec<f64> = Vec::with_capacity((n+1)*n/2); // results
self.iter().enumerate().for_each(|(i,v)|
// its dependencies up to and including the diagonal
self.iter().take(i+1).for_each(|vj| {
codp.push(f(v,vj));
}));
codp
}
/// acentroid = simple multidimensional arithmetic mean
fn acentroid(self) -> Vec<f64> {
let mut centre = vec![0_f64; self[0].len()];
for v in self { centre.mutvadd(v) }
centre.mutsmultf64(1.0 / (self.len() as f64));
centre
}
/// gcentroid = multidimensional geometric mean
fn gcentroid(self) -> Vec<f64> {
let nf = self.len() as f64; // number of points
let d = self[0].len(); // dimensions
let mut centre = vec![0_f64; d];
let mut lnvec = vec![0_f64; d];
for v in self {
v.iter().zip(&mut lnvec).for_each(|(&vi,lni)| *lni = f64::from(vi).ln() );
centre.mutvaddf64(&lnvec)
}
centre.iter().map(|comp| (comp/nf).exp()).collect()
}
/// hcentroid = multidimensional harmonic mean
fn hcentroid(self) -> Vec<f64> {
let mut centre = vec![0_f64; self[0].len()];
for v in self { centre.mutvaddf64(&v.vinverse().unwrap()) }
centre.smultf64(1.0/(self.len() as f64)).vinverse().unwrap()
}
/// For each member point, gives its sum of distances to all other points and their MinMax
fn distsums(self) -> Vec<f64> {
let n = self.len();
let mut dists = vec![0_f64; n]; // distances accumulator for all points
// examine all unique pairings (lower triangular part of symmetric flat matrix)
self.iter().enumerate().for_each(|(i,thisp)|
self.iter().take(i).enumerate().for_each(|(j,thatp)| {
let d = thisp.vdist(thatp); // calculate each distance relation just once
dists[i] += d;
dists[j] += d; // but add it to both points' sums
}));
dists
}
/// The sum of distances from one member point, given by its `indx`, to all the other points in self.
/// For all the points, use more efficient `distsums`.
fn distsuminset(self, indx: usize) -> f64 {
let thisp = &self[indx];
self.iter().enumerate().map(|(i,thatp)|
if i == indx { 0.0 } else {thisp.vdist(thatp)}).sum()
}
/// Medoid and Outlier (Medout)
/// Medoid is the member point (point belonging to the set of points `self`),
/// which has the least sum of distances to all other points.
/// Outlier is the point with the greatest sum of distances.
/// In other words, they are the members nearest and furthest from the geometric median.
/// Returns struct MinMax{min,minindex,max,maxindex}
fn medout(self) -> MinMax<f64> {
minmax(&self.distsums())
}
/// Finds approximate vectors from each member point towards the geometric median.
/// Twice as fast using symmetry, as doing them individually.
fn eccentricities(self) -> Vec<Vec<f64>> {
let n = self.len();
// allocate vectors for the results
let mut eccs = vec![vec![0_f64; self[0].len()]; n];
let mut recips = vec![0_f64; n];
// ecentricities vectors accumulator for all points
// examine all unique pairings (lower triangular part of symmetric flat matrix)
for i in 1..n {
let thisp = &self[i];
for j in 0..i {
// calculate each unit vector between any pair of points just once
let dvmag = self[j].vdist(thisp);
if !dvmag.is_normal() { continue }
let rec = 1.0_f64/dvmag;
eccs[i].mutvaddf64(&self[j].smultf64(rec));
recips[i] += rec;
// mind the vector's opposite orientations w.r.t. to the two points!
eccs[j].mutvsubf64(&self[j].smultf64(rec));
recips[j] += rec; // but scalar distances are the same
}
}
for i in 0..n {
eccs[i].mutsmultf64(1.0/recips[i]);
eccs[i].mutvsub(&self[i])
}
eccs
}
/// Exact eccentricity vectors from all member points by first finding the Geometric Median.
/// Usually faster than the approximate `eccentricities` above, especially when there are many points.
fn exacteccs(self, eps:f64) -> Vec<Vec<f64>> {
let mut eccs = Vec::with_capacity(self.len()); // Vectors for the results
let gm:Vec<f64> = self.gmedian(eps);
for v in self {
eccs.push(gm.vsub(v))
}
eccs
}
/// Estimated (computed) eccentricity vector for a member point.
/// It points towards the geometric median.
/// The true geometric median is as yet unknown.
/// The true geometric median would return zero vector.
/// The member point in question is specified by its index `indx`.
/// This function is suitable for a single member point.
/// When eccentricities of all the points are wanted, use `exacteccs` above.
fn eccmember(self, indx: usize) -> Vec<f64> {
self.nxmember(indx).vsub(&self[indx])
}
/// Estimated (computed) eccentricity vector for a non member point.
/// The true geometric median is as yet unknown.
/// Returns the eccentricity vector.
/// The true geometric median would return zero vector.
/// This function is suitable for a single non-member point.
fn eccnonmember(self, p:&[f64]) -> Vec<f64> {
self.nxnonmember(p).vsubf64(p)
}
/// Mean and Std (in MStats struct), Median and quartiles (in Med struct), Median and Outlier (in MinMax struct)
/// of scalar eccentricities of points in self.
/// These are new robust measures of a cloud of multidimensional points (or multivariate sample).
fn eccinfo(self, eps: f64) -> (MStats, Med, MinMax<f64>) where Vec<f64>:FromIterator<f64> {
let gm = self.gmedian(eps);
let eccs:Vec<f64> = self.iter().map(|v| gm.vdist(v)).collect();
(eccs.ameanstd().unwrap(),eccs.median().unwrap(),minmax(&eccs))
}
/// MADn multidimensional median of distances from gm: data spread estimator that is more stable than variances
fn madn(self, eps: f64) -> f64 {
let gm = self.gmedian(eps);
let eccs:Vec<f64> = self.iter().map(|v| gm.vdist(v)).collect();
let Med{median,..} = eccs.median().unwrap_or_else(|_| panic!("{},median failed\n",here!()));
median
}
/// GM and sorted eccentricities magnitudes.
/// Describing a set of points `self` in n dimensions
fn sortedeccs(self, ascending:bool, gm:&[f64]) -> Vec<f64> {
let mut eccs = Vec::with_capacity(self.len());
// collect raw ecentricities magnitudes
for v in self { eccs.push(v.vdistf64(gm)) }
sortm(&eccs,ascending)
}
/// Eccentricities of Medoid and Outlier.
/// Same as just the third element of a tuple returned by eccinfo
fn emedoid(self, eps: f64) -> MinMax<f64> where Vec<f64>:FromIterator<f64> {
let gm:Vec<f64> = self.gmedian(eps);
let eccs:Vec<f64> = self.iter().map(|v| gm.vdist(v)).collect();
minmax(&eccs)
}
/// Initial (first) point for geometric medians.
/// Same as eccnonmember('origin') but saving the subtractions of zeroes.
fn firstpoint(self) -> Vec<f64> {
let mut rsum = 0_f64;
let mut vsum = vec![0_f64; self[0].len()];
for p in self {
let mag = p.vmag();
if mag.is_normal() { // skip if p is at the origin
let rec = 1.0_f64/mag;
// the sum of reciprocals of magnitudes for the final scaling
rsum += rec;
// so not using simply .unitv
vsum.mutvaddf64(&p.smultf64(rec)) // add all unit vectors
}
}
vsum.mutsmultf64(1.0/rsum); // scale by the sum of reciprocals
vsum
}
/// Next approximate gm computed from a member point
/// specified by its index `indx` to self.
fn nxmember(self, indx: usize) -> Vec<f64> {
let mut vsum = vec![0_f64; self[0].len()];
let p = &tof64(&self[indx]);
let mut recip = 0_f64;
for (i,x) in self.iter().enumerate() {
if i != indx { // not point p
let mag:f64 = x.iter().zip(p).map(|(&xi,&pi)|(f64::from(xi)-pi).powi(2)).sum::<f64>().sqrt();
if mag.is_normal() { // ignore this point should distance be zero
let rec = 1.0_f64/mag;
vsum.iter_mut().zip(x).for_each(|(vi,xi)| *vi += rec*f64::from(*xi));
recip += rec // add separately the reciprocals
}
}
};
vsum.iter_mut().for_each(|vi| *vi /= recip);
vsum
}
/// Next approximate gm computed from a non-member point p
fn nxnonmember(self, p:&[f64]) -> Vec<f64> {
let mut vsum = vec![0_f64; self[0].len()];
let mut recip = 0_f64;
for x in self {
let mag:f64 = x.iter().zip(p).map(|(&xi,&pi)|(f64::from(xi)-pi).powi(2)).sum::<f64>().sqrt();
if mag.is_normal() { // ignore this point should distance be zero
let rec = 1.0_f64/mag;
vsum.iter_mut().zip(x).for_each(|(vi,xi)| *vi += rec*f64::from(*xi));
recip += rec // add separately the reciprocals
}
}
vsum.iter_mut().for_each(|vi| *vi /= recip);
vsum
}
/// Geometric Median (gm) is the point that minimises the sum of distances to a given set of points.
/// It has (provably) only vector iterative solutions.
/// Search methods are slow and difficult in highly dimensional space.
/// Weiszfeld's fixed point iteration formula had known problems with sometimes failing to converge.
/// Especially, when the points are dense in the close proximity of the gm, or gm coincides with one of them.
/// However, these problems are fixed in my new algorithm here.
/// There will eventually be a multithreaded version.
fn gmedian(self, eps: f64) -> Vec<f64> {
let eps2 = eps.powi(2);
let mut p = self.acentroid(); // start iterating from the Centre
loop { // vector iteration till accuracy eps2 is reached
let nextp = self.nxnonmember(&p);
if nextp.iter().zip(p).map(|(&xi,pi)|(xi-pi).powi(2)).sum::<f64>() < eps2 { return nextp }; // termination
p = nextp;
};
}
/// Same a gmedian but returns also the number of iterations
fn igmedian(self, eps: f64) -> ( Vec<f64>, usize ) {
let eps2 = eps.powi(2);
let mut p = self.acentroid();
let mut iterations = 1_usize;
loop {
let nextp = self.nxnonmember(&p);
if nextp.vdistsqf64(&p) < eps2 { return (nextp,iterations) }; // termination
iterations +=1;
p = nextp;
};
}
}