"""Проверка табличной карты и восстановление энергии из той же интерполяции."""
import math
from numerics import FluxMap,reciprocity_error

def means(plan,rows):
    result={}
    for point in plan['points']:
        selected=sorted([r for r in rows if r['point']==point['id']],key=lambda r:r['theta_deg'])
        n=point['angles']
        if len(selected)!=n or any(abs(r['theta_deg']-360*k/n)>1e-8 for k,r in enumerate(selected)):raise ValueError('Неполная угловая серия: '+point['id'])
        if any(r['currents_dqf_A']!=point['currents_A'] for r in selected):raise ValueError('Изменён ток точки карты')
        average=lambda values:sum(values)/len(values)
        result[point['id']]={'current_A':point['currents_A'], 'lambda_Wb':[average([r['lambda_dqf_Wb'][j] for r in selected]) for j in range(3)],
            'energy_J':average([r['energy_J'] for r in selected]),'coenergy_J':average([r['coenergy_J'] for r in selected]),
            'torque_Nm':average([r['force_torque'][2] for r in selected]),
            'half_samples_lambda_Wb':[average([r['lambda_dqf_Wb'][j] for r in selected[::2]]) for j in range(3)]}
    return result

def coenergy_increment(magnetic_map,start,end):
    """Точный интеграл трилинейной λ вдоль отрезка (2-point Gauss по ячейкам)."""
    delta=[end[j]-start[j] for j in range(3)];cuts=[0.,1.];metric=[1.5,1.5,1.]
    for axis,a,d in zip(magnetic_map.axes,start,delta):
        if d:
            cuts.extend((x-a)/d for x in axis[1:-1] if 0<(x-a)/d<1)
    cuts=sorted(set(cuts));total=0.
    for lo,hi in zip(cuts,cuts[1:]):
        for sign in (-1,1):
            t=(hi+lo)/2+sign*(hi-lo)/(2*math.sqrt(3));point=[start[j]+t*delta[j] for j in range(3)]
            flux,_=magnetic_map.evaluate(point)
            total+=(hi-lo)/2*sum(metric[j]*delta[j]*flux[j] for j in range(3))
    return total

def magnetic_energy(magnetic_map,current,anchor,anchor_energy):
    if hasattr(magnetic_map,'energy'):return magnetic_map.energy(current)
    weights=[1.5,1.5,1.];flux,_=magnetic_map.evaluate(current);base,_=magnetic_map.evaluate(anchor)
    coenergy0=sum(weights[j]*anchor[j]*base[j] for j in range(3))-anchor_energy
    return sum(weights[j]*current[j]*flux[j] for j in range(3))-coenergy0-coenergy_increment(magnetic_map,anchor,current)

def process(plan,rows,interpolation='component-trilinear'):
    averaged=means(plan,rows);axes=plan['axes_A']
    values=[[[averaged['node-%d%d%d'%(i,j,k)]['lambda_Wb'] for k in range(len(axes[2]))] for j in range(len(axes[1]))] for i in range(len(axes[0]))]
    coenergy=[[[averaged['node-%d%d%d'%(i,j,k)]['coenergy_J'] for k in range(len(axes[2]))] for j in range(len(axes[1]))] for i in range(len(axes[0]))]
    if interpolation=='coenergy-tricubic-hermite':
        from energy_map import EnergyFluxMap
        magnetic_map=EnergyFluxMap(axes,values,coenergy)
    elif interpolation=='component-trilinear':magnetic_map=FluxMap(axes,values)
    else:raise ValueError('Неизвестная интерполяция')
    criteria=plan['criteria'];checks=[]
    anchor=averaged['open-circuit']
    for point in plan['points']:
        if point['purpose']!='independent-control':continue
        reference=averaged[point['id']];flux,jac=magnetic_map.evaluate(point['currents_A'])
        error=[abs(a-b)/max(abs(b),criteria['flux_absolute_scale_Wb']) for a,b in zip(flux,reference['lambda_Wb'])]
        angular=[abs(a-b)/max(abs(b),criteria['flux_absolute_scale_Wb']) for a,b in zip(reference['half_samples_lambda_Wb'],reference['lambda_Wb'])]
        reciprocity=reciprocity_error(jac,criteria['reciprocity_absolute_scale_H'])
        energy=magnetic_energy(magnetic_map,point['currents_A'],anchor['current_A'],anchor['energy_J'])
        energy_error=abs(energy-reference['energy_J'])/max(abs(reference['energy_J']),1000.)
        checks.append({'point':point['id'],'flux_relative_errors':error,'angular_relative_errors':angular,'reciprocity_error':reciprocity,
                       'energy_J_from_map':energy,'energy_relative_error':energy_error,
                       'passed':max(error)<=criteria['flux_relative'] and max(angular)<=criteria['angular_relative'] and reciprocity<=criteria['reciprocity_relative'] and energy_error<=criteria['energy_balance_relative']})
    # Necessary local invertibility/energy convexity check throughout each cell.
    # Sampling is a numerical check, not an analytic global convexity proof.
    samples=[sorted(set(a+(b-a)*k/4 for a,b in zip(axis,axis[1:]) for k in range(5))) for axis in axes]
    minimum_pivot=1.;count=0;positive=True
    for d in samples[0]:
        for q in samples[1]:
            for f in samples[2]:
                _,jac=magnetic_map.evaluate([d,q,f]);weights=(1.5,1.5,1.)
                h=[[weights[j]*jac[j][k] for k in range(3)] for j in range(3)]
                scale=max(abs(h[j][j]) for j in range(3));lower=[[0.]*3 for _ in range(3)]
                for j in range(3):
                    pivot=h[j][j]-sum(lower[j][k]**2 for k in range(j))
                    minimum_pivot=min(minimum_pivot,pivot/max(scale,1e-30))
                    if pivot<=0:positive=False;break
                    lower[j][j]=math.sqrt(pivot)
                    for k in range(j+1,3):lower[k][j]=(h[k][j]-sum(lower[k][m]*lower[j][m] for m in range(j)))/lower[j][j]
                count+=1
    convexity={'points':count,'positive_at_all_samples':positive,'minimum_scaled_cholesky_pivot':minimum_pivot,'scope':'sampled, not a global proof'}
    return {'axes_A':axes,'values_Wb':values,'coenergy_J':coenergy,'interpolation':interpolation,'means':averaged,'checks':checks,'convexity':convexity,'validated':all(row['passed'] for row in checks) and positive,
            'criteria':criteria,'anchor_current_A':anchor['current_A'],'anchor_energy_J':anchor['energy_J'],
            'scope':plan['scope'],'extrapolation':'forbidden'}