#!/usr/bin/env python3
"""Reviewed independent small-strain nonlinear axial-network reference.

Usage: python3 Lingo-Nonlinear28-worker.py exact-job.json receipt.json
Requires NumPy/SciPy; the OpenSees profile additionally needs OpenSeesPy. Does not run
commands, contact a service, control hardware,
or interpret arbitrary Python from a job. Source and complete job bytes are checked.
"""
from __future__ import annotations
import sys, json, math, hashlib, pathlib, datetime, os, tempfile
VERSION='lingo-nonlinear-reference28-v1'
LIMIT=1048576

def digest(text):
    return hashlib.sha256(text.encode('utf-8')).hexdigest()
def object_pairs(pairs):
    out={}
    for key,value in pairs:
        if key in out or key in ('__proto__','prototype','constructor'):
            raise ValueError('Duplicate or reserved JSON field')
        out[key]=value
    return out
def number(text):
    value=float(text)
    if not math.isfinite(value) or value==0 and any(c in '123456789' for c in text.split('e')[0].split('E')[0]):
        raise ValueError('Nonfinite or underflowed JSON number')
    return value
def integer(text):
    value=int(text)
    if abs(value)>9007199254740991:
        raise ValueError('Integer exceeds exact cross-runtime range')
    return value
def parse(text):
    return json.loads(text,object_pairs_hook=object_pairs,parse_float=number,parse_int=integer,parse_constant=lambda _: (_ for _ in ()).throw(ValueError('Nonfinite JSON')))
def dump(value):
    return json.dumps(value,sort_keys=True,separators=(',',':'),ensure_ascii=False,allow_nan=False)
def exact(v,keys):
    if not isinstance(v,dict) or set(v)!=set(keys):raise ValueError('Unsupported or missing job fields')
def finite(value,lo=-1e30,hi=1e30):
    if isinstance(value,bool) or not isinstance(value,(int,float)) or not math.isfinite(value) or not lo<=value<=hi:raise ValueError('Number outside supported range')
    return float(value)
def validate(model,cases,settings):
    exact(model,['format','units','nodes','members','drives','unsupportedFeatures'])
    if model['format']!='LingoAxialNetwork28' or model['units']!={'length':'mm','force':'N','stress':'MPa'} or model['unsupportedFeatures']:raise ValueError('Unsupported formulation, units or features')
    if not 2<=len(model['nodes'])<=24 or not 1<=len(model['members'])<=48 or not 1<=len(model['drives'])<=12 or not 1<=len(cases)<=8:raise ValueError('Model/case limits')
    import re
    valid=lambda x:isinstance(x,str) and re.fullmatch(r'[A-Za-z][A-Za-z0-9_-]{0,63}',x)
    ids=[]
    for n in model['nodes']:
        exact(n,['id','xyzMm','fixed'])
        if not valid(n['id']) or n['id'] in ids or len(n['xyzMm'])!=3 or len(n['fixed'])!=3 or any(type(v)!=bool for v in n['fixed']):raise ValueError('Node identity/constraints')
        for x in n['xyzMm']:finite(x,-1e6,1e6)
        ids.append(n['id'])
    mids=[];touched=set();adj={n:[] for n in ids}
    for m in model['members']:
        exact(m,['id','start','end','areaMm2','youngMPa','yieldMPa','hardeningMPa','source'])
        if not valid(m['id']) or m['id'] in mids or m['start'] not in ids or m['end'] not in ids or m['start']==m['end'] or not isinstance(m['source'],str) or not 1<=len(m['source'])<=500:raise ValueError('Member source/identity')
        for key,lo,hi in [('areaMm2',1e-9,1e9),('youngMPa',1e-9,1e6),('yieldMPa',1e-9,1e5),('hardeningMPa',1e-6,1e6)]:finite(m[key],lo,hi)
        a=model['nodes'][ids.index(m['start'])]['xyzMm'];b=model['nodes'][ids.index(m['end'])]['xyzMm']
        if math.dist(a,b)<.001:raise ValueError('Degenerate member')
        mids.append(m['id']);touched.update([m['start'],m['end']]);adj[m['start']].append(m['end']);adj[m['end']].append(m['start'])
    visited=set();todo=[ids[0]]
    while todo:
        n=todo.pop()
        if n in visited:continue
        visited.add(n);todo.extend(adj[n])
    if touched!=set(ids) or visited!=set(ids):raise ValueError('Disconnected nodes/components')
    drive_keys=set();steps=None
    for d in model['drives']:
        exact(d,['node','dof','valuesMm'])
        if d['node'] not in ids or type(d['dof'])!=int or d['dof'] not in (0,1,2) or model['nodes'][ids.index(d['node'])]['fixed'][d['dof']] or (d['node'],d['dof']) in drive_keys:raise ValueError('Invalid drive')
        if not 2<=len(d['valuesMm'])<=129 or d['valuesMm'][0]!=0 or steps is not None and steps!=len(d['valuesMm']):raise ValueError('Incomplete displacement path')
        steps=len(d['valuesMm']);drive_keys.add((d['node'],d['dof']))
        for v in d['valuesMm']:finite(v,-10000,10000)
    cids=set()
    for c in cases:
        exact(c,['id','amplitude','youngScale','yieldScale'])
        if not valid(c['id']) or c['id'] in cids:raise ValueError('Case identity')
        cids.add(c['id'])
        for k in ('amplitude','youngScale','yieldScale'):finite(c[k],.01,10)
    exact(settings,['forceToleranceN','maxIterations','relativeTolerance'])
    finite(settings['forceToleranceN'],1e-8,1);finite(settings['relativeTolerance'],1e-9,.001)
    if type(settings['maxIterations'])!=int or not 2<=settings['maxIterations']<=80:raise ValueError('Iteration bound')
    if len(cases)*steps*(len(ids)*12+len(mids)*10+30)>35000:raise ValueError('Full history output exceeds planning budget')

def material(strain,old,E,fy,H):
    trial=E*(strain-old[0]);relative=trial-old[1];f=abs(relative)-fy
    if f<=1e-10*max(1,fy):return trial,E,old.copy()
    sign=1. if relative>=0 else -1.;gamma=f/(E+H)
    return trial-E*gamma*sign,E*H/(E+H),[old[0]+gamma*sign,old[1]+H*gamma*sign,old[2]+gamma]

def run_case(model,c,settings):
    import numpy as np
    from scipy.linalg import solve as scipy_solve
    ids={n['id']:i for i,n in enumerate(model['nodes'])};size=len(ids)*3
    members=[]
    for e in model['members']:
        a=ids[e['start']];b=ids[e['end']];delta=np.array(model['nodes'][b]['xyzMm'],dtype=float)-model['nodes'][a]['xyzMm'];L=float(np.linalg.norm(delta));direction=delta/L
        members.append((e,L,np.concatenate([-direction/L,direction/L]),np.array([3*a,3*a+1,3*a+2,3*b,3*b+1,3*b+2])))
    controlled={3*ids[d['node']]+d['dof'] for d in model['drives']}
    free=np.array([3*i+j for i,n in enumerate(model['nodes']) for j,fixed in enumerate(n['fixed']) if not fixed and 3*i+j not in controlled],dtype=int)
    states=[[0.,0.,0.] for _ in members];u=np.zeros(size);last=np.zeros(size);work=0.;steps=[]
    def evaluate(v):
        K=np.zeros((size,size));internal=np.zeros(size);values=[]
        for i,(e,L,B,index) in enumerate(members):
            strain=float(B@v[index]);E=e['youngMPa']*c['youngScale'];fy=e['yieldMPa']*c['yieldScale'];stress,tangent,state=material(strain,states[i],E,fy,e['hardeningMPa']);volume=e['areaMm2']*L
            internal[index]+=B*stress*volume;K[np.ix_(index,index)]+=np.outer(B,B)*tangent*volume
            values.append((e,L,strain,E,fy,stress,state))
        return K,internal,values
    try:
        for step in range(len(model['drives'][0]['valuesMm'])):
            previous=u.copy()
            for d in model['drives']:u[3*ids[d['node']]+d['dof']]=d['valuesMm'][step]*c['amplitude']
            K,internal,values=evaluate(u);iteration=0
            if step==0 and free.size and np.linalg.cond(K[np.ix_(free,free)])>1e13:raise ValueError('The free tangent is singular; inspect supports and axial mechanisms.')
            while True:
                residual=float(np.max(np.abs(internal[free]),initial=0));scale=max(1.,float(np.max(np.abs(internal),initial=0)))
                if residual<=settings['forceToleranceN']+settings['relativeTolerance']*scale:break
                if iteration>=settings['maxIterations']:raise ValueError('Nonlinear equilibrium did not converge; accepted prior frames remain inspectable.')
                iteration+=1
                A=K[np.ix_(free,free)]
                if not np.isfinite(A).all() or np.linalg.cond(A)>1e13:raise ValueError('The free tangent is singular; inspect supports and axial mechanisms.')
                delta=scipy_solve(A,-internal[free],assume_a='sym',check_finite=True)
                accepted=False
                for cut in range(12):
                    candidate=u.copy();factor=2.**-cut;candidate[free]+=delta*factor;nK,nF,nV=evaluate(candidate);r=float(np.max(np.abs(nF[free]),initial=0))
                    if r<residual*(1-1e-5*factor) or r<=settings['forceToleranceN']:
                        u=candidate;K,internal,values=nK,nF,nV;accepted=True;break
                if not accepted:raise ValueError('Line search failed to reduce the equilibrium residual.')
            if max(abs(v[2]) for v in values)>.1:raise ValueError('Strain exceeds the supported small-strain network bound of 0.1.')
            elastic=sum(.5*stress**2/E*e['areaMm2']*L for e,L,strain,E,fy,stress,state in values)
            hardening=sum(.5*state[1]**2/e['hardeningMPa']*e['areaMm2']*L for e,L,strain,E,fy,stress,state in values)
            plastic=sum(fy*state[2]*e['areaMm2']*L for e,L,strain,E,fy,stress,state in values)
            work+=float(.5*(internal+last)@(u-previous));last=internal.copy();states=[v[6] for v in values]
            steps.append({'index':step,'nodes':[{'id':n['id'],'uMm':u[3*i:3*i+3].tolist(),'reactionN':internal[3*i:3*i+3].tolist()} for i,n in enumerate(model['nodes'])],
                'members':[{'id':e['id'],'strain':strain,'stressMPa':stress,'axialForceN':stress*e['areaMm2'],'plasticStrain':state[0],'backStressMPa':state[1],'plasticDissipationNmm':fy*state[2]*e['areaMm2']*L} for e,L,strain,E,fy,stress,state in values],
                'iterations':iteration,'maxFreeResidualN':float(np.max(np.abs(internal[free]),initial=0)),'elasticEnergyNmm':elastic,'hardeningEnergyNmm':hardening,'plasticDissipationNmm':plastic,'externalWorkTrapezoidNmm':work,'energyQuadratureDefectNmm':work-elastic-hardening-plastic})
        return {'id':c['id'],'status':'complete','reason':'Every supplied displacement frame converged with committed bilinear kinematic states. External work uses the declared path trapezoid; its energy quadrature defect is reported, not hidden.','steps':steps}
    except Exception as error:
        return {'id':c['id'],'status':'failed','reason':str(error),'steps':steps}

def run_opensees_case(model,c,settings,log_entries):
    import numpy as np
    import openseespy.opensees as ops
    tags={n['id']:i+1 for i,n in enumerate(model['nodes'])};indices={n['id']:i for i,n in enumerate(model['nodes'])};size=len(tags)*3
    controlled={3*indices[d['node']]+d['dof'] for d in model['drives']};free=[3*i+j for i,n in enumerate(model['nodes']) for j,fixed in enumerate(n['fixed']) if not fixed and 3*i+j not in controlled]
    if not free:return {'id':c['id'],'status':'failed','reason':'The installed OpenSees network adapter requires at least one free equilibrium degree of freedom. Fully prescribed paths can use the SciPy profile.','steps':[]}
    initial=np.zeros((size,size))
    for e in model['members']:
        a=indices[e['start']];b=indices[e['end']];delta=np.array(model['nodes'][b]['xyzMm'])-model['nodes'][a]['xyzMm'];L=float(np.linalg.norm(delta));direction=delta/L;B=np.concatenate([-direction/L,direction/L]);index=[3*a,3*a+1,3*a+2,3*b,3*b+1,3*b+2];initial[np.ix_(index,index)]+=np.outer(B,B)*e['youngMPa']*c['youngScale']*e['areaMm2']*L
    if np.linalg.cond(initial[np.ix_(free,free)])>1e13:return {'id':c['id'],'status':'failed','reason':'The free tangent is singular; inspect supports and axial mechanisms.','steps':[]}
    steps=[];old_plastic=[0. for _ in model['members']];accumulated=[0. for _ in model['members']];last_u=np.zeros(size);last_reaction=np.zeros(size);work=0.
    try:
        ops.wipe();ops.model('basic','-ndm',3,'-ndf',3)
        for n in model['nodes']:
            ops.node(tags[n['id']],*n['xyzMm']);ops.fix(tags[n['id']],*[int(v) for v in n['fixed']])
        for i,e in enumerate(model['members'],1):
            E=e['youngMPa']*c['youngScale'];H=e['hardeningMPa'];ops.uniaxialMaterial('Steel01',i,e['yieldMPa']*c['yieldScale'],E,H/(E+H));ops.element('Truss',i,tags[e['start']],tags[e['end']],e['areaMm2'],i)
        for i,d in enumerate(model['drives'],1):
            ops.timeSeries('Path',i,'-time',*list(range(len(d['valuesMm']))),'-values',*[v*c['amplitude'] for v in d['valuesMm']],'-useLast');ops.pattern('Plain',i,i);ops.sp(tags[d['node']],d['dof']+1,1.)
        ops.constraints('Transformation');ops.numberer('RCM');ops.system('BandGeneral');ops.test('NormUnbalance',settings['forceToleranceN'],settings['maxIterations']);ops.algorithm('NewtonLineSearch');ops.integrator('LoadControl',1.);ops.analysis('Static')
        for step in range(len(model['drives'][0]['valuesMm'])):
            if step==0:
                log_entries.append('case '+c['id']+' frame 0: initial zero state; no analyze call')
            else:
                increment=1.;cutbacks=0;accepted_substeps=0
                while float(ops.getTime())<step-1e-10:
                    previous_time=float(ops.getTime());dt=min(increment,step-previous_time);ops.integrator('LoadControl',dt);code=int(ops.analyze(1))
                    if code!=0:
                        if abs(float(ops.getTime())-previous_time)>1e-10:raise ValueError('Failed OpenSees step did not restore its committed time; response withheld')
                        increment=dt/2;cutbacks+=1
                        if increment<1/4096 or cutbacks>24:raise ValueError('Installed OpenSees adaptive equilibrium failed at original frame '+str(step))
                    else:
                        accepted_substeps+=1
                        if accepted_substeps>4096:raise ValueError('Installed OpenSees adaptive substep bound exceeded')
                        increment=min(1.,increment*2)
                log_entries.append('case '+c['id']+' frame '+str(step)+' analyze(1)=0; accepted '+str(accepted_substeps)+' substeps; '+str(cutbacks)+' cutbacks')
            ops.reactions();u=np.array([v for n in model['nodes'] for v in ops.nodeDisp(tags[n['id']])]);reaction=np.array([v for n in model['nodes'] for v in ops.nodeReaction(tags[n['id']])]);rows=[];elastic=hardening=plastic=0.
            for i,e in enumerate(model['members']):
                values=ops.eleResponse(i+1,'material','stressStrain')
                if len(values)!=2:raise ValueError('Installed OpenSees material response omitted stress or strain')
                stress,strain=[float(x) for x in values]
                if abs(strain)>.1:raise ValueError('Strain exceeds the supported small-strain network bound of 0.1.')
                E=e['youngMPa']*c['youngScale'];H=e['hardeningMPa'];fy=e['yieldMPa']*c['yieldScale'];ep=strain-stress/E;back=H*ep;accumulated[i]+=abs(ep-old_plastic[i]);old_plastic[i]=ep;L=math.dist(model['nodes'][indices[e['start']]]['xyzMm'],model['nodes'][indices[e['end']]]['xyzMm']);volume=e['areaMm2']*L;dissipation=fy*accumulated[i]*volume
                elastic+=.5*stress**2/E*volume;hardening+=.5*back**2/H*volume;plastic+=dissipation
                rows.append({'id':e['id'],'strain':strain,'stressMPa':stress,'axialForceN':stress*e['areaMm2'],'plasticStrain':ep,'backStressMPa':back,'plasticDissipationNmm':dissipation})
            work+=float(.5*(reaction+last_reaction)@(u-last_u));last_u=u.copy();last_reaction=reaction.copy();residual=max([abs(float(reaction[i])) for i in free]+[0.])
            if residual>settings['forceToleranceN']+settings['relativeTolerance']*max(1.,float(np.max(np.abs(reaction)))):raise ValueError('Installed OpenSees result violates the declared free force tolerance')
            steps.append({'index':step,'nodes':[{'id':n['id'],'uMm':u[3*i:3*i+3].tolist(),'reactionN':reaction[3*i:3*i+3].tolist()} for i,n in enumerate(model['nodes'])],'members':rows,'iterations':0 if step==0 else int(ops.testIter()),'maxFreeResidualN':residual,'elasticEnergyNmm':elastic,'hardeningEnergyNmm':hardening,'plasticDissipationNmm':plastic,'externalWorkTrapezoidNmm':work,'energyQuadratureDefectNmm':work-elastic-hardening-plastic})
        return {'id':c['id'],'status':'complete','reason':'Every supplied displacement frame converged with committed bilinear kinematic states. External work uses the declared path trapezoid; its energy quadrature defect is reported, not hidden.','steps':steps}
    except Exception as error:
        return {'id':c['id'],'status':'failed','reason':str(error),'steps':steps}
    finally:
        ops.wipe()

def main():
    if len(sys.argv)==2 and sys.argv[1] in ('--help','-h'):
        print(__doc__);return 0
    if len(sys.argv)!=3:raise ValueError('Use exact job.json and distinct receipt.json paths')
    source=pathlib.Path(sys.argv[1]);destination=pathlib.Path(sys.argv[2])
    if source.resolve()==destination.resolve() or destination.exists():raise ValueError('Output must be a new path, distinct from the original job')
    data=source.read_bytes()
    if len(data)>LIMIT:raise ValueError('Job exceeds 1 MiB')
    job=parse(data.decode('utf-8'));exact(job,['format','version','bodyJSON','bodySHA256'])
    if job['format']!='LingoNonlinearJob28' or job['version']!=VERSION or digest(job['bodyJSON'])!=job['bodySHA256']:raise ValueError('Job version or exact source hash changed')
    body=parse(job['bodyJSON']);exact(body,['format','version','adapter','modelJSON','casesJSON','settingsJSON','sourceRefs','workerSourceSHA256','consent'])
    worker_hash=hashlib.sha256(pathlib.Path(__file__).read_bytes()).hexdigest()
    if body['format']!='LingoNonlinearBody28' or body['version']!=1 or body['adapter'] not in ('scipy-reference','openseespy') or body['consent'] is not True or body['workerSourceSHA256']!=worker_hash:raise ValueError('Wrong reviewed worker or missing exact consent')
    if not isinstance(body['sourceRefs'],list) or len(body['sourceRefs'])>32:raise ValueError('Source reference bound')
    import re
    seen=set()
    for ref in body['sourceRefs']:
        exact(ref,['id','sha256'])
        if not re.fullmatch(r'[a-f0-9]{8}-[a-f0-9]{4}-[1-8][a-f0-9]{3}-[89ab][a-f0-9]{3}-[a-f0-9]{12}',ref['id']) or not re.fullmatch(r'[a-f0-9]{64}',ref['sha256']) or ref['id'] in seen:raise ValueError('Malformed exact source references')
        seen.add(ref['id'])
    model=parse(body['modelJSON']);cases=parse(body['casesJSON']);settings=parse(body['settingsJSON']);validate(model,cases,settings)
    import numpy,scipy
    started=datetime.datetime.now(datetime.timezone.utc).isoformat();log_entries=[]
    if body['adapter']=='openseespy':
        import openseespy.opensees as ops
        engine='OpenSeesPy '+str(ops.version())+' Truss/Steel01 axial network'
        output={'cases':[run_opensees_case(model,c,settings,log_entries) for c in cases]}
    else:
        engine='SciPy '+scipy.__version__+' independent axial network'
        output={'cases':[run_case(model,c,settings) for c in cases]}
    out=dump(output)
    if len(out.encode())>1000000:raise ValueError('Complete output exceeds receipt bound; reduce selected model/cases')
    status='complete' if all(c['status']=='complete' for c in output['cases']) else 'failed'
    log='Reviewed '+engine+'; NumPy '+numpy.__version__+'; SciPy '+scipy.__version__+'; '+str(len(cases))+' cases; '+str(sum(len(c['steps']) for c in output['cases']))+' accepted frames; no external network calls.'+'\n'+'\n'.join(log_entries)
    receipt={'format':'LingoNonlinearReceipt28','version':VERSION,'jobSHA256':job['bodySHA256'],'workerSourceSHA256':worker_hash,'engine':engine,'startedAt':started,'finishedAt':datetime.datetime.now(datetime.timezone.utc).isoformat(),'status':status,'outputJSON':out,'outputSHA256':digest(out),'log':log,'logSHA256':digest(log)}
    with tempfile.NamedTemporaryFile('w',encoding='utf-8',dir=str(destination.parent),delete=False) as f:
        temporary=f.name;f.write(dump(receipt));f.flush();os.fsync(f.fileno())
    try:
        os.link(temporary,destination)
    finally:
        os.unlink(temporary)
    print('Wrote exact-source '+status+' nonlinear receipt: '+str(destination))
    return 0 if status=='complete' else 2
if __name__=='__main__':
    try:sys.exit(main())
    except Exception as error:
        print('Nonlinear reference refused: '+str(error),file=sys.stderr);sys.exit(1)