SR-RCD

Apply Refining-Cover-by-Defect algorithm to solve Sound Ranging problem in time-dependent-metric (and, in particular, (quasi-)metric) spaces.

Part of "Decaennead: ruminations about sound ranging in abstract space(time)s and related themes" project.

Sound ranging

There is an unknown point s, source, which at unknown instant t₀ of time emits the sound. The sound wave propagates through space and reaches known points r_i, sensors, at known instants t_i. The sound ranging (SR) problem, also called source localization or time-difference-of-arrival (TDOA) problem, is to determine s and t₀ from {r_i, t_i}.

Supposedly the earliest records of a SR problem date from the World War I period; inter alia, in 1915, it was briefly mentioned in the correspondence between E. Borel and H. Lebesgue.

By replacing sound with electromagnetic waves or some kind of "propagation", e.g. in a graph, one moves from acoustical to optical and other contexts.

Timetrics

F-timetric. Let X be a non-empty set and τ be a nonnegative function defined on ℝ×X×X that satisfies the following axioms:

1. Forward reflexivity: for any t ∊ ℝ: τ_t(x, y) = 0 if and only if x = y,

2. Forward timed triangle inequality: for any t ∊ ℝ and any x, y, z ∊ X: τ_t(x, y) ≤ τ_t(x, z) + τ_{t + τt(x, z)}(z, y).

Then (X, τ) is called an f-timetric space.

In SR context, τ_t(x, y) is interpreted as the time needed for the sound wave that was emitted in x at the instant t to propagate from x to y. This wave reaches y at the instant t + τ_t(x, y).

B-timetric. Let X be a non-empty set and τ̅ be a nonnegative function defined on ℝ×X×X that satisfies

1. Backward reflexivity: for any t ∊ ℝ: τ̅_t(x, y) = 0 if and only if x = y,

2. Backward timed triangle inequality: for any t ∊ ℝ and any x, y, z ∊ X: τ̅_t(x, y) ≤ τ̅_t(z, y) + τ̅_{t - τ̅t(z, y)}(x, z).

Then (X, τ̅) is called a b-timetric space.

In SR context, τ̅_t(x, y) is interpreted as the time needed for the sound wave that was emitted in x and reached y at the instant t to propagate from x to y. This wave was emitted in x at the instant t - τ̅_t(x, y).

Timed triangle inequalities represent "a wave propagates through a medium as soon as possible" notion. Note that there is τ_{t + τt(x, z)}(z, y), not just τ_t(z, y), in e.g. forward TTI, because a medium changes and waves propagate through it simultaneously.

The connection between f- and b-timetrics:

τ̅_t(x, y) = τ_{t - τ̅t(x, y)}(x, y)

τ_t(x, y) = τ̅_{t + τt(x, y)}(x, y)

If τ does not depend on t, then we have a quasi-metric:

1. Reflexivity: τ(x, y) = 0 if and only if x = y,

2. Oriented triangle inequality: for any x, y, z ∊ X: τ(x, y) ≤ τ(x, z) + τ(z, y).

If a quasi-metric τ satisfies

3. Symmetry: for any x, y ∊ X: τ(x, y) = τ(y, x)

then (X, τ) is a usual metric space.

For a fixed t, in general case, a timetric τ_t(·) does not have to be a metric or even a quasi-metric.

Quick start

extern crate sr_rcd;

use std::time::Instant;

use sr_rcd::{    
    Point,
    TimetricSpace,
    rcd,
    RMPAtmo,
    Sensor
};

fn main() {
    let space = RMPAtmo::new(3, 3.14, &Point::new(vec![3.4, -5.6, 7.8]), 0.03, 1e-4);

    let sensorium = vec![
        Sensor::make(&Point::new(vec![19.894, -437.750, -455.413]), -244.940),
        Sensor::make(&Point::new(vec![216.980, -481.563, -45.759]), -153.382),
        Sensor::make(&Point::new(vec![85.351, -196.650, -422.263]), -300.745),
        Sensor::make(&Point::new(vec![61.788, 287.007, 384.448]), 383.750),
        Sensor::make(&Point::new(vec![-467.058, -97.973, -241.471]), 246.294),
        Sensor::make(&Point::new(vec![-234.397, -150.994, -429.143]), 10.128),
        Sensor::make(&Point::new(vec![-86.850, -215.597, 248.349]), 175.276),
        Sensor::make(&Point::new(vec![50.284, 460.754, -16.951]), 315.329),
        Sensor::make(&Point::new(vec![-294.371, -61.068, -427.287]), 82.082),
        Sensor::make(&Point::new(vec![475.401, 353.986, 62.538]), 242.594),
        Sensor::make(&Point::new(vec![263.802, -321.952, -408.792]), -480.919),
        Sensor::make(&Point::new(vec![-112.307, 493.069, -168.420]), 343.031)
    ];

    let start = Instant::now();
    match rcd(&space, &sensorium, &Point::new_default(space.dim()), 1000.0, 0.1, 0.001, false) {
        Ok(z) => {
            println!("Time: {:.3} sec", start.elapsed().as_secs_f64());
            println!("RCD-approximated source: {:.4?}", &z);
        },
        Err(s) => println!("ERROR: {}", s)
    }
}

The execution result looks like

Iteration 1: 75 coverands, d = 3342.0306
Iteration 2: 428 coverands, d = 4013.1741
Iteration 3: 501 coverands, d = 2568.8301
Iteration 4: 87 coverands, d = 848.0291
Iteration 5: 59 coverands, d = 357.5475
Iteration 6: 44 coverands, d = 178.7737
Iteration 7: 20 coverands, d = 40.2004
Iteration 8: 13 coverands, d = 20.1002
Iteration 9: 15 coverands, d = 10.0501
Iteration 10: 15 coverands, d = 3.8286
Iteration 11: 17 coverands, d = 2.1544
Iteration 12: 16 coverands, d = 0.9355
Iteration 13: 17 coverands, d = 0.4677
Iteration 14: 16 coverands, d = 0.2339
Iteration 15: 17 coverands, d = 0.1365
Iteration 16: 15 coverands, d = 0.0585
Time: 4.348 sec
RCD-approximated source: [262.0173, -319.4046, -392.4598]

Time here was obtained with release profile; dev one is nearly 5 times slower.

See bin/random.rs and bin/manual.rs.

RCD

The algorithm is considered in PDF preprint; see also preceding another PDF preprint and doi:10.12697/ACUTM.2020.24.14 for simpler versions regarding quasi-metric and proper metric spaces respectively. (Or maybe this method was conceived long ago, in some paper from 1920–50?..) Basically, part of space is covered by balls, then the cover is iteratively refined by replacing each ball with its cover by balls of halved radius and discarding the ones such that certain "spread of pasts" at their centers is larger than their doubled radius, until the cover becomes small enough. It is a particular case of the so-called "exclude & enclose" approach.

Its implementation is in rcd.rs.

Custom spaces

In principle, the algorithm works in any b-timetric space under some assumptions, here there is only a normed space ℝ^m_p with periodically shifting atmosphere or constant wind, — see rmp.rs. As soon as TimetricSpace trait is implemented likewise for AnotherSpace, RCD can work with SR problem in AnotherSpace. Also, RCD itself does not require ftau() method of that trait, hence, if SR data is already present, this method can be stubbed.

Robustness

To simulate noise in measurements, add N_0,σ²-distributed error to each t_i by setting err_dev parameter of make_random_sr_problem() to σ. When other parameters take default values, in particular, the approximation precision δ = 0.1, for σ = δ/10 the algorithm fails in nearly 30% of runs as the cover becomes empty, and for σ = δ/5 the failure rate is nearly 90%.

Higher dimensions and stronger winds

In 4D (moreover 5D etc.) or with angular frequency of atmosphere shift increased from 0.03 to 0.07, the process almost always stucks for a long, on a general purpose computer of nowadays, time at 2nd iteration, with the number of coverands after 1st iteration in hundreds. Extensively, this can be circumvented by using a fast enough computer, but there may be other optimizations.

So, for higher dimensions or stronger winds this is mostly in theoretical realm for now...