"""Direct analytic compilation. No learning, behavioral feedback or trajectories.
Run with NumPy from the directory containing public_models.
"""
from pathlib import Path
import json
import xml.etree.ElementTree as ET
import numpy as np
ROOT=Path('public_models')
PUBLIC=json.loads((ROOT/'franka_public.json').read_text())
HOME=np.array(PUBLIC['home'],dtype=float)

def numbers(s,default):
    return np.fromstring(s,sep=' ') if s is not None else np.array(default,dtype=float)
def quat_rotation(q):
    w,x,y,z=q/np.linalg.norm(q)
    return np.array([[1-2*(y*y+z*z),2*(x*y-z*w),2*(x*z+y*w)],
                     [2*(x*y+z*w),1-2*(x*x+z*z),2*(y*z-x*w)],
                     [2*(x*z-y*w),2*(y*z+x*w),1-2*(x*x+y*y)]])
def rz(q):
    c,s=np.cos(q),np.sin(q)
    return np.array([[c,-s,0],[s,c,0],[0,0,1.]])
def fixed_transform(body):
    t=np.eye(4)
    t[:3,:3]=quat_rotation(numbers(body.get('quat'),[1,0,0,0]))
    t[:3,3]=numbers(body.get('pos'),[0,0,0])
    return t
def inertia_matrix(e):
    if e.get('fullinertia') is not None:
        xx,yy,zz,xy,xz,yz=numbers(e.get('fullinertia'),[])
        return np.array([[xx,xy,xz],[xy,yy,yz],[xz,yz,zz]])
    R=quat_rotation(numbers(e.get('quat'),[1,0,0,0]))
    return R@np.diag(numbers(e.get('diaginertia'),[]))@R.T

def static_home_model(xmlroot):
    """Single static pose, inertial sums only. No dynamics or integration."""
    M=np.eye(7)*.1; g=np.zeros(7); found={}
    def visit(body,parent,axes,origins):
        nonlocal M,g
        T=parent@fixed_transform(body)
        axes,origins=dict(axes),dict(origins)
        joint=body.find('joint')
        if joint is not None:
            j=int(joint.get('name').replace('joint',''))-1
            axes[j]=T[:3,:3]@numbers(joint.get('axis'),[0,0,1])
            origins[j]=T[:3,3].copy()
            R=np.eye(4); R[:3,:3]=rz(HOME[j]); T=T@R
        def jac(p):
            V,A=np.zeros((3,7)),np.zeros((3,7))
            for j in axes:
                A[:,j]=axes[j]; V[:,j]=np.cross(axes[j],p-origins[j])
            return V,A
        inertial=body.find('inertial')
        if inertial is not None:
            mass=float(inertial.get('mass'))
            com=(T@np.r_[numbers(inertial.get('pos'),[0,0,0]),1])[:3]
            V,A=jac(com)
            I=T[:3,:3]@inertia_matrix(inertial)@T[:3,:3].T
            M+=mass*(V.T@V)+A.T@I@A
            g+=mass*9.81*V[2]
        for site in body.findall('site'):
            if site.get('name')=='push_tip':
                p=(T@np.r_[numbers(site.get('pos'),[0,0,0]),1])[:3]
                V,A=jac(p); found.update(tip=p,J=V,Jr=A,R=T[:3,:3])
        for child in body.findall('body'): visit(child,T,axes,origins)
    visit(xmlroot.find("worldbody/body[@name='link0']"),np.eye(4),{},{})
    found.update(M=M,g=g)
    return found

# Exact Fourier-polynomial matrix algebra on the q2/q4/q6 plane.
# All seven derivatives are carried to obtain physical gravity coefficients.
ZERO=(0,0,0)
def clean(p): return {k:v for k,v in p.items() if np.max(np.abs(v))>2e-12}
def padd(a,b):
    out={k:v.copy() for k,v in a.items()}
    for k,v in b.items(): out[k]=out.get(k,np.zeros_like(v))+v
    return clean(out)
def pmatmul(a,b):
    out={}
    for ka,va in a.items():
        for kb,vb in b.items():
            k=tuple(x+y for x,y in zip(ka,kb)); v=va@vb
            out[k]=out.get(k,np.zeros_like(v))+v
    return clean(out)
def rotation_poly(j):
    if j not in (1,3,5):
        R=np.eye(4); R[:3,:3]=rz(HOME[j])
        D=np.zeros((4,4)); D[:3,:3]=rz(HOME[j])@np.array([[0,-1,0],[1,0,0],[0,0,0.]])
        return {ZERO:R.astype(complex)},{ZERO:D.astype(complex)}
    axis=(1,3,5).index(j)
    kp=tuple(int(i==axis) for i in range(3)); km=tuple(-x for x in kp)
    A=np.zeros((4,4),dtype=complex); A[:2,:2]=[[.5,.5j],[-.5j,.5]]
    C=np.zeros((4,4),dtype=complex); C[2,2]=C[3,3]=1
    return {ZERO:C,kp:A,km:A.conj()},{kp:1j*A,km:-1j*A.conj()}
def planar_fourier(xmlroot):
    gravity={}; tip={}
    def visit(body,parent,parent_d):
        nonlocal gravity,tip
        F={ZERO:fixed_transform(body).astype(complex)}
        T=pmatmul(parent,F); D=[pmatmul(d,F) for d in parent_d]
        joint=body.find('joint')
        if joint is not None:
            j=int(joint.get('name').replace('joint',''))-1
            R,DR=rotation_poly(j)
            D=[padd(pmatmul(d,R),pmatmul(T,DR) if i==j else {}) for i,d in enumerate(D)]
            T=pmatmul(T,R)
        e=body.find('inertial')
        if e is not None:
            com=np.r_[numbers(e.get('pos'),[0,0,0]),1]; mg=float(e.get('mass'))*9.81
            for j,d in enumerate(D):
                for k,v in d.items():
                    if k not in gravity: gravity[k]=np.zeros(7,dtype=complex)
                    gravity[k][j]+=mg*(v[2]@com)
        for site in body.findall('site'):
            if site.get('name')=='push_tip':
                loc=np.r_[numbers(site.get('pos'),[0,0,0]),1]
                tip={k:(v@loc)[:3] for k,v in T.items()}
        for child in body.findall('body'): visit(child,T,D)
    visit(xmlroot.find("worldbody/body[@name='link0']"),{ZERO:np.eye(4,dtype=complex)},[{} for _ in range(7)])
    return clean(tip),clean(gravity)
def peval(poly,angles):
    return sum(v*np.exp(1j*np.dot(k,angles)) for k,v in poly.items()).real

class Expr:
    def __init__(self,terms=None,bias=0.):
        self.terms={} if terms is None else {k:float(v) for k,v in terms.items() if v!=0}
        self.bias=float(bias)
    def __add__(self,other):
        other=as_expr(other); terms=dict(self.terms)
        for k,v in other.terms.items(): terms[k]=terms.get(k,0.)+v
        return Expr(terms,self.bias+other.bias)
    __radd__=__add__
    def __neg__(self): return self*-1
    def __sub__(self,other): return self+-as_expr(other)
    def __rsub__(self,other): return as_expr(other)+-self
    def __mul__(self,a):
        if not np.isscalar(a): raise TypeError('Affine scalar only')
        return Expr({k:v*a for k,v in self.terms.items()},self.bias*a)
    __rmul__=__mul__
    def __truediv__(self,a): return self*(1/a)
def as_expr(x): return x if isinstance(x,Expr) else Expr(bias=x)
def var(name): return Expr({name:1})
def dot(coefs,names): return sum((float(a)*var(n) for a,n in zip(coefs,names)),Expr())

class MLP:
    def __init__(self,n):
        self.names=['o'+str(i) for i in range(n)]; self.input_dim=n; self.layers=[]
    def add_layer(self,W,b,act):
        W,b=np.asarray(W,dtype=np.float64),np.asarray(b,dtype=np.float64)
        assert W.ndim==2 and b.shape==(W.shape[0],)
        assert W.shape[1]==(self.layers[-1][0].shape[0] if self.layers else self.input_dim)
        assert np.isfinite(W).all() and np.isfinite(b).all()
        self.layers.append((W,b,act))
    def affine(self,outputs,keep=()):
        exprs={n:var(n) for n in keep}; exprs.update(outputs)
        W=np.zeros((len(exprs),len(self.names))); b=np.zeros(len(exprs))
        idx={n:i for i,n in enumerate(self.names)}
        for i,e in enumerate(exprs.values()):
            e=as_expr(e); b[i]=e.bias
            for n,a in e.terms.items(): W[i,idx[n]]+=a
        self.add_layer(W,b,'linear'); self.names=list(exprs)
    def stage(self): return Stage(self)
    def bounded_torques(self,expressions,limits):
        # Exact stable clipping, avoiding subtraction of two huge ReLUs.
        idx={n:i for i,n in enumerate(self.names)}
        W=np.zeros((len(expressions),len(self.names))); b=np.zeros(len(expressions))
        for i,(e,L) in enumerate(zip(expressions,limits)):
            e=as_expr(e); b[i]=1-e.bias/L
            for n,a in e.terms.items(): W[i,idx[n]]-=a/L
        self.add_layer(W,b,'relu')
        self.add_layer(-np.eye(len(expressions)),np.full(len(expressions),2.),'relu')
        self.add_layer(np.eye(len(expressions)),np.full(len(expressions),-1.),'linear')
        self.names=['output'+str(i) for i in range(len(expressions))]
    def save(self,path):
        arrays={'activations':np.array([a for _,_,a in self.layers],dtype='<U6')}
        for i,(W,b,_) in enumerate(self.layers): arrays['w'+str(i)]=W; arrays['b'+str(i)]=b
        assert 1<=len(self.layers)<=64
        assert max(max(W.shape) for W,_,_ in self.layers)<=4096
        assert sum(a.nbytes for a in arrays.values())<128*1024**2
        np.savez(path,**arrays)
    def evaluate(self,x):
        x=np.asarray(x,dtype=float)
        for W,b,a in self.layers:
            x=x@W.T+b
            if a=='relu': x=np.maximum(x,0)
            elif a=='tanh': x=np.tanh(x)
        return x
    def architecture(self): return [self.input_dim]+[W.shape[0] for W,_,_ in self.layers]

class Stage:
    def __init__(self,mlp): self.mlp=mlp; self.features=[]; self.cache={}
    def relu(self,e):
        e=as_expr(e); assert all(n in self.mlp.names for n in e.terms)
        key=(tuple(sorted(e.terms.items())),e.bias)
        if key not in self.cache:
            self.cache[key]='@'+str(len(self.features)); self.features.append(e)
        return var(self.cache[key])
    def clip(self,e,lo,hi): return lo+self.relu(e-lo)-self.relu(e-hi)
    def pwl(self,e,knots,values):
        # Hinge expansion of an analytic scalar linear interpolant, not fitting.
        t,y=np.asarray(knots),np.asarray(values); m=np.diff(y)/np.diff(t)
        coeff=np.r_[m[0],np.diff(m),-m[-1]]
        return float(y[0])+sum((float(a)*self.relu(e-float(k)) for a,k in zip(coeff,t) if abs(a)>1e-16),Expr())
    def product(self,a,b,bound_a,bound_b,segments=256):
        u,v=as_expr(a)/bound_a,as_expr(b)/bound_b; knots=np.linspace(-2.,2.,segments+1)
        return (bound_a*bound_b/4)*(self.pwl(u+v,knots,knots**2)-self.pwl(u-v,knots,knots**2))
    def finish(self,outputs,keep=()):
        mlp=self.mlp; n=len(mlp.names); idx={name:i for i,name in enumerate(mlp.names)}
        W=np.zeros((2*n+len(self.features),n)); b=np.zeros(W.shape[0])
        W[:n]=np.eye(n); W[n:2*n]=-np.eye(n)
        for i,e in enumerate(self.features,2*n):
            b[i]=e.bias
            for name,a in e.terms.items(): W[i,idx[name]]+=a
        mlp.add_layer(W,b,'relu')
        exprs={name:var(name) for name in keep}; exprs.update(outputs)
        V=np.zeros((len(exprs),W.shape[0])); c=np.zeros(len(exprs))
        for i,e in enumerate(exprs.values()):
            e=as_expr(e); c[i]=e.bias
            for name,a in e.terms.items():
                if name.startswith('@'): V[i,2*n+int(name[1:])]+=a
                else: V[i,idx[name]]+=a; V[i,n+idx[name]]-=a
        mlp.add_layer(V,c,'linear'); mlp.names=list(exprs)


def compile_swingup():
    net=MLP(3)
    st=net.stage()
    st.finish({'omega':st.clip(8*var('o2'),-12.,12.),
               'absomega':st.relu(8*var('o2'))+st.relu(-8*var('o2'))},keep=['o0','o1','o2'])
    st=net.stage(); knots=np.linspace(-12,12,513)
    st.finish({'deficit':4.905*(1+var('o1'))-st.pwl(var('omega'),knots,knots**2)/6,
               'position_gate':st.clip((-var('o1')-.85)/.08,0.,1.),
               'speed_gate':st.clip((2.1-var('absomega'))/.7,0.,1.),
               'seed_position':st.clip(2*var('o1')-1,0.,1.),
               'seed_speed':st.clip(1-2*var('absomega'),0.,1.)},keep=['o0','o1','omega'])
    st=net.stage()
    st.finish({'energy_product':st.product(var('deficit'),var('omega'),24.,12.,512),
               'gate':var('position_gate')-st.relu(var('position_gate')-var('speed_gate')),
               'seed':var('seed_position')-st.relu(var('seed_position')-var('seed_speed'))},keep=['o0','omega'])
    st=net.stage()
    st.finish({'pump':st.clip(.75*var('energy_product')+.05*var('omega')+.10*var('seed'),-2.,2.),
               'balance':st.clip(12.*var('o0')-3.*var('omega'),-2.,2.)},keep=['gate'])
    st=net.stage()
    st.finish({'torque':var('pump')+st.product(var('gate'),var('balance')-var('pump'),1.,4.,512)})
    net.bounded_torques([var('torque')],[2.])
    return net

def positive_key(k): return next((x for x in k if x),0)>0
def fourier_expr(poly,index):
    e=Expr(bias=poly[ZERO][index].real if ZERO in poly else 0.)
    for k,v in poly.items():
        if positive_key(k):
            tag='_'.join(map(str,k))
            e+=2*v[index].real*var('cos_'+tag)-2*v[index].imag*var('sin_'+tag)
    return e

def compile_franka(xmlroot,static,tip_poly,gravity_poly):
    net=MLP(29)
    qnames=['o'+str(i) for i in range(7)]; vnames=['o'+str(i) for i in range(7,14)]
    keys=sorted({k for p in (tip_poly,gravity_poly) for k in p if positive_key(k)})
    joint_bounds=np.array([[-1.7628,1.7628],[-3.0718,-.0698],[-.0175,3.7525]])
    st=net.stage(); trig={}
    for k in keys:
        a=dot(k,[qnames[i] for i in (1,3,5)])+float(np.dot(k,HOME[[1,3,5]]))
        lo=sum(min(c*x,c*y) for c,(x,y) in zip(k,joint_bounds))-.10
        hi=sum(max(c*x,c*y) for c,(x,y) in zip(k,joint_bounds))+.10
        knots=np.linspace(lo,hi,int(np.ceil((hi-lo)/.02))+1)
        tag='_'.join(map(str,k))
        trig['sin_'+tag]=st.pwl(a,knots,np.sin(knots))
        trig['cos_'+tag]=st.pwl(a,knots,np.cos(knots))
    st.finish(trig,keep=qnames+vnames+['o15','o27','o28'])
    out={'px':fourier_expr(tip_poly,0),'pz':fourier_expr(tip_poly,2)}
    for j in range(7): out['g'+str(j)]=fourier_expr(gravity_poly,j)
    for m,j in enumerate((1,3,5)):
        dpoly={k:1j*k[m]*v for k,v in tip_poly.items()}
        out['jx'+str(j)]=fourier_expr(dpoly,0); out['jz'+str(j)]=fourier_expr(dpoly,2)
    net.affine(out,keep=qnames+vnames+['o15','o27','o28'])
    st=net.stage()
    out={'fx':st.clip(220.*(float(PUBLIC['tip_home'][0])+.2*var('o27')-.038-var('px')),-4.,4.),
         'fz':st.clip(800.*(float(PUBLIC['tip_home'][2])-var('pz')),-30.,30.),
         'pitch_moment':st.clip(-35.*(var('o1')-var('o3')-var('o5')),-8.,8.),
         'lateral_goal':st.clip(.2*(1.15*var('o15')-.15*var('o28')),-.04,.04)}
    inv_grid=np.linspace(.28,.85,115)
    out['inverse_radius']=st.pwl(var('px'),inv_grid,1/inv_grid)
    keep=qnames+vnames+['g'+str(j) for j in range(7)]+['jx'+str(j) for j in (1,3,5)]+['jz'+str(j) for j in (1,3,5)]
    st.finish(out,keep=keep)
    st=net.stage(); out={}
    for j in (1,3,5):
        out['tx'+str(j)]=st.product(var('jx'+str(j)),var('fx'),1.1,4.,256)
        out['tz'+str(j)]=st.product(var('jz'+str(j)),var('fz'),1.1,30.,256)
    out['yaw_goal']=st.product(var('lateral_goal'),var('inverse_radius'),.04,3.6,256)
    st.finish(out,keep=qnames+vnames+['g'+str(j) for j in range(7)]+['pitch_moment'])
    # Modal damping from physical static inertia and designed stiffness only.
    J,Jr,M=static['J'],static['Jr'],static['M']; K=np.zeros((7,7))
    for j,k in ((0,70.),(2,90.),(4,30.)): K[j,j]+=k
    yaw_vector=np.zeros(7); yaw_vector[6]=1.; yaw_vector[0]=-1.
    K+=20.*np.outer(yaw_vector,yaw_vector)
    even=[1,3,5]; Jplanar=np.array([J[0,even],J[2,even],Jr[1,even]])
    K[np.ix_(even,even)]+=Jplanar.T@np.diag([220.,800.,35.])@Jplanar
    eig,U=np.linalg.eigh(M); Mh=(U*np.sqrt(eig))@U.T; Mi=(U*(1/np.sqrt(eig)))@U.T
    freq2,V=np.linalg.eigh(Mi@K@Mi)
    Dtotal=1.8*Mh@((V*np.sqrt(np.maximum(freq2,0)))@V.T)@Mh
    d,Q=np.linalg.eigh(Dtotal-np.eye(7)); Dfeedback=(Q*np.maximum(d,.2))@Q.T
    G=np.array(PUBLIC['gravity_gradient']); odd=[0,2,4,6]; torques=[]
    for j in range(7):
        t=var('g'+str(j))+dot(G[j,odd],[qnames[i] for i in odd])-3.*dot(Dfeedback[j],vnames)
        if j in even:
            t+=var('tx'+str(j))+var('tz'+str(j))
            t+=(1. if j==1 else -1.)*var('pitch_moment')
        elif j==0: t+=70.*(var('yaw_goal')-var('o0'))+20.*(var('o6')-var('o0'))
        elif j==2: t-=90.*var('o2')
        elif j==4: t-=30.*var('o4')
        elif j==6: t-=20.*(var('o6')-var('o0'))
        torques.append(t)
    net.bounded_torques(torques,PUBLIC['torque_limits'])
    design={'fourier_keys':keys,'home_mass_eigenvalues':eig.tolist(),
            'undamped_home_modal_frequencies':np.sqrt(freq2).tolist(),
            'feedback_damping':Dfeedback.tolist(),
            'max_sine_cosine_interpolation_error_bound':.02**2/8}
    return net,design

def validate_artifact(path,inputs):
    # Independent synthetic input checks; no transitions, times or rewards.
    with np.load(path,allow_pickle=False) as data:
        acts=data['activations']; assert acts.dtype.kind=='U'
        assert set(data.files)=={'activations'}|{p+str(i) for i in range(len(acts)) for p in ('w','b')}
        x=inputs.copy()
        for i,a in enumerate(acts):
            W,b=data['w'+str(i)],data['b'+str(i)]
            assert W.ndim==2 and b.shape==(W.shape[0],)
            assert W.dtype.kind=='f' and b.dtype.kind=='f'
            assert np.isfinite(W).all() and np.isfinite(b).all()
            x=x@W.T+b
            if a=='relu': x=np.maximum(x,0)
            elif a=='tanh': x=np.tanh(x)
            else: assert a=='linear'
        assert np.isfinite(x).all()
        assert np.all(x>=-1.) and np.all(x<=1.)
    return {'layers':len(acts),'bytes':Path(path).stat().st_size,
            'output_min':float(x.min()),'output_max':float(x.max())}

def main():
    ET.parse(ROOT/'swingup.xml')
    xmlroot=ET.parse(ROOT/'franka_push.xml').getroot(); static=static_home_model(xmlroot)
    np.testing.assert_allclose(static['tip'],PUBLIC['tip_home'],atol=2e-12)
    np.testing.assert_allclose(static['J'],PUBLIC['jacobian'],atol=2e-12)
    np.testing.assert_allclose(static['Jr'],PUBLIC['rotation_jacobian'],atol=2e-12)
    np.testing.assert_allclose(static['g'],PUBLIC['gravity_home'],atol=2e-10)
    tip_poly,g_poly=planar_fourier(xmlroot)
    np.testing.assert_allclose(peval(tip_poly,HOME[[1,3,5]]),static['tip'],atol=2e-12)
    np.testing.assert_allclose(peval(g_poly,HOME[[1,3,5]]),static['g'],atol=2e-10)
    swing=compile_swingup(); franka,design=compile_franka(xmlroot,static,tip_poly,g_poly)
    swing.save('swingup.npz'); franka.save('franka_push.npz')
    rng=np.random.default_rng(20261105)
    q=rng.uniform(-np.pi,np.pi,128); v=rng.uniform(-30,30,128)
    si=np.column_stack([np.sin(q),np.cos(q),v/8])
    si=np.vstack([si,[[0,1,0],[0,-1,0],[0,-1,1e4],[0,1,-1e4]]])
    fi=rng.uniform(-1,1,(128,29))
    for row in fi: row[17:21]/=np.linalg.norm(row[17:21])
    fi=np.vstack([fi,np.zeros((1,29)),np.full((1,29),1e4),np.full((1,29),-1e4)])
    checks={'swingup':validate_artifact('swingup.npz',si),'franka_push':validate_artifact('franka_push.npz',fi)}
    for name,net in [('swingup',swing),('franka_push',franka)]:
        checks[name]['architecture']=net.architecture()
        checks[name]['parameters']=sum(W.size+b.size for W,b,_ in net.layers)
        checks[name]['activations']=[a for _,_,a in net.layers]
    report={'checks':checks,'static_design':design,'trajectory_steps':0,'optimization_steps':0,'behavioral_evaluations':0}
    Path('construction_checks.json').write_text(json.dumps(report,indent=2)+'\n')
    rationale="""ZERO-BEHAVIORAL-FEEDBACK DIRECT CONSTRUCTION

Both artifacts are genuine float64 affine/ReLU feed-forward MLPs. There are no
custom layers, hidden states, runtime controllers, simulator calls or downloads.
All model arithmetic below is encoded in numerical matrices and biases.
The final three layers compute -1 + ReLU(2 - ReLU(1 - torque/limit)), giving an
exact bounded [-1,1] output without attenuating unsaturated gravity support.

SWINGUP
Energy deficit: 4.905*(1+cos(q)) - omega^2/6. Pump torque is the saturated
0.75*deficit*omega + 0.05*omega, plus a 0.10 Nm low-speed/bottom startup seed.
A continuous position-and-speed gate blends this with saturated upright PD:
12*sin(q) - 3*omega. The gate is zero for cos(q)>=-0.85 or |omega|>=2.1,
and one for cos(q)<=-0.93 and |omega|<=1.4. Internal energy speed is clipped
at 12 rad/s, with braking outside that range. Squares and products use analytic
piecewise-linear scalar identities, not fitted controllers.

FRANKA PUSH
The XML is used for exact algebraic Fourier kinematics and gravity on the
planar q2/q4/q6 slice, including gravity torques on all seven joints. Sin/cos
are analytic linear interpolants with angle spacing at most 0.02 rad. Odd
joint departures receive the public first-order gravity correction. The
planar Jacobian is differentiated algebraically. A Cartesian spring holds
the home tip height (800 N/m) and home pitch (35 Nm/rad), while advancing the
tip toward target_x - 0.038 m with a 220 N/m spring capped at 4 N. This offset
accounts approximately for the cube half-width and finger front surface.
Joint 1 centers the pusher near the cube y, biased slightly to the opposite
side of the target's lateral error; joints 3 and 5 stay near home, while joint
7 cancels base yaw. The lateral command is capped at +/-0.04 m. Damping is an
analytic 0.9-ratio modal construction from the static XML mass matrix and
home stiffness, accounting for existing joint damping. All online products,
kinematics, gravity support, damping and saturation are inside the MLP.

VALIDATION AND LIMITATIONS
Only static model identities and network schema/finite-output/bound checks
on independent synthetic observations were performed. No time integration,
trajectories, reward queries, optimization, training or teacher fitting.
Pendulum capture is not formally certified. Franka uses a planar model with
first-order off-plane support; Coriolis terms are omitted. Contact, cube yaw,
lateral drift, finite-control-rate transients, collisions, and sample-and-hold
stability have not been behaviorally evaluated. The small initial workspace
and target footprint motivate this design, but neither task's success is
claimed without evaluation.

ARCHITECTURES (input then affine output widths)
"""
    for name in ['swingup','franka_push']:
        c=checks[name]; rationale+=name+': '+str(c['architecture'])+'\n'
        rationale+='activations: '+', '.join(c['activations'])+'\n'
        rationale+=str(c['parameters'])+' parameters; '+str(c['bytes'])+' file bytes.\n'
    Path('rationale.txt').write_text(rationale)
    print(json.dumps(report,indent=2))

if __name__=='__main__': main()