scirs2_optimize/parallel/
mod.rs1use ndarray::{Array1, ArrayView1};
25use scirs2_core::parallel_ops::*;
26
27#[cfg(feature = "parallel")]
29use scirs2_core::parallel_ops::{IntoParallelIterator, ParallelIterator};
30
31#[derive(Debug, Clone)]
33pub struct ParallelOptions {
34 pub num_workers: Option<usize>,
36
37 pub min_parallel_size: usize,
39
40 pub chunk_size: usize,
42
43 pub parallel_evaluations: bool,
45
46 pub parallel_gradient: bool,
48}
49
50impl Default for ParallelOptions {
51 fn default() -> Self {
52 ParallelOptions {
53 num_workers: None,
54 min_parallel_size: 10,
55 chunk_size: 1,
56 parallel_evaluations: true,
57 parallel_gradient: true,
58 }
59 }
60}
61
62#[allow(dead_code)]
67pub fn parallel_finite_diff_gradient<F>(
68 f: F,
69 x: ArrayView1<f64>,
70 options: &ParallelOptions,
71) -> Array1<f64>
72where
73 F: Fn(&ArrayView1<f64>) -> f64 + Sync,
74{
75 let n = x.len();
76 let eps = 1e-8;
77
78 if n < options.min_parallel_size || !options.parallel_gradient {
80 return sequential_finite_diff_gradient(f, x, eps);
81 }
82
83 compute_parallel_gradient(&f, x, eps)
85}
86
87#[allow(dead_code)]
89fn sequential_finite_diff_gradient<F>(f: F, x: ArrayView1<f64>, eps: f64) -> Array1<f64>
90where
91 F: Fn(&ArrayView1<f64>) -> f64,
92{
93 let n = x.len();
94 let mut grad = Array1::zeros(n);
95 let fx = f(&x);
96
97 for i in 0..n {
98 let mut x_plus = x.to_owned();
99 x_plus[i] += eps;
100 let fx_plus = f(&x_plus.view());
101 grad[i] = (fx_plus - fx) / eps;
102 }
103
104 grad
105}
106
107#[cfg(feature = "parallel")]
109#[allow(dead_code)]
110fn compute_parallel_gradient<F>(f: &F, x: ArrayView1<f64>, eps: f64) -> Array1<f64>
111where
112 F: Fn(&ArrayView1<f64>) -> f64 + Sync,
113{
114 let n = x.len();
115 let fx = f(&x);
116
117 let grad_vec: Vec<f64> = (0..n)
119 .into_par_iter()
120 .map(|i| {
121 let mut x_plus = x.to_owned();
122 x_plus[i] += eps;
123 let fx_plus = f(&x_plus.view());
124 (fx_plus - fx) / eps
125 })
126 .collect();
127
128 Array1::from_vec(grad_vec)
129}
130
131#[cfg(not(feature = "parallel"))]
133#[allow(dead_code)]
134fn compute_parallel_gradient<F>(f: &F, x: ArrayView1<f64>, eps: f64) -> Array1<f64>
135where
136 F: Fn(&ArrayView1<f64>) -> f64,
137{
138 sequential_finite_diff_gradient(|x| f(x), x, eps)
139}
140
141#[allow(dead_code)]
145pub fn parallel_evaluate_batch<F>(
146 f: F,
147 points: &[Array1<f64>],
148 options: &ParallelOptions,
149) -> Vec<f64>
150where
151 F: Fn(&ArrayView1<f64>) -> f64 + Sync,
152{
153 if points.len() < options.min_parallel_size || !options.parallel_evaluations {
154 points.iter().map(|x| f(&x.view())).collect()
156 } else {
157 #[cfg(feature = "parallel")]
159 {
160 points.par_iter().map(|x| f(&x.view())).collect()
161 }
162 #[cfg(not(feature = "parallel"))]
163 {
164 points.iter().map(|x| f(&x.view())).collect()
165 }
166 }
167}
168
169#[allow(dead_code)] pub struct ParallelMultiStart<F> {
174 objective: F,
175 bounds: Vec<(f64, f64)>,
176 options: ParallelOptions,
177}
178
179impl<F> ParallelMultiStart<F>
180where
181 F: Fn(&ArrayView1<f64>) -> f64 + Sync + Send,
182{
183 pub fn new(objective: F, bounds: Vec<(f64, f64)>, options: ParallelOptions) -> Self {
185 ParallelMultiStart {
186 objective,
187 bounds,
188 options,
189 }
190 }
191
192 pub fn run<O>(&self, starting_points: Vec<Array1<f64>>, optimizer: O) -> Vec<OptimizationResult>
194 where
195 O: Fn(&Array1<f64>, &F) -> OptimizationResult + Sync,
196 {
197 if starting_points.len() < self.options.min_parallel_size {
198 starting_points
200 .iter()
201 .map(|x0| optimizer(x0, &self.objective))
202 .collect()
203 } else {
204 #[cfg(feature = "parallel")]
206 {
207 starting_points
208 .par_iter()
209 .map(|x0| optimizer(x0, &self.objective))
210 .collect()
211 }
212 #[cfg(not(feature = "parallel"))]
213 {
214 starting_points
215 .iter()
216 .map(|x0| optimizer(x0, &self.objective))
217 .collect()
218 }
219 }
220 }
221
222 pub fn best_result(results: &[OptimizationResult]) -> Option<&OptimizationResult> {
224 results
225 .iter()
226 .min_by(|a, b| a.function_value.partial_cmp(&b.function_value).unwrap())
227 }
228}
229
230#[derive(Debug, Clone)]
232pub struct OptimizationResult {
233 pub x: Array1<f64>,
234 pub function_value: f64,
235 pub nit: usize,
236 pub success: bool,
237}
238
239#[allow(dead_code)]
241pub fn parallel_finite_diff_hessian<F>(
242 f: &F,
243 x: ArrayView1<f64>,
244 gradient: Option<&Array1<f64>>,
245 options: &ParallelOptions,
246) -> Array1<f64>
247where
249 F: Fn(&ArrayView1<f64>) -> f64 + Sync,
250{
251 let n = x.len();
252 let eps = 1e-8;
253
254 if n < options.min_parallel_size {
256 return sequential_finite_diff_hessian(f, x, gradient, eps);
257 }
258
259 let _grad = match gradient {
261 Some(g) => g.clone(),
262 None => parallel_finite_diff_gradient(f, x, options),
263 };
264
265 let num_elements = n * (n + 1) / 2;
267
268 #[cfg(feature = "parallel")]
270 let hessian_elements: Vec<f64> = (0..num_elements)
271 .into_par_iter()
272 .map(|idx| {
273 let (i, j) = index_to_upper_triangle(idx, n);
275
276 if i == j {
277 let mut x_plus = x.to_owned();
279 let mut x_minus = x.to_owned();
280 x_plus[i] += eps;
281 x_minus[i] -= eps;
282
283 let f_plus = f(&x_plus.view());
284 let f_minus = f(&x_minus.view());
285 let f_center = f(&x);
286
287 (f_plus - 2.0 * f_center + f_minus) / (eps * eps)
288 } else {
289 let mut x_pp = x.to_owned();
291 let mut x_pm = x.to_owned();
292 let mut x_mp = x.to_owned();
293 let mut x_mm = x.to_owned();
294
295 x_pp[i] += eps;
296 x_pp[j] += eps;
297 x_pm[i] += eps;
298 x_pm[j] -= eps;
299 x_mp[i] -= eps;
300 x_mp[j] += eps;
301 x_mm[i] -= eps;
302 x_mm[j] -= eps;
303
304 let f_pp = f(&x_pp.view());
305 let f_pm = f(&x_pm.view());
306 let f_mp = f(&x_mp.view());
307 let f_mm = f(&x_mm.view());
308
309 (f_pp - f_pm - f_mp + f_mm) / (4.0 * eps * eps)
310 }
311 })
312 .collect();
313
314 #[cfg(not(feature = "parallel"))]
315 let hessian_elements: Vec<f64> = (0..num_elements)
316 .map(|idx| {
317 let (i, j) = index_to_upper_triangle(idx, n);
319
320 if i == j {
321 let mut x_plus = x.to_owned();
323 let mut x_minus = x.to_owned();
324 x_plus[i] += eps;
325 x_minus[i] -= eps;
326
327 let f_plus = f(&x_plus.view());
328 let f_minus = f(&x_minus.view());
329 let f_center = f(&x);
330
331 (f_plus - 2.0 * f_center + f_minus) / (eps * eps)
332 } else {
333 let mut x_pp = x.to_owned();
335 let mut x_pm = x.to_owned();
336 let mut x_mp = x.to_owned();
337 let mut x_mm = x.to_owned();
338
339 x_pp[i] += eps;
340 x_pp[j] += eps;
341 x_pm[i] += eps;
342 x_pm[j] -= eps;
343 x_mp[i] -= eps;
344 x_mp[j] += eps;
345 x_mm[i] -= eps;
346 x_mm[j] -= eps;
347
348 let f_pp = f(&x_pp.view());
349 let f_pm = f(&x_pm.view());
350 let f_mp = f(&x_mp.view());
351 let f_mm = f(&x_mm.view());
352
353 (f_pp - f_pm - f_mp + f_mm) / (4.0 * eps * eps)
354 }
355 })
356 .collect();
357
358 Array1::from_vec(hessian_elements)
359}
360
361#[allow(dead_code)]
363fn sequential_finite_diff_hessian<F>(
364 f: &F,
365 x: ArrayView1<f64>,
366 _gradient: Option<&Array1<f64>>,
367 eps: f64,
368) -> Array1<f64>
369where
370 F: Fn(&ArrayView1<f64>) -> f64,
371{
372 let n = x.len();
373 let num_elements = n * (n + 1) / 2;
374 let mut hessian = Array1::zeros(num_elements);
375
376 for idx in 0..num_elements {
377 let (i, j) = index_to_upper_triangle(idx, n);
378
379 if i == j {
380 let mut x_plus = x.to_owned();
382 let mut x_minus = x.to_owned();
383 x_plus[i] += eps;
384 x_minus[i] -= eps;
385
386 let f_plus = f(&x_plus.view());
387 let f_minus = f(&x_minus.view());
388 let f_center = f(&x);
389
390 hessian[idx] = (f_plus - 2.0 * f_center + f_minus) / (eps * eps);
391 } else {
392 let mut x_pp = x.to_owned();
394 let mut x_pm = x.to_owned();
395 let mut x_mp = x.to_owned();
396 let mut x_mm = x.to_owned();
397
398 x_pp[i] += eps;
399 x_pp[j] += eps;
400 x_pm[i] += eps;
401 x_pm[j] -= eps;
402 x_mp[i] -= eps;
403 x_mp[j] += eps;
404 x_mm[i] -= eps;
405 x_mm[j] -= eps;
406
407 let f_pp = f(&x_pp.view());
408 let f_pm = f(&x_pm.view());
409 let f_mp = f(&x_mp.view());
410 let f_mm = f(&x_mm.view());
411
412 hessian[idx] = (f_pp - f_pm - f_mp + f_mm) / (4.0 * eps * eps);
413 }
414 }
415
416 hessian
417}
418
419#[allow(dead_code)]
421fn index_to_upper_triangle(_idx: usize, n: usize) -> (usize, usize) {
422 let mut i = 0;
424 let mut cumsum = 0;
425
426 while cumsum + i < _idx {
427 i += 1;
428 cumsum += i;
429 }
430
431 let j = _idx - cumsum;
432 (j, i)
433}
434
435#[allow(dead_code)]
439pub fn parallel_line_search<F>(
440 f: F,
441 x: &Array1<f64>,
442 direction: &Array1<f64>,
443 step_sizes: &[f64],
444 options: &ParallelOptions,
445) -> (f64, f64)
446where
448 F: Fn(&ArrayView1<f64>) -> f64 + Sync,
449{
450 let evaluations: Vec<(f64, f64)> = if step_sizes.len() < options.min_parallel_size {
451 step_sizes
453 .iter()
454 .map(|&alpha| {
455 let x_new = x + alpha * direction;
456 let value = f(&x_new.view());
457 (alpha, value)
458 })
459 .collect()
460 } else {
461 #[cfg(feature = "parallel")]
463 {
464 step_sizes
465 .par_iter()
466 .map(|&alpha| {
467 let x_new = x + alpha * direction;
468 let value = f(&x_new.view());
469 (alpha, value)
470 })
471 .collect()
472 }
473 #[cfg(not(feature = "parallel"))]
474 {
475 step_sizes
476 .iter()
477 .map(|&alpha| {
478 let x_new = x + alpha * direction;
479 let value = f(&x_new.view());
480 (alpha, value)
481 })
482 .collect()
483 }
484 };
485
486 evaluations
488 .into_iter()
489 .min_by(|(_, v1), (_, v2)| v1.partial_cmp(v2).unwrap())
490 .unwrap_or((0.0, f64::INFINITY))
491}
492
493#[cfg(test)]
494mod tests {
495 use super::*;
496 use ndarray::array;
497
498 fn quadratic(x: &ArrayView1<f64>) -> f64 {
499 x.iter().map(|&xi| xi * xi).sum()
500 }
501
502 #[test]
503 fn test_parallel_gradient() {
504 let x = array![1.0, 2.0, 3.0, 4.0, 5.0];
505 let options = ParallelOptions::default();
506
507 let grad_parallel = parallel_finite_diff_gradient(quadratic, x.view(), &options);
508 let grad_sequential = sequential_finite_diff_gradient(quadratic, x.view(), 1e-8);
509
510 for i in 0..x.len() {
512 assert!((grad_parallel[i] - grad_sequential[i]).abs() < 1e-10);
513 assert!((grad_parallel[i] - 2.0 * x[i]).abs() < 1e-6);
515 }
516 }
517
518 #[test]
519 fn test_parallel_batch_evaluation() {
520 let points = vec![
521 array![1.0, 2.0],
522 array![3.0, 4.0],
523 array![5.0, 6.0],
524 array![7.0, 8.0],
525 ];
526
527 let options = ParallelOptions {
528 min_parallel_size: 2,
529 ..Default::default()
530 };
531
532 let values = parallel_evaluate_batch(quadratic, &points, &options);
533
534 for (i, point) in points.iter().enumerate() {
536 let expected = quadratic(&point.view());
537 assert_eq!(values[i], expected);
538 }
539 }
540
541 #[test]
542 fn test_parallel_line_search() {
543 let x = array![1.0, 1.0];
544 let direction = array![-1.0, -1.0];
545 let step_sizes: Vec<f64> = (0..20).map(|i| i as f64 * 0.1).collect();
546
547 let options = ParallelOptions::default();
548 let (best_step, best_value) =
549 parallel_line_search(quadratic, &x, &direction, &step_sizes, &options);
550
551 assert!((best_step - 1.0).abs() < 0.1);
554 assert!(best_value < 0.1);
555 }
556}