"""Profile the current native BDF solver in an isolated diagnostic source copy. Example: python tests/manual/profile_native_solver.py --run Builds BOTH UD00 modes before timing; one excluded warmup and three serial pairs per mode by default. Production sources and cached workers are not instrumented. Requires the production GCC toolchain and SUNDIALS 7.4 public interfaces. """ from __future__ import annotations import argparse import hashlib import json import math import os from pathlib import Path import shutil import statistics import subprocess import sys import time ROOT = Path(__file__).resolve().parents[2] sys.path.insert(0, str(ROOT)) from app.main import compile_system_xml_network from app.simulation.native_codegen.build import ( build_native, toolchain, platform_build_inputs, link_library_arguments, _command, RUNTIME_SOURCES, ) from app.simulation.native_codegen.compiler import compile_native_program from app.simulation.native_codegen.input import load_input from app.simulation.native_codegen.modules import component_modules from tests.manual.native_compute_profile import replace_once, scope_function HERE = Path(__file__).resolve().parent APIS = { 'STEPS': 'CVodeGetNumSteps', 'RHS': 'CVodeGetNumRhsEvals', 'LINEAR_RHS': 'CVodeGetNumLinRhsEvals', 'JACOBIANS': 'CVodeGetNumJacEvals', 'SETUPS': 'CVodeGetNumLinSolvSetups', 'ERROR_FAILS': 'CVodeGetNumErrTestFails', 'NEWTON_ITERS': 'CVodeGetNumNonlinSolvIters', 'NEWTON_FAILS': 'CVodeGetNumNonlinSolvConvFails', 'STEP_SOLVE_FAILS': 'CVodeGetNumStepSolveFails', } # Every exclusive scope occurs exactly once in this disjoint partition. PARTITION = { 'residual_evaluation': ['rhs', 'nonlinear_residual'], 'jacobian_assembly': ['jacobian', 'jacobian_probe'], 'numerical_factorization': ['numerical_factorization'], 'linear_solve': ['linear_solve'], 'linear_system_setup_other': ['linear_setup'], 'newton_overhead': ['newton'], 'cvode_controller_including_error_estimation': ['cvode_controller'], 'event_detection': ['event_detection'], 'accepted_property_checks': ['accepted_property_checks'], 'sampling_and_storage': ['sampling', 'append', 'flush'], 'output_replay': ['output_replay'], 'result_encoding': ['result_encoding'], 'poll_and_progress': ['poll'], 'framework_other': ['total', 'integration'], } def digest(path): with path.open('rb') as stream: return hashlib.file_digest(stream, 'sha256').hexdigest() def write(path, data): path.write_text(json.dumps(data, ensure_ascii=False, indent=2) + '\n', encoding='utf-8') def scope(text, function, category): return scope_function(text, function, category).replace('PROFILE_' + category.upper(), 'P_' + category.upper()) def instrument(native): for name in ('solver_profile.h', 'solver_profile_hooks.h'): shutil.copy2(HERE / name, native / 'include' / name) shutil.copy2(HERE / 'solver_profile.c', native / 'runtime/solver_profile.c') for filename, functions in { 'common.c': {'native_poll': 'poll', 'native_accept': 'event', 'native_append': 'append', 'check_property_temperatures': 'property'}, 'cvode_solver.c': {'native_bdf': 'integration', 'cv_rhs': 'rhs', 'cv_jacobian': 'jacobian', 'jac_rhs': 'jacobian_probe'}, 'sample_storage.c': {'native_samples_flush': 'flush', 'native_samples_outputs': 'replay'}, 'main.c': {'main': 'total', 'write_result': 'encoding'}, }.items(): path = native / 'runtime' / filename text = '#include "solver_profile.h"\n' + path.read_text(encoding='utf-8') for name, category in functions.items(): text = scope(text, name, category) if filename == 'common.c': text = replace_once(text, ' if (r->options.record_samples) {', ' if (r->options.record_samples) {\n PROFILE_SCOPE(P_SAMPLING);') text = replace_once(text, ' return model_friction_drives(t,state,drives);', ' profile_add(C_FRICTION_EVALS,0,1);\n return model_friction_drives(t,state,drives);') elif filename == 'cvode_solver.c': text = replace_once(text, 'typedef struct { void *solver;', '#include "solver_profile_hooks.h"\ntypedef struct { void *solver;') text = replace_once(text, ' if(!context->jacobian_workspace)return jac_rhs(context,t,y,f);', ' if(!context->jacobian_workspace)return jac_rhs(context,t,y,f);\n PROFILE_SCOPE(P_JACOBIAN_PROBE);') text = replace_once(text, ' int success=0, initialized=0;', ' SUNNonlinearSolver profile_nls=NULL;\n int success=0, initialized=0;') text = replace_once(text, ' if (!linear) goto allocation_failure;', ' if (!linear) goto allocation_failure;\n' ' original_lu=linear->ops->setup;original_solve=linear->ops->solve;\n' ' linear->ops->setup=profiled_lu;linear->ops->solve=profiled_solve;') text = replace_once(text, ' CV_CHECK(CVodeSetMaxStep(solver,r->options.max_step));', ' CV_CHECK(CVodeSetMaxStep(solver,r->options.max_step));\n' ' profile_nls=profile_attach_newton(y,ctx);\n' ' allocation="SUNNonlinSol_Newton (profile)";\n' ' if(!profile_nls)goto allocation_failure;\n' ' CV_CHECK(CVodeSetNonlinearSolver(solver,profile_nls));') text = replace_once(text, ' if (linear) SUNLinSolFree(linear);', ' if(profile_nls)SUNNonlinSolFree(profile_nls);\n if (linear) SUNLinSolFree(linear);') text = replace_once(text, 'int flag=CVode(solver,end,y,&next,CV_ONE_STEP);', 'int flag=profiled_cvode(solver,end,y,&next,CV_ONE_STEP);') text = replace_once(text, ' if (toptions.stop) {', ' if (toptions.stop) {\n profile_add(C_BOUNDARIES,0,1);') snapshots = '\n profile_add(C_SEGMENTS,0,1);\n' for tag, api in APIS.items(): snapshots += f' {{long int n=0;int code={api}(solver,&n);profile_add(C_{tag},code,n);}}\n' text = replace_once(text, 'static void counters(NativeRun *r, void *solver) {', 'static void counters(NativeRun *r, void *solver) {' + snapshots) elif filename == 'main.c': text = replace_once(text, ' native_run_free(&r); return code;', ' native_run_free(&r); profile_end(&profile_scope);profile_dump();return code;') path.write_text(text, encoding='utf-8', newline='\n') def prepare(args): out = args.output.resolve() if not out.is_relative_to(ROOT / 'test') or out.exists(): raise ValueError('Use a NEW output directory beneath the project test/ directory') if args.repeats < 1 or args.warmups < 0: raise ValueError('repeats must be positive and warmups nonnegative') out.mkdir(parents=True) compiler, sundials, compiler_version = toolchain() flags, libraries, dlls, executable = platform_build_inputs(sundials) project = json.loads(args.input.read_text(encoding='utf-8')) variants = [] for cyclic in (False, True): case = out / ('cyclic' if cyclic else 'noncyclic') case.mkdir() ud00 = [n for n in project['nodes'] if n['data'].get('modelType') == 'amesim_ud00'] if not ud00: raise ValueError('UD00 comparison requires at least one UD00 signal') for node in ud00: node['data']['parameters']['iscyclic'] = int(cyclic) project['simulation']['t_stop'] = args.stop write(case / 'input.json', project) _, document = load_input(case / 'input.json') program = compile_native_program(compile_system_xml_network(document)) # The assembly timer observes our callback; CVODE's private dense-DQ # fallback has no separate public timer. Reject instead of mislabelling. if '#define MODEL_JACOBIAN_COLORED 1' not in program.header: raise ValueError('This diagnostic currently requires generated colored Jacobian support') build = build_native(program) try: shutil.copytree(build.executable.parent, case / 'control') manifest = build.manifest finally: build.close() native = case / 'native' shutil.copytree(ROOT / 'native', native) target = case / 'profiled' target.mkdir() for name in ('model.c', 'model.h', 'THIRD_PARTY_NOTICES.txt'): shutil.copy2(case / 'control' / name, target / name) for dll in dlls: shutil.copy2(dll, target / dll.name) instrument(native) sources = [native / name for name in RUNTIME_SOURCES] sources += [native / 'components/modules' / f'{name}.c' for name in component_modules(program.source + '\n' + program.header)] sources.append(native / 'runtime/solver_profile.c') command = [compiler, *flags, '-I', str(target), '-I', str(native / 'include'), '-I', str(sundials / 'include'), str(target / 'model.c'), *map(str, sources), *link_library_arguments(libraries), '-lm', '-o', str(target / executable)] log = [] try: _command(command, log=log, timeout=180) finally: (case / 'build.log').write_text('\n'.join(log), encoding='utf-8') assert digest(target / 'model.c') == digest(case / 'control/model.c') variants.append(dict(name=case.name, ud00Count=len(ud00), buildKey=manifest['buildKey'], executable=executable, inputSha256=digest(case / 'input.json'), controlSha256=digest(case / 'control' / executable), profiledSha256=digest(target / executable), nativeSourceHashes={p.relative_to(ROOT / 'native').as_posix(): digest(p) for p in (ROOT / 'native').rglob('*') if p.is_file()}, instrumentedHashes={p.relative_to(native).as_posix(): digest(p) for p in native.rglob('*') if p.is_file()})) print(f'Prepared {case.name}: {manifest["buildKey"]}', flush=True) result = dict(variants=variants, stop=args.stop, repeats=args.repeats, warmups=args.warmups, compiler=compiler_version, platform=sys.platform, revision=subprocess.check_output(['git', 'rev-parse', 'HEAD'], cwd=ROOT, text=True).strip(), settings=dict(method='BDF', start=0, sample_step=.01, rtol=1e-8, max_step=1e30), scope='Native worker: model initialization through result encoding and cleanup. ' 'Excludes build, process launch, API and browser. Error estimation remains inside CVODE controller. ' 'Dense LU has no symbolic factorization. Explicit Newton uses the same library implementation.') write(out / 'prepared.json', result) return result def analyze(profile, result): assert profile['counterErrors'] == 0 s, c = profile['scopes'], profile['counters'] assert set(s) == {k for group in PARTITION.values() for k in group} times = {name: sum(s[k]['exclusiveSeconds'] for k in scopes) for name, scopes in PARTITION.items()} total = s['total']['inclusiveSeconds'] assert math.isclose(sum(times.values()), total, abs_tol=1e-8, rel_tol=1e-10) assert all(t >= 0 for t in times.values()) for name, key in (('residual_evaluations', 'cvodeRhsCalls'), ('linear_rhs_evaluations', 'cvodeLinearRhsCalls'), ('jacobian_evaluations', 'njev'), ('linear_setups', 'nlu'), ('error_test_failures', 'rejectedSteps'), ('counter_segments', 'solverStarts')): assert c[name] == result[key], (name, c[name], result[key]) assert s['jacobian']['calls'] == c['jacobian_evaluations'] assert s['numerical_factorization']['calls'] == c['linear_setups'] assert c['jacobian_reuses'] + c['jacobian_refresh_setups'] + c['linear_setup_failures'] == c['linear_setups'] assert s['rhs']['calls'] == c['residual_evaluations'] + c['linear_rhs_evaluations'] assert result['nfev'] == s['rhs']['calls'] + result['jacobianRhsCalls'] + c['event_model_evaluations'] c.update(rejected_steps=c['error_test_failures'] + c['nonlinear_step_failures'], application_accepted_steps=result['acceptedSteps'], same_time_returns=result['solverControl']['sameTimeReturns'], max_same_time_streak=result['solverControl']['maxSameTimeStreak'], nonlinear_residual_calls=s['nonlinear_residual']['calls'], jacobian_probe_evaluations=result['jacobianRhsCalls'], LU_factorizations=s['numerical_factorization']['calls'], LU_solves=s['linear_solve']['calls'], newton_solve_calls=s['newton']['calls'], event_count=result['stateTransitions'], event_detection_calls=s['event_detection']['calls'], solver_starts=result['solverStarts'], all_counted_model_evaluations=result['nfev'], accepted_property_checks=s['accepted_property_checks']['calls'], sample_count=len(result['series']['time']), jacobian_kernel_reuse=result['jacobianReuse']) return dict(totalSeconds=total, exclusiveSeconds=times, counters=c) def run(args, prepared): rows = [] for variant in prepared['variants']: case = args.output.resolve() / variant['name'] baseline = None for index in range(-args.warmups, args.repeats): label = f'warmup-{index + args.warmups}' if index < 0 else f'run-{index}' for kind in ('control', 'profiled'): dest = case / f'{kind}-{label}' dest.mkdir() command = [str(case / kind / variant['executable']), '--method', 'BDF', '--start', '0', '--stop', str(args.stop), '--sample-step', '.01', '--rtol', '1e-8', '--max-step', '1e30', '--timeout', '300', '--output', str(dest / 'result.json'), '--sample-file', str(dest / 'states.bin'), '--output-block-file', str(dest / 'outputs.bin')] env = dict(os.environ) env.pop('NATIVE_SOLVER_PROFILE', None) if kind == 'profiled': env['NATIVE_SOLVER_PROFILE'] = str(dest / 'profile.json') start = time.perf_counter() with (dest / 'stderr.log').open('wb') as err: process = subprocess.run(command, stdout=subprocess.DEVNULL, stderr=err, env=env, timeout=330) elapsed = time.perf_counter() - start result = json.loads((dest / 'result.json').read_bytes()) assert process.returncode == 0 and result['success'], str(dest) comparable = {k: v for k, v in result.items() if k not in ('solveSeconds', 'solveCpuSeconds')} comparable.update(statesSha256=digest(dest / 'states.bin'), outputsSha256=digest(dest / 'outputs.bin')) if baseline is None: baseline = comparable differences = [k for k in baseline.keys() | comparable.keys() if baseline.get(k) != comparable.get(k)] assert not differences, (str(dest), differences) row = dict(case=case.name, variant=kind, label=label, warmup=index < 0, processWallSeconds=elapsed, solveSeconds=result['solveSeconds'], solveCpuSeconds=result['solveCpuSeconds'], fullParity=True, statesSha256=comparable['statesSha256'], outputsSha256=comparable['outputsSha256']) if kind == 'profiled': row['profile'] = analyze(json.loads((dest / 'profile.json').read_bytes()), result) write(dest / 'measurement.json', row) rows.append(row) print(f'{case.name}/{kind}/{label}: solve {result["solveSeconds"]:.4f}s; parity OK', flush=True) # Release the large decoded series before the next case. del baseline write(args.output / 'measurements.json', rows) return summarize(args.output, prepared, rows) def summarize(output, prepared, rows): summary = dict(scope=prepared['scope'], cases=[]) for variant in prepared['variants']: selected = [r for r in rows if r['case'] == variant['name'] and not r['warmup']] profiled = [r for r in selected if r['variant'] == 'profiled'] control = [r for r in selected if r['variant'] == 'control'] c = profiled[0]['profile']['counters'] assert all(r['profile']['counters'] == c for r in profiled) raw_profiles = [json.loads((output / variant['name'] / f'profiled-{r["label"]}' / 'profile.json').read_bytes()) for r in profiled] times = {key: statistics.mean(r['profile']['exclusiveSeconds'][key] for r in profiled) for key in PARTITION} total = statistics.mean(r['profile']['totalSeconds'] for r in profiled) medians = {name: {key: statistics.median(r[key] for r in data) for key in ('solveSeconds', 'solveCpuSeconds', 'processWallSeconds')} for name, data in (('control', control), ('profiled', profiled))} summary['cases'].append(dict(name=variant['name'], meanTotalSeconds=total, phases={key: dict(seconds=value, percent=100 * value / total) for key, value in times.items()}, meanScopes={key: {metric: statistics.mean(r['scopes'][key][metric] for r in raw_profiles) for metric in ('inclusiveSeconds', 'exclusiveSeconds')} for key in raw_profiles[0]['scopes']}, symbolic_factorization=dict(seconds=0, status='not applicable: dense LU'), error_estimation=dict(seconds=None, status='included in cvode_controller_including_error_estimation'), counters=c, medians=medians, observedOverheadPercent={key: 100 * (medians['profiled'][key] / medians['control'][key] - 1) for key in medians['control']}, allRunsFullParity=True, samples=len(profiled))) write(output / 'summary.json', summary) print(json.dumps(summary, ensure_ascii=False, indent=2), flush=True) return summary def report_existing(output): prepared = json.loads((output / 'prepared.json').read_bytes()) rows = json.loads((output / 'measurements.json').read_bytes()) for row in rows: if row['variant'] != 'profiled': continue directory = output / row['case'] / f'profiled-{row["label"]}' row['profile'] = analyze(json.loads((directory / 'profile.json').read_bytes()), json.loads((directory / 'result.json').read_bytes())) write(directory / 'measurement.json', row) write(output / 'measurements.json', rows) return summarize(output, prepared, rows) if __name__ == '__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', type=Path, default=ROOT / 'test/solver-profile-20260916') parser.add_argument('--stop', type=float, default=21.7) parser.add_argument('--repeats', type=int, default=3) parser.add_argument('--warmups', type=int, default=1) mode = parser.add_mutually_exclusive_group() mode.add_argument('--run', action='store_true', help='also run after all builds complete') mode.add_argument('--report-only', action='store_true', help='rebuild summary from completed runs; no build or simulation') options = parser.parse_args() if options.report_only: report_existing(options.output) else: prepared = prepare(options) if options.run: run(options, prepared)