From e5deaad772c2ed0a38f62284254333e116eec8e4 Mon Sep 17 00:00:00 2001 From: Nick LaSorte <37697540+nicklasorte@users.noreply.github.com> Date: Thu, 26 Mar 2026 12:13:50 -0400 Subject: [PATCH] Add two-level MC chunk planning and grouped-save rev9 wrapper --- dynamic_mc_chunks_rev2.m | 132 +++++++++++++++++++ parfor_randchunk_aggcheck_rev9_claude.m | 92 ++++++++++++++ validate_grouped_chunk_strategy_rev1.m | 162 ++++++++++++++++++++++++ 3 files changed, 386 insertions(+) create mode 100644 dynamic_mc_chunks_rev2.m create mode 100644 parfor_randchunk_aggcheck_rev9_claude.m create mode 100644 validate_grouped_chunk_strategy_rev1.m diff --git a/dynamic_mc_chunks_rev2.m b/dynamic_mc_chunks_rev2.m new file mode 100644 index 0000000..1182110 --- /dev/null +++ b/dynamic_mc_chunks_rev2.m @@ -0,0 +1,132 @@ +function [chunk_plan]=dynamic_mc_chunks_rev2(app,mc_size,num_bs,num_sim_azi,worker_memory_mb,target_memory_utilization,max_saved_chunks) +%DYNAMIC_MC_CHUNKS_REV2 Two-level Monte Carlo chunk planner. +% This revision separates compute chunking (memory safety) from save +% chunking (file-count control). +% +% chunk_plan fields (primary): +% - num_save_chunks +% - save_chunk_idx_ranges [num_save_chunks x 2] +% - save_chunk_to_compute_subchunks {num_save_chunks x 1} +% - save_chunk_rand_order randomized save-chunk order +% - compute_chunk_size +% - save_chunk_size_mc +% - cell_compute_chunk_idx per-compute-chunk MC indices +% +% Inputs: +% mc_size total Monte Carlo iterations +% num_bs number of BS rows in hot path +% num_sim_azi number of simulated azimuths +% worker_memory_mb memory per worker (default 2048) +% target_memory_utilization fraction of worker memory for chunk working +% set (default 0.55) +% max_saved_chunks hard cap on grouped save files (default 128) + +if nargin < 4 + error('dynamic_mc_chunks_rev2:NotEnoughInputs', ... + 'Need at least app, mc_size, num_bs, num_sim_azi.'); +end + +if nargin < 5 || isempty(worker_memory_mb) + worker_memory_mb = 2048; +end +if nargin < 6 || isempty(target_memory_utilization) + target_memory_utilization = 0.55; +end +if nargin < 7 || isempty(max_saved_chunks) + max_saved_chunks = 128; +end + +worker_memory_mb = max(1, worker_memory_mb); +target_memory_utilization = min(max(target_memory_utilization,0.10),0.90); +max_saved_chunks = max(1, floor(max_saved_chunks)); + +if mc_size <= 0 + error('dynamic_mc_chunks_rev2:InvalidMcSize','mc_size must be positive.'); +end +if num_bs <= 0 || num_sim_azi <= 0 + error('dynamic_mc_chunks_rev2:InvalidDimensions','num_bs and num_sim_azi must be positive.'); +end + +% Conservative memory model: +% - per-MC raw arrays include [num_bs x 1] PR/EIRP/clutter and one +% [num_bs x num_sim_azi] style intermediate footprint. +% - multiply by a safety factor for MATLAB temporaries / conversion buffers. +bytes_per_double = 8; +raw_bytes_per_mc = double(num_bs) * double(num_sim_azi + 3) * bytes_per_double; +memory_safety_multiplier = 4.0; % conservative practical mid-point (3x to 6x) +estimated_bytes_per_mc = raw_bytes_per_mc * memory_safety_multiplier; + +worker_budget_bytes = double(worker_memory_mb) * 1024 * 1024 * target_memory_utilization; +compute_chunk_size = floor(worker_budget_bytes / max(1, estimated_bytes_per_mc)); +compute_chunk_size = max(1, min(mc_size, compute_chunk_size)); + +num_compute_chunks = ceil(mc_size / compute_chunk_size); +cell_compute_chunk_idx = cell(num_compute_chunks,1); +compute_chunk_idx_ranges = zeros(num_compute_chunks,2); + +for compute_subchunk_idx = 1:num_compute_chunks + start_idx = (compute_subchunk_idx-1) * compute_chunk_size + 1; + stop_idx = min(compute_subchunk_idx * compute_chunk_size, mc_size); + idx = start_idx:stop_idx; + cell_compute_chunk_idx{compute_subchunk_idx} = idx; + compute_chunk_idx_ranges(compute_subchunk_idx,:) = [start_idx, stop_idx]; +end + +preferred_saved_chunks = min(max_saved_chunks, 128); +num_save_chunks = min(num_compute_chunks, preferred_saved_chunks); +compute_chunks_per_save = ceil(num_compute_chunks / num_save_chunks); + +save_chunk_to_compute_subchunks = cell(num_save_chunks,1); +save_chunk_idx_ranges = zeros(num_save_chunks,2); + +for save_chunk_idx = 1:num_save_chunks + comp_start = (save_chunk_idx-1) * compute_chunks_per_save + 1; + comp_stop = min(save_chunk_idx * compute_chunks_per_save, num_compute_chunks); + subchunks = comp_start:comp_stop; + save_chunk_to_compute_subchunks{save_chunk_idx} = subchunks; + + mc_start = compute_chunk_idx_ranges(comp_start,1); + mc_stop = compute_chunk_idx_ranges(comp_stop,2); + save_chunk_idx_ranges(save_chunk_idx,:) = [mc_start, mc_stop]; +end + +save_chunk_size_mc = ceil(mc_size / num_save_chunks); + +[tf_ml_toolbox]=check_ml_toolbox(app); +if tf_ml_toolbox==1 + save_chunk_rand_order = randsample(num_save_chunks,num_save_chunks,false); +else + save_chunk_rand_order = randperm(num_save_chunks); +end + +all_compute_idx = horzcat(cell_compute_chunk_idx{:}); +if length(all_compute_idx) ~= mc_size || ~isequal(sort(all_compute_idx),1:mc_size) + error('dynamic_mc_chunks_rev2:CoverageError', ... + 'Compute chunks do not cover MC indices exactly once.'); +end + +chunk_plan = struct(); +chunk_plan.mc_size = mc_size; +chunk_plan.num_bs = num_bs; +chunk_plan.num_sim_azi = num_sim_azi; +chunk_plan.worker_memory_mb = worker_memory_mb; +chunk_plan.target_memory_utilization = target_memory_utilization; +chunk_plan.worker_budget_mb = worker_budget_bytes/(1024*1024); +chunk_plan.max_saved_chunks = max_saved_chunks; +chunk_plan.memory_safety_multiplier = memory_safety_multiplier; +chunk_plan.estimated_bytes_per_mc = estimated_bytes_per_mc; +chunk_plan.estimated_mb_per_mc = estimated_bytes_per_mc/(1024*1024); +chunk_plan.compute_chunk_size = compute_chunk_size; +chunk_plan.num_compute_chunks = num_compute_chunks; +chunk_plan.compute_chunk_idx_ranges = compute_chunk_idx_ranges; +chunk_plan.cell_compute_chunk_idx = cell_compute_chunk_idx; +chunk_plan.num_save_chunks = num_save_chunks; +chunk_plan.save_chunk_size_mc = save_chunk_size_mc; +chunk_plan.compute_chunks_per_save = compute_chunks_per_save; +chunk_plan.save_chunk_idx_ranges = save_chunk_idx_ranges; +chunk_plan.save_chunk_to_compute_subchunks = save_chunk_to_compute_subchunks; +chunk_plan.save_chunk_rand_order = save_chunk_rand_order; + +chunk_plan.saved_chunk_cap_respected = (num_save_chunks <= max_saved_chunks); +chunk_plan.saved_chunk_preference_respected = (num_save_chunks <= 128); +end diff --git a/parfor_randchunk_aggcheck_rev9_claude.m b/parfor_randchunk_aggcheck_rev9_claude.m new file mode 100644 index 0000000..9ccab79 --- /dev/null +++ b/parfor_randchunk_aggcheck_rev9_claude.m @@ -0,0 +1,92 @@ +function [save_chunk_agg_check_mc_dBm]=parfor_randchunk_aggcheck_rev9_claude(app,agg_check_file_name,agg_dist_file_name,chunk_plan,save_chunk_loop_idx,point_idx,sim_number,data_label1,cell_aas_dist_data,array_bs_azi_data,radar_beamwidth,min_azimuth,max_azimuth,base_protection_pts,on_list_bs,rand_seed1,agg_check_reliability,on_full_Pr_dBm,clutter_loss,custom_antenna_pattern,parallel_flag,single_search_dist) +%PARFOR_RANDCHUNK_AGGCHECK_REV9_CLAUDE Grouped-save Monte Carlo wrapper. +% Rev9 processes one save chunk at a time, where each save chunk contains +% multiple memory-safe compute subchunks. It writes one MAT file per save +% chunk (instead of one file per compute subchunk). + +save_chunk_agg_check_mc_dBm = NaN(1,1); + +[var_exist1]=persistent_var_exist_with_corruption(app,agg_check_file_name); +[var_exist2]=persistent_var_exist_with_corruption(app,agg_dist_file_name); +if var_exist1==2 && var_exist2==2 + % Final files already exist; no further work needed. + return +end + +if ~isstruct(chunk_plan) || ~isfield(chunk_plan,'save_chunk_rand_order') + error('parfor_randchunk_aggcheck_rev9_claude:InvalidPlan', ... + 'chunk_plan must be a struct created by dynamic_mc_chunks_rev2.'); +end + +save_chunk_idx = chunk_plan.save_chunk_rand_order(save_chunk_loop_idx); +mc_range = chunk_plan.save_chunk_idx_ranges(save_chunk_idx,:); + +file_name_agg_check_chunk = strcat( ... + 'subg_',num2str(save_chunk_idx), ... + '_array_agg_check_mc_dBm_',num2str(point_idx),'_',num2str(sim_number),'_', ... + data_label1,'_',num2str(single_search_dist),'km_', ... + 'mc',num2str(mc_range(1)),'_',num2str(mc_range(2)),'.mat'); + +[var_exist_group]=persistent_var_exist_with_corruption(app,file_name_agg_check_chunk); +if var_exist_group==2 && parallel_flag==0 + retry_load=1; + while(retry_load==1) + try + load(file_name_agg_check_chunk,'save_chunk_agg_check_mc_dBm'); + temp_data=save_chunk_agg_check_mc_dBm; + clear save_chunk_agg_check_mc_dBm; + save_chunk_agg_check_mc_dBm=temp_data; + clear temp_data; + retry_load=0; + catch + retry_load=1; + pause(1); + end + end + return +elseif var_exist_group==2 && parallel_flag==1 + % In parallel mode, avoid loading if already present. + return +end + +compute_subchunk_idx_list = chunk_plan.save_chunk_to_compute_subchunks{save_chunk_idx}; +num_compute_subchunks = length(compute_subchunk_idx_list); +sub_results = cell(num_compute_subchunks,1); + +for local_idx = 1:num_compute_subchunks + compute_subchunk_idx = compute_subchunk_idx_list(local_idx); + sub_results{local_idx} = subchunk_agg_check_rev8( ... + app,cell_aas_dist_data,array_bs_azi_data,radar_beamwidth,min_azimuth, ... + max_azimuth,base_protection_pts,point_idx,on_list_bs, ... + chunk_plan.cell_compute_chunk_idx,rand_seed1,agg_check_reliability, ... + on_full_Pr_dBm,clutter_loss,custom_antenna_pattern,compute_subchunk_idx); +end + +save_chunk_agg_check_mc_dBm = vertcat(sub_results{:}); +expected_rows = mc_range(2) - mc_range(1) + 1; +if size(save_chunk_agg_check_mc_dBm,1) ~= expected_rows + error('parfor_randchunk_aggcheck_rev9_claude:SizeMismatch', ... + 'Grouped save chunk rows (%d) != expected MC rows (%d).', ... + size(save_chunk_agg_check_mc_dBm,1), expected_rows); +end + +[var_exist_after_compute]=persistent_var_exist_with_corruption(app,file_name_agg_check_chunk); +if var_exist_after_compute==0 + temp_file_name = strcat(file_name_agg_check_chunk,'.tmp_',num2str(feature('getpid')),'_',num2str(randi(1e9)),'.mat'); + retry_save=1; + while(retry_save==1) + try + save(temp_file_name,'save_chunk_agg_check_mc_dBm','save_chunk_idx','mc_range','compute_subchunk_idx_list','-v7.3'); + movefile(temp_file_name,file_name_agg_check_chunk,'f'); + retry_save=0; + catch + retry_save=1; + if exist(temp_file_name,'file')==2 + delete(temp_file_name); + end + pause(1); + end + end +end + +end diff --git a/validate_grouped_chunk_strategy_rev1.m b/validate_grouped_chunk_strategy_rev1.m new file mode 100644 index 0000000..9cfc01f --- /dev/null +++ b/validate_grouped_chunk_strategy_rev1.m @@ -0,0 +1,162 @@ +function report = validate_grouped_chunk_strategy_rev1(app,chunk_plan,varargin) +%VALIDATE_GROUPED_CHUNK_STRATEGY_REV1 Validation checks for rev2 + rev9 design. +% This helper validates chunk-plan integrity and can optionally run a +% numerical A/B comparison between direct per-compute assembly and grouped +% save-chunk assembly using subchunk_agg_check_rev8. +% +% Required: +% app, chunk_plan +% +% Optional name/value: +% 'RunNumericalCheck' (default false) +% 'NumericalInputs' struct with fields needed by subchunk_agg_check_rev8 +% 'UseRandomSaveOrder' (default true) + +p = inputParser; +p.addParameter('RunNumericalCheck',false,@(x)islogical(x)&&isscalar(x)); +p.addParameter('NumericalInputs',struct(),@isstruct); +p.addParameter('UseRandomSaveOrder',true,@(x)islogical(x)&&isscalar(x)); +p.parse(varargin{:}); +opts = p.Results; + +report = struct(); + +%% 1) Planner coverage, no gaps, no duplicates +all_idx = horzcat(chunk_plan.cell_compute_chunk_idx{:}); +report.coverage_exact_once = (numel(all_idx)==chunk_plan.mc_size) && isequal(sort(all_idx),1:chunk_plan.mc_size); +report.no_gaps = isequal(unique(all_idx),1:chunk_plan.mc_size); +report.no_duplicates = (numel(unique(all_idx))==chunk_plan.mc_size); + +%% 2) Save chunks nested correctly +nested_ok = true; +save_idx_union = []; +for save_chunk_idx = 1:chunk_plan.num_save_chunks + subchunks = chunk_plan.save_chunk_to_compute_subchunks{save_chunk_idx}; + from_subchunks = []; + for k = 1:length(subchunks) + from_subchunks = [from_subchunks, chunk_plan.cell_compute_chunk_idx{subchunks(k)}]; %#ok + end + range_idx = chunk_plan.save_chunk_idx_ranges(save_chunk_idx,1):chunk_plan.save_chunk_idx_ranges(save_chunk_idx,2); + if ~isequal(from_subchunks,range_idx) + nested_ok = false; + break + end + save_idx_union = [save_idx_union, range_idx]; %#ok +end +report.save_nested_correct = nested_ok; +report.save_union_exact_once = (numel(save_idx_union)==chunk_plan.mc_size) && isequal(sort(save_idx_union),1:chunk_plan.mc_size); + +%% 3) Basic operational summary values +report.num_save_chunks = chunk_plan.num_save_chunks; +report.compute_chunk_size = chunk_plan.compute_chunk_size; +report.worker_budget_mb = chunk_plan.worker_budget_mb; +report.saved_chunk_cap_respected = chunk_plan.saved_chunk_cap_respected; +report.saved_chunk_preference_respected = chunk_plan.saved_chunk_preference_respected; + +%% 4) Optional numerical and randomized-order invariance checks +report.size_equal = NaN; +report.max_abs_diff = NaN; +report.nan_pattern_equal = NaN; +report.randomized_order_invariant = NaN; +report.restart_skip_detected = NaN; + +if opts.RunNumericalCheck + req = {'cell_aas_dist_data','array_bs_azi_data','radar_beamwidth','min_azimuth', ... + 'max_azimuth','base_protection_pts','point_idx','on_list_bs','rand_seed1', ... + 'agg_check_reliability','on_full_Pr_dBm','clutter_loss','custom_antenna_pattern'}; + for i = 1:length(req) + if ~isfield(opts.NumericalInputs,req{i}) + error('validate_grouped_chunk_strategy_rev1:MissingInput', ... + 'Missing NumericalInputs.%s', req{i}); + end + end + + ni = opts.NumericalInputs; + + % Baseline assembly: compute chunks in numeric order + baseline_cells = cell(chunk_plan.num_compute_chunks,1); + for compute_subchunk_idx = 1:chunk_plan.num_compute_chunks + baseline_cells{compute_subchunk_idx} = subchunk_agg_check_rev8( ... + app,ni.cell_aas_dist_data,ni.array_bs_azi_data,ni.radar_beamwidth, ... + ni.min_azimuth,ni.max_azimuth,ni.base_protection_pts,ni.point_idx, ... + ni.on_list_bs,chunk_plan.cell_compute_chunk_idx,ni.rand_seed1, ... + ni.agg_check_reliability,ni.on_full_Pr_dBm,ni.clutter_loss, ... + ni.custom_antenna_pattern,compute_subchunk_idx); + end + baseline = vertcat(baseline_cells{:}); + + % Grouped assembly: save chunk order can be randomized and should match. + grouped = NaN(size(baseline)); + if opts.UseRandomSaveOrder + save_order = chunk_plan.save_chunk_rand_order; + else + save_order = 1:chunk_plan.num_save_chunks; + end + + for t = 1:length(save_order) + save_chunk_idx = save_order(t); + compute_subchunks = chunk_plan.save_chunk_to_compute_subchunks{save_chunk_idx}; + local_cells = cell(length(compute_subchunks),1); + for k = 1:length(compute_subchunks) + cidx = compute_subchunks(k); + local_cells{k} = subchunk_agg_check_rev8( ... + app,ni.cell_aas_dist_data,ni.array_bs_azi_data,ni.radar_beamwidth, ... + ni.min_azimuth,ni.max_azimuth,ni.base_protection_pts,ni.point_idx, ... + ni.on_list_bs,chunk_plan.cell_compute_chunk_idx,ni.rand_seed1, ... + ni.agg_check_reliability,ni.on_full_Pr_dBm,ni.clutter_loss, ... + ni.custom_antenna_pattern,cidx); + end + grouped_chunk = vertcat(local_cells{:}); + rng_idx = chunk_plan.save_chunk_idx_ranges(save_chunk_idx,1):chunk_plan.save_chunk_idx_ranges(save_chunk_idx,2); + grouped(rng_idx,:) = grouped_chunk; + end + + report.size_equal = isequal(size(baseline),size(grouped)); + report.nan_pattern_equal = isequal(isnan(baseline),isnan(grouped)); + + valid_mask = ~isnan(baseline) & ~isnan(grouped); + if any(valid_mask,'all') + report.max_abs_diff = max(abs(baseline(valid_mask)-grouped(valid_mask)),[],'all'); + else + report.max_abs_diff = NaN; + end + + % Compare randomized order vs deterministic order reconstruction. + grouped_det = NaN(size(baseline)); + for save_chunk_idx = 1:chunk_plan.num_save_chunks + compute_subchunks = chunk_plan.save_chunk_to_compute_subchunks{save_chunk_idx}; + local_cells = cell(length(compute_subchunks),1); + for k = 1:length(compute_subchunks) + cidx = compute_subchunks(k); + local_cells{k} = subchunk_agg_check_rev8( ... + app,ni.cell_aas_dist_data,ni.array_bs_azi_data,ni.radar_beamwidth, ... + ni.min_azimuth,ni.max_azimuth,ni.base_protection_pts,ni.point_idx, ... + ni.on_list_bs,chunk_plan.cell_compute_chunk_idx,ni.rand_seed1, ... + ni.agg_check_reliability,ni.on_full_Pr_dBm,ni.clutter_loss, ... + ni.custom_antenna_pattern,cidx); + end + grouped_det(chunk_plan.save_chunk_idx_ranges(save_chunk_idx,1):chunk_plan.save_chunk_idx_ranges(save_chunk_idx,2),:) = vertcat(local_cells{:}); + end + report.randomized_order_invariant = isequaln(grouped,grouped_det); + + % Restart skip behavior check (exists -> skip). + tmp_dir = tempname; + mkdir(tmp_dir); + orig_dir = pwd; + cleanup_obj = onCleanup(@() cd(orig_dir)); %#ok + cd(tmp_dir); + + agg_check_file_name = "final_dummy_agg_check.mat"; + agg_dist_file_name = "final_dummy_agg_dist.mat"; + save(agg_check_file_name,'grouped'); + save(agg_dist_file_name,'grouped_det'); + + out = parfor_randchunk_aggcheck_rev9_claude(app,agg_check_file_name,agg_dist_file_name,chunk_plan,1, ... + ni.point_idx,1,'validation',ni.cell_aas_dist_data,ni.array_bs_azi_data, ... + ni.radar_beamwidth,ni.min_azimuth,ni.max_azimuth,ni.base_protection_pts, ... + ni.on_list_bs,ni.rand_seed1,ni.agg_check_reliability,ni.on_full_Pr_dBm, ... + ni.clutter_loss,ni.custom_antenna_pattern,0,1); + + report.restart_skip_detected = isequaln(out,NaN(1,1)); +end +end