terrustrial 0.1.0

A Rust library for geospatial statistics, variograms, and kriging.
Documentation 
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

use crate::geometry::support::Support;
use crate::variography::model_variograms::composite::CompositeVariogram;
use crate::variography::model_variograms::VariogramModel;

use faer::{Col, Mat};

use ultraviolet::{DVec3, DVec3x4};
use wide::f64x4;

pub struct SKGeneralBuilder;

impl SKGeneralBuilder {
    const NUM_LANES: usize = 4;
    pub fn build_cov_mat(
        cov_mat: &mut Mat<f64>,
        cond: &[Support],
        vgram: &CompositeVariogram,
        h_buffer: &mut Vec<DVec3>,
        pt_buffer: &mut Vec<DVec3>,
        var_buffer: &mut Vec<f64>,
        ind_buffer: &mut Vec<usize>,
    ) {
        // 1. Clear input buffers.
        h_buffer.clear();
        pt_buffer.clear();
        var_buffer.clear();
        ind_buffer.clear();

        ind_buffer.push(0);

        //1. Compute pair-wise distances for lower triangular pairs.
        for (i, p1) in cond.iter().enumerate() {
            for p2 in cond.iter().take(i + 1) {
                p1.dists_to_other(p2, h_buffer, pt_buffer);
                ind_buffer.push(h_buffer.len());
            }
        }

        // 2. Compute the var
        for chunk in h_buffer.chunks_exact(Self::NUM_LANES) {
            let &[h1, h2, h3, h4] = chunk else {
                unreachable!()
            };

            let vec_h = DVec3x4::new(
                f64x4::new([h1.x, h2.x, h3.x, h4.x]),
                f64x4::new([h1.y, h2.y, h3.y, h4.y]),
                f64x4::new([h1.z, h2.z, h3.z, h4.z]),
            );

            let var = vgram.covariogram_simd(vec_h);

            var_buffer.extend_from_slice(var.as_array_ref());
        }

        for h in h_buffer.chunks_exact(Self::NUM_LANES).remainder() {
            var_buffer.push(vgram.covariogram(*h));
        }

        // 3. Populate lower traingular elements
        let mut cnt = 0;
        for i in 0..cond.len() {
            for j in 0..i + 1 {
                let low_idx = ind_buffer[cnt];
                let high_idx = ind_buffer[cnt + 1];
                let vars = &var_buffer[low_idx..high_idx];

                let avg_var = vars.iter().copied().sum::<f64>() / vars.len() as f64;

                *cov_mat.get_mut(i, j) = avg_var;

                cnt += 1;
            }
        }
    }

    #[allow(clippy::too_many_arguments)]
    pub fn build_cov_vec(
        cov_vec: &mut Col<f64>,
        cond: &[Support],
        kriging_point: &Support,
        vgram: &CompositeVariogram,
        h_buffer: &mut Vec<DVec3>,
        pt_buffer: &mut Vec<DVec3>,
        var_buffer: &mut Vec<f64>,
        ind_buffer: &mut Vec<usize>,
    ) {
        // 1. Clear input buffers.
        h_buffer.clear();
        pt_buffer.clear();
        var_buffer.clear();
        ind_buffer.clear();

        ind_buffer.push(0);

        //1. Compute pair-wise distances for lower triangular pairs.
        for p1 in cond.iter() {
            kriging_point.dists_to_other(p1, h_buffer, pt_buffer);
            ind_buffer.push(h_buffer.len());
        }

        // 2. Compute the var
        for chunk in h_buffer.chunks_exact(Self::NUM_LANES) {
            let &[h1, h2, h3, h4] = chunk else {
                unreachable!()
            };

            let vec_h = DVec3x4::new(
                f64x4::new([h1.x, h2.x, h3.x, h4.x]),
                f64x4::new([h1.y, h2.y, h3.y, h4.y]),
                f64x4::new([h1.z, h2.z, h3.z, h4.z]),
            );

            let var = vgram.covariogram_simd(vec_h);

            var_buffer.extend_from_slice(var.as_array_ref());
        }

        for h in h_buffer.chunks_exact(Self::NUM_LANES).remainder() {
            var_buffer.push(vgram.covariogram(*h));
        }

        // 3. Populate lower traingular elements
        for i in 0..cond.len() {
            let low_idx = ind_buffer[i];
            let high_idx = ind_buffer[i + 1];
            let vars = &var_buffer[low_idx..high_idx];

            let avg_var = vars.iter().copied().sum::<f64>() / vars.len() as f64;

            *cov_vec.get_mut(i) = avg_var;
        }
    }
}