add inplace fft matrix by columns, phase correction and reimAbsolute

benchmark/reimMatrixFFTByColumnsComparison.ts

@@ -0,0 +1,108 @@
+/* eslint-disable no-console */
+import { reimMatrixFFT } from '../src/reimMatrix/reimMatrixFFT.ts';
+import { reimMatrixFFTByColumns } from '../src/reimMatrix/reimMatrixFFTByColumns.ts';
+import type { DataReImMatrix } from '../src/types/index.ts';
+
+/**
+ * Transpose a complex matrix (with separate re and im arrays)
+ * @param data
+ */
+function transposeComplexMatrix(data: DataReImMatrix): DataReImMatrix {
+  const { re, im } = data;
+  const numRows = re.length;
+  if (numRows === 0) return { re: [], im: [] };
+
+  const numCols = re[0].length;
+  const transRe: Float64Array[] = [];
+  const transIm: Float64Array[] = [];
+
+  for (let col = 0; col < numCols; col++) {
+    const newRow = new Float64Array(numRows);
+    const newIm = new Float64Array(numRows);
+    for (let row = 0; row < numRows; row++) {
+      newRow[row] = re[row][col];
+      newIm[row] = im[row][col];
+    }
+    transRe.push(newRow);
+    transIm.push(newIm);
+  }
+
+  return { re: transRe, im: transIm };
+}
+
+// Parameters
+const numRows = 512; // number of columns to transform
+const numCols = 4096; // length of each column (size of FFT)
+const iterations = 20;
+
+console.log(
+  `Matrix size: ${numRows} rows × ${numCols} columns (${numRows * numCols} complex numbers)`,
+);
+console.log(`Iterations: ${iterations}`);
+console.log('');
+
+// Build input data
+const inputData: DataReImMatrix = {
+  re: Array.from({ length: numRows }, () => {
+    const arr = new Float64Array(numCols);
+    for (let i = 0; i < numCols; i++) {
+      arr[i] = Math.random();
+    }
+    return arr;
+  }),
+  im: Array.from({ length: numRows }, () => {
+    const arr = new Float64Array(numCols);
+    for (let i = 0; i < numCols; i++) {
+      arr[i] = Math.random();
+    }
+    return arr;
+  }),
+};
+
+// Warmup runs
+reimMatrixFFTByColumns(inputData);
+const transposed = transposeComplexMatrix(inputData);
+reimMatrixFFT(transposed);
+
+console.log(
+  '--- Approach 1: reimMatrixFFTByColumns (Transform columns directly) ---',
+);
+{
+  let totalTime = 0;
+  for (let i = 0; i < iterations; i++) {
+    const start = performance.now();
+    reimMatrixFFTByColumns(inputData, { inPlace: true });
+    const elapsed = performance.now() - start;
+    totalTime += elapsed;
+  }
+  const avgTime = totalTime / iterations;
+  console.log(`Total time: ${totalTime.toFixed(2)} ms`);
+  console.log(`Average time per iteration: ${avgTime.toFixed(4)} ms`);
+  console.log(`Time per column: ${(avgTime / numRows).toFixed(6)} ms`);
+}
+
+console.log('');
+console.log(
+  '--- Approach 2: Transpose → reimMatrixFFT → Transpose (Transform rows via transposition) ---',
+);
+{
+  let totalTime = 0;
+  for (let i = 0; i < iterations; i++) {
+    const start = performance.now();
+    const transposed = transposeComplexMatrix(inputData);
+    const fftResult = reimMatrixFFT(transposed, { inPlace: true });
+    transposeComplexMatrix(fftResult);
+    const elapsed = performance.now() - start;
+    totalTime += elapsed;
+  }
+  const avgTime = totalTime / iterations;
+  console.log(`Total time: ${totalTime.toFixed(2)} ms`);
+  console.log(`Average time per iteration: ${avgTime.toFixed(4)} ms`);
+  console.log(`Time per column: ${(avgTime / numRows).toFixed(6)} ms`);
+}
+
+console.log('');
+console.log('Notes:');
+console.log('- Approach 1 is optimized for column-wise FFT transformation');
+console.log('- Approach 2 transposes twice, which may have overhead');
+console.log('- Both should produce mathematically equivalent results');

src/__tests__/__snapshots__/index.test.ts.snap

@@ -9,6 +9,8 @@ exports[`existence of exported functions 1`] = `
   "reimZeroFilling",
   "reimArrayFFT",
   "reimMatrixFFT",
+  "reimMatrixFFTByColumns",
+  "reimMatrixPhaseCorrection",
   "getOutputArray",
   "xAbsolute",
   "xAbsoluteMedian",

src/reim/__tests__/reimAbsolute.test.ts

@@ -9,3 +9,12 @@ test('reimAbsolute', () => {
 
   expect(result).toStrictEqual(Float64Array.from([0, 5, 10]));
 });
+
+test('reimAbsolute inplace on re', () => {
+  const re = [0, 3, 6];
+  const im = [0, 4, 8];
+  const result = reimAbsolute({ re, im }, { output: re });
+
+  expect(result).toStrictEqual([0, 5, 10]);
+  expect(re).toStrictEqual([0, 5, 10]);
+});

src/reim/__tests__/reimPhaseCorrection.test.ts

@@ -125,3 +125,105 @@ test('180 zero order phasing', () => {
 
   expect(diff[index]).toBeCloseTo(0, 4);
 });
+
+test('inPlace option: false (default) should create new arrays', () => {
+  const re = new Float64Array([0, 1, 2, 3]);
+  const im = new Float64Array([0, 1, 2, 3]);
+  const originalRe = Array.from(re);
+  const originalIm = Array.from(im);
+
+  const result = reimPhaseCorrection({ re, im }, Math.PI / 4, 0, {
+    inPlace: false,
+  });
+
+  // Input arrays should not be modified
+  expect(Array.from(re)).toStrictEqual(originalRe);
+  expect(Array.from(im)).toStrictEqual(originalIm);
+
+  // Result should be different objects
+  expect(result.re).not.toBe(re);
+  expect(result.im).not.toBe(im);
+
+  // Result should have correct values
+  expect(result.re).toHaveLength(4);
+  expect(result.im).toHaveLength(4);
+});
+
+test('inPlace option: true should modify input arrays', () => {
+  const re = new Float64Array([0, 1, 2, 3]);
+  const im = new Float64Array([0, 1, 2, 3]);
+
+  const result = reimPhaseCorrection({ re, im }, Math.PI / 4, 0, {
+    inPlace: true,
+  });
+
+  // Result should reference the same arrays
+  expect(result.re).toBe(re);
+  expect(result.im).toBe(im);
+
+  // Input arrays should be modified
+  expect(Array.from(re)).not.toStrictEqual([0, 1, 2, 3]);
+  expect(Array.from(im)).not.toStrictEqual([0, 1, 2, 3]);
+});
+
+test('inPlace option: true and false should produce same values', () => {
+  const phi0 = Math.PI / 3;
+  const phi1 = Math.PI / 6;
+
+  const re1 = new Float64Array([1, 2, 3, 4, 5]);
+  const im1 = new Float64Array([1, 2, 3, 4, 5]);
+  const result1 = reimPhaseCorrection({ re: re1, im: im1 }, phi0, phi1, {
+    inPlace: false,
+  });
+
+  const re2 = new Float64Array([1, 2, 3, 4, 5]);
+  const im2 = new Float64Array([1, 2, 3, 4, 5]);
+  const result2 = reimPhaseCorrection({ re: re2, im: im2 }, phi0, phi1, {
+    inPlace: true,
+  });
+
+  // Values should be the same regardless of inPlace option
+  for (let i = 0; i < result1.re.length; i++) {
+    expect(result1.re[i]).toBeCloseTo(result2.re[i], 10);
+    expect(result1.im[i]).toBeCloseTo(result2.im[i], 10);
+  }
+});
+
+test('inPlace option with reverse flag', () => {
+  const re1 = new Float64Array([1, 2, 3, 4]);
+  const im1 = new Float64Array([1, 2, 3, 4]);
+  const result1 = reimPhaseCorrection({ re: re1, im: im1 }, Math.PI / 4, 0.5, {
+    inPlace: false,
+    reverse: true,
+  });
+
+  const re2 = new Float64Array([1, 2, 3, 4]);
+  const im2 = new Float64Array([1, 2, 3, 4]);
+  const result2 = reimPhaseCorrection({ re: re2, im: im2 }, Math.PI / 4, 0.5, {
+    inPlace: true,
+    reverse: true,
+  });
+
+  // Result should reference the same arrays when inPlace: true
+  expect(result2.re).toBe(re2);
+  expect(result2.im).toBe(im2);
+
+  // Values should be identical
+  for (let i = 0; i < result1.re.length; i++) {
+    expect(result1.re[i]).toBeCloseTo(result2.re[i], 10);
+    expect(result1.im[i]).toBeCloseTo(result2.im[i], 10);
+  }
+});
+
+test('inPlace option with zero phase correction', () => {
+  const re = new Float64Array([1, 2, 3]);
+  const im = new Float64Array([1, 2, 3]);
+
+  const result = reimPhaseCorrection({ re, im }, 0, 0, { inPlace: true });
+
+  // With zero phase correction, values should remain unchanged
+  expect(result.re).toBe(re);
+  expect(result.im).toBe(im);
+  expect(Array.from(result.re)).toStrictEqual([1, 2, 3]);
+  expect(Array.from(result.im)).toStrictEqual([1, 2, 3]);
+});

src/reim/reimAbsolute.ts

@@ -1,18 +1,25 @@
+import type { DoubleArray } from 'cheminfo-types';
+
 import type { DataReIm } from '../types/index.ts';
 
 /**
  * Calculates reimAbsolute value of a complex spectrum.
  * @param data - complex spectrum
+ * @param options
+ * @param options.output
  * @returns - reimAbsolute value
  */
-export function reimAbsolute(data: DataReIm): Float64Array<ArrayBuffer> {
+export function reimAbsolute(
+  data: DataReIm,
+  options: { output?: DoubleArray } = {},
+): DoubleArray {
   const length = data.re.length;
+  const { output = new Float64Array(length) } = options;
   const re = data.re;
   const im = data.im;
-  const newArray = new Float64Array(length);
   for (let i = 0; i < length; i++) {
-    newArray[i] = Math.hypot(re[i], im[i]);
+    output[i] = Math.hypot(re[i], im[i]);
   }
 
-  return newArray;
+  return output;
 }

src/reim/reimAutoPhaseCorrection.ts

@@ -36,6 +36,11 @@ export interface AutoPhaseCorrectionOptions {
    * @default false
    */
   reverse?: boolean;
+  /**
+   * Apply the phase correction directly in the input data
+   * @default false
+   */
+  inPlace?: boolean;
 }
 
 /**
@@ -48,13 +53,14 @@ export interface AutoPhaseCorrectionOptions {
 export function reimAutoPhaseCorrection(
   data: DataReIm,
   options: AutoPhaseCorrectionOptions = {},
-): { data: DataReIm<Float64Array>; ph0: number; ph1: number } {
+): { data: DataReIm; ph0: number; ph1: number } {
   const {
     magnitudeMode = true,
     minRegSize = 30,
     factorNoise = 3,
     maxDistanceToJoin = 256,
     reverse = false,
+    inPlace = false,
   } = options;
 
   const finalPeaks = detectBaselineRegions(data, {
@@ -96,7 +102,7 @@ export function reimAutoPhaseCorrection(
     { re, im },
     toRadians(ph0),
     toRadians(ph1),
-    { reverse },
+    { reverse, inPlace },
   );
 
   return { data: phased, ph0, ph1 };

src/reim/reimPhaseCorrection.ts

@@ -1,7 +1,16 @@
 import type { DataReIm } from '../types/index.ts';
 
 export interface ReimPhaseCorrectionOptions {
+  /**
+   * Apply the phase correction from the last element of the array
+   * @default false
+   */
   reverse?: boolean;
+  /**
+   * Apply the phase correction directly in the input data
+   * @default false
+   */
+  inPlace?: boolean;
 }
 
 /**
@@ -10,46 +19,56 @@ export interface ReimPhaseCorrectionOptions {
  * @param phi0 - Angle in radians for zero order phase correction
  * @param phi1 - Angle in radians for first order phase correction
  * @param options
- * @returns - returns a new object {re:[], im:[]}
+ * @returns - returns a new object {re:[], im:[]} unless inPlace=true
  */
 export function reimPhaseCorrection(
   data: DataReIm,
   phi0 = 0,
   phi1 = 0,
   options: ReimPhaseCorrectionOptions = {},
-): DataReIm<Float64Array> {
-  const { reverse = false } = options;
+): DataReIm {
+  const { reverse = false, inPlace = false } = options;
 
   phi0 = Number.isFinite(phi0) ? phi0 : 0;
   phi1 = Number.isFinite(phi1) ? phi1 : 0;
 
+  const length = data.re.length;
+
+  // Decide target arrays
   const re = data.re;
   const im = data.im;
-  const length = data.re.length;
+
+  const outRe = inPlace ? re : new Float64Array(length);
+  const outIm = inPlace ? im : new Float64Array(length);
 
   let firstAngle = phi0;
   let delta = phi1 / length;
+
   if (reverse) {
     delta *= -1;
     firstAngle += phi1;
   }
 
   const alpha = 2 * Math.sin(delta / 2) ** 2;
   const beta = Math.sin(delta);
+
   let cosTheta = Math.cos(firstAngle);
   let sinTheta = Math.sin(firstAngle);
 
-  const newRe = new Float64Array(length);
-  const newIm = new Float64Array(length);
   for (let i = 0; i < length; i++) {
-    newRe[i] = re[i] * cosTheta - im[i] * sinTheta;
-    newIm[i] = im[i] * cosTheta + re[i] * sinTheta;
-    // calculate angles i+1 from i
+    const r = re[i];
+    const ii = im[i];
+
+    outRe[i] = r * cosTheta - ii * sinTheta;
+    outIm[i] = ii * cosTheta + r * sinTheta;
+
+    // Recursive angle update (stable incremental rotation)
     const newCosTheta = cosTheta - (alpha * cosTheta + beta * sinTheta);
     const newSinTheta = sinTheta - (alpha * sinTheta - beta * cosTheta);
+
     cosTheta = newCosTheta;
     sinTheta = newSinTheta;
   }
 
-  return { re: newRe, im: newIm };
+  return { re: outRe, im: outIm };
 }