dfa.js

/* 
    Detrended Fluctuation Analysis (DFA)
    Javascript port
    available on https://github.com/deemeetree/dfa

    based on the approach described on https://noduslabs.com/featured/fractal-variability-feedback-system/
    based on the python script https://github.com/dokato/dfa


    ~ Concept 

    DFA is used to measure the behaviour of a time series. 
    
    An obtained alpha component (closely related to the Hurst exponent) will indicate 
    the presence of correlations in the time series. 

    It is based on the relationship between the length of an observation and cumulated variability.
    

    ~ Algorithm:

    1) represent a time series as a one-dimensional vector
    2) transform it into cumulative sum of variances from the mean
    3) generate the scales of different window sizes for the time series — each containing the number of observations
    4) split the cumulative variance vector into those chunks
    5) for each chunk: generate polynomial in the range in a window (to detrend it)
    6) calculate RMS of the fluctuations of the original from the fit
    7) get the average RMS value for each window in the scale
    8) take the average for each squared RMS for each scale
    9) if they align along in a straight line in loglog plot, there is power law relation
    i.e. as observations increase in length, the amplitude of average fluctuations increases as well 
    that is, there is an exponential growth of fluctuations when there's an exponential growth of scale (2^n)
    on the smaller scales the deviations are smaller. on the bigger scales they are much bigger.

    ~ Simple use

    let time_series = [8, 10, 6, 9, 7, 5, 5, 11, 11, 8, 6, 7, 9, 10, 7, 9]
    let dfa = new DFA(time_series)
    let alpha_component = dfa.compute()
    console.log(alpha_component)

    */

let DFA = (function () {
	// Helper functions
	function mean(x) {
		return x.reduce((a, b) => a + b, 0) / x.length;
	}

	function variancePopulation(x) {
		const m = mean(x);
		let s = 0;
		for (let i = 0; i < x.length; i++) {
			const d = x[i] - m;
			s += d * d;
		}
		return s / x.length; // population variance
	}

	function sdnn(x) {
		return Math.sqrt(variancePopulation(x));
	}

	function rmssd(x) {
		if (x.length < 2) return 0;
		let s = 0;
		for (let i = 1; i < x.length; i++) {
			const d = x[i] - x[i - 1];
			s += d * d;
		}
		return Math.sqrt(s / (x.length - 1));
	}

	function pnn50(x) {
		if (x.length < 2) return 0;
		let c = 0;
		for (let i = 1; i < x.length; i++) {
			if (Math.abs(x[i] - x[i - 1]) > 50) c++;
		}
		return (c / (x.length - 1)) * 100;
	}

	function integratedProfile(rr) {
		const m = mean(rr);
		const y = new Array(rr.length);
		let acc = 0;
		for (let i = 0; i < rr.length; i++) {
			acc += rr[i] - m;
			y[i] = acc;
		}
		return { profile: y, meanValue: m };
	}

	function generateScales(
		N,
		{
			minWindow = 4,
			step = 2, // linear step for small scales (and short series)
			expStep = 0.25, // log2 spacing for geometric progression (≈√2)
			shortMax = 16,
			longMin = 16,
			longMaxFraction = 0.25,
		} = {}
	) {
		// Use the full data length (N) for global alpha, no hard cap
		const maxReli = N;
		const S = new Set();

		// α1 anchors: minWindow..shortMax inclusive, in +step increments
		for (let s = minWindow; s <= Math.min(shortMax, maxReli); s += step)
			S.add(s);

		// α2 range: add linear steps from 16 to N/4 with step 2
		const alpha2Start = 16;
		const alpha2Max = Math.floor(N * longMaxFraction);
		for (let s = alpha2Start; s <= Math.min(alpha2Max, maxReli); s += 2) {
			S.add(s);
		}

		// Add geometric progression for scales beyond alpha2Max (for global alpha only)
		// Start geometric progression after the alpha2 range ends
		let g = [];
		if (alpha2Max < maxReli) {
			let s = alpha2Max;
			const mult = Math.pow(2, expStep);
			// Start from next geometric step after alpha2Max
			s = s * mult;
			while (s <= maxReli) {
				g.push(Math.round(s));
				s *= mult;
			}
		}

		// Add all geometric scales (they'll be used for global alpha)
		for (const s of g) S.add(s);

		// Ensure at least 3 scales overall (global fit needs ≥3 points)
		let scales = Array.from(S).filter((s) => s >= minWindow && s <= maxReli);
		scales.sort((a, b) => a - b);
		if (scales.length < 3) {
			// minimal fallback
			scales = [
				minWindow,
				Math.max(minWindow + step, Math.floor((minWindow + maxReli) / 2)),
				maxReli,
			]
				.filter((v, i, a) => a.indexOf(v) === i)
				.sort((a, b) => a - b);
		}
		return scales;
	}

	function linearRegression(x, y) {
		const n = x.length;
		if (n === 0) return { slope: 0, intercept: 0 };
		let sx = 0,
			sy = 0,
			sxx = 0,
			sxy = 0;
		for (let i = 0; i < n; i++) {
			sx += x[i];
			sy += y[i];
			sxx += x[i] * x[i];
			sxy += x[i] * y[i];
		}
		const denom = n * sxx - sx * sx;
		if (denom === 0) return { slope: 0, intercept: sy / n };
		const slope = (n * sxy - sx * sy) / denom;
		const intercept = (sy - slope * sx) / n;
		return { slope, intercept };
	}

	function fitQuadratic(x, y) {
		// Least-squares fit y = a + b*x + c*x^2 via 3x3 normal equations.
		const n = x.length;
		if (n < 3) return null;
		let S0 = n,
			S1 = 0,
			S2 = 0,
			S3 = 0,
			S4 = 0,
			T0 = 0,
			T1 = 0,
			T2 = 0;
		for (let i = 0; i < n; i++) {
			const xi = x[i],
				yi = y[i],
				x2 = xi * xi;
			S1 += xi;
			S2 += x2;
			S3 += x2 * xi;
			S4 += x2 * x2;
			T0 += yi;
			T1 += xi * yi;
			T2 += x2 * yi;
		}
		// Augmented matrix for Gaussian elimination with partial pivoting.
		const m = [
			[S0, S1, S2, T0],
			[S1, S2, S3, T1],
			[S2, S3, S4, T2],
		];
		for (let col = 0; col < 3; col++) {
			let piv = col;
			for (let r = col + 1; r < 3; r++)
				if (Math.abs(m[r][col]) > Math.abs(m[piv][col])) piv = r;
			if (Math.abs(m[piv][col]) < 1e-12) return null; // singular
			const tmp = m[col];
			m[col] = m[piv];
			m[piv] = tmp;
			for (let r = 0; r < 3; r++) {
				if (r === col) continue;
				const f = m[r][col] / m[col][col];
				for (let k = col; k < 4; k++) m[r][k] -= f * m[col][k];
			}
		}
		return { a: m[0][3] / m[0][0], b: m[1][3] / m[1][1], c: m[2][3] / m[2][2] };
	}

	// Characterise the shape of an h(q) curve: overall direction, curvature
	// (U vs inverted-U), deviation from a straight line, and where h is lowest.
	function analyzeHCurve(qValues, hq) {
		const xs = [],
			ys = [];
		for (let i = 0; i < hq.length; i++) {
			if (hq[i] !== null && hq[i] !== undefined) {
				xs.push(qValues[i]);
				ys.push(hq[i]);
			}
		}
		const empty = {
			slope: null,
			curvature: null,
			nonlinearity: null,
			hMinLocation: null,
		};
		if (xs.length < 3) return empty;

		// General direction: slope of the best-fit straight line through h(q).
		const { slope, intercept } = linearRegression(xs, ys);

		// Curvature: 2c from the quadratic fit (the 2nd derivative).
		// > 0 → U-shape (convex), < 0 → inverted-U (concave).
		const quad = fitQuadratic(xs, ys);
		const curvature = quad ? 2 * quad.c : null;

		// Nonlinearity: RMS of residuals around the straight-line fit.
		let ss = 0;
		for (let i = 0; i < xs.length; i++) {
			const r = ys[i] - (intercept + slope * xs[i]);
			ss += r * r;
		}
		const nonlinearity = Math.sqrt(ss / xs.length);

		// hMinLocation: the q where h(q) is smallest (least persistent).
		let minH = ys[0],
			hMinLocation = xs[0];
		for (let i = 1; i < xs.length; i++) {
			if (ys[i] < minH) {
				minH = ys[i];
				hMinLocation = xs[i];
			}
		}

		return { slope, curvature, nonlinearity, hMinLocation };
	}

	function detrendLinear(ySeg) {
		// returns fitted values for the segment
		const n = ySeg.length;
		const x = new Array(n);
		for (let i = 0; i < n; i++) x[i] = i;
		const { slope, intercept } = linearRegression(x, ySeg);
		const trend = new Array(n);
		for (let i = 0; i < n; i++) trend[i] = slope * i + intercept;
		return trend;
	}

	function fluctuationForScale(profile, s, minWindow = 4, q = 2) {
		const n = profile.length;
		const numSegments = Math.floor(n / s);
		if (numSegments <= 0) return { F: 0, segments: 0 };

		const rms = [];

		// forward segments
		for (let seg = 0; seg < numSegments; seg++) {
			const start = seg * s,
				end = start + s;
			const segData = profile.slice(start, end);
			const trend = detrendLinear(segData);
			let ss = 0;
			for (let i = 0; i < s; i++) {
				const r = segData[i] - trend[i];
				ss += r * r;
			}
			rms.push(Math.sqrt(ss / s));
		}

		// backward segments if remainder is sufficiently long
		const remainder = n % s;
		if (remainder >= minWindow) {
			for (let seg = 0; seg < numSegments; seg++) {
				const end = n - seg * s;
				const start = end - s;
				if (start < 0) break;
				const segData = profile.slice(start, end);
				const trend = detrendLinear(segData);
				let ss = 0;
				for (let i = 0; i < s; i++) {
					const r = segData[i] - trend[i];
					ss += r * r;
				}
				rms.push(Math.sqrt(ss / s));
			}
		}

		if (rms.length === 0) return { F: 0, segments: numSegments };

		// q-th order fluctuation function (MFDFA generalization).
		// q === 2 reproduces canonical DFA: F(s) = sqrt( mean( RMS^2 ) ).
		// q === 0 uses the logarithmic-averaging special case (1/q is undefined).
		let F;
		if (Math.abs(q) < 1e-10) {
			// F_0(s) = exp( (1/2N_s) * sum( ln( RMS^2 ) ) )
			let sumLog = 0;
			for (let i = 0; i < rms.length; i++) {
				const sq = rms[i] * rms[i];
				// guard against log(0) for degenerate (perfectly-fit) segments
				sumLog += Math.log(sq > 0 ? sq : Number.MIN_VALUE);
			}
			F = Math.exp(sumLog / (2 * rms.length));
		} else {
			// F_q(s) = ( (1/N_s) * sum( (RMS^2)^(q/2) ) )^(1/q)
			let sumPow = 0;
			for (let i = 0; i < rms.length; i++) {
				sumPow += Math.pow(rms[i] * rms[i], q / 2);
			}
			F = Math.pow(sumPow / rms.length, 1 / q);
		}

		return { F, segments: numSegments };
	}

	function fitAlphaInRange(
		scales,
		logs,
		range,
		segments = null,
		minSegments = 0
	) {
		const xs = [],
			ys = [],
			used = [];
		for (let i = 0; i < scales.length; i++) {
			const s = scales[i];
			const ylog = logs.fluct[i];
			if (s < range[0] || s > range[1]) continue;
			if (ylog === null) continue; // filtered out (e.g., F==0)
			if (segments && segments[i] < minSegments) continue;
			xs.push(logs.scale[i]);
			ys.push(ylog);
			used.push(s);
		}
		if (xs.length < 3) return { alpha: null, usedRange: null, usedScales: [] };
		const { slope } = linearRegression(xs, ys);
		return {
			alpha: slope,
			usedRange: [Math.min(...used), Math.max(...used)],
			usedScales: used,
		};
	}

	function alphaScoreNumeric(alpha) {
		const dev = Math.abs(alpha - 1.0);
		return Math.max(0, (1.0 - dev) * 100);
	}

	function getThresholds(level) {
		if (level === "relaxed") return [0.65, 0.9, 1.1];
		if (level === "strict") return [0.5, 0.98, 1.02];
		return [0.55, 0.95, 1.05]; // moderate (default)
	}

	function alphaToLabel(alpha, level = "moderate") {
		const t = getThresholds(level);
		if (alpha <= t[0]) return "random";
		if (alpha < t[1]) return "regular";
		if (alpha <= t[2]) return "fractal";
		return "complex";
	}

	function getDfaThresholds(level) {
		if (level === "relaxed") return [0.40, 0.60, 0.90, 1.10, 1.40];
		if (level === "strict") return [0.48, 0.52, 0.98, 1.02, 1.48];
		return [0.45, 0.55, 0.95, 1.05, 1.45]; // moderate (default)
	}

	function alphaToDfaLabel(alpha, level = "moderate") {
		const t = getDfaThresholds(level);
		if (alpha < t[0]) return "anti-correlated";
		if (alpha <= t[1]) return "white noise";
		if (alpha < t[2]) return "correlated";
		if (alpha <= t[3]) return "1/f noise";
		if (alpha < t[4]) return "strongly correlated";
		return "Brownian motion";
	}

	// DFA class definition
	let DFA = function (x) {
		if (!(this instanceof DFA)) {
			return new DFA(x);
		}

		if (
			!(
				x instanceof Array ||
				x instanceof Float32Array ||
				x instanceof Float64Array
			)
		) {
			throw new Error("x must be an array");
		}

		this.x = x;
	};

	// Legacy methods kept for backward compatibility
	DFA.prototype.cumsumVariance = function (x, mean) {
		let cumulativeSumVariance = (
			(sum) => (value) =>
				(sum += value - mean)
		)(0);
		return x.map(cumulativeSumVariance);
	};

	DFA.prototype.averageVariance = function (x, mean) {
		// Note: This is actually mean absolute deviation, kept for compatibility
		let totalDeviation = 0;
		for (let i = 0; i < x.length; i++) {
			totalDeviation += Math.abs(x[i] - mean);
		}
		return totalDeviation / x.length;
	};

	DFA.prototype.deviationsFromMean = function (x, mean) {
		const deviationsVector = x.map((value) => value - mean);
		return deviationsVector;
	};

	DFA.prototype.squareVector = function (x) {
		return x.map((value) => Math.pow(value, 2));
	};

	DFA.prototype.SDNN = function (x) {
		return sdnn(x);
	};

	DFA.prototype.averageDifferences = function (x) {
		let sumOfDifferences = 0;
		for (let i = 0; i < x.length - 1; i++) {
			sumOfDifferences += Math.abs(x[i + 1] - x[i]);
		}
		return sumOfDifferences / (x.length - 1);
	};

	DFA.prototype.RMSSD = function (x) {
		return rmssd(x);
	};

	DFA.prototype.naturalLog = function (value) {
		return Math.log(value);
	};

	DFA.prototype.PNN50 = function (rrIntervals) {
		return pnn50(rrIntervals);
	};

	DFA.prototype.meanOfVector = function (x) {
		let y = x || this.x;
		return mean(y);
	};

	DFA.prototype.meanOfVectorNoBias = function (x) {
		let y = x || this.x;
		let meanNoBias = (arr) => arr.reduce((p, c) => p + c, 0) / (arr.length - 1);
		return meanNoBias(y);
	};

	DFA.prototype.alphaScore = function (alpha, level = "moderate") {
		const t = getThresholds(level);
		if (alpha <= t[0]) return "recovering";
		if (alpha < t[1]) return "regular";
		if (alpha <= t[2]) return "resilient";
		return "tension";
	};

	DFA.prototype.compute = function (
		minWindow = 4,
		expStep = 0.25, // geometric log2 step for longer series (≈√2)
		step = 2, // linear increment for α1 and for short series fallback
		shortMax = 16,
		longMin = 16,
		longMaxFraction = 0.25,
		level = "moderate" // threshold sensitivity: "relaxed", "moderate", or "strict"
	) {
		const rr = this.x;
		const N = rr.length;
		if (N <= minWindow) {
			return {
				averageVariance: 0,
				meanValue: 0,
				lengthOfData: N,
				SDNN: 0,
				RMSSD: 0,
				lnRMSSD: 0,
				PNN50: 0,
				averageDifferences: 0,
				scales: [],
				segments: [],
				fluctuations: [],
				scalesLog: [],
				fluctuationsLog: [],
				scales_log: [], // Legacy alias
				fluctuations_log: [], // Legacy alias
				coefficients: { slope: 0, intercept: 0 },
				alpha: 1.0,
				alphaScore: "resilient",
				alphaScoreNumeric: 100,
				alphaLabel: "fractal",
				dfaLabel: "1/f noise",
				alpha1: null,
				alpha1Score: null,
				alpha1ScoreNumeric: null,
				alpha1Label: null,
				dfa1Label: null,
				alpha2: null,
				alpha2Score: null,
				alpha2ScoreNumeric: null,
				alpha2Label: null,
				dfa2Label: null,
				alpha1Range: null,
				alpha2Range: null,
			};
		}

		// Basic stats (parity with Swift)
		const meanValue = mean(rr);
		const avgVar = variancePopulation(rr);
		const SDNN = Math.sqrt(avgVar);
		const RMSSD = rmssd(rr);
		const lnRMSSD = RMSSD > 0 ? Math.log(RMSSD) : 0;
		const PNN50 = pnn50(rr);
		// Keep averageDifferences to match legacy JS result
		let sumAbs = 0;
		for (let i = 1; i < N; i++) sumAbs += Math.abs(rr[i] - rr[i - 1]);
		const averageDifferences = N > 1 ? sumAbs / (N - 1) : 0;

		// Profile
		const { profile } = integratedProfile(rr);

		// Scales
		const scales = generateScales(N, {
			minWindow,
			step,
			expStep,
			shortMax,
			longMin,
			longMaxFraction,
		});
		const segments = new Array(scales.length);
		const fluctuations = new Array(scales.length);

		// F(s) per scale (forward+backward); record forward segments only
		for (let i = 0; i < scales.length; i++) {
			const s = scales[i];
			const { F, segments: segs } = fluctuationForScale(profile, s, minWindow);
			fluctuations[i] = F;
			segments[i] = segs; // parity with Swift (forward only)
		}

		// Filter zeros before logs
		const scalesLog = new Array(scales.length);
		const fluctuationsLog = new Array(scales.length).fill(null);
		for (let i = 0; i < scales.length; i++) {
			scalesLog[i] = Math.log(scales[i]); // natural log
			if (fluctuations[i] > 0) fluctuationsLog[i] = Math.log(fluctuations[i]);
		}

		// Global α on valid (non-null) points
		const xAll = [],
			yAll = [];
		for (let i = 0; i < scales.length; i++) {
			if (fluctuationsLog[i] !== null) {
				xAll.push(scalesLog[i]);
				yAll.push(fluctuationsLog[i]);
			}
		}
		const coefficients = linearRegression(xAll, yAll);
		const alpha = coefficients.slope;
		const alphaScoreVal = alphaScoreNumeric(alpha);
		const alphaScoreCat = this.alphaScore(alpha, level);

		// α1 (minWindow..shortMax)
		const {
			alpha: alpha1,
			usedRange: alpha1Range,
			usedScales: scalesAlpha1,
		} = fitAlphaInRange(
			scales,
			{ scale: scalesLog, fluct: fluctuationsLog },
			[minWindow, Math.min(shortMax, N)],
			null,
			0
		);

		// α2 (16..maxReliable), require ≥4 forward segments
		// Start from 16 and cap at N/4 for alpha2 reliability
		const alpha2Start = 16;
		const maxReliableAlpha2 = Math.floor(N * 0.25);
		const {
			alpha: alpha2,
			usedRange: alpha2Range,
			usedScales: scalesAlpha2,
		} = fitAlphaInRange(
			scales,
			{ scale: scalesLog, fluct: fluctuationsLog },
			[alpha2Start, maxReliableAlpha2],
			segments,
			4
		);

		return {
			averageVariance: avgVar,
			meanValue,
			lengthOfData: N,
			SDNN,
			RMSSD,
			lnRMSSD,
			PNN50,
			averageDifferences,
			scales,
			segments,
			scalesAlpha1,
			scalesAlpha2,
			fluctuations,
			scalesLog,
			fluctuationsLog,
			scales_log: scalesLog, // Legacy alias
			fluctuations_log: fluctuationsLog, // Legacy alias
			coefficients,
			alpha,
			alphaScore: alphaScoreCat,
			alphaScoreNumeric: alphaScoreVal,
			alphaLabel: alphaToLabel(alpha, level),
			dfaLabel: alphaToDfaLabel(alpha, level),
			alpha1,
			alpha1Score: alpha1 !== null ? this.alphaScore(alpha1, level) : null,
			alpha1ScoreNumeric: alpha1 !== null ? alphaScoreNumeric(alpha1) : null,
			alpha1Label: alpha1 !== null ? alphaToLabel(alpha1, level) : null,
			dfa1Label: alpha1 !== null ? alphaToDfaLabel(alpha1, level) : null,
			alpha2,
			alpha2Score: alpha2 !== null ? this.alphaScore(alpha2, level) : null,
			alpha2ScoreNumeric: alpha2 !== null ? alphaScoreNumeric(alpha2) : null,
			alpha2Label: alpha2 !== null ? alphaToLabel(alpha2, level) : null,
			dfa2Label: alpha2 !== null ? alphaToDfaLabel(alpha2, level) : null,
			alpha1Range,
			alpha2Range,
		};
	};

	// Multifractal DFA (MFDFA).
	// Sweeps the q-order parameter to obtain the generalized Hurst exponent h(q),
	// the mass exponent tau(q), and the singularity spectrum f(alpha).
	// q === 2 reproduces the monofractal `compute()` alpha (returned as hq2).
	DFA.prototype.computeMultifractal = function ({
		qMin = -5,
		qMax = 5,
		qStep = 0.5,
		minWindow = 4,
		expStep = 0.25,
		step = 2,
		shortMax = 16,
		longMin = 16,
		longMaxFraction = 0.25,
	} = {}) {
		const rr = this.x;
		const N = rr.length;

		// Build the q list (inclusive of qMax up to floating-point tolerance).
		const qValues = [];
		for (let q = qMin; q <= qMax + 1e-9; q += qStep) {
			// snap tiny floating-point drift so e.g. q=0 is exactly 0
			qValues.push(Math.abs(q) < 1e-10 ? 0 : Math.round(q * 1e6) / 1e6);
		}

		if (N <= minWindow) {
			return {
				q: qValues,
				hq: qValues.map(() => null),
				hq1: qValues.map(() => null),
				hq2: qValues.map(() => null),
				tau: qValues.map(() => null),
				alpha: qValues.map(() => null),
				falpha: qValues.map(() => null),
				hMin: null,
				hMax: null,
				multifractalWidth: null,
				width1: null,
				width2: null,
				hCurveSlope: null,
				hCurveCurvature: null,
				hCurveNonlinearity: null,
				hMinLocation: null,
				hCurveSlope1: null,
				hCurveCurvature1: null,
				hCurveNonlinearity1: null,
				hMinLocation1: null,
				hCurveSlope2: null,
				hCurveCurvature2: null,
				hCurveNonlinearity2: null,
				hMinLocation2: null,
				monofractal: { alpha: null, alpha1: null, alpha2: null },
				ranges: { global: null, alpha1: null, alpha2: null },
				scales: [],
				fluctuationsByQ: qValues.map(() => []),
				lengthOfData: N,
			};
		}

		// Same integrated profile and scales as the monofractal path.
		const { profile } = integratedProfile(rr);
		const scales = generateScales(N, {
			minWindow,
			step,
			expStep,
			shortMax,
			longMin,
			longMaxFraction,
		});
		const scalesLog = scales.map((s) => Math.log(s));

		// Fit ranges mirror compute(): global (all scales), α₁ (short), α₂ (long).
		const alpha1Range = [minWindow, Math.min(shortMax, N)];
		const alpha2Range = [16, Math.floor(N * 0.25)];

		const hq = []; // global h(q): fitted across all valid scales
		const hq1 = []; // α₁ h(q): fitted over short scales (4..shortMax)
		const hq2 = []; // α₂ h(q): fitted over long scales (16..N/4)
		const fluctuationsByQ = [];

		for (let qi = 0; qi < qValues.length; qi++) {
			const q = qValues[qi];
			const fluct = new Array(scales.length);
			const fluctLog = new Array(scales.length).fill(null);
			const segs = new Array(scales.length);
			const xs = [];
			const ys = [];
			for (let i = 0; i < scales.length; i++) {
				const { F, segments } = fluctuationForScale(
					profile,
					scales[i],
					minWindow,
					q
				);
				fluct[i] = F;
				segs[i] = segments;
				if (F > 0) {
					fluctLog[i] = Math.log(F);
					xs.push(scalesLog[i]);
					ys.push(fluctLog[i]);
				}
			}
			fluctuationsByQ.push(fluct);

			// global fit
			hq.push(xs.length >= 3 ? linearRegression(xs, ys).slope : null);

			// α₁ and α₂ fits via the same range-restricted helper compute() uses
			const logs = { scale: scalesLog, fluct: fluctLog };
			hq1.push(fitAlphaInRange(scales, logs, alpha1Range, null, 0).alpha);
			hq2.push(fitAlphaInRange(scales, logs, alpha2Range, segs, 4).alpha);
		}

		// tau(q) = q * h(q) - 1
		const tau = qValues.map((q, i) => (hq[i] !== null ? q * hq[i] - 1 : null));

		// Singularity spectrum via numerical Legendre transform of tau(q):
		//   alpha = d(tau)/dq ,  f(alpha) = q * alpha - tau(q)
		const alpha = new Array(qValues.length).fill(null);
		const falpha = new Array(qValues.length).fill(null);
		for (let i = 0; i < qValues.length; i++) {
			if (tau[i] === null) continue;
			// central difference on q where neighbours are available, else one-sided
			let lo = i,
				hi = i;
			if (i > 0 && tau[i - 1] !== null) lo = i - 1;
			if (i < qValues.length - 1 && tau[i + 1] !== null) hi = i + 1;
			if (hi === lo) continue; // no usable neighbour → can't differentiate
			const dAlpha =
				(tau[hi] - tau[lo]) / (qValues[hi] - qValues[lo]);
			alpha[i] = dAlpha;
			falpha[i] = qValues[i] * dAlpha - tau[i];
		}

		// Multifractal width = spread of h(q). h(q) decreases with q, so for the
		// global curve hMax sits at qMin and hMin at qMax. ≈0 → monofractal.
		const widthOf = (arr) => {
			const v = arr.filter((h) => h !== null);
			return v.length ? Math.max(...v) - Math.min(...v) : null;
		};
		const validH = hq.filter((h) => h !== null);
		const hMax = validH.length ? Math.max(...validH) : null;
		const hMin = validH.length ? Math.min(...validH) : null;
		const multifractalWidth = widthOf(hq);
		const width1 = widthOf(hq1);
		const width2 = widthOf(hq2);

		// Shape of each h(q) curve: direction, curvature, nonlinearity,
		// and the q at which persistence is lowest.
		const shape = analyzeHCurve(qValues, hq);
		const shape1 = analyzeHCurve(qValues, hq1);
		const shape2 = analyzeHCurve(qValues, hq2);

		// h(2) for each curve — should match compute()'s alpha / alpha1 / alpha2
		// (a built-in sanity check that the multifractal path is consistent).
		const i2 = qValues.indexOf(2);
		const monofractal = {
			alpha: i2 >= 0 ? hq[i2] : null,
			alpha1: i2 >= 0 ? hq1[i2] : null,
			alpha2: i2 >= 0 ? hq2[i2] : null,
		};

		return {
			q: qValues,
			hq, // global h(q), per q
			hq1, // α₁-range h(q) (short scales), per q
			hq2, // α₂-range h(q) (long scales), per q — null where insufficient
			tau, // mass exponent τ(q) from the global curve
			alpha, // Hölder exponent α from the global curve
			falpha, // singularity spectrum f(α) from the global curve
			hMin,
			hMax,
			multifractalWidth, // width of the global h(q) curve (extremes)
			width1, // width of the α₁ h(q) curve
			width2, // width of the α₂ h(q) curve
			// shape of the global h(q) curve
			hCurveSlope: shape.slope, // general direction of h(q) vs q
			hCurveCurvature: shape.curvature, // >0 U-shape, <0 inverted-U
			hCurveNonlinearity: shape.nonlinearity, // RMS deviation from a straight line
			hMinLocation: shape.hMinLocation, // q where h(q) (persistence) is lowest
			// shape of the α₁ (short-scale) h(q) curve
			hCurveSlope1: shape1.slope,
			hCurveCurvature1: shape1.curvature,
			hCurveNonlinearity1: shape1.nonlinearity,
			hMinLocation1: shape1.hMinLocation,
			// shape of the α₂ (long-scale) h(q) curve — null where insufficient
			hCurveSlope2: shape2.slope,
			hCurveCurvature2: shape2.curvature,
			hCurveNonlinearity2: shape2.nonlinearity,
			hMinLocation2: shape2.hMinLocation,
			monofractal, // h(2) per range; matches compute() alpha/alpha1/alpha2
			ranges: {
				global: [scales[0], scales[scales.length - 1]],
				alpha1: alpha1Range,
				alpha2: alpha2Range,
			},
			scales,
			fluctuationsByQ,
			lengthOfData: N,
		};
	};

	// Deprecated/unused methods kept for backward compatibility
	DFA.prototype.generateRange = function (array, startAt = 0, step = 1) {
		let size = array.length;
		let result = [];
		for (let i = 0; i < size; i++) {
			result.push(i * step + startAt);
		}
		return result;
	};

	DFA.prototype.generateZeroArray = function (size) {
		let result = [];
		for (let i = 0; i < size; i++) {
			result.push(0);
		}
		return result;
	};

	DFA.prototype.generateScaleRange = function (
		array,
		startAt = 4,
		step = 0.25
	) {
		// Legacy method - deprecated but kept for compatibility
		let startPow = Math.sqrt(startAt);
		let size = Math.floor((Math.log2(array.length) - startPow) / step + 1);
		let result = [];
		for (let i = 0; i < size; i++) {
			result.push(Math.floor(Math.pow(2, i * step + startPow)));
		}
		return result;
	};

	DFA.prototype.splitArrayIntoChunks = function (arr, size, strict = false) {
		// Legacy method - deprecated but kept for compatibility
		let chunkedArray = [];
		for (let i = 0; i < arr.length; i += size) {
			chunkedArray.push(arr.slice(i, i + size));
		}
		if (strict) {
			if (chunkedArray[chunkedArray.length - 1].length < size) {
				chunkedArray.pop();
			}
		}
		return chunkedArray;
	};

	DFA.prototype.log2Vector = function (x) {
		// Legacy method - deprecated but kept for compatibility
		if (!(x instanceof Array)) {
			throw new Error("log2Vector input must be an array");
		}
		let y = [];
		for (let i in x) {
			y[i] = Math.log2(x[i]);
		}
		return y;
	};

	DFA.prototype.RMS = function (cut, fit) {
		// Legacy method - deprecated but kept for compatibility
		if (cut.length != fit.length) {
			throw new Error(
				"for calculating the RMS of the two vectors they should be equal"
			);
		}
		let variance_vector = [];
		for (let i = 0; i < cut.length; i++) {
			variance_vector[i] = Math.pow(cut[i] - fit[i], 2);
		}
		return Math.sqrt(mean(variance_vector));
	};

	return DFA;
})();

module.exports = DFA;