sr_rcd/

lib.rs

/*!
# 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<sub>0</sub> of time emits the sound. The sound wave propagates through space and reaches known points **r**<sub>i</sub>, *sensors*, at known instants t<sub>i</sub>. The *sound ranging* (SR) problem, also called *source localization* or *time-difference-of-arrival* (TDOA) problem, is to determine **s** and t<sub>0</sub> from {**r**<sub>i</sub>, t<sub>i</sub>}.

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 ∊ ℝ: τ<sub>t</sub>(**x**, **y**) = 0 if and only if **x** = **y**,

2. *Forward timed triangle inequality:* for any t ∊ ℝ and any **x**, **y**, **z** ∊ X: τ<sub>t</sub>(**x**, **y**) ≤ τ<sub>t</sub>(**x**, **z**) + τ<sub>t + τ<sub>t</sub>(**x**, **z**)</sub>(**z**, **y**).

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

In SR context, τ<sub>t</sub>(**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 + τ<sub>t</sub>(**x**, **y**).

<hr>

**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 ∊ ℝ: τ̅<sub>t</sub>(**x**, **y**) = 0 if and only if **x** = **y**,

2. *Backward timed triangle inequality:* for any t ∊ ℝ and any **x**, **y**, **z** ∊ X: τ̅<sub>t</sub>(**x**, **y**) ≤ τ̅<sub>t</sub>(**z**, **y**) + τ̅<sub>t - τ̅<sub>t</sub>(**z**, **y**)</sub>(**x**, **z**).

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

In SR context, τ̅<sub>t</sub>(**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 - τ̅<sub>t</sub>(**x**, **y**).

<hr>

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

<hr>

The connection between f- and b-timetrics:

τ̅<sub>t</sub>(**x**, **y**) = τ<sub>t - τ̅<sub>t</sub>(**x**, **y**)</sub>(**x**, **y**)

τ<sub>t</sub>(**x**, **y**) = τ̅<sub>t + τ<sub>t</sub>(**x**, **y**)</sub>(**x**, **y**)

<hr>

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 τ<sub>t</sub>(·) does not have to be a metric or even a quasi-metric.

## Quick start

```ignore
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

```text
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](https://www.researchgate.net/publication/355344503_On_sound_ranging_in_time-dependent-metric_spaces); see also preceding [another PDF preprint](https://www.researchgate.net/publication/354236739_On_sound_ranging_in_quasi-metric_spaces) and [doi:10.12697/ACUTM.2020.24.14](https://doi.org/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 ℝ<sup>m</sup><sub>p</sub> 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<sub>0,σ<sup>2</sup></sub>-distributed error to each t<sub>i</sub> 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...
*/

mod point;
mod tmspace;
mod rcd;
mod rmp;

pub use point::Point;
pub use tmspace::TimetricSpace;
pub use rcd::{Sensor, rcd};
pub use rmp::{RMPAtmo, RMPWind};