scirs2_optimize/parallel/
mod.rs

1//! Parallel computation utilities for optimization algorithms
2//!
3//! This module provides parallel implementations of computationally intensive
4//! operations used in optimization algorithms, including:
5//! - Parallel function evaluations
6//! - Parallel gradient approximations
7//! - Multi-start optimization with parallel workers
8//!
9//! # Example
10//!
11//! ```
12//! use scirs2_optimize::parallel::{parallel_finite_diff_gradient, ParallelOptions};
13//! use ndarray::{array, ArrayView1};
14//!
15//! fn objective(x: &ArrayView1<f64>) -> f64 {
16//!     x.iter().map(|&xi| xi.powi(2)).sum()
17//! }
18//!
19//! let x = array![1.0, 2.0, 3.0];
20//! let options = ParallelOptions::default();
21//! let gradient = parallel_finite_diff_gradient(objective, x.view(), &options);
22//! ```
23
24use ndarray::{Array1, ArrayView1};
25use scirs2_core::parallel_ops::*;
26
27// Conditional imports for parallel operations
28#[cfg(feature = "parallel")]
29use scirs2_core::parallel_ops::{IntoParallelIterator, ParallelIterator};
30
31/// Options for parallel computation
32#[derive(Debug, Clone)]
33pub struct ParallelOptions {
34    /// Number of worker threads (None = use core parallel default)
35    pub num_workers: Option<usize>,
36
37    /// Minimum problem size to enable parallelization
38    pub min_parallel_size: usize,
39
40    /// Chunk size for parallel operations
41    pub chunk_size: usize,
42
43    /// Enable parallel function evaluations
44    pub parallel_evaluations: bool,
45
46    /// Enable parallel gradient computation
47    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/// Parallel finite difference gradient approximation
63///
64/// Computes the gradient using parallel function evaluations when
65/// the dimension is large enough to benefit from parallelization.
66#[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    // For small problems, use sequential computation
79    if n < options.min_parallel_size || !options.parallel_gradient {
80        return sequential_finite_diff_gradient(f, x, eps);
81    }
82
83    // Use core parallel abstractions
84    compute_parallel_gradient(&f, x, eps)
85}
86
87/// Sequential finite difference gradient (fallback)
88#[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/// Compute gradient in parallel
108#[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    // Parallel computation of gradient components
118    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/// Sequential fallback for gradient computation
132#[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/// Parallel evaluation of multiple points
142///
143/// Evaluates the objective function at multiple points in parallel.
144#[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        // Sequential evaluation for small batches
155        points.iter().map(|x| f(&x.view())).collect()
156    } else {
157        // Parallel evaluation
158        #[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/// Parallel multi-start optimization
170///
171/// Runs multiple optimization instances from different starting points in parallel.
172#[allow(dead_code)] // Bounds will be used in future implementation
173pub 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    /// Create a new parallel multi-start optimizer
184    pub fn new(objective: F, bounds: Vec<(f64, f64)>, options: ParallelOptions) -> Self {
185        ParallelMultiStart {
186            objective,
187            bounds,
188            options,
189        }
190    }
191
192    /// Run optimization from multiple starting points
193    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            // Sequential execution for small number of starts
199            starting_points
200                .iter()
201                .map(|x0| optimizer(x0, &self.objective))
202                .collect()
203        } else {
204            // Parallel execution
205            #[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    /// Find the best result among all runs
223    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/// Result from an optimization run
231#[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/// Parallel computation of Hessian matrix using finite differences
240#[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>
247// Returns upper triangle in row-major order
248where
249    F: Fn(&ArrayView1<f64>) -> f64 + Sync,
250{
251    let n = x.len();
252    let eps = 1e-8;
253
254    // For small problems, use sequential computation
255    if n < options.min_parallel_size {
256        return sequential_finite_diff_hessian(f, x, gradient, eps);
257    }
258
259    // Compute gradient if not provided
260    let _grad = match gradient {
261        Some(g) => g.clone(),
262        None => parallel_finite_diff_gradient(f, x, options),
263    };
264
265    // Number of elements in upper triangle (including diagonal)
266    let num_elements = n * (n + 1) / 2;
267
268    // Parallel computation of Hessian elements
269    #[cfg(feature = "parallel")]
270    let hessian_elements: Vec<f64> = (0..num_elements)
271        .into_par_iter()
272        .map(|idx| {
273            // Convert linear index to (i, j) coordinates
274            let (i, j) = index_to_upper_triangle(idx, n);
275
276            if i == j {
277                // Diagonal element: ∂²f/∂xᵢ²
278                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                // Off-diagonal element: ∂²f/∂xᵢ∂xⱼ
290                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            // Convert linear index to (i, j) coordinates
318            let (i, j) = index_to_upper_triangle(idx, n);
319
320            if i == j {
321                // Diagonal element: ∂²f/∂xᵢ²
322                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                // Off-diagonal element: ∂²f/∂xᵢ∂xⱼ
334                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/// Sequential Hessian computation (fallback)
362#[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            // Diagonal element
381            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            // Off-diagonal element
393            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/// Convert linear index to (i, j) coordinates in upper triangle
420#[allow(dead_code)]
421fn index_to_upper_triangle(_idx: usize, n: usize) -> (usize, usize) {
422    // Find row i such that i*(i+1)/2 <= _idx < (i+1)*(i+2)/2
423    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/// Parallel line search
436///
437/// Evaluates multiple step sizes in parallel to find the best one.
438#[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)
446// Returns (best_step, best_value)
447where
448    F: Fn(&ArrayView1<f64>) -> f64 + Sync,
449{
450    let evaluations: Vec<(f64, f64)> = if step_sizes.len() < options.min_parallel_size {
451        // Sequential evaluation
452        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        // Parallel evaluation
462        #[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    // Find the best step size
487    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        // Check that parallel and sequential give similar results
511        for i in 0..x.len() {
512            assert!((grad_parallel[i] - grad_sequential[i]).abs() < 1e-10);
513            // For quadratic function, gradient should be 2*x
514            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        // Verify results
535        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        // For quadratic function starting at (1,1) going in direction (-1,-1),
552        // the minimum should be at step size 1.0
553        assert!((best_step - 1.0).abs() < 0.1);
554        assert!(best_value < 0.1);
555    }
556}