"""Static C lowering for the complete built-in component catalog.

Python constructs the graph and eliminates constant linear constraints once.
All thermodynamics, flow/stream closure and derivatives execute in the EXE.
"""
from __future__ import annotations

import json
import math
import re
from importlib import import_module

from .compiler import NativeCapabilityError, NativeProgram, _Groups, _number as num
from .contracts import SUPPORTED_VERSIONS
from .schedule import Computation, EvaluationSchedule, references
from .tolerances import state_absolute_tolerance
from .jacobian import StateDependencies, expression_inputs
from app.simulation.components.amesim.semantics import contact_stiffness, validate_numerical_semantics


GAS_TYPES = {'amesim_pnch023', 'amesim_pnch012', 'amesim_pnl0001',
             'amesim_pnl0002', 'amesim_pnl0003', 'cylinder', 'tank'}
NODES = {'amesim_pn3node2', 'amesim_p4node2', 'tee'}
RESISTORS = {'amesim_pnor001', 'amesim_pnvo001_fixed', 'amesim_pnvo001',
             'amesim_pnl00r', 'pipe', 'orifice'}


def is_polytropic(c):
    return c.model_type in {'amesim_pnl0001', 'amesim_pnl0002', 'amesim_pnl0003'} and int(c.mode) == 1


def catalog_contracts():
    from app.simulation.components.amesim.library import LIBRARY as a
    from app.simulation.components.experimental.library import LIBRARY as e
    result = {}
    for entry in (*a.models, *e.models):
        module, name = entry.split(':')
        cls = getattr(import_module(module), name)
        result[cls.MODEL_TYPE] = cls
    return result


def linear_assignments(equations, unknowns, free_flows=()):
    """Eliminate constant coefficients; retain RHS expressions as C temporaries.

    Exact pivot elimination avoids a numeric pseudoinverse and its tiny spurious
    dependencies. Redundant rows are left to the nodal pressure closure.
    """
    rows = [[dict(a), {i: 1.0}] for i, (a, _) in enumerate(equations)]
    pivots = []
    for key in unknowns:
        found = next((i for i in range(len(pivots), len(rows)) if abs(rows[i][0].get(key, 0)) > 1e-14), None)
        if found is None and key in free_flows:
            # Ideal coupled pipe compliances have a redundant internal flow.
            # Their total m/U derivative is distributed by physical volume.
            index = len(equations)
            equations.append(({key: 1.0}, '0.0'))
            rows.append([{key: 1.0}, {index: 1.0}])
            found = len(rows)-1
        if found is None:
            raise NativeCapabilityError(f'Underdetermined native connection constraint: {key}')
        j = len(pivots)
        rows[j], rows[found] = rows[found], rows[j]
        a, b = rows[j]
        pivot = a[key]
        rows[j] = [{k: v/pivot for k, v in a.items()}, {k: v/pivot for k, v in b.items()}]
        for i, (a, b) in enumerate(rows):
            if i == j:
                continue
            scale = a.get(key, 0)
            if not scale:
                continue
            for target, source in zip((a, b), rows[j]):
                for k, v in source.items():
                    value = target.get(k, 0) - scale*v
                    if abs(value) < 1e-14:
                        target.pop(k, None)
                    else:
                        target[k] = value
        pivots.append(key)
    lines = [(f'b{i}', expr) for i, (_, expr) in enumerate(equations)
             if any(i in rows[j][1] for j in range(len(pivots)))]
    for j, key in enumerate(pivots):
        expr = ' + '.join(f'({num(v)})*b{i}' for i, v in rows[j][1].items()) or '0.0'
        lines.append((key, expr))
    return lines


def linear_schedule(equations, unknowns, free_flows=()):
    return [f'{"double " if key.startswith("b") else ""}{key} = {expr};'
            for key, expr in linear_assignments(equations, unknowns, free_flows)]


def compile_extended_program(network):
    network.validate_port_supplies()
    components = list(network.components.values())
    contracts = catalog_contracts()
    for c in components:
        if type(c) is not contracts.get(c.model_type):
            raise NativeCapabilityError(f'{c.name}: no native contract for {c.model_type} / {type(c).__name__}')
        if c.MODEL_VERSION != SUPPORTED_VERSIONS.get(c.model_type):
            raise NativeCapabilityError(f'{c.name}: native kernel version does not match the component contract')
        try:
            validate_numerical_semantics(c)
        except ValueError as exc:
            raise NativeCapabilityError(str(exc)) from exc
    variables = tuple(v for c in components for v in c.result_variable_metadata())
    if not variables:
        raise NativeCapabilityError('Simulation requires a runtime component; medium definitions alone have no outputs')
    slots = {v.key: i for i, v in enumerate(variables)}
    ports = {(c.name, p.name): p for c in components for p in c.active_port_definitions}
    groups = _Groups(ports)
    adjacent = {}
    for edge in network.connections:
        a, b = (p.key for p in edge.endpoints)
        if edge.kind == 'physical':
            adjacent[a], adjacent[b] = b, a
            groups.union(a, b)
        else:
            source, target = (a, b) if ports[a].nominal_role == 'output' else (b, a)
            adjacent[target] = source
    for c in components:
        for port in c.required_connection_ports:
            if (c.name, port) not in adjacent:
                raise NativeCapabilityError(f'{c.name}.{port}: unconnected required port')
        names = [p.name for p in c.active_port_definitions]
        if c.model_type in NODES | {'amesim_pnch023', 'amesim_pnch012', 'amesim_mecmas21', 'amesim_lmechn1'} or (c.model_type=='pipe' and c.lambda_darcy==0):
            for p in names[1:]:
                groups.union((c.name, names[0]), (c.name, p))
        if c.model_type == 'amesim_pnrp17':
            for a, b in [('port_2', 'port_5'), ('port_3', 'port_4')]:
                groups.union((c.name, a), (c.name, b))

    state_keys, initial, gas_initializers = [], [], []
    states = {}
    def add_states(c, fields, values):
        for field, value in zip(fields, values):
            states[c.name, field] = len(state_keys)
            state_keys.append(f'{c.name}.{field}')
            initial.append(float(value))
    mass_groups = {}
    for c in components:
        if c.model_type in GAS_TYPES:
            thermal = 'T' if is_polytropic(c) else 'U'
            fields = ('m1', thermal+'1', 'm2', thermal+'2') if c.model_type == 'amesim_pnl0003' else ('m', thermal)
            add_states(c, fields, [0.0]*len(fields))
        elif c.model_type == 'amesim_mecmas21':
            root = groups.find((c.name, 'port_1'))
            mass_groups.setdefault(root, []).append(c)
            if len(mass_groups[root]) == 1:
                add_states(c, ('v', 'x'), (c.v0, c.x0))
            else:
                ref = mass_groups[root][0]
                if abs(c.v0-ref.v0)>1e-10*max(abs(c.v0),1) or abs(c.x0-ref.x0)>1e-10*max(abs(c.x0),1):
                    raise NativeCapabilityError('Rigidly connected masses require consistent initial x/v')
                for field in ('v', 'x'):
                    states[c.name, field] = states[ref.name, field]
    friction_groups = []
    for group in mass_groups.values():
        dry = [c for c in group if c.use_friction and int(c.stoptype) != 3 and c.fstick > 0]
        if dry:
            ref = group[0]
            threshold = min(c.dvel if int(c.frictionType) == 2 else 0 for c in dry)
            mode = (1 if ref.v0 > 0 else -1) if abs(ref.v0) > threshold else 0
            if mode == 0:
                initial[states[ref.name, 'v']] = 0.0
            add_states(ref, ('_friction_mode',), (mode,))
            friction_groups.append((group, dry, threshold))
    if len(state_keys)>1024 or len(variables)>16384:
        raise NativeCapabilityError('Native model exceeds the 1024-state / 16384-output resource limit')
    # Algebraic models use an internal constant state; the public state map stays empty.
    nstates = max(len(state_keys), 1)
    if not initial:
        initial = [0.0]
    dependencies = StateDependencies(len(state_keys))
    lines, declarations, breaks = [], [], []
    assigned = set()
    def w(c, field):
        name = c if isinstance(c, str) else c.name
        return f'w[{slots[name+"."+field]}]'
    def put(c, field, expr, dest=None):
        (lines if dest is None else dest).append(f'{w(c,field)} = {expr};')
        dependencies.expression(w(c,field), expr)
        assigned.add((c if isinstance(c,str) else c.name)+'.'+field)
    def y(c, field):
        return f'y[{states[c.name,field]}]'
    def ep(c, port):
        return c.name, port
    def pnames(c):
        return [p.name for p in c.active_port_definitions if p.domain == 'pneumatic']
    def mnames(c):
        return [p.name for p in c.active_port_definitions if p.domain == 'mechanical']
    for c in components:
        if c.model_type == 'amesim_mecmas21':
            for f in ('v', 'x'):
                put(c, f, y(c, f))
        for name in mnames(c):
            root = groups.find(ep(c, name))
            if root not in mass_groups:
                raise NativeCapabilityError(f'{c.name}.{name}: mechanical group has no inertia anchor')
            mass = mass_groups[root][0]
            for f in ('v', 'x'):
                put(c, name+'.'+f, y(mass, f))
        if c.model_type == 'amesim_step0':
            put(c, 'y', f't < {num(c.time)} ? {num(c.initial)} : {num(c.final)}')
            breaks.append(f'if(t < {num(c.time)}) result=fmin(result,{num(c.time)});')
        elif c.model_type == 'amesim_ud00':
            ident = 'signal_'+str(len(declarations))
            declarations.append(f'static const double {ident}[24] = {{'+','.join(num(v) for vs in (c.starts,c.ends,c.durations) for v in vs)+'};')
            args=f'{num(c.tstart)}, {c.nstages}, {int(c.iscyclic)}, {ident}'
            put(c, 'y', f'native_signal(t, {args})')
            breaks.append(f'result=fmin(result,native_signal_break(t,end,{args}));')
        else:
            continue
        put(c, 'out.signal', w(c, 'y'))
    for c in components:
        for port in c.active_port_definitions:
            if port.kind == 'signal' and port.nominal_role == 'input':
                target = adjacent.get(ep(c, port.name))
                if target:
                    expr = w(target[0], target[1]+'.signal')
                elif c.model_type == 'amesim_pnvo001':
                    expr = num(c.opening0)
                else:
                    raise NativeCapabilityError(f'{c.name}.{port.name}: signal input missing')
                put(c, port.name+'.signal', expr)

    media = {}
    def medium(c):
        m = c.medium
        key = (getattr(m, 'SUBSTANCE_ID', None), getattr(m, 'PROPERTY_METHOD_ID', None),
               m.R_gas, m.cp_ref, m.T_ref, m.cp_slope, m.viscosity_ref, m.viscosity_T_ref, m.sutherland_constant)
        if key not in media:
            if key[1] not in (None, 'ideal_gas', 'peng_robinson') or (key[1]=='peng_robinson' and key[0]!='helium'):
                raise NativeCapabilityError(f'{c.name}: unsupported medium contract {key[:2]}')
            ident = f'medium_{len(media)}'
            declarations.append(f'static const NativeMedium {ident} = {{'+str(int(key[1]=='peng_robinson'))+','+','.join(num(v) for v in key[2:])+'};')
            media[key]=ident
        return '&'+media[key]

    pistons = [c for c in components if c.model_type == 'amesim_pnrp17']
    for c in pistons:
        put(c, 'length', f'{num(c.x0)}+{w(c,"port_5.x")}-{w(c,"port_4.x")}')
        put(c, 'volume', f'{num(c.effective_area)}*{w(c,"length")}')
        put(c, 'volume_flow', f'{num(c.effective_area)}*({w(c,"port_5.v")}-{w(c,"port_4.v")})')

    pneu = [endpoint for endpoint, port in ports.items() if port.domain == 'pneumatic']
    pi = {endpoint: i for i, endpoint in enumerate(pneu)}
    pgroups = list(dict.fromkeys(groups.find(endpoint) for endpoint in pneu))
    pgi = {root: i for i, root in enumerate(pgroups)}
    def p(c, name): return f'p[{pgi[groups.find(ep(c,name))]}]'
    def q(c, name): return f'q[{pi[ep(c,name)]}]'
    def h(c, name): return f'h[{pi[ep(c,name)]}]'
    def hin(c, name, temperature=False):
        other = adjacent[ep(c,name)]
        obj = network.components[other[0]]
        if temperature and obj.model_type in NODES-{'tee'}:
            other = adjacent[(obj.name, 'port_2')]
        return f'h[{pi[other]}]'
    gases, anchor, port_gas, volume_rate = {}, {}, {}, {}
    anchor_partitions = {}
    gas_count = 0
    for c in components:
        if c.model_type not in GAS_TYPES:
            continue
        kind = c.model_type
        if kind == 'amesim_pnch023':
            volume, rate = num(c.cvol), '0.0'
        elif kind == 'amesim_pnch012':
            attached = [d for d in pistons if adjacent[ep(d,'port_1')][0] == c.name]
            volume = num(c.cvol0+sum(c.external_volumes.values()))+''.join('+'+w(d,'volume') for d in attached)
            put(c, 'vol', f'fmax({num(c.cvol0/100)}, {volume})')
            rate = num(sum(c.external_volume_rates.values()))+''.join('+'+w(d,'volume_flow') for d in attached)
            put(c, 'dvol', f'{w(c,"vol")} <= {num(c.cvol0/100)} ? 0.0 : ({rate})')
            volume, rate = w(c,'vol'), w(c,'dvol')
        elif kind in ('cylinder', 'tank'):
            volume, rate = num(c.V), '0.0'
        else:
            volume, rate = num(c.volume), '0.0'
        volume_rate[c.name] = rate
        halves = (1,2) if kind == 'amesim_pnl0003' else (0,)
        for half in halves:
            suffix = str(half) if half else ''
            gas = f'g[{gas_count}]';gas_count += 1
            gases[c.name,half] = gas
            V = num(c.compliance_volume) if half else volume
            if half:
                p0,T0=c.parameter_values[f'p{half}_0'],c.parameter_values[f'T{half}_0']
                initial_volume=c.compliance_volume
            else:
                p0,T0=c.parameter_values['p0'],c.parameter_values['T0']
                initial_volume=(max(c.cvol0+sum(c.external_volumes.values()),c.cvol0/100)
                                if kind=='amesim_pnch012' else
                                c.cvol if kind=='amesim_pnch023' else c.V if kind in ('cylinder','tank') else c.volume)
                if kind == 'amesim_pnch012':
                    def initial_x(d, port):
                        return mass_groups[groups.find(ep(d, port))][0].x0
                    initial_volume = max(c.cvol0 + sum(c.external_volumes.values()) + sum(
                        d.effective_area * (d.x0 + initial_x(d, 'port_5') - initial_x(d, 'port_4'))
                        for d in attached), c.cvol0 / 100)
            thermal = ('T' if is_polytropic(c) else 'U') + suffix
            if is_polytropic(c):
                gas_initializers += [f'{y(c,"m"+suffix)}={num(p0*initial_volume/(c.medium.R_gas*T0))};', f'{y(c,thermal)}={num(T0)};']
                gas_function = 'native_polytropic_gas_context'
            else:
                gas_initializers.append(f'if(!native_medium_init({medium(c)},{num(p0)},{num(T0)},{num(initial_volume)},{int(kind in ("cylinder","tank"))},&y[{states[c.name,"m"+suffix]}])) return 0;')
                gas_function = 'native_medium_gas_context'
            lines.append(f'if(!{gas_function}(gas_properties,{medium(c)}, {y(c,"m"+suffix)}, {y(c,thermal)}, {V}, &{gas})) return 0;')
            gas_inputs = expression_inputs(f'{y(c,"m"+suffix)}+{y(c,thermal)}+({V})')
            for field in ('p','T','rho','u','h'):
                dependencies.assign(gas+'.'+field, gas_inputs)
            put(c, 'm'+suffix, y(c, 'm'+suffix))
            put(c, 'U'+suffix, f'{y(c,"m"+suffix)}*{gas}.u' if is_polytropic(c) else y(c,'U'+suffix))
            for field in ('p','T','rho','u','h'):
                put(c, field+suffix, gas+'.'+field)
            anchored = pnames(c)
            if kind == 'amesim_pnl0001': anchored=['port_2']
            if kind == 'amesim_pnl0002': anchored=[]
            if half: anchored=['port_'+suffix]
            for name in anchored:
                root=groups.find(ep(c,name))
                anchor[root]=gas
                partition=(c.name,half)
                anchor_partitions.setdefault(root, {})[partition]=(c, float(c.compliance_volume if half else c.volume) if kind in ('amesim_pnl0001','amesim_pnl0003') else None)
            for name in (['port_'+suffix] if half else pnames(c)):
                port_gas[ep(c,name)]=gas
    coupled=[]
    project=[]
    for root,partitions in anchor_partitions.items():
        if len(partitions)<2: continue
        if any(volume is None for _,volume in partitions.values()) or len({id(c.medium) for c,_ in partitions.values()})!=1:
            raise NativeCapabilityError('Direct gas-storage coupling requires compatible fixed pipe compliances; insert a resistance between independent chambers')
        laws = {(is_polytropic(c), c.k if is_polytropic(c) else None) for c,_ in partitions.values()}
        if len(laws) != 1:
            raise NativeCapabilityError('Directly coupled pipe compliances require the same thermodynamic mode and polytropic exponent')
        poly = next(iter(laws))[0]
        offsets=[];volumes=[]
        for (name,half),(c,volume) in partitions.items():
            offsets.append(states[name,'m'+(str(half) if half else '')]);volumes.append(volume)
        for prefix in ('p','T'):
            values=[c.parameter_values[f'{prefix}{half}_0'] if half else getattr(c,prefix+'0')
                    for (_,half),(c,_) in partitions.items()]
            if max(values)-min(values)>1e-9*max(1,*map(abs,values)):
                raise NativeCapabilityError('Ideally coupled pipe compliances require consistent initial pressure and temperature')
        total=sum(volumes)
        thermal_sum = '+'.join(f'y[{i+1}]*{num(v/total)}' if poly else f'y[{i+1}]' for i,v in zip(offsets,volumes))
        project.append('{ double mass='+ '+'.join(f'y[{i}]' for i in offsets)+',thermal='+thermal_sum+';')
        for i,volume in zip(offsets,volumes):
            project += [f'projected[{i}]=mass*{num(volume/total)};projected[{i+1}]=thermal*{num(1 if poly else volume/total)};']
        project.append('}')
        coupled.append((root,offsets,volumes,poly))
        dependencies.project_states(offsets)
        dependencies.project_states([i+1 for i in offsets])
    for root, gas in anchor.items():
        lines.append(f'p[{pgi[root]}]={gas}.p;')
        dependencies.assign(f'p[{pgi[root]}]', (gas+'.p',))
    h_initial = {f'h[{pi[endpoint]}]': gas+'.h' for endpoint,gas in port_gas.items()}

    operations, flow_known, flow_eq = [], set(), []
    pipe_cache_count = 0
    def pipe_flow(expression):
        nonlocal pipe_cache_count
        result = f'native_pipe_flow_cached_context(properties,&pipe_cache[{pipe_cache_count}],{expression})'
        pipe_cache_count += 1
        return result
    def flow(c, name, expr):
        target=q(c,name)
        operations.append(Computation.assignment(f'flow:{c.name}.{name}',target,expr))
        flow_known.add(target)
    def flow_equation(terms):
        flow_eq.append(({q(c,name):coef for c,name,coef in terms},'0.0'))
    for c in components:
        kind=c.model_type
        names=pnames(c)
        if kind in RESISTORS or kind in NODES:
            flow_equation([(c,name,1) for name in names])
        if kind in {'amesim_pnpl01','amesim_pnrp17'}:
            flow(c,'port_1','0.0')
        if kind in RESISTORS:
            a,b=names
            pa,pb=p(c,a),p(c,b)
            if kind=='pipe' and c.lambda_darcy==0:
                # A zero-loss pipe is an ideal connection, whose flow is
                # determined by the neighbouring constitutive equations.
                continue
            if kind=='orifice':
                flow(c,a,f'{num(c.K_eff)}*copysign(sqrt(fabs({pa}-{pb})),{pa}-{pb})')
            elif kind=='pipe':
                resistance=c.lambda_darcy*c.L/c.D
                flow(c,a,f'copysign(sqrt(fabs({pa}-{pb})*2*fmax(native_density_context(properties,{medium(c)},fmax(.5*({pa}+{pb}),1),{num(c.T0)}),1e-12)*{num(c.area*c.area/resistance)}),{pa}-{pb})')
            elif kind=='amesim_pnl00r':
                T=f'native_temperature_ph_context(properties,{medium(c)},fmax(fmax({pa},{pb}),1),{pa}>={pb}?{hin(c,a,True)}:{hin(c,b,True)})'
                flow(c,a,pipe_flow(f'{medium(c)},{pa},{pb},{T},{num(c.diam)},{num(c.le)},{num(c.rr)},0'))
            else:
                opening = w(c,'xv') if kind!='amesim_pnor001' else '1.0'
                if kind!='amesim_pnor001':
                    value = f'fmax(0,fmin(1,{w(c,"res.signal")}))' if kind=='amesim_pnvo001' else num(c.opening)
                    put(c,'xv',value)
                area = c.effective_cq*(c.effective_area if kind=='amesim_pnor001' else c.maximum_area)
                inputs=f'{medium(c)},{pa},{pb},{hin(c,a)},{hin(c,b)},{num(area)},{opening}'
                outputs=(q(c,a),w(c,'cm'),w(c,'gasvel'))
                code=f'if(!native_medium_orifice_context(properties,{inputs},&{outputs[0]},&{outputs[1]},&{outputs[2]})) return 0;'
                operations.append(Computation(f'flow:{c.name}.{a}',outputs,references(inputs),(code,)))
                flow_known.add(q(c,a));assigned.update((c.name+'.cm',c.name+'.gasvel'))
            flow(c,b,f'-{q(c,a)}')
        elif kind in ('amesim_pnl0001','amesim_pnl0002'):
            gas=gases[c.name,0]
            for name in (['port_1'] if kind.endswith('1') else names):
                # Resistance uses the gas arriving from the upstream side,
                # including inflow into a PNL0001 storage volume.
                T=f'({p(c,name)}>{gas}.p?native_temperature_ph_context(properties,{medium(c)},fmax({p(c,name)},1),{hin(c,name,True)}):{gas}.T)'
                flow(c,name,pipe_flow(f'{medium(c)},{p(c,name)},{gas}.p,{T},{num(c.diam)},{num(c.le/(2 if kind.endswith("2") else 1))},{num(c.rr)},1'))
    for edge in network.connections:
        if edge.domain=='pneumatic':
            flow_eq.append(({f'q[{pi[e.key]}]':1 for e in edge.endpoints},'0.0'))
    unknownq=[f'q[{i}]' for i in range(len(pneu)) if f'q[{i}]' not in flow_known]
    reduced=[]
    for terms,rhs in flow_eq:
        known=''.join(f'-({num(v)})*{k}' for k,v in terms.items() if k in flow_known)
        reduced.append(({k:v for k,v in terms.items() if k not in flow_known},rhs+known))
    free_flows={f'q[{pi[e]}]' for root,_,_,_ in coupled for e in pneu if groups.find(e)==root}
    flow_bindings=linear_assignments(reduced,unknownq,free_flows)
    flow_temporary_count=max((int(key[1:])+1 for key,_ in flow_bindings if key.startswith('b')),default=0)
    for key,expr in flow_bindings:
        target=re.sub(r'\bb(\d+)\b',r'fb[\1]',key)
        value=re.sub(r'\bb(\d+)\b',r'fb[\1]',expr)
        operations.append(Computation.assignment('connection:'+key,target,value,'linear'))

    unknownp=[root for root in pgroups if root not in anchor]
    for root in unknownp:
        terms=[f'q[{pi[e]}]' for e in pneu if groups.find(e)==root and f'q[{pi[e]}]' in flow_known]
        if not terms:
            raise NativeCapabilityError('Unanchored pneumatic pressure group has no constitutive flow relation')
        expr=' + '.join(terms)
        operations.append(Computation('pressure:'+str(root),(f'p[{pgi[root]}]',),references(expr),(),'pressure',expr))

    for c in components:
        names=pnames(c);kind=c.model_type
        if kind in {'amesim_pnpl01','amesim_pnrp17'}:
            operations.append(Computation.assignment(f'alias:{c.name}.port_1',h(c,'port_1'),hin(c,'port_1'),'alias'))
        elif kind in RESISTORS:
            a,b=names
            for name,other in ((a,b),(b,a)):
                operations.append(Computation.assignment(f'alias:{c.name}.{name}',h(c,name),hin(c,other),'alias'))
        elif kind in NODES:
            if kind!='tee':
                # Reference copies have no flow dependency. Keep them separate
                # from the reference port's returned energy/mixing calculation.
                for name in names:
                    if name!='port_2':
                        operations.append(Computation.assignment(f'alias:{c.name}.{name}',h(c,name),hin(c,'port_2'),'alias'))
            stream_lines = ['{ double total=0,energy=0,average=0;', *[f'if({q(c,name)}>1e-12) {{total+={q(c,name)};energy+={q(c,name)}*{hin(c,name)};}} average+={hin(c,name)};' for name in names]]
            if kind=='tee':
                stream_lines += [f'double mixed=total>1e-12?energy/total:average/{len(names)};', *[f'{h(c,name)}=mixed;' for name in names]]
            else:
                stream_lines += [f'double ref={hin(c,"port_2")},mixed=total>1e-12?energy/total:ref;', f'{h(c,"port_2")}=mixed;',f'if({q(c,"port_2")}<0) {{ double e=0,scale=0;',
                                 *[f'e+={q(c,name)}*({q(c,name)}>1e-12?{hin(c,name)}:ref);scale+=fabs({q(c,name)});' for name in names if name!='port_2'],
                                 f'double flow={q(c,"port_2")},transition=fmax(.05*scale,1e-12);',
                                 'double inv=-flow>=transition?1/flow:flow*(2*transition*transition-flow*flow)/pow(transition,4);',
                                 f'{h(c,"port_2")}=mixed-(e+flow*mixed)*inv;','}']
                # Nodes expose their reference temperature as an output too.
                if c.name+'.T' in slots:
                    raise NativeCapabilityError('Unexpected node temperature output contract')
            stream_lines += ['}']
            outputs=tuple(h(c,name) for name in (names if kind=='tee' else ['port_2']))
            inputs=frozenset(q(c,name) for name in names)|frozenset(hin(c,name) for name in names)
            operations.append(Computation(f'stream:{c.name}',outputs,inputs,tuple(stream_lines),'stream'))

    # Current state/parameter/signal values are available before this phase.
    # Default h guesses are deliberately absent from the known-source map.
    labels={f'w[{i}]': key for key,i in slots.items()}
    known={f'w[{slots[key]}]': 'prepared:'+key for key in assigned}
    gas_origins={gas: 'state:'+','.join(f'{name}.{field}{str(half) if half else ""}' for field in ('m','U'))
                 for (name,half),gas in gases.items()}
    for gas,origin in gas_origins.items():
        for field in ('p','T','rho','u','h'):
            known[gas+'.'+field]=origin
    for root,gas in anchor.items():
        known[f'p[{pgi[root]}]']=gas_origins[gas]
    for endpoint,gas in port_gas.items():
        known[f'h[{pi[endpoint]}]']=gas_origins[gas]
    for endpoint,index in pi.items():
        labels[f'h[{index}]']='.'.join(endpoint)+'.h_outflow'
        labels[f'q[{index}]']='.'.join(endpoint)+'.m_flow'
    for root,index in pgi.items():
        labels[f'p[{index}]']='pressure:'+','.join('.'.join(ep) for ep in pneu if groups.find(ep)==root)
    schedule=EvaluationSchedule(operations,known,labels)
    for operation in operations:
        dependencies.computation(operation)
    for target,expr in h_initial.items():
        if target not in schedule.producers:
            lines.append(f'{target}={expr};')
            dependencies.expression(target, expr)
    schedule_helpers, scheduled_lines=schedule.emit()
    lines += scheduled_lines
    for c in components:
        for name in pnames(c):
            for field,expr in [('p',p(c,name)),('m_flow',q(c,name)),('h_outflow',h(c,name))]:
                put(c,name+'.'+field,expr)

    # Mechanical force balance includes shared accelerations for rigid groups.
    feq=[]
    for edge in network.connections:
        if edge.domain=='mechanical':
            feq.append(({w(e.component,e.port+'.f'):1 for e in edge.endpoints},'0.0'))
    for c in components:
        kind=c.model_type
        if kind=='amesim_forc':
            put(c,'force',f'{num(c.direction)}*{w(c,"res.signal")}')
            feq.append(({w(c,'port_2.f'):1},'-'+w(c,'force')))
        elif kind=='amesim_f000':
            feq.append(({w(c,'port_1.f'):1},'0.0'))
        elif kind=='amesim_lstp00a':
            put(c,'gap',f'{num(c.gap0)}+{w(c,"port_2.x")}-{w(c,"port_1.x")}')
            put(c,'penetration',f'fmax(-{w(c,"gap")},0)')
            put(c,'force',f'native_contact({w(c,"penetration")},{w(c,"port_1.v")}-{w(c,"port_2.v")},{num(contact_stiffness(c))},{num(c.rcont)},{num(c.Pdis)},{int(c.discContactOption)})')
            feq += [({w(c,'port_1.f'):1},w(c,'force')),({w(c,'port_2.f'):1},'-'+w(c,'force'))]
        elif kind=='amesim_pnrp17':
            put(c,'pressure_force',f'({w(c,"port_1.p")}-101300)*{num(c.effective_area)}')
            feq += [({w(c,'port_2.f'):1,w(c,'port_5.f'):1},'-'+w(c,'pressure_force')),({w(c,'port_3.f'):1,w(c,'port_4.f'):1},w(c,'pressure_force'))]
        elif kind=='amesim_lmechn1':
            feq.append(({w(c,name+'.f'):1 for name in mnames(c)},'0.0'))
        elif kind=='amesim_mecmas21':
            v,x=w(c,'v'),w(c,'x')
            put(c,'Fvisc',f'-{num(c.rvisc)}*{v}' if c.use_friction else '0.0')
            put(c,'Ffric','0.0')
            for field,penetration,velocity,suffix in [('Fmin',num(c.xmin)+'-'+x,'-'+v,'min'),('Fmax',x+'-'+num(c.xmax),v,'max')]:
                expr=f'native_limit_force({penetration},{velocity},{num(getattr(c,"Kb"+suffix))},{num(getattr(c,"Db"+suffix))},{num(getattr(c,"Pd"+suffix))},{int(c.discContactOption)})' if int(c.stoptype)==2 else '0.0'
                put(c,field,expr)
            ref=mass_groups[groups.find(ep(c,'port_1'))][0]
            extra=f'{w(c,"Fvisc")}+{w(c,"Ffric")}+{w(c,"Fmin")}-{w(c,"Fmax")}'
            if c.use_friction: extra+=f'-{num(c.wind)}*{v}*fabs({v})'
            extra += f'+({num(c.mass * 9.80665 * math.sin(math.radians(c.theta)))})'
            feq.append(({w(c,'port_1.f'):1,w(c,'port_2.f'):1,w(ref,'a'):-c.mass},f'-({extra})'))
    unknownf=[w(c,name+'.f') for c in components for name in mnames(c)]+[w(group[0],'a') for group in mass_groups.values()]
    # Keep the same elimination/order while registering its structured bindings.
    mechanical_assignments = linear_assignments(feq, unknownf)
    for target, expression in mechanical_assignments:
        lines.append(f'{"double " if target.startswith("b") else ""}{target} = {expression};')
        dependencies.expression(target, expression)
    for index, (group, dry, threshold) in enumerate(friction_groups):
        ref = group[0]
        drive = f'friction_drive_{index}'
        expression = f'{num(sum(c.mass for c in group))}*{w(ref,"a")}'
        lines.append(f'double {drive}={expression};')
        dependencies.expression(drive, expression)
        # Amesim ideal stops take priority over dry friction. A blocked mass
        # carries its load through the stop, not through an extra friction force.
        blocked = f'friction_blocked_{index}'
        ideal = [c for c in group if int(c.stoptype) == 1]
        conditions = []
        for c in ideal:
            x, v = y(ref, 'x'), y(ref, 'v')
            conditions += [f'({x}<={num(c.xmin+1e-12*max(abs(c.xmin),1))} && {v}<=1e-12 && {drive}<=0)',
                           f'({x}>={num(c.xmax-1e-12*max(abs(c.xmax),1))} && {v}>=-1e-12 && {drive}>=0)']
        blocked_expression = ' || '.join(conditions) or '0'
        lines.append(f'double {blocked}=({blocked_expression});')
        dependencies.expression(blocked, blocked_expression)
        lines.append(f'{drive}={blocked}?0:{drive};')
        dependencies.expression(drive, f'{blocked}?0:{drive}')
        lines.append(f'if(friction_drives) friction_drives[{index}]={drive};')
        capacity = sum(c.fstick for c in dry)
        for c in dry:
            put(c, 'Ffric', f'{blocked}?0:native_dry_friction({y(ref,"v")},{drive}*{num(c.fstick/capacity)},'
                f'{num(c.fcoul)},{num(c.fstick)},{num(c.astrib)},'
                f'{int(int(c.frictionType)==2 and int(c.strib)==2)},{y(ref,"_friction_mode")})')
        lines.append(f'dy[{states[ref.name,"_friction_mode"]}]=0;')
        dependencies.expression(f'dy[{states[ref.name,"_friction_mode"]}]', '0')
    if friction_groups:
        # The first linear solve supplies the drive without dry friction. The
        # second redistributes actual friction forces through the same network.
        # No pressure/flow iteration or mutable RHS mode is introduced.
        for target, expression in mechanical_assignments:
            lines.append(f'{target}={expression};')
            dependencies.expression(target, expression)
    for c in components:
        for name in mnames(c): assigned.add(c.name+'.'+name+'.f')
        if c.model_type=='amesim_lmechn1':
            put(c,'tforce',' + '.join(w(c,name+'.f') for name in mnames(c)[:-1]))
    stops=[]
    for group in mass_groups.values():
        ref=group[0];vi=states[ref.name,'v'];xi=states[ref.name,'x']
        limits=[c for c in group if int(c.stoptype) in (1,3)]
        lines += [f'dy[{vi}]={w(ref,"a")};dy[{xi}]={y(ref,"v")};']
        dependencies.expression(f'dy[{vi}]', w(ref,'a'))
        dependencies.expression(f'dy[{xi}]', y(ref,'v'))
        if limits:
            lower=max(c.xmin for c in limits);upper=min(c.xmax for c in limits)
            if lower>upper or ref.x0<lower-1e-12 or ref.x0>upper+1e-12:
                raise NativeCapabilityError('Inconsistent discrete endstop bounds or initial position')
            restitution, thresholds = [], []
            for parameter, bound in [('xmin',lower),('xmax',upper)]:
                active=[c for c in limits if abs(getattr(c,parameter)-bound)<=1e-12*max(abs(bound),1)]
                restitution.append(0 if any(int(c.stoptype)==1 for c in active) else min(c.restcoeff for c in active))
                thresholds.append(max((c.restdvel for c in active if int(c.stoptype)==3),default=0))
            stops.append((vi,lower,upper,*restitution,*thresholds))
            lines.append(f'native_stop_motion({y(ref,"x")},{y(ref,"v")},{num(lower)},{num(upper)},&dy[{vi}],&dy[{xi}]);')
            dependencies.stop_motion(vi, xi)
        for c in group: put(c,'a',f'dy[{vi}]')

    for c in components:
        kind=c.model_type
        if kind not in GAS_TYPES and kind!='amesim_pnl00r': continue
        if kind in GAS_TYPES:
            halves=(1,2) if kind=='amesim_pnl0003' else (0,)
            center='0.0'
            if kind=='amesim_pnl0003':
                a,b=gases[c.name,1],gases[c.name,2]
                center=w(c,'dmctr')
                put(c,'dmctr',f'native_pipe_flow_context(properties,{medium(c)},{a}.p,{b}.p,{a}.p>={b}.p?{a}.T:{b}.T,{num(c.diam)},{num(c.le)},{num(c.rr)},3)')
            for half in halves:
                gas=gases[c.name,half];suffix=str(half) if half else ''
                names=['port_'+suffix] if half else pnames(c)
                mass=' + '.join(q(c,name) for name in names)
                energy=' + '.join(f'{q(c,name)}*({q(c,name)}>0?{hin(c,name)}:{gas}.h)' for name in names)
                if half:
                    sign='-' if half==1 else '+'
                    mass+=sign+center
                    energy+=f'{sign}{center}*({center}>0?{gases[c.name,1]}.h:{gases[c.name,2]}.h)'
                heat='0.0'
                if kind.startswith('amesim_pnch'):
                    heat=f'{num(c.kth*c.sth)}*({num(c.extemp)}-{gas}.T)-{gas}.p*({volume_rate[c.name]})'
                elif kind.startswith('amesim_pnl') and int(c.mode)!=1:
                    temp=gas+'.T'
                    heat=f'{num(c.kth*c.exchange_area/(2 if half else 1))}*({num(c.extemp)}-({temp}))'
                thermal = ('T' if is_polytropic(c) else 'U') + suffix
                thermal_rhs = f'{num(c.k-1)}*{gas}.T/({y(c,"m"+suffix)})*({mass})' if is_polytropic(c) else f'{energy}+({heat})'
                lines += [f'dy[{states[c.name,"m"+suffix]}]={mass};',f'dy[{states[c.name,thermal]}]={thermal_rhs};']
                dependencies.expression(f'dy[{states[c.name,"m"+suffix]}]', mass)
                dependencies.expression(f'dy[{states[c.name,thermal]}]', thermal_rhs)
        if kind.startswith('amesim_pnl'):
            diag=[]
            if kind=='amesim_pnl0002':
                gas=gases[c.name,0]
                for name in pnames(c):
                    flow=q(c,name);pp=f'({flow}>=0?fmax({p(c,name)},1):fmax({gas}.p,1))'
                    temp=f'({flow}>=0?fmax(native_temperature_ph_context(properties,{medium(c)},{pp},{hin(c,name,True)}),1):{gas}.T)'
                    diag.append((flow,pp,temp,c.le/2,0))
            elif kind=='amesim_pnl0003':
                a,b=gases[c.name,1],gases[c.name,2];flow=w(c,'dmctr')
                diag=[(flow,f'fmax(fmax({a}.p,{b}.p),1)',f'({flow}>=0?{a}.T:{b}.T)',c.le,1)]
            else:
                pa,pb=p(c,'port_1'),p(c,'port_2');pp=f'fmax(fmax({pa},{pb}),1)'
                if kind=='amesim_pnl0001':
                    gas=gases[c.name,0]
                    temp=f'({pa}>{gas}.p?fmax(native_temperature_ph_context(properties,{medium(c)},{pp},{hin(c,"port_1",True)}),1):{gas}.T)'
                else:
                    temp=f'fmax(native_temperature_ph_context(properties,{medium(c)},{pp},{pa}>={pb}?{hin(c,"port_1",True)}:{hin(c,"port_2",True)}),1)'
                diag=[(q(c,'port_1'),pp,temp,c.le,0)]
            lines.append('{ double d[4],acc[4]={0};')
            for flow,pp,temp,length,diagnostic in diag:
                lines.append(f'native_pipe_diagnostics_context(properties,{medium(c)},{flow},{pp},{temp},{num(c.diam)},{num(length)},{num(c.rr)},{diagnostic},d);')
                if len(diag)>1: lines.append('d[2]=fabs(d[2]);')
                lines.append('for(int i=0;i<4;i++) acc[i]+=d[i];')
            for i,field in enumerate(('re','cm','v','ff')):
                expr=f'acc[{i}]/{len(diag)}'
                put(c,field,f'fmin({expr},64000000)' if field=='ff' else expr)
            lines.append('}')
    for _,offsets,volumes,poly in coupled:
        thermal_sum = '+'.join(f'dy[{i+1}]*{num(v/sum(volumes))}' if poly else f'dy[{i+1}]' for i,v in zip(offsets,volumes))
        lines.append('{ double mass='+ '+'.join(f'dy[{i}]' for i in offsets)+',thermal='+thermal_sum+';')
        for i,volume in zip(offsets,volumes):
            lines.append(f'dy[{i}]=mass*{num(volume/sum(volumes))};dy[{i+1}]=thermal*{num(1 if poly else volume/sum(volumes))};')
            dependencies.assign(f'dy[{i}]', (f'dy[{j}]' for j in offsets))
            dependencies.assign(f'dy[{i+1}]', (f'dy[{j+1}]' for j in offsets))
        lines.append('}')
    missing=set(slots)-assigned
    if missing: raise NativeCapabilityError(f'Native output mapping incomplete: {sorted(missing)}')
    if not state_keys: lines.append('dy[0]=0;')
    jacobian = dependencies.build()
    np,ng,nq=max(1,len(pgroups)),max(1,gas_count),max(1,len(pneu))
    # Advanced-friction reference runs at 5 N and 7 N confirm Amesim 2404's
    # 0.1% breakaway hysteresis. Simple friction uses the exact static threshold.
    source='\n'.join(['#include "model.h"','#include <math.h>',*declarations,
        f'const NativeFriction model_frictions[{max(1,len(friction_groups))}] = {{'+(','.join(
            '{'+str(states[group[0].name,'v'])+','+str(states[group[0].name,'_friction_mode'])+','+num(threshold)+','+
            num(sum(c.fstick*(1.001 if int(c.frictionType)==2 else 1) for c in dry))+'}'
            for group,dry,threshold in friction_groups) or '{0,0,0,0}')+'};',
        f'const NativeStop model_stops[{max(1,len(stops))}] = {{'+(','.join('{'+str(s[0])+','+','.join(num(v) for v in s[1:])+'}' for s in stops) or '{0,0,0,0,0,0,0}')+'};',
        'const double model_atol[NSTATES] = {'+','.join(map(state_absolute_tolerance, state_keys or ['dummy']))+'};',
        'const char *const model_output_keys[NOUTPUTS] = {'+(','.join(json.dumps(v.key,ensure_ascii=True) for v in variables) or '""')+'};',
        *jacobian.source_lines(),
        *schedule_helpers,
        'int model_init(double *y) {',*[f'y[{i}]={num(v)};' for i,v in enumerate(initial)],*gas_initializers,'return 1;}',
        'static int model_eval_internal(double t,const double *y,double *dy,double *w,int canonical,double *friction_drives) {',
        '(void)friction_drives;',
        f'NativePropertyState property_states[{min(256,max(16,4*gas_count+2*len(components)))}];',
        'NativePropertyCache property_cache, *properties=&property_cache;',
        'native_properties_init(properties,property_states,sizeof(property_states)/sizeof(property_states[0]));',
        'NativePropertyCache *gas_properties=canonical?NULL:properties;(void)gas_properties;',
        *(['double projected[NSTATES];for(int i=0;i<NSTATES;i++) projected[i]=y[i];',*project,'y=projected;'] if project else []),
        f'double p[{np}]={{0}},h[{nq}]={{0}},q[{nq}]={{0}};NativeGas g[{ng}];',
        f'double fb[{max(1,flow_temporary_count)}]={{0}};(void)fb;',
        *([f'NativePipeCache pipe_cache[{max(1,pipe_cache_count)}]={{0}};(void)pipe_cache;'] if pneu else []),
        '(void)t;(void)y;(void)w;(void)p;(void)h;(void)q;(void)g;',*lines,
        'for(int i=0;i<NSTATES;i++) if(!isfinite(dy[i])) return 0;',
        'for(int i=0;i<NOUTPUTS;i++) if(!isfinite(w[i])) return 0;','return 1;}',
        'int model_eval(double t,const double *y,double *dy,double *w) {return model_eval_internal(t,y,dy,w,0,NULL);}',
        'int model_eval_jacobian(double t,const double *y,double *dy,double *w) {return model_eval_internal(t,y,dy,w,1,NULL);}',
        'int model_friction_drives(double t,const double *y,double *drives) {double dy[NSTATES],w[NOUTPUTS];return model_eval_internal(t,y,dy,w,0,drives);}',
        'double model_next_break(double t,double end) { double result=end;(void)t;',*breaks,'return result;}',''])
    header=f'''#ifndef GENERATED_NATIVE_MODEL_H
#define GENERATED_NATIVE_MODEL_H
#include "kernels.h"
#define NSTATES {nstates}
#define NOUTPUTS {max(1,len(variables))}
#define NSTOPS {len(stops)}
#define NFRICTIONS {len(friction_groups)}
extern const NativeFriction model_frictions[{max(1,len(friction_groups))}];
int model_friction_drives(double t,const double *y,double *drives);
extern const NativeStop model_stops[{max(1,len(stops))}];
extern const double model_atol[NSTATES];
extern const char *const model_output_keys[NOUTPUTS];
{chr(10).join(jacobian.header_lines(canonical_rhs=True))}
int model_init(double *y);
int model_eval(double t,const double *y,double *dy,double *w);
double model_next_break(double t,double end);
#endif
'''
    return NativeProgram(source,header,tuple(state_keys),variables,tuple(sorted({c.model_type for c in components})),schedule.report(),jacobian.manifest(canonical_rhs=True))