"""Малые численные процедуры без NumPy; встроенный Python FEMM 5.

Эти процедуры сами не создают магнитную карту и не заменяют проверку FEMM.
"""
import math


def solve_linear(matrix, rhs):
    n=len(rhs)
    a=[list(row)+[rhs[k]] for k,row in enumerate(matrix)]
    if len(a)!=n or any(len(row)!=n+1 for row in a):raise ValueError('Размеры матрицы')
    norm=max(abs(v) for row in matrix for v in row)
    if not norm:raise ValueError('Нулевая матрица')
    for k in range(n):
        pivot=max(range(k,n),key=lambda j:abs(a[j][k]))
        if abs(a[pivot][k])<1e-12*norm:raise ValueError('Вырожденный якобиан')
        a[k],a[pivot]=a[pivot],a[k]
        divisor=a[k][k]
        for j in range(k,n+1):a[k][j]/=divisor
        for i in range(n):
            if i==k:continue
            multiplier=a[i][k]
            for j in range(k,n+1):a[i][j]-=multiplier*a[k][j]
    return [row[-1] for row in a]


def newton_bounded(evaluate, initial, bounds, scales, tolerance=1e-3, maximum=30, step_fraction=.2):
    """evaluate(x) -> (невязка, якобиан). Журнал содержит каждую принятую точку.

    Масштабы невязок положительны; шаг ограничен долей ширины области.
    Невозможность уменьшить невязку возвращается как ошибка, не как режим.
    """
    n=len(initial)
    if len(bounds)!=n or len(scales)!=n or any(s<=0 for s in scales):raise ValueError('Размеры и масштабы')
    if not n or tolerance<=0 or maximum<0:raise ValueError('Некорректный критерий остановки')
    if any(not math.isfinite(v) for v in list(initial)+list(scales)+[x for pair in bounds for x in pair]):
        raise ValueError('Нечисловые входы')
    if not 0<step_fraction<=1:raise ValueError('Ограничение шага')
    if any(not lo<hi for lo,hi in bounds):raise ValueError('Пустые границы')
    x=list(initial); log=[]
    if any(x[k]<bounds[k][0] or x[k]>bounds[k][1] for k in range(n)):
        raise ValueError('Начальная точка вне области')
    for iteration in range(maximum+1):
        residual,jacobian=evaluate(x)
        if len(residual)!=n or len(jacobian)!=n or any(len(row)!=n for row in jacobian):raise ValueError('Размеры отклика')
        if any(not math.isfinite(v) for v in residual+[v for row in jacobian for v in row]):raise ValueError('Нечисловой отклик')
        norm=max(abs(residual[k])/scales[k] for k in range(n))
        log.append({'iteration':iteration,'x':x[:],'residual':residual[:],'norm':norm})
        if norm<=tolerance:return {'converged':True,'x':x,'log':log}
        if iteration==maximum:break
        try:
            step=solve_linear([[jacobian[i][j]/scales[i] for j in range(n)] for i in range(n)],
                              [-residual[i]/scales[i] for i in range(n)])
        except ValueError as error:
            return {'converged':False,'x':x,'log':log,'reason':str(error)}
        factor=min([1.0]+[step_fraction*(bounds[k][1]-bounds[k][0])/abs(step[k]) for k in range(n) if step[k]])
        accepted=False
        for attempt in range(24):
            trial=[x[k]+factor*step[k] for k in range(n)]
            if all(bounds[k][0]<=trial[k]<=bounds[k][1] for k in range(n)):
                values,_=evaluate(trial)
                if max(abs(values[k])/scales[k] for k in range(n))<norm:
                    x=trial;accepted=True;break
            factor*=.5
        if not accepted:return {'converged':False,'x':x,'log':log,'reason':'Шаг не уменьшает невязку внутри области'}
    return {'converged':False,'x':x,'log':log,'reason':'Исчерпано число итераций'}


class FluxMap:
    """Трилинейная карта λd,λq,λf на прямоугольной сетке id,iq,If.

    values[i][j][k][component], все токи А, потокосцепления Вб·виток.
    Якобиан — точная производная этой же интерполяции внутри ячейки.
    На границе ячеек выбирается правая ячейка, на внешней границе — левая.
    """
    def __init__(self, axes, values):
        if len(axes)!=3:raise ValueError('Нужны три оси')
        for axis in axes:
            if len(axis)<2 or any(not math.isfinite(x) for x in axis) or any(b<=a for a,b in zip(axis,axis[1:])):raise ValueError('Оси должны возрастать')
        if len(values)!=len(axes[0]):raise ValueError('Размер карты d')
        for plane in values:
            if len(plane)!=len(axes[1]):raise ValueError('Размер карты q')
            for row in plane:
                if len(row)!=len(axes[2]) or any(len(v)!=3 for v in row):raise ValueError('Размер карты f')
                if any(not math.isfinite(x) for v in row for x in v):raise ValueError('Нечисловое значение карты')
        self.axes=axes;self.values=values

    def evaluate(self, point):
        if len(point)!=3:raise ValueError('Нужны id,iq,If')
        indices=[];weights=[];widths=[]
        for axis,x in zip(self.axes,point):
            if not math.isfinite(x) or x<axis[0] or x>axis[-1]:raise ValueError('Экстраполяция карты запрещена')
            k=0
            while k<len(axis)-2 and x>=axis[k+1]:k+=1
            width=axis[k+1]-axis[k];indices.append(k);widths.append(width);weights.append((x-axis[k])/width)
        flux=[0.0]*3;jac=[[0.0]*3 for _ in range(3)]
        for a in (0,1):
            for b in (0,1):
                for c in (0,1):
                    bits=(a,b,c);w=[weights[d] if bits[d] else 1-weights[d] for d in range(3)]
                    value=self.values[indices[0]+a][indices[1]+b][indices[2]+c]
                    for component in range(3):
                        flux[component]+=value[component]*w[0]*w[1]*w[2]
                        for d in range(3):
                            jac[component][d]+=value[component]*(1 if bits[d] else -1)/widths[d]*w[(d+1)%3]*w[(d+2)%3]
        return flux,jac

    def inverse(self, flux, initial, scales, tolerance=1e-8):
        def residual(current):
            value,jac=self.evaluate(current)
            return [value[k]-flux[k] for k in range(3)],jac
        result=newton_bounded(residual,initial,[(a[0],a[-1]) for a in self.axes],scales,tolerance,maximum=60)
        if not result['converged']:raise ValueError('Обращение карты: '+result['reason'])
        return result


def reciprocity_error(jacobian, absolute_scale):
    """Проверяется симметрия diag(3/2,3/2,1) * dλ/di, не просто dλ/di."""
    weights=(1.5,1.5,1.0)
    if absolute_scale<=0:raise ValueError('Нужна абсолютная шкала индуктивности')
    return max(abs(weights[i]*jacobian[i][j]-weights[j]*jacobian[j][i]) /
        max(absolute_scale,abs(weights[i]*jacobian[i][j]),abs(weights[j]*jacobian[j][i]))
        for i in range(3) for j in range(i))


def rk4(rhs, time, state, step):
    """Один шаг; выход за область карты распространяется как ошибка."""
    if step<=0:raise ValueError('Шаг должен быть положительным')
    n=len(state);k1=rhs(time,state)
    k2=rhs(time+step/2,[state[j]+step*k1[j]/2 for j in range(n)])
    k3=rhs(time+step/2,[state[j]+step*k2[j]/2 for j in range(n)])
    k4=rhs(time+step,[state[j]+step*k3[j] for j in range(n)])
    return [state[j]+step*(k1[j]+2*k2[j]+2*k3[j]+k4[j])/6 for j in range(n)]