sgrust 0.8.6

A sparse grid library written in Rust.
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
135

use crate::const_generic::{iterators::grid_iterator::{GridIteratorT, HashMapGridIterator}, storage::SparseGridData};

use super::refinement::RefinementFunctor;

/// Expands the kept-node set to include the level-zero boundary closure required by
/// boundary hierarchization sweeps.
///
/// For every kept node and every dimension, both left/right level-zero projections in that
/// dimension must exist on the active fiber. Applying this recursively yields the minimal
/// boundary skeleton, e.g. `2^D` corners for a constant boundary grid.
fn add_required_boundary_closure<const D: usize>(iterator: &mut HashMapGridIterator<D>, storage: &SparseGridData<D>, keep: &mut [bool])
{
    let mut changed = true;
    while changed
    {
        changed = false;
        for seq in 0..storage.len()
        {
            if !keep[seq]
            {
                continue;
            }
            let point = storage.nodes()[seq];
            for dim in 0..D
            {
                iterator.set_index(point);
                if iterator.reset_to_left_level_zero(dim)
                {
                    let boundary_seq = iterator.index().unwrap();
                    if !keep[boundary_seq]
                    {
                        keep[boundary_seq] = true;
                        changed = true;
                    }
                }
                iterator.set_index(point);
                if iterator.reset_to_right_level_zero(dim)
                {
                    let boundary_seq = iterator.index().unwrap();
                    if !keep[boundary_seq]
                    {
                        keep[boundary_seq] = true;
                        changed = true;
                    }
                }
            }
        }
    }
}

fn keep_mask_to_indices(keep: &[bool]) -> Vec<usize>
{
    let mut kept_points = Vec::with_capacity(keep.iter().filter(|&&keep| keep).count());
    for (seq, &is_kept) in keep.iter().enumerate()
    {
        if is_kept
        {
            kept_points.push(seq);
        }
    }
    kept_points
}

/// Coarsen the grid by removing leaf nodes with low error estimates.
/// 
/// # Arguments
/// * `storage` - The sparse grid storage to coarsen
/// * `functor` - Refinement functor that computes error estimates
/// * `alpha` - Hierarchical surplus coefficients
/// * `values` - Function values at grid points
/// * `threshold` - Nodes with error below this are candidates for removal
/// * `remove_boundary` - If true, also remove boundary nodes that are no longer needed
///
/// # Returns
/// Vec containing the indices of points that were kept, in original storage order.
pub(crate) fn coarsen<const D: usize, const DIM_OUT: usize>( storage: &mut SparseGridData<D>, functor: &dyn RefinementFunctor<D, DIM_OUT>, alpha: &[[f64; DIM_OUT]], values: &[[f64; DIM_OUT]], threshold: f64, remove_boundary: bool) -> Vec<usize>
{
    let mut iterator = HashMapGridIterator::new(storage);
    let mut keep = vec![false; storage.len()];
    iterator.reset_to_level_zero();
    let zero_index = iterator.index().unwrap();
    let values = functor.eval(storage.points(), alpha, values);
    
    // First pass: determine which inner nodes to keep
    for (seq, (point, r)) in storage.nodes().iter().zip(values.iter()).enumerate()
    {
        if point.is_inner_point()
        {
            if point.is_leaf() && seq != zero_index
            {
                if *r >= threshold
                {
                    keep[seq] = true;
                }
            }
            else
            {
                // Keep non-leaf inner nodes and zero index
                keep[seq] = true;
            }
        }
    }
    
    if remove_boundary
    {
        // Keep structurally significant boundary nodes before adding the required closure.
        for (seq, (point, r)) in storage.nodes().iter().zip(values.iter()).enumerate()
        {
            if !point.is_inner_point()
            {
                // Always keep zero index, if not a leaf (has children that are kept), or boundary leaf nodes with significant surplus
                if seq == zero_index || !point.is_leaf() || *r >= threshold
                {
                    keep[seq] = true;
                }
                // Otherwise, this boundary node can be removed
            }
        }
        add_required_boundary_closure(&mut iterator, storage, &mut keep);
    }
    else
    {
        // Keep all boundary nodes
        for (seq, point) in storage.nodes().iter().enumerate()
        {
            if !point.is_inner_point()
            {
                keep[seq] = true;
            }
        }
    }
    let kept_points = keep_mask_to_indices(&keep);
    storage.remove(&kept_points);
    kept_points
}