scirs2_core/memory_efficient/
zerocopy.rs

1//! Zero-copy operations for memory-mapped arrays.
2//!
3//! This module provides a set of operations that can be performed on memory-mapped
4//! arrays without loading the entire array into memory or making unnecessary copies.
5//! These operations maintain memory-mapping where possible and only load the minimum
6//! required data.
7
8use super::chunked::ChunkingStrategy;
9use super::memmap::{AccessMode, MemoryMappedArray};
10use crate::error::{CoreError, CoreResult, ErrorContext, ErrorLocation};
11// Use crate::ndarray for SciRS2 POLICY compliance (supports both ndarray 0.16 and 0.17)
12use crate::ndarray;
13use num_traits::Zero;
14use std::ops::{Add, Div, Mul, Sub};
15
16/// Trait for zero-copy operations on memory-mapped arrays.
17///
18/// This trait provides methods for performing operations on memory-mapped arrays
19/// without unnecessary memory allocations or copies. The operations are designed
20/// to work efficiently with large datasets by processing data in chunks and
21/// maintaining memory-mapping where possible.
22pub trait ZeroCopyOps<A: Clone + Copy + 'static + Send + Sync> {
23    /// Maps a function over each element of the array without loading the entire array.
24    ///
25    /// This is similar to the `map` function in functional programming, but implemented
26    /// to work efficiently with memory-mapped arrays by processing chunks.
27    ///
28    /// # Arguments
29    ///
30    /// * `f` - A function that takes an element of type `A` and returns a new element of the same type
31    ///
32    /// # Returns
33    ///
34    /// A new memory-mapped array containing the mapped values
35    ///
36    /// # Example
37    ///
38    /// ```no_run
39    /// # use scirs2_core::memory_efficient::{MemoryMappedArray, ZeroCopyOps};
40    /// # let mmap: MemoryMappedArray<f64> = unimplemented!();
41    /// // Double each element
42    /// let doubled = mmap.map_zero_copy(|x| x * 2.0);
43    /// ```
44    fn map_zero_copy<F>(&self, f: F) -> CoreResult<MemoryMappedArray<A>>
45    where
46        F: Fn(A) -> A + Send + Sync;
47
48    /// Reduces the array to a single value by applying a binary operation.
49    ///
50    /// This is similar to the `fold` or `reduce` function in functional programming,
51    /// but implemented to work efficiently with memory-mapped arrays by processing chunks.
52    ///
53    /// # Arguments
54    ///
55    /// * `init` - The initial value for the reduction
56    /// * `f` - A function that takes two values of type `A` and combines them into one
57    ///
58    /// # Returns
59    ///
60    /// The reduced value
61    ///
62    /// # Example
63    ///
64    /// ```no_run
65    /// # use scirs2_core::memory_efficient::{MemoryMappedArray, ZeroCopyOps};
66    /// # let mmap: MemoryMappedArray<f64> = unimplemented!();
67    /// // Sum all elements
68    /// let sum = mmap.reduce_zero_copy(0.0, |acc, x| acc + x);
69    /// ```
70    fn reduce_zero_copy<F>(&self, init: A, f: F) -> CoreResult<A>
71    where
72        F: Fn(A, A) -> A + Send + Sync;
73
74    /// Performs a binary operation between two memory-mapped arrays element-wise.
75    ///
76    /// This allows for operations like addition, subtraction, etc. between two arrays
77    /// without loading both arrays entirely into memory.
78    ///
79    /// # Arguments
80    ///
81    /// * `other` - Another memory-mapped array with the same shape
82    /// * `f` - A function that takes two elements (one from each array) and returns a new element
83    ///
84    /// # Returns
85    ///
86    /// A new memory-mapped array containing the result of the binary operation
87    ///
88    /// # Example
89    ///
90    /// ```no_run
91    /// # use scirs2_core::memory_efficient::{MemoryMappedArray, ZeroCopyOps};
92    /// # let mmap1: MemoryMappedArray<f64> = unimplemented!();
93    /// # let mmap2: MemoryMappedArray<f64> = unimplemented!();
94    /// // Add two arrays element-wise
95    /// let sum_array = mmap1.combine_zero_copy(&mmap2, |a, b| a + b);
96    /// ```
97    fn combine_zero_copy<F>(&self, other: &Self, f: F) -> CoreResult<MemoryMappedArray<A>>
98    where
99        F: Fn(A, A) -> A + Send + Sync;
100
101    /// Filters elements based on a predicate function.
102    ///
103    /// Returns a new array containing only the elements that satisfy the predicate.
104    ///
105    /// # Arguments
106    ///
107    /// * `predicate` - A function that takes an element and returns a boolean
108    ///
109    /// # Returns
110    ///
111    /// A new array containing only the elements that satisfy the predicate
112    ///
113    /// # Example
114    ///
115    /// ```no_run
116    /// # use scirs2_core::memory_efficient::{MemoryMappedArray, ZeroCopyOps};
117    /// # let mmap: MemoryMappedArray<f64> = unimplemented!();
118    /// // Get only positive elements
119    /// let positives = mmap.filter_zero_copy(|&x| x > 0.0);
120    /// ```
121    fn filter_zero_copy<F>(&self, predicate: F) -> CoreResult<Vec<A>>
122    where
123        F: Fn(&A) -> bool + Send + Sync;
124
125    /// Returns the maximum element in the array.
126    ///
127    /// # Returns
128    ///
129    /// The maximum element, or an error if the array is empty
130    ///
131    /// # Example
132    ///
133    /// ```no_run
134    /// # use scirs2_core::memory_efficient::{MemoryMappedArray, ZeroCopyOps};
135    /// # let mmap: MemoryMappedArray<f64> = unimplemented!();
136    /// let max_value = mmap.max_zero_copy();
137    /// ```
138    fn max_zero_copy(&self) -> CoreResult<A>
139    where
140        A: PartialOrd;
141
142    /// Returns the minimum element in the array.
143    ///
144    /// # Returns
145    ///
146    /// The minimum element, or an error if the array is empty
147    ///
148    /// # Example
149    ///
150    /// ```no_run
151    /// # use scirs2_core::memory_efficient::{MemoryMappedArray, ZeroCopyOps};
152    /// # let mmap: MemoryMappedArray<f64> = unimplemented!();
153    /// let min_value = mmap.min_zero_copy();
154    /// ```
155    fn min_zero_copy(&self) -> CoreResult<A>
156    where
157        A: PartialOrd;
158
159    /// Calculates the sum of all elements in the array.
160    ///
161    /// # Returns
162    ///
163    /// The sum of all elements
164    ///
165    /// # Example
166    ///
167    /// ```no_run
168    /// # use scirs2_core::memory_efficient::{MemoryMappedArray, ZeroCopyOps};
169    /// # let mmap: MemoryMappedArray<f64> = unimplemented!();
170    /// let total = mmap.sum_zero_copy();
171    /// ```
172    fn sum_zero_copy(&self) -> CoreResult<A>
173    where
174        A: Add<Output = A> + From<u8>;
175
176    /// Calculates the product of all elements in the array.
177    ///
178    /// # Returns
179    ///
180    /// The product of all elements
181    ///
182    /// # Example
183    ///
184    /// ```no_run
185    /// # use scirs2_core::memory_efficient::{MemoryMappedArray, ZeroCopyOps};
186    /// # let mmap: MemoryMappedArray<f64> = unimplemented!();
187    /// let product = mmap.product_zero_copy();
188    /// ```
189    fn product_zero_copy(&self) -> CoreResult<A>
190    where
191        A: Mul<Output = A> + From<u8>;
192
193    /// Calculates the mean of all elements in the array.
194    ///
195    /// # Returns
196    ///
197    /// The mean of all elements
198    ///
199    /// # Example
200    ///
201    /// ```ignore
202    /// # use scirs2_core::memory_efficient::{MemoryMappedArray, ZeroCopyOps};
203    /// # let mmap: MemoryMappedArray<f64> = unimplemented!();
204    /// let avg = mmap.mean_zero_copy();
205    /// ```
206    fn mean_zero_copy(&self) -> CoreResult<A>
207    where
208        A: Add<Output = A> + Div<Output = A> + From<u8> + From<usize>;
209}
210
211impl<A: Clone + Copy + 'static + Send + Sync + Send + Sync + Zero> ZeroCopyOps<A>
212    for MemoryMappedArray<A>
213{
214    fn map_zero_copy<F>(&self, f: F) -> CoreResult<MemoryMappedArray<A>>
215    where
216        F: Fn(A) -> A + Send + Sync,
217    {
218        // Create a temporary file for the result
219        let temp_file = tempfile::NamedTempFile::new()?;
220        let temp_path = temp_file.path().to_path_buf();
221
222        // Create an output memory-mapped array with the same shape
223        let shape = &self.shape;
224        let element_size = std::mem::size_of::<A>();
225        let file_size = self.size * element_size;
226
227        // Resize the temp file to the required size
228        temp_file.as_file().set_len(file_size as u64)?;
229        drop(temp_file); // Close the file before memory-mapping it
230
231        // Create the output memory-mapped array with zeros to initialize the header
232        // First create with Write mode to initialize
233        let zeros = crate::ndarray::ArrayD::zeros(crate::ndarray::IxDyn(&self.shape));
234        {
235            let _ = MemoryMappedArray::<A>::new::<
236                crate::ndarray::OwnedRepr<A>,
237                crate::ndarray::IxDyn,
238            >(Some(&zeros), &temp_path, AccessMode::Write, 0)?;
239        }
240
241        // Now reopen in ReadWrite mode to allow modifications
242        let mut output = MemoryMappedArray::<A>::new::<
243            crate::ndarray::OwnedRepr<A>,
244            crate::ndarray::IxDyn,
245        >(None, &temp_path, AccessMode::ReadWrite, 0)?;
246
247        // Process the input array in chunks
248        #[cfg(feature = "parallel")]
249        {
250            // Use rayon directly since process_chunks_parallel has trait bounds
251            // that we can't satisfy here (Send + Sync)
252
253            let chunk_size = (self.size / rayon::current_num_threads()).max(1024);
254
255            // Calculate the number of chunks
256            let num_chunks = self.size.div_ceil(chunk_size);
257
258            // Process each chunk sequentially to avoid mutable borrow issues
259            let array = self.as_array::<crate::ndarray::IxDyn>()?;
260            let mut out_array = output.as_array_mut::<crate::ndarray::IxDyn>()?;
261
262            for chunk_idx in 0..num_chunks {
263                // Calculate chunk bounds
264                let start = chunk_idx * chunk_size;
265                let end = (start + chunk_size).min(self.size);
266
267                // Get the data for this chunk
268                let slice = array.as_slice().ok_or_else(|| {
269                    CoreError::ValidationError(
270                        ErrorContext::new("Array is not contiguous in memory".to_string())
271                            .with_location(ErrorLocation::new(file!(), line!())),
272                    )
273                })?;
274                let chunk = &slice[start..end];
275
276                // Apply the mapping function to each element in the chunk
277                let mapped_chunk: Vec<A> = chunk.iter().map(|&x| f(x)).collect();
278
279                // Copy the mapped chunk to the output at the same position
280                let out_slice_full = out_array.as_slice_mut().ok_or_else(|| {
281                    CoreError::ValidationError(
282                        ErrorContext::new("Output array is not contiguous in memory".to_string())
283                            .with_location(ErrorLocation::new(file!(), line!())),
284                    )
285                })?;
286                let out_slice = &mut out_slice_full[start..end];
287
288                // Copy the mapped chunk to the output
289                out_slice.copy_from_slice(&mapped_chunk);
290            }
291        }
292
293        #[cfg(not(feature = "parallel"))]
294        {
295            // Use sequential processing
296
297            let chunk_size = 1024 * 1024; // 1M elements
298            let strategy = ChunkingStrategy::Fixed(chunk_size);
299
300            // Manually process chunks instead of using process_chunks_mut
301            for chunk_idx in 0..self.size.div_ceil(chunk_size) {
302                // Calculate chunk bounds
303                let start = chunk_idx * chunk_size;
304                let end = (start + chunk_size).min(self.size);
305
306                // Get the data for this chunk
307                let array = self.as_array::<crate::ndarray::IxDyn>()?;
308                let slice = array.as_slice().ok_or_else(|| {
309                    CoreError::ValidationError(
310                        ErrorContext::new("Array is not contiguous in memory".to_string())
311                            .with_location(ErrorLocation::new(file!(), line!())),
312                    )
313                })?;
314                let chunk = &slice[start..end];
315
316                // Apply the mapping function to each element in the chunk
317                let mapped_chunk: Vec<A> = chunk.iter().map(|&x| f(x)).collect();
318
319                // Copy the mapped chunk to the output at the same position
320                // Get a mutable view of the output array
321                let mut out_array = output.as_array_mut::<crate::ndarray::IxDyn>()?;
322                let out_slice_full = out_array.as_slice_mut().ok_or_else(|| {
323                    CoreError::ValidationError(
324                        ErrorContext::new("Output array is not contiguous in memory".to_string())
325                            .with_location(ErrorLocation::new(file!(), line!())),
326                    )
327                })?;
328                let out_slice = &mut out_slice_full[start..end];
329
330                // Copy the mapped chunk to the output
331                out_slice.copy_from_slice(&mapped_chunk);
332            }
333        }
334
335        Ok(output)
336    }
337
338    fn reduce_zero_copy<F>(&self, init: A, f: F) -> CoreResult<A>
339    where
340        F: Fn(A, A) -> A + Send + Sync,
341    {
342        // Process the input array in chunks
343        let chunk_size = 1024 * 1024; // 1M elements
344        let strategy = ChunkingStrategy::Fixed(chunk_size);
345
346        // Since we can't use process_chunks directly, we'll implement manually
347        let num_chunks = self.size.div_ceil(chunk_size);
348        let mut chunk_results = Vec::with_capacity(num_chunks);
349
350        // Process each chunk
351        for chunk_idx in 0..num_chunks {
352            // Calculate chunk bounds
353            let start = chunk_idx * chunk_size;
354            let end = (start + chunk_size).min(self.size);
355
356            // Load the array
357            let array = self.as_array::<crate::ndarray::IxDyn>()?;
358            let slice = array.as_slice().ok_or_else(|| {
359                CoreError::ValidationError(
360                    ErrorContext::new("Array is not contiguous in memory".to_string())
361                        .with_location(ErrorLocation::new(file!(), line!())),
362                )
363            })?;
364            let chunk = &slice[start..end];
365
366            // Reduce the chunk
367            let chunk_result = chunk.iter().fold(init, |acc, &x| f(acc, x));
368            chunk_results.push(chunk_result);
369        }
370
371        // Combine chunk results
372        let final_result = chunk_results.into_iter().fold(init, f);
373
374        Ok(final_result)
375    }
376
377    fn combine_zero_copy<F>(&self, other: &Self, f: F) -> CoreResult<MemoryMappedArray<A>>
378    where
379        F: Fn(A, A) -> A + Send + Sync,
380    {
381        // Check that the arrays have the same shape
382        if self.shape != other.shape {
383            return Err(CoreError::ShapeError(ErrorContext::new(format!(
384                "Arrays have different shapes: {:?} vs {:?}",
385                self.shape, other.shape
386            ))));
387        }
388
389        // Create a temporary file for the result
390        let temp_file = tempfile::NamedTempFile::new()?;
391        let temp_path = temp_file.path().to_path_buf();
392
393        // Create an output memory-mapped array with the same shape
394        let shape = &self.shape;
395        let element_size = std::mem::size_of::<A>();
396        let file_size = self.size * element_size;
397
398        // Resize the temp file to the required size
399        temp_file.as_file().set_len(file_size as u64)?;
400        drop(temp_file); // Close the file before memory-mapping it
401
402        // Create the output memory-mapped array with zeros to initialize the header
403        // First create with Write mode to initialize
404        let zeros = crate::ndarray::ArrayD::zeros(crate::ndarray::IxDyn(&self.shape));
405        {
406            let _ = MemoryMappedArray::<A>::new::<
407                crate::ndarray::OwnedRepr<A>,
408                crate::ndarray::IxDyn,
409            >(Some(&zeros), &temp_path, AccessMode::Write, 0)?;
410        }
411
412        // Now reopen in ReadWrite mode to allow modifications
413        let mut output = MemoryMappedArray::<A>::new::<
414            crate::ndarray::OwnedRepr<A>,
415            crate::ndarray::IxDyn,
416        >(None, &temp_path, AccessMode::ReadWrite, 0)?;
417
418        // Process the arrays in chunks
419        let chunk_size = 1024 * 1024; // 1M elements
420        let strategy = ChunkingStrategy::Fixed(chunk_size);
421
422        // Calculate the number of chunks
423        let num_chunks = self.size.div_ceil(chunk_size);
424
425        // Process each chunk
426        for chunk_idx in 0..num_chunks {
427            // Calculate chunk bounds
428            let start = chunk_idx * chunk_size;
429            let end = (start + chunk_size).min(self.size);
430            let len = end - start;
431
432            // Load chunks from both arrays
433            let self_array = self.as_array::<crate::ndarray::IxDyn>()?;
434            let other_array = other.as_array::<crate::ndarray::IxDyn>()?;
435
436            let self_slice = self_array.as_slice().ok_or_else(|| {
437                CoreError::ValidationError(
438                    ErrorContext::new("Self array is not contiguous in memory".to_string())
439                        .with_location(ErrorLocation::new(file!(), line!())),
440                )
441            })?;
442            let other_slice = other_array.as_slice().ok_or_else(|| {
443                CoreError::ValidationError(
444                    ErrorContext::new("Other array is not contiguous in memory".to_string())
445                        .with_location(ErrorLocation::new(file!(), line!())),
446                )
447            })?;
448            let self_chunk = &self_slice[start..end];
449            let other_chunk = &other_slice[start..end];
450
451            // Apply the binary operation
452            let mut result_chunk = Vec::with_capacity(len);
453            for i in 0..len {
454                result_chunk.push(f(self_chunk[i], other_chunk[i]));
455            }
456
457            // Write the result to the output array
458            let mut out_array = output.as_array_mut::<crate::ndarray::IxDyn>()?;
459            let out_slice_full = out_array.as_slice_mut().ok_or_else(|| {
460                CoreError::ValidationError(
461                    ErrorContext::new("Output array is not contiguous in memory".to_string())
462                        .with_location(ErrorLocation::new(file!(), line!())),
463                )
464            })?;
465            let out_slice = &mut out_slice_full[start..end];
466            out_slice.copy_from_slice(&result_chunk);
467        }
468
469        Ok(output)
470    }
471
472    fn filter_zero_copy<F>(&self, predicate: F) -> CoreResult<Vec<A>>
473    where
474        F: Fn(&A) -> bool + Send + Sync,
475    {
476        // Process the input array in chunks manually
477        let chunk_size = 1024 * 1024; // 1M elements
478        let num_chunks = self.size.div_ceil(chunk_size);
479        let mut result = Vec::new();
480
481        // Process each chunk
482        for chunk_idx in 0..num_chunks {
483            // Calculate chunk bounds
484            let start = chunk_idx * chunk_size;
485            let end = (start + chunk_size).min(self.size);
486
487            // Load the array
488            let array = self.as_array::<crate::ndarray::IxDyn>()?;
489            let array_slice = array.as_slice().ok_or_else(|| {
490                CoreError::ValidationError(
491                    ErrorContext::new("Array is not contiguous in memory".to_string())
492                        .with_location(ErrorLocation::new(file!(), line!())),
493                )
494            })?;
495            let slice = &array_slice[start..end];
496
497            // Filter the chunk
498            let filtered_chunk = slice
499                .iter()
500                .filter(|&x| predicate(x))
501                .cloned()
502                .collect::<Vec<A>>();
503
504            // Add filtered elements to the result
505            result.extend(filtered_chunk);
506        }
507
508        Ok(result)
509    }
510
511    fn max_zero_copy(&self) -> CoreResult<A>
512    where
513        A: PartialOrd,
514    {
515        // Handle empty array
516        if self.size == 0 {
517            return Err(CoreError::ValueError(ErrorContext::new(
518                "Array is empty".to_string(),
519            )));
520        }
521
522        // Read the first element to initialize
523        let first_element = {
524            let array = self.as_array::<crate::ndarray::IxDyn>()?;
525            let slice = array.as_slice().ok_or_else(|| {
526                CoreError::ValidationError(
527                    ErrorContext::new("Array is not contiguous in memory".to_string())
528                        .with_location(ErrorLocation::new(file!(), line!())),
529                )
530            })?;
531            slice[0]
532        };
533
534        // Reduce the array to find the maximum
535        self.reduce_zero_copy(first_element, |acc, x| if x > acc { x } else { acc })
536    }
537
538    fn min_zero_copy(&self) -> CoreResult<A>
539    where
540        A: PartialOrd,
541    {
542        // Handle empty array
543        if self.size == 0 {
544            return Err(CoreError::ValueError(ErrorContext::new(
545                "Array is empty".to_string(),
546            )));
547        }
548
549        // Read the first element to initialize
550        let first_element = {
551            let array = self.as_array::<crate::ndarray::IxDyn>()?;
552            let slice = array.as_slice().ok_or_else(|| {
553                CoreError::ValidationError(
554                    ErrorContext::new("Array is not contiguous in memory".to_string())
555                        .with_location(ErrorLocation::new(file!(), line!())),
556                )
557            })?;
558            slice[0]
559        };
560
561        // Reduce the array to find the minimum
562        self.reduce_zero_copy(first_element, |acc, x| if x < acc { x } else { acc })
563    }
564
565    fn sum_zero_copy(&self) -> CoreResult<A>
566    where
567        A: Add<Output = A> + From<u8>,
568    {
569        // Initialize with zero
570        let zero = A::from(0u8);
571
572        // Sum all elements
573        self.reduce_zero_copy(zero, |acc, x| acc + x)
574    }
575
576    fn product_zero_copy(&self) -> CoreResult<A>
577    where
578        A: Mul<Output = A> + From<u8>,
579    {
580        // Handle empty array
581        if self.size == 0 {
582            return Err(CoreError::ValueError(ErrorContext::new(
583                "Array is empty".to_string(),
584            )));
585        }
586
587        // Initialize with one
588        let one = A::from(1u8);
589
590        // Multiply all elements
591        self.reduce_zero_copy(one, |acc, x| acc * x)
592    }
593
594    fn mean_zero_copy(&self) -> CoreResult<A>
595    where
596        A: Add<Output = A> + Div<Output = A> + From<u8> + From<usize>,
597    {
598        // Handle empty array
599        if self.size == 0 {
600            return Err(CoreError::ValueError(ErrorContext::new(
601                "Array is empty".to_string(),
602            )));
603        }
604
605        // Calculate sum
606        let sum = self.sum_zero_copy()?;
607
608        // Divide by count
609        let count = A::from(self.size);
610
611        Ok(sum / count)
612    }
613}
614
615/// Trait for broadcasting operations between memory-mapped arrays of different shapes.
616///
617/// This trait provides methods for performing broadcasting operations between
618/// memory-mapped arrays without unnecessary memory allocations or copies.
619pub trait BroadcastOps<A: Clone + Copy + 'static + Send + Sync> {
620    /// Broadcasts an operation between two arrays of compatible shapes.
621    ///
622    /// Follows the `NumPy` broadcasting rules:
623    /// 1. If arrays don't have the same rank, prepend shape with 1s
624    /// 2. Two dimensions are compatible if:
625    ///    - They are equal, or
626    ///    - One of them is 1
627    ///
628    /// # Arguments
629    ///
630    /// * `other` - Another memory-mapped array with a compatible shape
631    /// * `f` - A function that takes two elements (one from each array) and returns a new element
632    ///
633    /// # Returns
634    ///
635    /// A new memory-mapped array containing the result of the broadcasted operation
636    ///
637    /// # Example
638    ///
639    /// ```no_run
640    /// # use scirs2_core::memory_efficient::{MemoryMappedArray, BroadcastOps};
641    /// # let mmap1: MemoryMappedArray<f64> = unimplemented!(); // Shape [3, 4]
642    /// # let mmap2: MemoryMappedArray<f64> = unimplemented!(); // Shape [4]
643    /// // Broadcast and multiply
644    /// let result = mmap1.broadcast_op(&mmap2, |a, b| a * b);
645    /// ```
646    fn broadcast_op<F>(&self, other: &Self, f: F) -> CoreResult<MemoryMappedArray<A>>
647    where
648        F: Fn(A, A) -> A + Send + Sync;
649}
650
651impl<A: Clone + Copy + 'static + Send + Sync + Send + Sync + Zero> BroadcastOps<A>
652    for MemoryMappedArray<A>
653{
654    fn broadcast_op<F>(&self, other: &Self, f: F) -> CoreResult<MemoryMappedArray<A>>
655    where
656        F: Fn(A, A) -> A + Send + Sync,
657    {
658        // Check shape compatibility for broadcasting
659        let selfshape = &self.shape;
660        let othershape = &other.shape;
661
662        // Get the dimensions
663        let self_ndim = selfshape.len();
664        let other_ndim = othershape.len();
665        let output_ndim = std::cmp::max(self_ndim, other_ndim);
666
667        // Convert shapes to vectors with leading 1s as needed
668        let mut self_dims = Vec::with_capacity(output_ndim);
669        let mut other_dims = Vec::with_capacity(output_ndim);
670
671        // Prepend 1s to the shape with fewer dimensions
672        self_dims.resize(output_ndim - self_ndim, 1);
673        for dim in selfshape.iter() {
674            self_dims.push(*dim);
675        }
676
677        other_dims.resize(output_ndim - other_ndim, 1);
678        for dim in othershape.iter() {
679            other_dims.push(*dim);
680        }
681
682        // Determine the output shape
683        let mut outputshape = Vec::with_capacity(output_ndim);
684        for i in 0..output_ndim {
685            #[allow(clippy::if_same_then_else)]
686            if self_dims[i] == 1 {
687                outputshape.push(other_dims[i]);
688            } else if other_dims[i] == 1 {
689                outputshape.push(self_dims[i]);
690            } else if self_dims[i] == other_dims[i] {
691                outputshape.push(self_dims[i]);
692            } else {
693                return Err(CoreError::ValueError(ErrorContext::new(format!(
694                    "Arrays cannot be broadcast together with shapes {selfshape:?} and {othershape:?}"
695                ))));
696            }
697        }
698
699        // Create a temporary file for the result
700        let temp_file = tempfile::NamedTempFile::new()?;
701        let temp_path = temp_file.path().to_path_buf();
702
703        // Calculate the output array size
704        let output_size = outputshape.iter().product::<usize>();
705        let element_size = std::mem::size_of::<A>();
706        let file_size = output_size * element_size;
707
708        // Resize the temp file to the required size
709        temp_file.as_file().set_len(file_size as u64)?;
710        drop(temp_file); // Close the file before memory-mapping it
711
712        // Create the output memory-mapped array with zeros to initialize the header
713        // Use the calculated output shape instead of self.shape
714        let zeros = crate::ndarray::ArrayD::zeros(crate::ndarray::IxDyn(&outputshape));
715        {
716            let _ = MemoryMappedArray::<A>::new::<
717                crate::ndarray::OwnedRepr<A>,
718                crate::ndarray::IxDyn,
719            >(Some(&zeros), &temp_path, AccessMode::Write, 0)?;
720        }
721
722        // Now reopen in ReadWrite mode to allow modifications
723        let mut output = MemoryMappedArray::<A>::new::<
724            crate::ndarray::OwnedRepr<A>,
725            crate::ndarray::IxDyn,
726        >(None, &temp_path, AccessMode::ReadWrite, 0)?;
727
728        // Load both arrays into memory (for broadcasting, we need random access)
729        let self_array = self.as_array::<crate::ndarray::IxDyn>()?;
730        let other_array = other.as_array::<crate::ndarray::IxDyn>()?;
731
732        // Create ndarray views for easier broadcasting
733        let self_view = self_array.view();
734        let other_view = other_array.view();
735
736        // Perform the broadcasted operation
737        let mut output_array = output.as_array_mut::<crate::ndarray::IxDyn>()?;
738
739        // Use ndarray's broadcasting capability
740        crate::ndarray::Zip::from(&mut output_array)
741            .and_broadcast(&self_view)
742            .and_broadcast(&other_view)
743            .for_each(|out, &a, &b| {
744                *out = f(a, b);
745            });
746
747        Ok(output)
748    }
749}
750
751/// Extension trait for standard arithmetic operations on memory-mapped arrays.
752///
753/// This trait provides implementations of standard arithmetic operations
754/// (addition, subtraction, multiplication, division) for memory-mapped arrays
755/// using the zero-copy infrastructure.
756pub trait ArithmeticOps<A: Clone + Copy + 'static + Send + Sync> {
757    /// Adds two arrays element-wise.
758    ///
759    /// # Arguments
760    ///
761    /// * `other` - Another memory-mapped array with the same shape
762    ///
763    /// # Returns
764    ///
765    /// A new memory-mapped array containing the sum
766    fn add(&self, other: &Self) -> CoreResult<MemoryMappedArray<A>>
767    where
768        A: Add<Output = A>;
769
770    /// Subtracts another array from this one element-wise.
771    ///
772    /// # Arguments
773    ///
774    /// * `other` - Another memory-mapped array with the same shape
775    ///
776    /// # Returns
777    ///
778    /// A new memory-mapped array containing the difference
779    fn sub(&self, other: &Self) -> CoreResult<MemoryMappedArray<A>>
780    where
781        A: Sub<Output = A>;
782
783    /// Multiplies two arrays element-wise.
784    ///
785    /// # Arguments
786    ///
787    /// * `other` - Another memory-mapped array with the same shape
788    ///
789    /// # Returns
790    ///
791    /// A new memory-mapped array containing the product
792    fn mul(&self, other: &Self) -> CoreResult<MemoryMappedArray<A>>
793    where
794        A: Mul<Output = A>;
795
796    /// Divides this array by another element-wise.
797    ///
798    /// # Arguments
799    ///
800    /// * `other` - Another memory-mapped array with the same shape
801    ///
802    /// # Returns
803    ///
804    /// A new memory-mapped array containing the quotient
805    fn div(&self, other: &Self) -> CoreResult<MemoryMappedArray<A>>
806    where
807        A: Div<Output = A>;
808}
809
810impl<A: Clone + Copy + 'static + Send + Sync + Send + Sync + Zero> ArithmeticOps<A>
811    for MemoryMappedArray<A>
812{
813    fn add(&self, other: &Self) -> CoreResult<MemoryMappedArray<A>>
814    where
815        A: Add<Output = A>,
816    {
817        self.combine_zero_copy(other, |a, b| a + b)
818    }
819
820    fn sub(&self, other: &Self) -> CoreResult<MemoryMappedArray<A>>
821    where
822        A: Sub<Output = A>,
823    {
824        self.combine_zero_copy(other, |a, b| a - b)
825    }
826
827    fn mul(&self, other: &Self) -> CoreResult<MemoryMappedArray<A>>
828    where
829        A: Mul<Output = A>,
830    {
831        self.combine_zero_copy(other, |a, b| a * b)
832    }
833
834    fn div(&self, other: &Self) -> CoreResult<MemoryMappedArray<A>>
835    where
836        A: Div<Output = A>,
837    {
838        self.combine_zero_copy(other, |a, b| a / b)
839    }
840}
841
842#[cfg(test)]
843mod tests {
844    use super::*;
845    use ::ndarray::Array2;
846    use std::fs::File;
847    use std::io::Write;
848    use tempfile::tempdir;
849
850    #[test]
851    fn test_map_zero_copy() {
852        // Create a temporary directory for our test files
853        let dir = tempdir().expect("Operation failed");
854        let file_path = dir.path().join("test_map.bin");
855
856        // Create a test array and save it with proper header using save_array
857        let data = crate::ndarray::Array1::from_vec((0..1000).map(|i| i as f64).collect());
858        MemoryMappedArray::<f64>::save_array(&data, &file_path, None).expect("Operation failed");
859
860        // Open the file for zero-copy operations
861        let mmap = MemoryMappedArray::<f64>::open_zero_copy(&file_path, AccessMode::ReadOnly)
862            .expect("Operation failed");
863
864        // Map operation: double each element
865        let result = mmap.map_zero_copy(|x| x * 2.0).expect("Operation failed");
866
867        // Verify the result
868        let result_array = result
869            .readonlyarray::<crate::ndarray::Ix1>()
870            .expect("Operation failed");
871        for i in 0..1000 {
872            assert_eq!(result_array[i], (i as f64) * 2.0);
873        }
874    }
875
876    #[test]
877    fn test_reduce_zero_copy() {
878        // Create a temporary directory for our test files
879        let dir = tempdir().expect("Operation failed");
880        let file_path = dir.path().join("test_reduce.bin");
881
882        // Create a test array and save it to a file
883        let data: Vec<f64> = (0..1000).map(|i| i as f64).collect();
884        let mut file = File::create(&file_path).expect("Operation failed");
885        for val in &data {
886            file.write_all(&val.to_ne_bytes())
887                .expect("Operation failed");
888        }
889        drop(file);
890
891        // Create a memory-mapped array
892        let mmap = MemoryMappedArray::<f64>::path(&file_path, &[1000]).expect("Operation failed");
893
894        // Reduce operation: sum all elements
895        let sum = mmap
896            .reduce_zero_copy(0.0, |acc, x| acc + x)
897            .expect("Operation failed");
898
899        // Verify the result (sum of 0..999 = 499500)
900        assert_eq!(sum, 499500.0);
901    }
902
903    #[test]
904    fn test_combine_zero_copy() {
905        // Create a temporary directory for our test files
906        let dir = tempdir().expect("Operation failed");
907        let file_path1 = dir.path().join("test_combine1.bin");
908        let file_path2 = dir.path().join("test_combine2.bin");
909
910        // Create two test arrays and save them with proper headers using save_array
911        let data1 = crate::ndarray::Array1::from_vec((0..1000).map(|i| i as f64).collect());
912        let data2 = crate::ndarray::Array1::from_vec((0..1000).map(|i| (i * 2) as f64).collect());
913
914        MemoryMappedArray::<f64>::save_array(&data1, &file_path1, None).expect("Operation failed");
915        MemoryMappedArray::<f64>::save_array(&data2, &file_path2, None).expect("Operation failed");
916
917        // Open the files for zero-copy operations
918        let mmap1 = MemoryMappedArray::<f64>::open_zero_copy(&file_path1, AccessMode::ReadOnly)
919            .expect("Operation failed");
920        let mmap2 = MemoryMappedArray::<f64>::open_zero_copy(&file_path2, AccessMode::ReadOnly)
921            .expect("Operation failed");
922
923        // Combine operation: add the arrays
924        let result = mmap1
925            .combine_zero_copy(&mmap2, |a, b| a + b)
926            .expect("Operation failed");
927
928        // Verify the result (each element should be 3*i)
929        let result_array = result
930            .readonlyarray::<crate::ndarray::Ix1>()
931            .expect("Operation failed");
932        for i in 0..1000 {
933            assert_eq!(result_array[i], (i as f64) * 3.0);
934        }
935    }
936
937    #[test]
938    fn test_filter_zero_copy() {
939        // Create a temporary directory for our test files
940        let dir = tempdir().expect("Operation failed");
941        let file_path = dir.path().join("test_filter.bin");
942
943        // Create a test array and save it to a file
944        let data: Vec<f64> = (0..1000).map(|i| i as f64).collect();
945        let mut file = File::create(&file_path).expect("Operation failed");
946        for val in &data {
947            file.write_all(&val.to_ne_bytes())
948                .expect("Operation failed");
949        }
950        drop(file);
951
952        // Create a memory-mapped array
953        let mmap = MemoryMappedArray::<f64>::path(&file_path, &[1000]).expect("Operation failed");
954
955        // Filter operation: keep only even numbers
956        let even_numbers = mmap
957            .filter_zero_copy(|&x| (x as usize) % 2 == 0)
958            .expect("Operation failed");
959
960        // Verify the result (should be 0, 2, 4, ..., 998)
961        assert_eq!(even_numbers.len(), 500);
962        for (i, val) in even_numbers.iter().enumerate() {
963            assert_eq!(*val, (i * 2) as f64);
964        }
965    }
966
967    #[test]
968    fn test_arithmetic_ops() {
969        // Create a temporary directory for our test files
970        let dir = tempdir().expect("Operation failed");
971        let file_path1 = dir.path().join("test_arithmetic1.bin");
972        let file_path2 = dir.path().join("test_arithmetic2.bin");
973
974        // Create two test arrays and save them with proper headers using save_array
975        let data1 = crate::ndarray::Array1::from_vec((0..100).map(|i| i as f64).collect());
976        let data2 = crate::ndarray::Array1::from_vec((0..100).map(|i| (i + 5) as f64).collect());
977
978        MemoryMappedArray::<f64>::save_array(&data1, &file_path1, None).expect("Operation failed");
979        MemoryMappedArray::<f64>::save_array(&data2, &file_path2, None).expect("Operation failed");
980
981        // Open the files for zero-copy operations
982        let mmap1 = MemoryMappedArray::<f64>::open_zero_copy(&file_path1, AccessMode::ReadOnly)
983            .expect("Operation failed");
984        let mmap2 = MemoryMappedArray::<f64>::open_zero_copy(&file_path2, AccessMode::ReadOnly)
985            .expect("Operation failed");
986
987        // Test addition
988        let add_result = mmap1.add(&mmap2).expect("Operation failed");
989        let add_array = add_result
990            .readonlyarray::<crate::ndarray::Ix1>()
991            .expect("Operation failed");
992        for i in 0..100 {
993            assert_eq!(add_array[i], (i as f64) + ((i + 5) as f64));
994        }
995
996        // Test subtraction
997        let sub_result = mmap1.sub(&mmap2).expect("Operation failed");
998        let sub_array = sub_result
999            .readonlyarray::<crate::ndarray::Ix1>()
1000            .expect("Operation failed");
1001        for i in 0..100 {
1002            assert_eq!(sub_array[i], (i as f64) - ((i + 5) as f64));
1003        }
1004
1005        // Test multiplication
1006        let mul_result = mmap1.mul(&mmap2).expect("Operation failed");
1007        let mul_array = mul_result
1008            .readonlyarray::<crate::ndarray::Ix1>()
1009            .expect("Operation failed");
1010        for i in 0..100 {
1011            assert_eq!(mul_array[i], (i as f64) * ((i + 5) as f64));
1012        }
1013
1014        // Test division (avoid division by zero)
1015        let div_result = mmap2
1016            .div(&mmap1.map_zero_copy(|x| x + 1.0).expect("Operation failed"))
1017            .expect("Test: operation failed");
1018        let div_array = div_result
1019            .readonlyarray::<crate::ndarray::Ix1>()
1020            .expect("Operation failed");
1021        for i in 0..100 {
1022            assert_eq!(div_array[i], ((i + 5) as f64) / ((i + 1) as f64));
1023        }
1024    }
1025
1026    #[test]
1027    fn test_broadcast_op() {
1028        // Create a temporary directory for our test files
1029        let dir = tempdir().expect("Operation failed");
1030        let file_path1 = dir.path().join("test_broadcast1.bin");
1031        let file_path2 = dir.path().join("test_broadcast2.bin");
1032
1033        // Create a 2D array (3x4) and a 1D array (4)
1034        let data1 = Array2::<f64>::from_shape_fn((3, 4), |(i, j)| (i * 4 + j) as f64);
1035        let data2 = crate::ndarray::Array1::from_vec((0..4).map(|i| (i + 1) as f64).collect());
1036
1037        // Save the arrays with proper headers using save_array
1038        MemoryMappedArray::<f64>::save_array(&data1, &file_path1, None).expect("Operation failed");
1039        MemoryMappedArray::<f64>::save_array(&data2, &file_path2, None).expect("Operation failed");
1040
1041        // Open the files for zero-copy operations
1042        let mmap1 = MemoryMappedArray::<f64>::open_zero_copy(&file_path1, AccessMode::ReadOnly)
1043            .expect("Operation failed");
1044        let mmap2 = MemoryMappedArray::<f64>::open_zero_copy(&file_path2, AccessMode::ReadOnly)
1045            .expect("Operation failed");
1046
1047        // Test broadcasting
1048        let result = mmap1
1049            .broadcast_op(&mmap2, |a, b| a * b)
1050            .expect("Operation failed");
1051
1052        // Verify the result
1053        let result_array = result
1054            .readonlyarray::<crate::ndarray::Ix2>()
1055            .expect("Operation failed");
1056        assert_eq!(result_array.shape(), &[3, 4]);
1057
1058        for i in 0..3 {
1059            for j in 0..4 {
1060                let expected = (i * 4 + j) as f64 * (j + 1) as f64;
1061                assert_eq!(result_array[[i, j]], expected);
1062            }
1063        }
1064    }
1065}