Manufacturing-MoE / src /hole_physics.py

Sync from rasitcesur/Manufacturing-MoE via hub-sync

547ac30 verified 23 days ago

8.12 kB

"""Canonical hole-cycle physics approximation, independent of training targets.
One Hole Processing Expert for G17/Z and G19/X; physical axis retained only for
machine-energy accounting. New synthetic metadata is NOT suitable for G-code-only
production deployments unless corresponding tool/stock data are available.
This is NOT a LinuxCNC runtime simulator.
"""
from __future__ import annotations
import math,re,json
import numpy as np
from hole_constants import MAT,TOOL,move_time
from hole_energy_hardware import components_by_row
RE=re.compile(r'([A-Z])\s*([-+]?(?:\d+(?:\.\d*)?|\.\d+))',re.I)
OPS=['drilling','boring','reaming','tapping']
CODE={'G81':46,'G82':47,'G83':48,'G84':49,'G85':50,'G89':51,'G80':52,'G0':1,'G4':5,'M3':40,'M5':41,'M8':42,'M9':43,'M2':44,'STATE':45}

def simulate_program(rows, feature_dim=42, accel_nominal=590, rapid_nominal=5700):
    profile=json.loads(rows[0]['initial_profile']);op=profile.get('operation',rows[0].get('operation','drilling'))
    if op not in OPS: raise ValueError(op)
    orient='horizontal' if any(re.search(r'\bG19\b',r['gcode']) for r in rows) else 'vertical'
    # Geometric stock/tool descriptors from simulated machine metadata. These are
    # optional input data (not tool-wear sensor feedback); no target is read here.
    diameter=float(profile.get('hidden_diameter_mm',0.0));pilot=float(profile.get('hidden_pilot_diameter_mm',0.0))
    depth=float(profile.get('hole_depth_mm',0.0));pitch=float(profile.get('pitch_mm',0.0))
    form=bool(profile.get('tap_form',False));material=profile.get('hidden_work_material','AL6061')
    tool=profile.get('hidden_tool_material','CARBIDE');density,se,_=MAT.get(material,MAT['AL6061'])
    tool_power=TOOL.get(tool,TOOL['CARBIDE'])[0]
    area=math.pi/4*max(0,diameter*diameter-pilot*pilot)
    if op=='tapping':
        area=0 if form else math.pi*diameter*pitch*.11
    area=max(0,area)
    # Energy teacher uses operation-specific friction/contact. No hidden wear or
    # hidden motor gain is supplied to this baseline; they remain unpredictable.
    contact_area=math.pi*diameter*pitch*.25 if op=='tapping' else math.pi/4*max(0,diameter**2-pilot**2)
    contact_geometry={'drilling':1.,'boring':.78,'reaming':1.42,'tapping':2.4 if form else 1.7}[op]
    feed=600.;rpm=0.;on=False;cool=False;axispos=5.;crosspos=(12.,12.);cycle=None;cycle_peck=4.5;cycle_dwell=0.3
    out_codes=[];out_feats=[];out_prior=[];out_mass=[];total_energy=0.
    running_total_distance=0.;running_cutvol=0.;running_motion=0.
    running_power_components=np.zeros(8,np.float64)
    ax='X' if orient=='horizontal' else 'Z';cross=('Y','Z') if orient=='horizontal' else ('X','Y')
    for ir,row in enumerate(rows):
        gc=row['gcode'].strip().upper(); pairs=[(k,float(v)) for k,v in RE.findall(gc)]; w=dict(pairs)
        gv=[int(v) for k,v in pairs if k=='G'];mv=[int(v) for k,v in pairs if k=='M']
        previous_cycle=cycle
        for g in gv:
            if g in (81,82,83,84,85,89):cycle='G'+str(g)
            if g==80:cycle=None
        if 'F' in w:feed=w['F']
        if cycle and any(g in (81,82,83,84,85,89) for g in gv):
            cycle_peck=float(w.get('Q',cycle_peck));cycle_dwell=float(w.get('P',cycle_dwell))
        if 'S' in w:rpm=w['S']
        for m in mv:
            if m in (3,4):on=True
            elif m==5:on=False
            elif m==8:cool=True
            elif m==9:cool=False
        code=('G'+str(gv[0])) if gv else ('M'+str(mv[0]) if mv else ('STATE'))
        dt=0.;volume=0.;axial_distance=0.;lateral_distance=0.;is_cycle=False;dwell=0.;rapid_distance=0.
        if 4 in gv:dt=max(0.,w.get('P',0.))
        if 0 in gv:
            new_axis=w.get(ax,axispos)
            new_cross=(w.get(cross[0],crosspos[0]),w.get(cross[1],crosspos[1]))
            axial_distance=abs(new_axis-axispos);lateral_distance=math.dist(crosspos,new_cross)
            rapid_distance=math.hypot(axial_distance,lateral_distance)
            dt=move_time(rapid_distance,rapid_nominal/60,accel_nominal)
            axispos=new_axis;crosspos=new_cross
        # A canned-cycle command OR a positional-only line while a cycle is modal.
        if cycle and (int(cycle[1:]) in gv or (not gv and not mv and any(k in w for k in cross))):
            nc=(w.get(cross[0],crosspos[0]),w.get(cross[1],crosspos[1]))
            lateral_distance=math.dist(crosspos,nc);crosspos=nc
            axial_distance=depth+2.
            speed=max(feed/60,0.01)
            dt=move_time(lateral_distance,rapid_nominal/60,accel_nominal)+move_time(axial_distance,speed,accel_nominal)
            dt+=move_time(axial_distance,speed if cycle in ('G85','G89','G84') else rapid_nominal/60,accel_nominal)
            if cycle=='G83':
                q=max(.1,cycle_peck)
                dt+=.55*(max(1,math.ceil(depth/q))-1)*depth/max(speed,.1)
            if cycle in ('G82','G84','G89'):
                dwell=cycle_dwell
                dt+=dwell
            if cycle=='G84': dt+=.27+.035*depth
            volume=area*depth;is_cycle=True;code=cycle
        energy=0.
        if dt>0:
            base_w=58.+18*(orient=='horizontal')+22.
            spindle_w=(48+rpm*.022+rpm*rpm*1e-6) if on else 0.
            aux_w=(82. if cool else 0.)
            motion_dist=lateral_distance+axial_distance if is_cycle else rapid_distance
            speed_motion=motion_dist/max(dt,1e-8)
            motion_w=(.4*speed_motion+.012*speed_motion**2) if motion_dist else 0.
            # Form tapping consumes contact work EVEN when zero chips are produced.
            vol_for_contact=contact_area*depth if (is_cycle and op=='tapping') else volume
            cutting_J=se*tool_power*vol_for_contact*contact_geometry if is_cycle else 0.
            power_comp=np.asarray([base_w*dt,0.,motion_w*dt*.5,motion_w*dt*.5,0.,spindle_w*dt,cutting_J,aux_w*dt],float)/3600.
            running_power_components+=power_comp
            energy=float(power_comp.sum())
            running_motion+=dt if (rapid_distance or is_cycle) else 0.
            running_total_distance+=motion_dist
        running_cutvol+=volume
        vec=np.zeros(feature_dim,np.float32)
        # Canonical axial/lateral motion ONLY. No raw X/Y/Z positions, signs,
        # physical orientation, or direction indicators leak into Hole Expert.
        fields=[math.log1p(axial_distance+lateral_distance),math.log1p(max(feed/60,0)),math.log1p(max(rpm,0)),
                math.log1p(max(diameter,0)),math.log1p(max(depth,0)),math.log1p(axial_distance),
                math.log1p(lateral_distance),math.log1p(abs(depth)),float(is_cycle),
                float(code in ('G81','G82','G83')),float(code in ('G84','G85','G89')),
                float(code=='G0'),math.log1p(max(volume,0)),float(on),float(cool),float(is_cycle),
                float(code=='G0'),float(code in ('G81','G82','G83')),float(code=='G83'),
                float(op=='boring'),float(op=='reaming'),math.log1p(max(feed,0)),math.log1p(dt),
                float(dwell>0),math.log1p(len(rows)),math.log1p(max(pitch,0)),
                float(form),math.log1p(max(pilot,0)),math.log1p(max(area,0)),math.log1p(max(contact_area,0))]
        vec[:len(fields)]=fields
        vec[37]=1.
        # These 2 are density and machine speed, which prior architecture uses.
        # They are not raw X/Z coordinates. For the hole canonical invariance
        # test, physical orientation affects only the hardware energy baseline.
        vec[40]=density*1000
        vec[41]=math.log1p(rapid_nominal)
        out_codes.append(CODE.get(code, CODE['STATE']))
        out_feats.append(vec)
        out_prior.append([dt,0.,volume*density]);out_mass.append(volume*density)
    # Corrected job-level, twelve-component electromechanical energy prior.
    # Here out_prior durations are G-code-only approximate time predictions.
    # The inference-stage prior is recomputed with the trained time prediction.
    _,raw=components_by_row(rows,np.asarray(out_prior)[:,0],teacher=False)
    return (np.asarray(out_codes,np.int64),np.asarray(out_feats,np.float32),np.asarray(out_prior,np.float32),
            np.asarray(out_mass,np.float64),raw)

	"""Canonical hole-cycle physics approximation, independent of training targets.
	One Hole Processing Expert for G17/Z and G19/X; physical axis retained only for
	machine-energy accounting. New synthetic metadata is NOT suitable for G-code-only
	production deployments unless corresponding tool/stock data are available.
	This is NOT a LinuxCNC runtime simulator.
	"""
	from __future__ import annotations
	import math,re,json
	import numpy as np
	from hole_constants import MAT,TOOL,move_time
	from hole_energy_hardware import components_by_row
	RE=re.compile(r'([A-Z])\s([-+]?(?:\d+(?:\.\d)?\|\.\d+))',re.I)
	OPS=['drilling','boring','reaming','tapping']
	CODE={'G81':46,'G82':47,'G83':48,'G84':49,'G85':50,'G89':51,'G80':52,'G0':1,'G4':5,'M3':40,'M5':41,'M8':42,'M9':43,'M2':44,'STATE':45}

	def simulate_program(rows, feature_dim=42, accel_nominal=590, rapid_nominal=5700):
	profile=json.loads(rows[0]['initial_profile']);op=profile.get('operation',rows[0].get('operation','drilling'))
	if op not in OPS: raise ValueError(op)
	orient='horizontal' if any(re.search(r'\bG19\b',r['gcode']) for r in rows) else 'vertical'
	# Geometric stock/tool descriptors from simulated machine metadata. These are
	# optional input data (not tool-wear sensor feedback); no target is read here.
	diameter=float(profile.get('hidden_diameter_mm',0.0));pilot=float(profile.get('hidden_pilot_diameter_mm',0.0))
	depth=float(profile.get('hole_depth_mm',0.0));pitch=float(profile.get('pitch_mm',0.0))
	form=bool(profile.get('tap_form',False));material=profile.get('hidden_work_material','AL6061')
	tool=profile.get('hidden_tool_material','CARBIDE');density,se,_=MAT.get(material,MAT['AL6061'])
	tool_power=TOOL.get(tool,TOOL['CARBIDE'])[0]
	area=math.pi/4max(0,diameterdiameter-pilot*pilot)
	if op=='tapping':
	area=0 if form else math.pidiameterpitch*.11
	area=max(0,area)
	# Energy teacher uses operation-specific friction/contact. No hidden wear or
	# hidden motor gain is supplied to this baseline; they remain unpredictable.
	contact_area=math.pidiameterpitch.25 if op=='tapping' else math.pi/4max(0,diameter2-pilot2)
	contact_geometry={'drilling':1.,'boring':.78,'reaming':1.42,'tapping':2.4 if form else 1.7}[op]
	feed=600.;rpm=0.;on=False;cool=False;axispos=5.;crosspos=(12.,12.);cycle=None;cycle_peck=4.5;cycle_dwell=0.3
	out_codes=[];out_feats=[];out_prior=[];out_mass=[];total_energy=0.
	running_total_distance=0.;running_cutvol=0.;running_motion=0.
	running_power_components=np.zeros(8,np.float64)
	ax='X' if orient=='horizontal' else 'Z';cross=('Y','Z') if orient=='horizontal' else ('X','Y')
	for ir,row in enumerate(rows):
	gc=row['gcode'].strip().upper(); pairs=[(k,float(v)) for k,v in RE.findall(gc)]; w=dict(pairs)
	gv=[int(v) for k,v in pairs if k=='G'];mv=[int(v) for k,v in pairs if k=='M']
	previous_cycle=cycle
	for g in gv:
	if g in (81,82,83,84,85,89):cycle='G'+str(g)
	if g==80:cycle=None
	if 'F' in w:feed=w['F']
	if cycle and any(g in (81,82,83,84,85,89) for g in gv):
	cycle_peck=float(w.get('Q',cycle_peck));cycle_dwell=float(w.get('P',cycle_dwell))
	if 'S' in w:rpm=w['S']
	for m in mv:
	if m in (3,4):on=True
	elif m==5:on=False
	elif m==8:cool=True
	elif m==9:cool=False
	code=('G'+str(gv[0])) if gv else ('M'+str(mv[0]) if mv else ('STATE'))
	dt=0.;volume=0.;axial_distance=0.;lateral_distance=0.;is_cycle=False;dwell=0.;rapid_distance=0.
	if 4 in gv:dt=max(0.,w.get('P',0.))
	if 0 in gv:
	new_axis=w.get(ax,axispos)
	new_cross=(w.get(cross[0],crosspos[0]),w.get(cross[1],crosspos[1]))
	axial_distance=abs(new_axis-axispos);lateral_distance=math.dist(crosspos,new_cross)
	rapid_distance=math.hypot(axial_distance,lateral_distance)
	dt=move_time(rapid_distance,rapid_nominal/60,accel_nominal)
	axispos=new_axis;crosspos=new_cross
	# A canned-cycle command OR a positional-only line while a cycle is modal.
	if cycle and (int(cycle[1:]) in gv or (not gv and not mv and any(k in w for k in cross))):
	nc=(w.get(cross[0],crosspos[0]),w.get(cross[1],crosspos[1]))
	lateral_distance=math.dist(crosspos,nc);crosspos=nc
	axial_distance=depth+2.
	speed=max(feed/60,0.01)
	dt=move_time(lateral_distance,rapid_nominal/60,accel_nominal)+move_time(axial_distance,speed,accel_nominal)
	dt+=move_time(axial_distance,speed if cycle in ('G85','G89','G84') else rapid_nominal/60,accel_nominal)
	if cycle=='G83':
	q=max(.1,cycle_peck)
	dt+=.55(max(1,math.ceil(depth/q))-1)depth/max(speed,.1)
	if cycle in ('G82','G84','G89'):
	dwell=cycle_dwell
	dt+=dwell
	if cycle=='G84': dt+=.27+.035*depth
	volume=area*depth;is_cycle=True;code=cycle
	energy=0.
	if dt>0:
	base_w=58.+18*(orient=='horizontal')+22.
	spindle_w=(48+rpm.022+rpmrpm*1e-6) if on else 0.
	aux_w=(82. if cool else 0.)
	motion_dist=lateral_distance+axial_distance if is_cycle else rapid_distance
	speed_motion=motion_dist/max(dt,1e-8)
	motion_w=(.4speed_motion+.012speed_motion**2) if motion_dist else 0.
	# Form tapping consumes contact work EVEN when zero chips are produced.
	vol_for_contact=contact_area*depth if (is_cycle and op=='tapping') else volume
	cutting_J=setool_powervol_for_contact*contact_geometry if is_cycle else 0.
	power_comp=np.asarray([base_wdt,0.,motion_wdt.5,motion_wdt.5,0.,spindle_wdt,cutting_J,aux_w*dt],float)/3600.
	running_power_components+=power_comp
	energy=float(power_comp.sum())
	running_motion+=dt if (rapid_distance or is_cycle) else 0.
	running_total_distance+=motion_dist
	running_cutvol+=volume
	vec=np.zeros(feature_dim,np.float32)
	# Canonical axial/lateral motion ONLY. No raw X/Y/Z positions, signs,
	# physical orientation, or direction indicators leak into Hole Expert.
	fields=[math.log1p(axial_distance+lateral_distance),math.log1p(max(feed/60,0)),math.log1p(max(rpm,0)),
	math.log1p(max(diameter,0)),math.log1p(max(depth,0)),math.log1p(axial_distance),
	math.log1p(lateral_distance),math.log1p(abs(depth)),float(is_cycle),
	float(code in ('G81','G82','G83')),float(code in ('G84','G85','G89')),
	float(code=='G0'),math.log1p(max(volume,0)),float(on),float(cool),float(is_cycle),
	float(code=='G0'),float(code in ('G81','G82','G83')),float(code=='G83'),
	float(op=='boring'),float(op=='reaming'),math.log1p(max(feed,0)),math.log1p(dt),
	float(dwell>0),math.log1p(len(rows)),math.log1p(max(pitch,0)),
	float(form),math.log1p(max(pilot,0)),math.log1p(max(area,0)),math.log1p(max(contact_area,0))]
	vec[:len(fields)]=fields
	vec[37]=1.
	# These 2 are density and machine speed, which prior architecture uses.
	# They are not raw X/Z coordinates. For the hole canonical invariance
	# test, physical orientation affects only the hardware energy baseline.
	vec[40]=density*1000
	vec[41]=math.log1p(rapid_nominal)
	out_codes.append(CODE.get(code, CODE['STATE']))
	out_feats.append(vec)
	out_prior.append([dt,0.,volumedensity]);out_mass.append(volumedensity)
	# Corrected job-level, twelve-component electromechanical energy prior.
	# Here out_prior durations are G-code-only approximate time predictions.
	# The inference-stage prior is recomputed with the trained time prediction.
	_,raw=components_by_row(rows,np.asarray(out_prior)[:,0],teacher=False)
	return (np.asarray(out_codes,np.int64),np.asarray(out_feats,np.float32),np.asarray(out_prior,np.float32),
	np.asarray(out_mass,np.float64),raw)