Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

What are typical problem sizes? #6

Merged
merged 5 commits into from
Oct 12, 2021
Merged

What are typical problem sizes? #6

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
4c7e115
Generate larger problems for benchmarks, presumably to increase reali…
adamreichold Sep 30, 2021
e913fdd
Do not parallelize the second loop level of the summator_incompr func…
adamreichold Sep 30, 2021
c688b82
Parallelize the second loop level of the krige functions
adamreichold Sep 30, 2021
8f58226
Switch loop nesting for summator incompr implementation to avoid repe…
adamreichold Oct 5, 2021
931faf3
Redo parallel implementation of summator_incompr
adamreichold Oct 8, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
51 changes: 12 additions & 39 deletions benches/gen_benchmark_inputs.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,38 +8,21 @@


def gen_field_summate(path, seed):
pos_no = 1000
mode_no = 100
pos_no = 100_000
mode_no = 1_000
x = np.linspace(0.0, 10.0, pos_no)
y = np.linspace(-5.0, 5.0, pos_no)
z = np.linspace(-6.0, 8.0, pos_no)

model = gs.Gaussian(dim=3, var=1.0, len_scale=1.0)

rm = RandMeth(model, mode_no, seed)
np.savetxt(path / 'field_bench_x.txt', x)
np.savetxt(path / 'field_bench_y.txt', y)
np.savetxt(path / 'field_bench_z.txt', z)
np.savetxt(path / 'field_bench_cov_samples.txt', rm._cov_sample)
np.savetxt(path / 'field_bench_z_1.txt', rm._z_1)
np.savetxt(path / 'field_bench_z_2.txt', rm._z_2)

pos_no = 100000
mode_no = 1000
x = np.linspace(0.0, 10.0, pos_no)
y = np.linspace(-5.0, 5.0, pos_no)
z = np.linspace(-6.0, 8.0, pos_no)

model = gs.Gaussian(dim=3, var=1.0, len_scale=1.0)

rm = RandMeth(model, mode_no, seed)

np.savetxt(path / 'field_vs_x.txt', x)
np.savetxt(path / 'field_vs_y.txt', y)
np.savetxt(path / 'field_vs_z.txt', z)
np.savetxt(path / 'field_vs_cov_samples.txt', rm._cov_sample)
np.savetxt(path / 'field_vs_z_1.txt', rm._z_1)
np.savetxt(path / 'field_vs_z_2.txt', rm._z_2)
np.savetxt(path / 'field_x.txt', x)
np.savetxt(path / 'field_y.txt', y)
np.savetxt(path / 'field_z.txt', z)
np.savetxt(path / 'field_cov_samples.txt', rm._cov_sample)
np.savetxt(path / 'field_z_1.txt', rm._z_1)
np.savetxt(path / 'field_z_2.txt', rm._z_2)

def gen_krige(path, seed):
def prepare_data(pos_no, cond_no):
Expand Down Expand Up @@ -67,21 +50,11 @@ def prepare_data(pos_no, cond_no):

return krige_mat, k_vec, krige_cond

pos_no = 10000
cond_no = 500
krige_mat, k_vec, krige_cond = prepare_data(pos_no, cond_no)

np.savetxt(path / 'krige_vs_krige_mat.txt', krige_mat)
np.savetxt(path / 'krige_vs_k_vec.txt', k_vec)
np.savetxt(path / 'krige_vs_krige_cond.txt', krige_cond)

pos_no = 1000
cond_no = 50
krige_mat, k_vec, krige_cond = prepare_data(pos_no, cond_no)
Comment on lines -70 to -80
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

See above comment

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The question about typical problem sizes is hard to answer, as you might have already thought.

I just checked a publication of mine, where I used a very old implementation in C++ of the summator_incompr function. I used a 2d grid of 4600 x 1800 ~ 8*10^6 cells, which I would consider a larger problem size (for a 2d problem). I also actually used 6400 modes (len of z1, z2), which is probably a bit much (nowadays I always use exactly 1000 modes). Running one calcuation took 48 minutes on JURECA.

I guess @MuellerSeb can make much better estimates on typical sizes of kriging problems, as he has much more experience with it.

And just for the record, in this publication, I estimated a variogram on an unstructured grid with 412 cells and a bin length of 28, which, computationally wise is rather small, but I think some variogram estimations are much smaller.
In this example, which uses non-synthetic data, I used 2000 data points and about 30 bins to estimate the variogram. But we actually only sampled a small part of the data due to computational costs.

But all in all, for small problem sizes, the run times are hardly noticeable, so I support the idea of using larger problem sizes and reducing the criterion sample sizes.

Regarding your final point to maybe only parallelize the outer most loops: as long as we are not working on 64 or more cores, that would probably at least be the more conservative solution and save some potential and surprising declines in performance.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

So

But all in all, for small problem sizes, the run times are hardly noticeable, so I support the idea of using larger problem sizes and reducing the criterion sample sizes.

would suggest, to just fix up the first commit to just drop the *_bench_* variants whereas

Regarding your final point to maybe only parallelize the outer most loops: as long as we are not working on 64 or more cores, that would probably at least be the more conservative solution and save some potential and surprising declines in performance.

would imply to drop the last commit, i.e. remove all parallelization except for the outermost loops.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

would imply to drop the last commit, i.e. remove all parallelization except for the outermost loops.

Will wait for @MuellerSeb to chime in before making that change.

krige_mat, k_vec, krige_cond = prepare_data(pos_no=10_000, cond_no=500)

np.savetxt(path / 'krige_bench_krige_mat.txt', krige_mat)
np.savetxt(path / 'krige_bench_k_vec.txt', k_vec)
np.savetxt(path / 'krige_bench_krige_cond.txt', krige_cond)
np.savetxt(path / 'krige_krige_mat.txt', krige_mat)
np.savetxt(path / 'krige_k_vec.txt', k_vec)
np.savetxt(path / 'krige_krige_cond.txt', krige_cond)


if __name__ == '__main__':
Expand Down
18 changes: 9 additions & 9 deletions benches/main.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -40,16 +40,16 @@ fn read_2d_from_file(file_path: &Path) -> Array2<f64> {
pub fn field_benchmark(c: &mut Criterion) {
let path = Path::new("benches/input");

let x = read_1d_from_file(&path.join("field_bench_x.txt"));
let y = read_1d_from_file(&path.join("field_bench_y.txt"));
let z = read_1d_from_file(&path.join("field_bench_z.txt"));
let x = read_1d_from_file(&path.join("field_x.txt"));
let y = read_1d_from_file(&path.join("field_y.txt"));
let z = read_1d_from_file(&path.join("field_z.txt"));

let pos = stack![Axis(0), x, y, z];

let cov_samples = read_2d_from_file(&path.join("field_bench_cov_samples.txt"));
let cov_samples = read_2d_from_file(&path.join("field_cov_samples.txt"));

let z_1 = read_1d_from_file(&path.join("field_bench_z_1.txt"));
let z_2 = read_1d_from_file(&path.join("field_bench_z_2.txt"));
let z_1 = read_1d_from_file(&path.join("field_z_1.txt"));
let z_2 = read_1d_from_file(&path.join("field_z_2.txt"));

c.bench_function("field summate", |b| {
b.iter(|| summator(cov_samples.view(), z_1.view(), z_2.view(), pos.view()))
Expand All @@ -63,9 +63,9 @@ pub fn field_benchmark(c: &mut Criterion) {
pub fn krige_benchmark(c: &mut Criterion) {
let path = Path::new("benches/input");

let krige_mat = read_2d_from_file(&path.join("krige_bench_krige_mat.txt"));
let k_vec = read_2d_from_file(&path.join("krige_bench_k_vec.txt"));
let krige_cond = read_1d_from_file(&path.join("krige_bench_krige_cond.txt"));
let krige_mat = read_2d_from_file(&path.join("krige_krige_mat.txt"));
let k_vec = read_2d_from_file(&path.join("krige_k_vec.txt"));
let krige_cond = read_1d_from_file(&path.join("krige_krige_cond.txt"));

c.bench_function("krige error", |b| {
b.iter(|| {
Expand Down
83 changes: 40 additions & 43 deletions src/field.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
use ndarray::{Array1, Array2, ArrayView1, ArrayView2, Zip};
use ndarray::{Array1, Array2, ArrayView1, ArrayView2, Axis, Zip};
use rayon::prelude::*;

pub fn summator(
Expand All @@ -7,11 +7,11 @@ pub fn summator(
z2: ArrayView1<'_, f64>,
pos: ArrayView2<'_, f64>,
) -> Array1<f64> {
assert!(cov_samples.shape()[0] == pos.shape()[0]);
assert!(cov_samples.shape()[1] == z1.shape()[0]);
assert!(z1.shape()[0] == z2.shape()[0]);
assert_eq!(cov_samples.dim().0, pos.dim().0);
assert_eq!(cov_samples.dim().1, z1.dim());
assert_eq!(cov_samples.dim().1, z2.dim());

let mut summed_modes = Array1::<f64>::zeros(pos.shape()[1]);
let mut summed_modes = Array1::<f64>::zeros(pos.dim().1);

Zip::from(&mut summed_modes)
.and(pos.columns())
Expand All @@ -35,54 +35,51 @@ pub fn summator_incompr(
z2: ArrayView1<'_, f64>,
pos: ArrayView2<'_, f64>,
) -> Array2<f64> {
assert!(cov_samples.shape()[0] == pos.shape()[0]);
assert!(cov_samples.shape()[1] == z1.shape()[0]);
assert!(z1.shape()[0] == z2.shape()[0]);

let dim = pos.shape()[0];

let mut summed_modes = Array2::<f64>::zeros(pos.dim());
assert_eq!(cov_samples.dim().0, pos.dim().0);
assert_eq!(cov_samples.dim().1, z1.dim());
assert_eq!(cov_samples.dim().1, z2.dim());

// unit vector in x dir.
let mut e1 = Array1::<f64>::zeros(dim);
let mut e1 = Array1::<f64>::zeros(pos.dim().0);
e1[0] = 1.0;

Zip::from(pos.columns())
.and(summed_modes.columns_mut())
.par_for_each(|pos, mut summed_modes| {
let sum = Zip::from(cov_samples.columns())
.and(z1)
.and(z2)
.into_par_iter()
.fold(
|| Array1::<f64>::zeros(dim),
|mut sum, (cov_samples, z1, z2)| {
let k_2 = cov_samples.dot(&cov_samples);
// FIXME: Use Zip::from(cov_samples.columns()).and(z1).and(z1) after
// https://github.com/rust-ndarray/ndarray/pull/1081 is merged.
cov_samples
.axis_iter(Axis(1))
.into_par_iter()
.zip(z1.axis_iter(Axis(0)))
.zip(z2.axis_iter(Axis(0)))
.with_min_len(100)
.fold(
|| Array2::<f64>::zeros(pos.dim()),
|mut summed_modes, ((cov_samples, z1), z2)| {
let k_2 = cov_samples[0] / cov_samples.dot(&cov_samples);
let z1 = z1.into_scalar();
let z2 = z2.into_scalar();

Zip::from(pos.columns())
.and(summed_modes.columns_mut())
.par_for_each(|pos, mut summed_modes| {
let phase = cov_samples.dot(&pos);
let z12 = z1 * phase.cos() + z2 * phase.sin();

Zip::from(&mut sum)
Zip::from(&mut summed_modes)
.and(&e1)
.and(cov_samples)
.for_each(|sum, e1, cs| {
let proj = *e1 - cs * cov_samples[0] / k_2;
*sum += proj * (z1 * phase.cos() + z2 * phase.sin());
*sum += (*e1 - cs * k_2) * z12;
});

sum
},
)
.reduce(
|| Array1::<f64>::zeros(dim),
|mut lhs, rhs| {
lhs += &rhs;
lhs
},
);

summed_modes.assign(&sum);
});

summed_modes
});

summed_modes
},
)
.reduce_with(|mut lhs, rhs| {
lhs += &rhs;
lhs
})
.unwrap()
}

#[cfg(test)]
Expand Down
41 changes: 31 additions & 10 deletions src/krige.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
use ndarray::{Array1, ArrayView1, ArrayView2, Zip};
use rayon::prelude::*;

pub fn calculator_field_krige_and_variance(
krig_mat: ArrayView2<'_, f64>,
Expand All @@ -16,14 +17,30 @@ pub fn calculator_field_krige_and_variance(
.and(error.view_mut())
.and(krig_vecs.columns())
.par_for_each(|f, e, v_col| {
Zip::from(cond)
let acc = Zip::from(cond)
.and(v_col)
.and(krig_mat.columns())
.for_each(|c, v, m_row| {
let krig_fac = m_row.dot(&v_col);
*f += c * krig_fac;
*e += v * krig_fac;
});
.into_par_iter()
.fold(
|| (0.0, 0.0),
|mut acc, (c, v, m_row)| {
let krig_fac = m_row.dot(&v_col);
acc.0 += c * krig_fac;
acc.1 += v * krig_fac;
acc
},
)
.reduce(
|| (0.0, 0.0),
|mut lhs, rhs| {
lhs.0 += rhs.0;
lhs.1 += rhs.1;
lhs
},
);

*f = acc.0;
*e = acc.1;
});

(field, error)
Expand All @@ -43,12 +60,16 @@ pub fn calculator_field_krige(
Zip::from(field.view_mut())
.and(krig_vecs.columns())
.par_for_each(|f, v_col| {
Zip::from(cond)
let acc = Zip::from(cond)
.and(krig_mat.columns())
.for_each(|c, m_row| {
.into_par_iter()
.map(|(c, m_row)| {
let krig_fac = m_row.dot(&v_col);
*f += c * krig_fac;
});
c * krig_fac
})
.sum();

*f = acc;
});

field
Expand Down