From e47162ee31621fa44830025e9e69dfca5ce59048 Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 27 Mar 2026 18:47:01 +0000 Subject: [PATCH] Add rev7 EIRP helper with corrected reshape and helper-level validation Rev6 had a reshape axis-swap bug: unmkpp stores coefs as [dim*pieces x order] with dim varying fastest, so the correct reshape is [dim, pieces, order]. Rev6 used [pieces, dim, order], scrambling BS/piece coefficient pairing and producing the observed semantic inversion (high query -> low EIRP, vice versa). Rev7 fixes only the reshape and downstream indexing, preserving the ~35x speedup from vectorized spline evaluation. Includes a helper-level validation script that checks absolute/relative error, correlation, and endpoints against rev5 before end-to-end integration. https://claude.ai/code/session_0189a75rK83NM29Ryy84NZ7X --- monte_carlo_super_bs_eirp_dist_rev7.m | 91 +++++++++++++ ...monte_carlo_super_bs_eirp_dist_rev5_rev7.m | 126 ++++++++++++++++++ 2 files changed, 217 insertions(+) create mode 100644 monte_carlo_super_bs_eirp_dist_rev7.m create mode 100644 validate_monte_carlo_super_bs_eirp_dist_rev5_rev7.m diff --git a/monte_carlo_super_bs_eirp_dist_rev7.m b/monte_carlo_super_bs_eirp_dist_rev7.m new file mode 100644 index 0000000..a63d4cd --- /dev/null +++ b/monte_carlo_super_bs_eirp_dist_rev7.m @@ -0,0 +1,91 @@ +function [rand_norm_eirp]=monte_carlo_super_bs_eirp_dist_rev7(app,super_array_bs_eirp_dist,reliability,rand_numbers) +%MONTE_CARLO_SUPER_BS_EIRP_DIST_REV7 Corrected fast RNG-free MC EIRP interpolation. +% Fix over rev6: correct reshape dimension order for unmkpp coefficient layout. +% unmkpp stores coefs as [dim*pieces x order] with dim consecutive rows per piece +% (dim varies fastest). Correct reshape is [dim, pieces, order], NOT [pieces, dim, order]. +% +% Preserves rev5 semantics exactly: per-row spline interpolation of EIRP vs reliability. +% Strategy: +% 1) Build one spline piecewise polynomial object for all BS rows at once. +% 2) Evaluate each BS at its own random reliability via direct pp coefficients. + +% app is intentionally unused (signature compatibility). + +[num_rows,num_cols]=size(super_array_bs_eirp_dist); + +if num_cols<=1 + rand_norm_eirp=zeros(num_rows,1); + return; +end + +rel_col=reliability(:); +if ~issorted(rel_col) + [rel_col,sort_idx]=sort(rel_col,'ascend'); + super_array_bs_eirp_dist=super_array_bs_eirp_dist(:,sort_idx); +end + +rel_min=rel_col(1); +rel_max=rel_col(end); +xi=min(max(rand_numbers(:),rel_min),rel_max); + +% Build spline PP for all rows in one call. +% spline(x,Y) with size(Y,2)==numel(x) treats each row of Y as a separate function. +pp=spline(rel_col,super_array_bs_eirp_dist); +[breaks,coefs,pieces,order,dim]=unmkpp(pp); + +if order~=4 + error('monte_carlo_super_bs_eirp_dist_rev7:UnexpectedPPOrder', ... + 'Expected cubic spline order 4, got order %d.',order); +end +if dim~=num_rows + error('monte_carlo_super_bs_eirp_dist_rev7:UnexpectedPPDim', ... + 'Expected PP dim %d, got %d.',num_rows,dim); +end + +% coefs layout from unmkpp: [dim*pieces x order] with dim consecutive rows per piece. +% Row (p-1)*dim + d = piece p, dimension (BS) d. +% Correct reshape: [dim, pieces, order] so coefs3(d, p, k) = coef for BS d, piece p, power k. +coefs3=reshape(coefs,[dim,pieces,order]); +a_all=coefs3(:,:,1); % [dim x pieces] = [num_rows x pieces] +b_all=coefs3(:,:,2); +c_all=coefs3(:,:,3); +d_all=coefs3(:,:,4); + +% Locate xi interval index (1..pieces), matching ppval boundary handling. +num_samples=numel(xi); +if num_samples~=num_rows + error('monte_carlo_super_bs_eirp_dist_rev7:SizeMismatch', ... + 'Expected rand_numbers length %d, got %d.',num_rows,num_samples); +end + +% Interval index via cumulative comparison. +idx=ones(num_rows,1); +for k=2:numel(breaks) + idx=idx + (xi>=breaks(k)); +end +idx=min(idx,pieces); + +base_break=breaks(idx); +dx=xi-base_break(:); + +% Linear index into [dim x pieces] arrays: row d, column p → (p-1)*dim + d. +row_idx=(1:num_rows).'; +lin_idx=row_idx + (idx-1)*dim; + +a=a_all(lin_idx); +b=b_all(lin_idx); +c=c_all(lin_idx); +d_coef=d_all(lin_idx); + +% Force column vectors. +a=a(:); +b=b(:); +c=c(:); +d_coef=d_coef(:); +dx=dx(:); + +% Horner evaluation of cubic: ((a*dx + b)*dx + c)*dx + d +rand_norm_eirp=((a.*dx+b).*dx+c).*dx+d_coef; +rand_norm_eirp=rand_norm_eirp(:); + +end diff --git a/validate_monte_carlo_super_bs_eirp_dist_rev5_rev7.m b/validate_monte_carlo_super_bs_eirp_dist_rev5_rev7.m new file mode 100644 index 0000000..4ae4a18 --- /dev/null +++ b/validate_monte_carlo_super_bs_eirp_dist_rev5_rev7.m @@ -0,0 +1,126 @@ +function results=validate_monte_carlo_super_bs_eirp_dist_rev5_rev7(app,super_array_bs_eirp_dist,reliability,rand_numbers) +%VALIDATE_MONTE_CARLO_SUPER_BS_EIRP_DIST_REV5_REV7 +% Helper-level validation: compare rev5 (golden) vs rev7 (fixed fast path). +% Must pass before any end-to-end integration. +% +% Checks: +% 1) Exact shape match of outputs +% 2) Max absolute error across all elements +% 3) Per-element relative error statistics +% 4) Known failure pattern detection (inversion check) +% 5) Endpoint / boundary behavior +% 6) Monotonicity preservation spot-check + +fprintf('\n=== HELPER VALIDATION: rev5 vs rev7 ===\n'); + +% --- Run both --- +out5=monte_carlo_super_bs_eirp_dist_rev5(app,super_array_bs_eirp_dist,reliability,rand_numbers); +out7=monte_carlo_super_bs_eirp_dist_rev7(app,super_array_bs_eirp_dist,reliability,rand_numbers); + +% --- Shape check --- +sz5=size(out5); +sz7=size(out7); +shape_ok=isequal(sz5,sz7); +fprintf('Shape rev5: [%s] rev7: [%s] match: %s\n', ... + num2str(sz5),num2str(sz7),yesno(shape_ok)); +if ~shape_ok + error('validate_rev5_rev7:ShapeMismatch','Output shapes differ.'); +end + +% --- Absolute error --- +abs_err=abs(out7-out5); +max_abs_err=max(abs_err); +mean_abs_err=mean(abs_err); +fprintf('Max absolute error: %.6e\n',max_abs_err); +fprintf('Mean absolute error: %.6e\n',mean_abs_err); + +% --- Relative error (relative to rev5 magnitude, guarded) --- +denom=max(abs(out5),1e-12); +rel_err=abs_err./denom; +max_rel_err=max(rel_err); +mean_rel_err=mean(rel_err); +fprintf('Max relative error: %.6e\n',max_rel_err); +fprintf('Mean relative error: %.6e\n',mean_rel_err); + +% --- Inversion detection (the rev6 failure signature) --- +% If sign(out7 - mean(out7)) is anti-correlated with sign(out5 - mean(out5)), +% we have the same axis-swap bug. +corr_val=corr(out5(:),out7(:)); +fprintf('Pearson correlation: %.8f\n',corr_val); +inversion_detected=corr_val<0.5; +if inversion_detected + fprintf('*** INVERSION DETECTED: correlation %.4f < 0.5 ***\n',corr_val); +end + +% --- Worst-case examples (for manual inspection) --- +[~,worst_idx]=sort(abs_err,'descend'); +n_show=min(5,numel(worst_idx)); +fprintf('\nWorst-case elements:\n'); +fprintf(' %-6s %-12s %-12s %-12s %-10s\n','idx','rev5','rev7','abs_err','query'); +for k=1:n_show + ii=worst_idx(k); + fprintf(' %-6d %12.6f %12.6f %12.6e %10.6f\n', ... + ii,out5(ii),out7(ii),abs_err(ii),rand_numbers(ii)); +end + +% --- Endpoint spot-check: query at rel_min and rel_max --- +fprintf('\nEndpoint spot-check:\n'); +rel_sorted=sort(reliability(:)); +rel_min=rel_sorted(1); +rel_max=rel_sorted(end); +% Test with first row's data at boundaries +rn_lo=rel_min*ones(size(rand_numbers)); +rn_hi=rel_max*ones(size(rand_numbers)); +out5_lo=monte_carlo_super_bs_eirp_dist_rev5(app,super_array_bs_eirp_dist,reliability,rn_lo); +out7_lo=monte_carlo_super_bs_eirp_dist_rev7(app,super_array_bs_eirp_dist,reliability,rn_lo); +out5_hi=monte_carlo_super_bs_eirp_dist_rev5(app,super_array_bs_eirp_dist,reliability,rn_hi); +out7_hi=monte_carlo_super_bs_eirp_dist_rev7(app,super_array_bs_eirp_dist,reliability,rn_hi); +lo_err=max(abs(out7_lo-out5_lo)); +hi_err=max(abs(out7_hi-out5_hi)); +fprintf(' At rel_min (%.6f): max abs error = %.6e\n',rel_min,lo_err); +fprintf(' At rel_max (%.6f): max abs error = %.6e\n',rel_max,hi_err); + +% --- Pass/fail thresholds --- +% Spline should be numerically identical (same algorithm, same coefficients). +% Allow for floating-point rounding only. +ABS_TOL=1e-10; +REL_TOL=1e-10; +pass_abs=max_abs_err<=ABS_TOL; +pass_rel=max_rel_err<=REL_TOL; +pass_corr=corr_val>0.999; +pass_endpoints=(lo_err<=ABS_TOL) && (hi_err<=ABS_TOL); +overall_pass=pass_abs && pass_rel && pass_corr && pass_endpoints; + +fprintf('\nPass/fail summary:\n'); +fprintf(' Absolute error <= %.1e: %s\n',ABS_TOL,passfail(pass_abs)); +fprintf(' Relative error <= %.1e: %s\n',REL_TOL,passfail(pass_rel)); +fprintf(' Correlation > 0.999: %s\n',passfail(pass_corr)); +fprintf(' Endpoints <= %.1e: %s\n',ABS_TOL,passfail(pass_endpoints)); +fprintf(' OVERALL: %s\n',passfail(overall_pass)); + +if ~overall_pass + error('validate_rev5_rev7:Failed', ... + 'Helper validation FAILED. Do NOT integrate rev7 into end-to-end.'); +end + +% --- Build results struct --- +results=struct(); +results.max_abs_err=max_abs_err; +results.mean_abs_err=mean_abs_err; +results.max_rel_err=max_rel_err; +results.mean_rel_err=mean_rel_err; +results.correlation=corr_val; +results.endpoint_lo_err=lo_err; +results.endpoint_hi_err=hi_err; +results.pass=overall_pass; +results.n_elements=numel(out5); + +end + +function s=yesno(tf) +if tf; s='YES'; else; s='NO'; end +end + +function s=passfail(tf) +if tf; s='PASS'; else; s='FAIL'; end +end