@@ -20,25 +20,25 @@ pub struct SctArgs {
20
20
pub num_max : usize ,
21
21
// FIXME: this doc comment can be improved
22
22
/// Radius in which OI will be reused. Unit: m
23
- pub inner_radius : f32 ,
23
+ pub inner_radius : f64 ,
24
24
/// Radius for computing OI and background. Unit: m
25
- pub outer_radius : f32 ,
25
+ pub outer_radius : f64 ,
26
26
/// The number of iterations of SCT to perform before returning.
27
27
pub num_iterations : u32 ,
28
28
/// Minimum number of observations to compute vertical profile.
29
29
pub num_min_prof : usize ,
30
30
/// Minimum elevation difference to compute vertical profile. Unit: m
31
- pub min_elev_diff : f32 ,
31
+ pub min_elev_diff : f64 ,
32
32
/// Minimum horizontal decorrelation length. Unit: m
33
- pub min_horizontal_scale : f32 ,
33
+ pub min_horizontal_scale : f64 ,
34
34
/// Vertical decorrelation length. Unit: m
35
- pub vertical_scale : f32 ,
35
+ pub vertical_scale : f64 ,
36
36
/// Positive deviation allowed. Unit: σ (standard deviations)
37
- pub pos : SingleOrVec < f32 > ,
37
+ pub pos : SingleOrVec < f64 > ,
38
38
/// Negative deviation allowed. Unit: σ (standard deviations)
39
- pub neg : SingleOrVec < f32 > ,
39
+ pub neg : SingleOrVec < f64 > ,
40
40
/// Ratio of observation error variance to background variance.
41
- pub eps2 : SingleOrVec < f32 > ,
41
+ pub eps2 : SingleOrVec < f64 > ,
42
42
}
43
43
44
44
fn subset < T : Copy > ( array : & [ T ] , indices : & [ usize ] ) -> Vec < T > {
@@ -53,16 +53,16 @@ fn subset<T: Copy>(array: &[T], indices: &[usize]) -> Vec<T> {
53
53
}
54
54
55
55
fn compute_vertical_profile_theil_sen (
56
- elevs : & [ f32 ] ,
57
- values : & [ f32 ] ,
56
+ elevs : & [ f64 ] ,
57
+ values : & [ f64 ] ,
58
58
num_min_prof : usize ,
59
- min_elev_diff : f32 ,
60
- ) -> Vec < f32 > {
59
+ min_elev_diff : f64 ,
60
+ ) -> Vec < f64 > {
61
61
let n = values. len ( ) ;
62
62
63
63
// Starting value guesses
64
- let gamma: f32 = -0.0065 ;
65
- let mean_t: f32 = values. iter ( ) . sum :: < f32 > ( ) / n as f32 ; // should this be f64?
64
+ let gamma: f64 = -0.0065 ;
65
+ let mean_t: f64 = values. iter ( ) . sum :: < f64 > ( ) / n as f64 ; // should this be f64?
66
66
67
67
// special case when all observations have the same elevation
68
68
if elevs. iter ( ) . min_by ( |a, b| a. total_cmp ( b) ) == elevs. iter ( ) . max_by ( |a, b| a. total_cmp ( b) ) {
@@ -81,7 +81,7 @@ fn compute_vertical_profile_theil_sen(
81
81
gamma
82
82
} else {
83
83
let nm = n * ( n - 1 ) / 2 ;
84
- let mut m: Vec < f32 > = Vec :: with_capacity ( nm) ;
84
+ let mut m: Vec < f64 > = Vec :: with_capacity ( nm) ;
85
85
for i in 0 ..( n - 1 ) {
86
86
for j in ( i + 1 ) ..n {
87
87
m. push ( if ( elevs[ i] - elevs[ j] ) . abs ( ) < 1. {
@@ -93,7 +93,7 @@ fn compute_vertical_profile_theil_sen(
93
93
}
94
94
compute_quantile ( 0.5 , & m)
95
95
} ;
96
- let q: Vec < f32 > = values
96
+ let q: Vec < f64 > = values
97
97
. iter ( )
98
98
. zip ( elevs)
99
99
. map ( |( val, elev) | val - m_median * elev)
@@ -107,8 +107,8 @@ fn compute_vertical_profile_theil_sen(
107
107
}
108
108
109
109
// TODO: replace assertions with errors or remove them
110
- fn compute_quantile ( quantile : f32 , array : & [ f32 ] ) -> f32 {
111
- let mut new_array: Vec < f32 > = array
110
+ fn compute_quantile ( quantile : f64 , array : & [ f64 ] ) -> f64 {
111
+ let mut new_array: Vec < f64 > = array
112
112
. iter ( )
113
113
. copied ( )
114
114
. filter ( |x| util:: is_valid ( * x) )
@@ -120,12 +120,12 @@ fn compute_quantile(quantile: f32, array: &[f32]) -> f32 {
120
120
assert ! ( n > 0 ) ;
121
121
122
122
// get the quantile from the sorted array
123
- let lower_index = ( quantile * ( n - 1 ) as f32 ) . floor ( ) as usize ;
124
- let upper_index = ( quantile * ( n - 1 ) as f32 ) . ceil ( ) as usize ;
123
+ let lower_index = ( quantile * ( n - 1 ) as f64 ) . floor ( ) as usize ;
124
+ let upper_index = ( quantile * ( n - 1 ) as f64 ) . ceil ( ) as usize ;
125
125
let lower_value = new_array[ lower_index] ;
126
126
let upper_value = new_array[ upper_index] ;
127
- let lower_quantile = lower_index as f32 / ( n - 1 ) as f32 ;
128
- let upper_quantile = upper_index as f32 / ( n - 1 ) as f32 ;
127
+ let lower_quantile = lower_index as f64 / ( n - 1 ) as f64 ;
128
+ let upper_quantile = upper_index as f64 / ( n - 1 ) as f64 ;
129
129
let exact_q = if lower_index == upper_index {
130
130
lower_value
131
131
} else {
@@ -142,16 +142,16 @@ fn compute_quantile(quantile: f32, array: &[f32]) -> f32 {
142
142
exact_q
143
143
}
144
144
145
- fn invert_matrix ( input : & Mat < f32 > ) -> Mat < f32 > {
145
+ fn invert_matrix ( input : & Mat < f64 > ) -> Mat < f64 > {
146
146
let lu = input. partial_piv_lu ( ) ;
147
147
lu. inverse ( )
148
148
}
149
149
150
150
fn remove_flagged < ' a > (
151
151
neighbours : Vec < & ' a SpatialPoint > ,
152
- distances : Vec < f32 > ,
152
+ distances : Vec < f64 > ,
153
153
flags : & [ Flag ] ,
154
- ) -> ( Vec < & ' a SpatialPoint > , Vec < f32 > ) {
154
+ ) -> ( Vec < & ' a SpatialPoint > , Vec < f64 > ) {
155
155
let vec_length = neighbours. len ( ) ;
156
156
let mut neighbours_new = Vec :: with_capacity ( vec_length) ;
157
157
let mut distances_new = Vec :: with_capacity ( vec_length) ;
@@ -214,7 +214,7 @@ fn remove_flagged<'a>(
214
214
/// \* optional, ou = Unit of the observation, σ = Standard deviations
215
215
#[ allow( clippy:: too_many_arguments) ]
216
216
pub fn sct (
217
- data : & [ Option < f32 > ] ,
217
+ data : & [ Option < f64 > ] ,
218
218
rtree : & SpatialTree ,
219
219
args : & SctArgs ,
220
220
obs_to_check : Option < & [ bool ] > ,
@@ -369,7 +369,7 @@ pub fn sct(
369
369
remove_flagged ( neighbours_unfiltered, distances_unfiltered, & flags) ;
370
370
371
371
if neighbours. len ( ) > args. num_max {
372
- let mut pairs: Vec < ( & SpatialPoint , f32 ) > =
372
+ let mut pairs: Vec < ( & SpatialPoint , f64 ) > =
373
373
neighbours. into_iter ( ) . zip ( distances. into_iter ( ) ) . collect ( ) ;
374
374
pairs. sort_by ( |a, b| a. 1 . partial_cmp ( & b. 1 ) . unwrap_or ( std:: cmp:: Ordering :: Equal ) ) ;
375
375
@@ -394,7 +394,7 @@ pub fn sct(
394
394
let values_box = subset ( data, & neighbour_indices)
395
395
. into_iter ( )
396
396
. map ( |v| v. unwrap ( ) )
397
- . collect :: < Vec < f32 > > ( ) ;
397
+ . collect :: < Vec < f64 > > ( ) ;
398
398
let eps2_box = match & args. eps2 {
399
399
SingleOrVec :: Single ( eps2_value) => SingleOrVec :: Single ( * eps2_value) ,
400
400
SingleOrVec :: Vec ( eps2_vec) => {
@@ -411,15 +411,15 @@ pub fn sct(
411
411
args. min_elev_diff ,
412
412
) ;
413
413
414
- let disth: Mat < f32 > = Mat :: from_fn ( box_size, box_size, |i, j| {
414
+ let disth: Mat < f64 > = Mat :: from_fn ( box_size, box_size, |i, j| {
415
415
// TODO: remove this unwrap
416
416
util:: calc_distance ( lats_box[ i] , lons_box[ i] , lats_box[ j] , lons_box[ j] ) . unwrap ( )
417
417
} ) ;
418
- let distz: Mat < f32 > = Mat :: from_fn ( box_size, box_size, |i, j| {
418
+ let distz: Mat < f64 > = Mat :: from_fn ( box_size, box_size, |i, j| {
419
419
( elevs_box[ i] - elevs_box[ j] ) . abs ( )
420
420
} ) ;
421
421
// TODO: remove dh, and just reduce straight into dh_mean?
422
- let dh: Vec < f32 > = ( 0 ..box_size)
422
+ let dh: Vec < f64 > = ( 0 ..box_size)
423
423
. map ( |i| {
424
424
let mut dh_vector = Vec :: with_capacity ( box_size - 1 ) ;
425
425
for j in 0 ..box_size {
@@ -431,13 +431,13 @@ pub fn sct(
431
431
} )
432
432
. collect ( ) ;
433
433
434
- let dh_mean: f32 = args
434
+ let dh_mean: f64 = args
435
435
. min_horizontal_scale
436
- . max ( dh. into_iter ( ) . sum :: < f32 > ( ) / box_size as f32 ) ;
436
+ . max ( dh. into_iter ( ) . sum :: < f64 > ( ) / box_size as f64 ) ;
437
437
438
- let mut s: Mat < f32 > = Mat :: from_fn ( box_size, box_size, |i, j| {
439
- let value = ( -0.5 * ( disth. get ( i, j) / dh_mean) . powi ( 2 )
440
- - 0.5 * ( distz. get ( i, j) / args. vertical_scale ) . powi ( 2 ) )
438
+ let mut s: Mat < f64 > = Mat :: from_fn ( box_size, box_size, |i, j| {
439
+ let value = ( -0.5 * ( disth. read ( i, j) / dh_mean) . powi ( 2 )
440
+ - 0.5 * ( distz. read ( i, j) / args. vertical_scale ) . powi ( 2 ) )
441
441
. exp ( ) ;
442
442
// weight the diagonal?? (0.5 default)
443
443
if i == j {
@@ -448,7 +448,7 @@ pub fn sct(
448
448
} ) ;
449
449
450
450
// difference between actual temp and temp from vertical profile
451
- let d: Vec < f32 > = ( 0 ..box_size)
451
+ let d: Vec < f64 > = ( 0 ..box_size)
452
452
. map ( |i| values_box[ i] - vertical_profile[ i] )
453
453
. collect ( ) ;
454
454
@@ -465,22 +465,22 @@ pub fn sct(
465
465
elem. sub_assign ( eps2_box. index ( i) ) ;
466
466
}
467
467
468
- let s_inv_d: Vec < f32 > = ( 0 ..box_size)
469
- . map ( |i| ( 0 ..box_size) . map ( |j| s_inv. get ( i, j) * d[ j] ) . sum ( ) )
468
+ let s_inv_d: Vec < f64 > = ( 0 ..box_size)
469
+ . map ( |i| ( 0 ..box_size) . map ( |j| s_inv. read ( i, j) * d[ j] ) . sum ( ) )
470
470
. collect ( ) ;
471
471
472
- let ares_temp: Vec < f32 > = ( 0 ..box_size)
473
- . map ( |i| ( 0 ..box_size) . map ( |j| s. get ( i, j) * s_inv_d[ j] ) . sum ( ) )
472
+ let ares_temp: Vec < f64 > = ( 0 ..box_size)
473
+ . map ( |i| ( 0 ..box_size) . map ( |j| s. read ( i, j) * s_inv_d[ j] ) . sum ( ) )
474
474
. collect ( ) ;
475
475
476
- let z_inv: Vec < f32 > = ( 0 ..box_size) . map ( |i| 1. / s_inv. get ( i, i) ) . collect ( ) ;
476
+ let z_inv: Vec < f64 > = ( 0 ..box_size) . map ( |i| 1. / s_inv. read ( i, i) ) . collect ( ) ;
477
477
478
- let ares: Vec < f32 > = ( 0 ..box_size) . map ( |i| ares_temp[ i] - d[ i] ) . collect ( ) ;
478
+ let ares: Vec < f64 > = ( 0 ..box_size) . map ( |i| ares_temp[ i] - d[ i] ) . collect ( ) ;
479
479
480
- let cvres: Vec < f32 > = ( 0 ..box_size) . map ( |i| -1. * z_inv[ i] * s_inv_d[ i] ) . collect ( ) ;
480
+ let cvres: Vec < f64 > = ( 0 ..box_size) . map ( |i| -1. * z_inv[ i] * s_inv_d[ i] ) . collect ( ) ;
481
481
482
- let sig2o = 0.01_f32
483
- . max ( ( 0 ..box_size) . map ( |i| d[ i] * -1. * ares[ i] ) . sum :: < f32 > ( ) / box_size as f32 ) ;
482
+ let sig2o = 0.01_f64
483
+ . max ( ( 0 ..box_size) . map ( |i| d[ i] * -1. * ares[ i] ) . sum :: < f64 > ( ) / box_size as f64 ) ;
484
484
485
485
let curr = i;
486
486
for i in 0 ..box_size {
@@ -493,7 +493,7 @@ pub fn sct(
493
493
}
494
494
let dist = distances[ i] ;
495
495
if dist <= args. inner_radius {
496
- let pog: f32 = cvres[ i] * ares[ i] / sig2o;
496
+ let pog: f64 = cvres[ i] * ares[ i] / sig2o;
497
497
assert ! ( util:: is_valid( pog) ) ;
498
498
prob_gross_error[ index] = pog. max ( prob_gross_error[ index] ) ;
499
499
if ( cvres[ i] < 0. && pog > * args. pos . index ( index) )
@@ -534,7 +534,7 @@ pub fn sct_cache(
534
534
535
535
for i in ( cache. num_leading_points as usize ) ..( series_len - cache. num_trailing_points as usize )
536
536
{
537
- let timeslice: Vec < Option < f32 > > = cache. data . iter ( ) . map ( |ts| ts. values [ i] ) . collect ( ) ;
537
+ let timeslice: Vec < Option < f64 > > = cache. data . iter ( ) . map ( |ts| ts. values [ i] ) . collect ( ) ;
538
538
let spatial_result = sct ( & timeslice, & cache. rtree , args, obs_to_check) ?;
539
539
540
540
for i in 0 ..spatial_result. len ( ) {
@@ -584,9 +584,9 @@ mod tests {
584
584
sct(
585
585
& vec![ Some ( 1. ) ; N ] ,
586
586
& SpatialTree :: from_latlons(
587
- ( 0 ..N ) . map( |i| ( ( i as f32 ) . powi( 2 ) * 0.001 ) % 1. ) . collect( ) ,
587
+ ( 0 ..N ) . map( |i| ( ( i as f64 ) . powi( 2 ) * 0.001 ) % 1. ) . collect( ) ,
588
588
( 0 ..N )
589
- . map( |i| ( ( i as f32 + 1. ) . powi( 2 ) * 0.001 ) % 1. )
589
+ . map( |i| ( ( i as f64 + 1. ) . powi( 2 ) * 0.001 ) % 1. )
590
590
. collect( ) ,
591
591
vec![ 1. ; N ] ,
592
592
) ,
0 commit comments