"""Instrument isolated pipe solvers; these runs are diagnostics, never benchmarks. Prepare: python tests/manual/profile_pipe_iterations.py --output-dir test/pipe-profile --prepare-only Run the prepared diagnostic programs: use the same command without --prepare-only. All mutations except this test helper stay below the ignored output directory. The replay contains EVERY resistance-law input from the guarded solver's RHS trajectory, including its scalar low-Re analytic branch. Cache hits, zero dp, and the separate PNL00R analytic law are counted but do not enter this replay. """ from __future__ import annotations import argparse from hashlib import sha256 import json import os from pathlib import Path import subprocess import sys ROOT = Path(__file__).resolve().parents[2] sys.path.insert(0, str(ROOT)) from app.main import compile_system_xml_network from app.simulation.backends import simulation_config from app.simulation.native_codegen import build as builder from app.simulation.native_codegen.compiler import compile_native_program from app.simulation.native_codegen.input import load_input from tests.manual.benchmark_native_pipe_solver import build_snapshot, pipe_source_path, working_native_paths SNAPSHOT_HELPER = ROOT / 'tests/manual/benchmark_native_pipe_solver.py' VARIANTS = ('guarded-newton', 'previous-newton', 'fixed-point') FIELDS = '''cache_requests cache_hits cache_misses flow_calls zero_pressure_calls pnl00r_analytic_calls resistance_calls scalar_analytic_calls iterative_calls iterations_total iterations_max reached_last_iteration exhausted_limit algorithm_converged failed_before_iteration finite_returns nonfinite_returns bisections residual_pass residual_fail residual_nonfinite residual_fail_after_algorithm_converged residual_pass_after_limit invalid_inputs upper_bracket_evaluations upper_bracket_exhausted bisection_after_poor_progress bisection_invalid_or_outside float_stagnation '''.split() PROFILE_PREFIX = r''' /* Test-only instrumentation, injected in an isolated source tree. */ #include #include #include #include #define PROFILE_VARIANT "@VARIANT@" #define PROFILE_FIXED @FIXED@ #define PROFILE_LIMIT(kind) @LIMIT@ #define PROFILE_FIELDS(X) @FIELDS@ typedef struct { #define PROFILE_DECLARE(name) unsigned long long name; PROFILE_FIELDS(PROFILE_DECLARE) #undef PROFILE_DECLARE unsigned long long histogram[129]; double maximum_relative_residual; } PipeProfile; static PipeProfile profile_stats[2][4]; static unsigned long long profile_rhs_calls,profile_capture_records; static int profile_in_rhs,profile_kind,profile_last_exhausted; static FILE *profile_capture; static PipeProfile *profile_bucket(int kind) { if(kind<0 || kind>3){fprintf(stderr,"Unexpected pipe kind %d\n",kind);exit(71);} return &profile_stats[profile_in_rhs?0:1][kind]; } void pipe_profile_rhs_enter(void){profile_in_rhs=1;profile_rhs_calls++;} void pipe_profile_rhs_leave(void){profile_in_rhs=0;} static void profile_add(PipeProfile *total,const PipeProfile *value) { #define PROFILE_ADD(name) total->name+=value->name; PROFILE_FIELDS(PROFILE_ADD) #undef PROFILE_ADD if(value->iterations_max>total->iterations_max)total->iterations_max=value->iterations_max; for(int i=0;i<129;i++)total->histogram[i]+=value->histogram[i]; if(value->maximum_relative_residual>total->maximum_relative_residual) total->maximum_relative_residual=value->maximum_relative_residual; } static void profile_write_bucket(FILE *f,const PipeProfile *value) { fprintf(f,"{"); #define PROFILE_WRITE(name) fprintf(f,"\"" #name "\":%llu,",value->name); PROFILE_FIELDS(PROFILE_WRITE) #undef PROFILE_WRITE fprintf(f,"\"maximum_relative_residual\":%.17g,\"iteration_histogram\":{",value->maximum_relative_residual); int comma=0; for(int i=0;i<129;i++)if(value->histogram[i]) { fprintf(f,"%s\"%d\":%llu",comma?",":"",i,value->histogram[i]);comma=1; } fprintf(f,"}}"); } static void profile_dump(void) { const char *path=getenv("PIPE_PROFILE_JSON"); if(profile_capture){if(fclose(profile_capture))exit(73);profile_capture=NULL;} if(!path)return; FILE *f=fopen(path,"wb");if(!f){perror(path);exit(73);} fprintf(f,"{\"variant\":\"%s\",\"rhs_calls\":%llu,\"capture_records\":%llu,",PROFILE_VARIANT,profile_rhs_calls,profile_capture_records); for(int scope=0;scope<2;scope++) { PipeProfile total={0}; for(int kind=0;kind<4;kind++)profile_add(&total,&profile_stats[scope][kind]); /* iterations_max is a maximum, unlike the additive counters. */ total.iterations_max=0; for(int kind=0;kind<4;kind++)if(profile_stats[scope][kind].iterations_max>total.iterations_max) total.iterations_max=profile_stats[scope][kind].iterations_max; fprintf(f,"%s\"%s\":{\"total\":",scope?",":"",scope?"non_rhs":"rhs"); profile_write_bucket(f,&total);fprintf(f,",\"by_kind\":{"); for(int kind=0;kind<4;kind++) { fprintf(f,"%s\"%d\":",kind?",":"",kind);profile_write_bucket(f,&profile_stats[scope][kind]); } fprintf(f,"}}"); } fprintf(f,"}\n");if(fclose(f))exit(73); } void pipe_profile_install(void) { const char *path=getenv("PIPE_PROFILE_CAPTURE"); if(path){profile_capture=fopen(path,"wb");if(!profile_capture){perror(path);exit(73);}} if(atexit(profile_dump)){fprintf(stderr,"Cannot register profile writer\n");exit(73);} } static void profile_save_input(double base,double d,double length,double rr,double den,int kind) { if(profile_in_rhs && profile_capture) { /* Six IEEE doubles, native endian; no sampling or deduplication. */ double input[]={base,d,length,rr,den,(double)kind}; if(fwrite(input,sizeof(input),1,profile_capture)!=1){perror("capture");exit(73);} profile_capture_records++; } } ''' PROFILE_SOLVE = r''' /* Independently factored Darcy law: no production slope or solver status is consulted. Long double reduces rounding noise in the returned-q residual. */ static long double profile_reference_friction(long double re,long double rr) { if(!(re>0))return NAN; long double laminar=64/re; if(re<=89.96829989L)return laminar; long double smooth=powl(-1.8L*log10l(6.9L/re),-2),turbulent=smooth; if(rr>0) { long double fully_rough=powl(-2*log10l(rr/3.7L),-2); long double weight=1/(1+powl(180/(re*rr),2)); turbulent=(1-weight)*smooth+weight*fully_rough; } long double blend=powl((re-89.96829989L)/2741.96700831L,8.37293695L); return (laminar+blend*turbulent)/(1+blend); } static void profile_returned_residual(PipeProfile *s,double q,double den,double K,double rr, int converged,int exhausted) { if(!isfinite(q)){s->nonfinite_returns++;s->residual_nonfinite++;return;} s->finite_returns++; long double re=fabsl((long double)q/(den/4)),residual; if(K==0 && q==0)residual=0; else residual=fabsl(re*re*profile_reference_friction(re,rr)/K-1); if(!isfinite(residual)){s->residual_nonfinite++;return;} if(residual>s->maximum_relative_residual)s->maximum_relative_residual=(double)residual; if(residual<=1e-9L){s->residual_pass++;if(exhausted)s->residual_pass_after_limit++;} else {s->residual_fail++;if(converged)s->residual_fail_after_algorithm_converged++;} } static double profile_solve(double base,double d,double length,double rr,double den,int kind) { PipeProfile *s=profile_bucket(kind); double K=pow(4*base/den,2)*d/length; profile_save_input(base,d,length,rr,den,kind); s->resistance_calls++; if(!(K>=0 && rr>=0 && den/4>0) || !isfinite(K) || !isfinite(rr) || !isfinite(den/4))s->invalid_inputs++; int iterations=0,converged=0,bisections=0,exhausted=0; double q; profile_kind=kind;profile_last_exhausted=0; #if PROFILE_FIXED double rough_limit=pipe_rough_limit(rr); q=sqrt(d/(length*.02))*base; for(int i=0;i128){fprintf(stderr,"Unexpected iteration count\n");exit(71);} s->histogram[iterations]++;s->iterations_total+=(unsigned)iterations;s->bisections+=(unsigned)bisections; if((unsigned)iterations>s->iterations_max)s->iterations_max=(unsigned)iterations; if(iterations)s->iterative_calls++; else if(converged)s->scalar_analytic_calls++; else s->failed_before_iteration++; if(iterations==PROFILE_LIMIT(kind))s->reached_last_iteration++; if(exhausted)s->exhausted_limit++; if(converged)s->algorithm_converged++; profile_returned_residual(s,q,den,K,rr,converged,exhausted); return q; } ''' REPLAY_MAIN = r''' #include "native/components/kernels.c" int main(int argc,char **argv) { if(argc!=2)return 64; pipe_profile_install();profile_in_rhs=1; FILE *f=fopen(argv[1],"rb");if(!f){perror(argv[1]);return 73;} double input[6];size_t count; while((count=fread(input,1,sizeof(input),f))==sizeof(input)) { profile_solve(input[0],input[1],input[2],input[3],input[4],(int)input[5]); } int failed=count || ferror(f);fclose(f);return failed?74:0; } ''' def replace_once(source: str, old: str, new: str) -> str: if source.count(old) != 1: raise ValueError(f'Expected exactly one audited source fragment: {old[:100]!r}') return source.replace(old, new, 1) def instrument(source: str, variant: str) -> str: fixed = variant == 'fixed-point' limit = '(kind==0?64:16)' if fixed else '80' if variant == 'previous-newton' else '128' prefix = PROFILE_PREFIX.replace('@VARIANT@', variant).replace('@FIXED@', str(int(fixed))) prefix = prefix.replace('@LIMIT@', limit).replace('@FIELDS@', ' '.join(f'X({name})' for name in FIELDS)) source = replace_once(source, '#include ', '#include \n' + prefix) # Exhaustion is marked at the actual loop fall-through, not inferred from # visiting the last allowed iteration (which can still converge). start = source.index('double native_pipe_resistance(') end = source.index('double native_pipe_flow(', start) resistance = source[start:end] ending = ' return NAN;\n}\n' if not resistance.endswith(ending): raise ValueError('Unexpected resistance function ending') resistance = resistance[:-len(ending)] + ' profile_last_exhausted=1;return NAN;\n}\n' if variant == 'guarded-newton': resistance = replace_once(resistance, ' double value=hi*hi*pipe_friction_prepared(hi,rr,rough);', ' profile_bucket(profile_kind)->upper_bracket_evaluations++;\n double value=hi*hi*pipe_friction_prepared(hi,rr,rough);') resistance = replace_once(resistance, ' if(!bracketed)return NAN;', ' if(!bracketed){profile_bucket(profile_kind)->upper_bracket_exhausted++;return NAN;}') resistance = replace_once(resistance, ' if(bisect) {', ''' if(bisect) { if(previous_newton && fabs(F)>.5*previous_residual)profile_bucket(profile_kind)->bisection_after_poor_progress++; if(!(slope>0) || !isfinite(slope) || !isfinite(next) || next<=lo || next>=hi) profile_bucket(profile_kind)->bisection_invalid_or_outside++;''') resistance = replace_once(resistance, ' if(next<=lo || next>=hi) {', ' if(next<=lo || next>=hi) {\n profile_bucket(profile_kind)->float_stagnation++;') elif variant == 'previous-newton': # Count each evaluation of the original loop condition without changing # its short-circuit behaviour or the original upper endpoint arithmetic. resistance = replace_once(resistance, 'i<128 && hi*hi*pipe_friction_prepared(hi,rr,rough)upper_bracket_evaluations++,hi*hi*pipe_friction_prepared(hi,rr,rough)bisections++;', 'status->bisections++;profile_bucket(profile_kind)->bisection_invalid_or_outside++;') source = source[:start] + resistance + PROFILE_SOLVE + source[end:] source = replace_once(source, ' if(fabs(p1-p2)<=1e-8) return 0;', ''' PipeProfile *profile=profile_bucket(kind);profile->flow_calls++; if(fabs(p1-p2)<=1e-8){profile->zero_pressure_calls++;return 0;}''') source = replace_once(source, ' if(4*lam/den<=1000) return sign*lam;', ' if(4*lam/den<=1000){profile->pnl00r_analytic_calls++;return sign*lam;}') source = replace_once(source, ' double base=area*p*cm/sqrt(T),K=pow(4*base/den,2)*d/length;\n return sign*native_pipe_resistance(K,rr,den/4,NULL)*den/4;', ' double base=area*p*cm/sqrt(T);\n return sign*profile_solve(base,d,length,rr,den,kind);') source = replace_once(source, ' if(cache->valid && cache->p1==p1 && cache->p2==p2 && cache->T==T &&', ' PipeProfile *profile=profile_bucket(kind);profile->cache_requests++;\n if(cache->valid && cache->p1==p1 && cache->p2==p2 && cache->T==T &&') source = replace_once(source, ' return cache->flow;\n double result=properties?', ' {profile->cache_hits++;return cache->flow;}\n profile->cache_misses++;\n double result=properties?') return source def instrument_common(source: str) -> str: """Count both ordinary and canonical Jacobian RHS evaluations exactly once.""" calls = [('model_eval', ' return model_eval(t,y,dy,w);')] if 'int native_jacobian_rhs(' in source: calls.append(('model_eval_jacobian', ' return model_eval_jacobian(t,y,dy,w);')) for function, fragment in calls: source = replace_once(source, fragment, ' extern void pipe_profile_rhs_enter(void),pipe_profile_rhs_leave(void);\n' f' pipe_profile_rhs_enter();int ok={function}(t,y,dy,w);pipe_profile_rhs_leave();return ok;') return source def git(*arguments: str) -> str: return subprocess.check_output(['git', *arguments], cwd=ROOT, text=True).strip() def prepare(args, out: Path) -> dict: if (out / 'prepared.json').exists(): metadata = json.loads((out / 'prepared.json').read_text()) if metadata['input_sha256'] != sha256(args.input.read_bytes()).hexdigest(): raise ValueError('Prepared model no longer matches the input file') if (metadata['script_sha256'] != sha256(Path(__file__).read_bytes()).hexdigest() or metadata.get('snapshot_helper_sha256') != sha256(SNAPSHOT_HELPER.read_bytes()).hexdigest()): raise ValueError('Diagnostic helper changed; choose a fresh output directory') return metadata if any((out / variant).exists() for variant in VARIANTS): raise ValueError('Partially prepared variant directories exist; choose a fresh output directory') previous = git('rev-parse', args.previous_ref) current = 'working-tree' if args.current_ref == 'working-tree' else git('rev-parse', args.current_ref) xml, doc = load_input(args.input) config = simulation_config(doc.simulation) if config.rtol != 1e-8 or config.t_start != 0 or config.t_stop != 10: raise ValueError('This experiment requires API default rtol=1e-8 and 0–10 s model settings') program = compile_native_program(compile_system_xml_network(doc)) out.mkdir(parents=True, exist_ok=True) (out / 'input.xml').write_bytes(xml) (out / 'input.json').write_bytes(args.input.read_bytes()) metadata = dict(input=str(args.input.resolve()), input_sha256=sha256(args.input.read_bytes()).hexdigest(), xml_sha256=sha256(xml).hexdigest(), script_sha256=sha256(Path(__file__).read_bytes()).hexdigest(), snapshot_helper_sha256=sha256(SNAPSHOT_HELPER.read_bytes()).hexdigest(), previous_revision=previous, current_revision=current, settings=vars(config), sample_step=doc.simulation.sample_step, variants={}, measurement_note='Instrumented times are diagnostic overhead and MUST NOT be used as production benchmark results.', counter_scope='rhs counts native_rhs/model_eval and, where present, native_jacobian_rhs/model_eval_jacobian calls including rejected trials and Jacobian differences; non_rhs includes initialization/output/probe evaluations.', replay_scope='All guarded RHS resistance calls, including scalar analytic low-Re cases; no sampling or deduplication. Cache hits, zero pressure difference and direct PNL00R analytic calls are excluded and counted separately.', capture_format='Native-endian IEEE-754 binary64 records: base, diameter, length, relative roughness, den=pi*d*mu, kind (six doubles, 48 bytes). Replay on the same host.', residual_test='At actual returned q, independently factored long-double f(Re) evaluates abs(Re^2*f(Re)/K-1)<=1e-9. Separate from each algorithm stopping rule.') original_native = builder.NATIVE prepared_builds = [] # Pin earlier variants while later variants are built. try: for variant in VARIANTS: directory = out / variant revision = current if variant == 'guarded-newton' else previous paths = (working_native_paths() if revision == 'working-tree' else git('ls-tree', '-r', '--name-only', revision, 'native').splitlines()) for name in paths: target = directory / name target.parent.mkdir(parents=True, exist_ok=True) data = ((ROOT / name).read_bytes() if revision == 'working-tree' else subprocess.check_output(['git', 'show', f'{revision}:{name}'], cwd=ROOT)) target.write_bytes(data) kernel = pipe_source_path(directory / 'native') original_hash = sha256(kernel.read_bytes()).hexdigest() kernel.write_text(instrument(kernel.read_text(), variant)) common = directory / 'native/runtime/common.c' common.write_text(instrument_common(common.read_text())) main = directory / 'native/runtime/main.c' main.write_text(replace_once(main.read_text(), 'int main(int argc, char **argv) {', 'int main(int argc, char **argv) {\n extern void pipe_profile_install(void);pipe_profile_install();')) build = build_snapshot(program, directory / 'native', out / 'cache', required_sources=(common, main)) prepared_builds.append(build) replay_source = directory / 'replay.c' replay_source.write_text(REPLAY_MAIN) replay = directory / ('replay.exe' if os.name == 'nt' else 'replay') compiler, _, _ = builder.toolchain() command = [compiler, '-std=c11', '-O3', '-Wall', '-Wextra', '-Werror', '-ffp-contract=off', '-fno-fast-math', '-I', str(directory / 'native/include'), str(replay_source), '-lm', '-o', str(replay)] compiled = subprocess.run(command, capture_output=True, text=True, timeout=60) (directory / 'replay-build.log').write_text(compiled.stdout + compiled.stderr) if compiled.returncode: raise RuntimeError(f'Replay compilation failed: {compiled.stderr}') metadata['variants'][variant] = dict(executable=str(build.executable), replay=str(replay), build_key=build.manifest['buildKey'], pipe_source=str(kernel.relative_to(directory)), compiled_source_hashes=build.manifest['sourceHashes'], original_kernel_sha256=original_hash, instrumented_kernel_sha256=sha256(kernel.read_bytes()).hexdigest()) print(f'Prepared {variant}', flush=True) finally: builder.NATIVE = original_native (out / 'prepared.json').write_text(json.dumps(metadata, ensure_ascii=False, indent=2) + '\n') return metadata def load_stats(path: Path) -> dict: stats = json.loads(path.read_text()) for scope in ('rhs', 'non_rhs'): for row in [stats[scope]['total'], *stats[scope]['by_kind'].values()]: assert row['cache_requests'] == row['cache_hits'] + row['cache_misses'] assert row['flow_calls'] == row['zero_pressure_calls'] + row['pnl00r_analytic_calls'] + row['resistance_calls'] or not row['flow_calls'] assert row['resistance_calls'] == row['iterative_calls'] + row['scalar_analytic_calls'] + row['failed_before_iteration'] assert row['resistance_calls'] == sum(row['iteration_histogram'].values()) assert row['iterations_total'] == sum(int(k) * v for k, v in row['iteration_histogram'].items()) assert row['resistance_calls'] == row['finite_returns'] + row['nonfinite_returns'] assert row['resistance_calls'] == row['residual_pass'] + row['residual_fail'] + row['residual_nonfinite'] assert row['exhausted_limit'] <= row['reached_last_iteration'] return stats def run(metadata: dict, out: Path, timeout: float): summary_path = out / 'summary.json' if summary_path.exists(): raise ValueError('Diagnostic results already exist; choose a fresh output directory') settings = metadata['settings'] capture = out / 'guarded-rhs-resistance-inputs.bin' rows = {} for variant in VARIANTS: directory = out / variant / 'trajectory' directory.mkdir(parents=True, exist_ok=False) result_path = directory / 'result.json' stats_path = directory / 'profile.json' env = os.environ.copy() env.pop('PIPE_PROFILE_CAPTURE', None) env['PIPE_PROFILE_JSON'] = str(stats_path) if variant == 'guarded-newton': env['PIPE_PROFILE_CAPTURE'] = str(capture) executable = metadata['variants'][variant]['executable'] command = [executable, '--method', settings['method'], '--start', str(settings['t_start']), '--stop', str(settings['t_stop']), '--sample-step', str(metadata['sample_step']), '--max-step', str(settings['max_step']), '--rtol', str(settings['rtol']), '--timeout', str(timeout), '--output', str(result_path)] print(f'Starting diagnostic trajectory: {variant}', flush=True) with (directory / 'worker.log').open('w') as log: try: process = subprocess.run(command, cwd=Path(executable).parent, env=env, stdout=subprocess.DEVNULL, stderr=log, timeout=timeout+15) exit_code = process.returncode except subprocess.TimeoutExpired: exit_code = 'external-timeout' data = json.loads(result_path.read_text()) if result_path.exists() else {} # Deliberately omit measured times from the cross-variant summary. trajectory = {key:data.get(key) for key in ('success', 'status', 'message', 'simulatedUntil', 'nfev', 'acceptedSteps', 'rejectedSteps', 'njev', 'nlu', 'stateTransitions')} trajectory['exit_code'] = exit_code rows[variant] = dict(trajectory=trajectory, profile=load_stats(stats_path) if stats_path.exists() else None) print(json.dumps(dict(variant=variant, **trajectory), ensure_ascii=False), flush=True) if capture.stat().st_size % 48: raise ValueError('Truncated replay capture') count = capture.stat().st_size // 48 guarded = rows['guarded-newton']['profile'] if guarded is None or count != guarded['rhs']['total']['resistance_calls'] or count != guarded['capture_records']: raise ValueError('Capture does not contain every guarded RHS resistance call') for variant in VARIANTS: directory = out / variant / 'replay-results' directory.mkdir(exist_ok=False) stats_path = directory / 'profile.json' env = os.environ.copy() env.pop('PIPE_PROFILE_CAPTURE', None) env['PIPE_PROFILE_JSON'] = str(stats_path) print(f'Replaying all {count} common inputs: {variant}', flush=True) with (directory / 'worker.log').open('w') as log: subprocess.run([metadata['variants'][variant]['replay'], str(capture)], env=env, stdout=subprocess.DEVNULL, stderr=log, timeout=300, check=True) replay = load_stats(stats_path) if replay['rhs']['total']['resistance_calls'] != count: raise ValueError('Replay count mismatch') rows[variant]['replay'] = replay # Identical input replay must reproduce every resistance-solve counter for # the guarded solver, excluding cache/flow bookkeeping performed upstream. exempt = {'cache_requests', 'cache_hits', 'cache_misses', 'flow_calls', 'zero_pressure_calls', 'pnl00r_analytic_calls'} for kind in ('total', '0', '1', '2', '3'): actual = guarded['rhs']['total'] if kind == 'total' else guarded['rhs']['by_kind'][kind] replay = rows['guarded-newton']['replay']['rhs']['total'] if kind == 'total' else rows['guarded-newton']['replay']['rhs']['by_kind'][kind] assert {k:v for k,v in actual.items() if k not in exempt} == {k:v for k,v in replay.items() if k not in exempt} summary = dict(metadata=metadata, capture=dict(path=str(capture), records=count, bytes=capture.stat().st_size, sha256=sha256(capture.read_bytes()).hexdigest(), sampling='none'), variants=rows) summary_path.write_text(json.dumps(summary, ensure_ascii=False, indent=2) + '\n') print(f'Diagnostic profile complete: {summary_path}', flush=True) def main(): parser = argparse.ArgumentParser(description=__doc__) parser.add_argument('--input', type=Path, default=ROOT / 'tests/data/test-mql-8-corrected.json') parser.add_argument('--output-dir', type=Path, required=True) parser.add_argument('--previous-ref', default='5d5a2e1') parser.add_argument('--current-ref', default='808c484', help='Audited guarded-solver revision, or working-tree for the current modular sources') parser.add_argument('--timeout', type=float, default=120) parser.add_argument('--prepare-only', action='store_true') args = parser.parse_args() out = args.output_dir.resolve() # Keep diagnostic native copies out of production sources and tracked data. if not out.is_relative_to(ROOT / 'test'): parser.error('--output-dir must be below the ignored repository test/ directory') metadata = prepare(args, out) if not args.prepare_only: run(metadata, out, args.timeout) if __name__ == '__main__': main()