<jupyter_start><jupyter_text><jupyter_code>import numpy as np

from keras.models import Sequential
from keras.layers import Dense, Flatten, Activation
from keras.layers import Dropout
from keras.layers.convolutional import Conv2D, MaxPooling2D
from keras.utils import np_utils
from tensorflow.keras.optimizers import SGD, Adam

import pandas as pd
import numpy as np
import scipy as sc
from scipy.interpolate import interp2d
from scipy.optimize import curve_fit
from scipy.integrate import quad, trapz
from scipy import optimize
from mpmath import *
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import seaborn as sns
import cv2
import imageio
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
import pandas as pd
from scipy.fftpack import fft2, ifft2
from numba import jit

import cv2

import numpy as np
import matplotlib.pyplot as plt
import abel

import plotly.express as px
from plotly.subplots import make_subplots

import tifffile as tif

import numpy as np
from scipy.special import gamma
from scipy.linalg import svd

from scipy import interpolate


def delt(m,m1):
    delta=0
    if m==m1:
        delta=1
    return delta

def E_f(m,n,m1,n1,D,r_0):
    a=float((n+n1-2*m)/2)
    a1=(-1)**a
    b=np.sqrt((n+1)*(n1+1))*delt(m,m1)
    c=gamma((n+n1-5/3)/2)*(D/r_0)**(5/3)
    d=gamma((n-n1+17/3)/2)*gamma((n1-n-17/3)/2)*gamma((n+n1+23/3)/2)
    E = a1*b*c/d
    #print(a1,b,c,d)
    return E



def rand_num(v,mu,N):
    v = v/0.3# variance
    sigma=np.sqrt(v)
    X = sigma * (np.random.rand(N)-0.5)*2 + mu
    return X

def turb(N,N_zern=25,D=18, r_0=8):
    C=np.zeros((N_zern-3,N_zern-3),dtype=np.complex)
    z_mn=np.zeros((2,N_zern-3))
    z_mn[:,0]=[0,2]
    z_mn[:,1]=[2,2]
    z_mn[:,2]=[-2,2]
    z_mn[:,3]=[1,3]
    z_mn[:,4]=[-1,3]
    z_mn[:,5]=[0,4]
    z_mn[:,6]=[3,3]
    z_mn[:,7]=[-3,3]

    z_mn[:,8]=[2,4]
    z_mn[:,9]=[-2,4]
    z_mn[:,10]=[1,5]
    z_mn[:,11]=[-1,5]
    z_mn[:,12]=[0,6]
    z_mn[:,13]=[4,4]
    z_mn[:,14]=[-4,4]

    z_mn[:,15]=[3,5]
    z_mn[:,16]=[-3,5]
    z_mn[:,17]=[2,6]
    z_mn[:,18]=[-2,6]
    z_mn[:,19]=[1,7]
    z_mn[:,20]=[-1,7]
    z_mn[:,21]=[0,8]


    for i in range(N_zern-3):
        for j in range(N_zern-3):
            C[i,j]=(E_f(z_mn[0,i],z_mn[1,i],z_mn[0,j],z_mn[1,j],D,r_0))

    U,S,V=svd(C)

    B=np.zeros((N,N_zern-3))
    for i in range(N_zern-3):
        B[:,i]=rand_num(S[i],0,N)

    A=B.dot(V.real) # A матрица столбцы-коэф полинома
    zern_coef=np.zeros((N,25))
    zern_coef[:,3:25]=A
    return zern_coef
def zern_act_reader(name,n_act,n_zern):
    n_zern=n_zern-2
    df = pd.read_csv(name,delimiter='\t',header=0,encoding='utf-8')
    zern_act=np.zeros((n_act,n_zern))
    for i in range(0,n_act):
        p=df.iloc[i]
        d=p[4:26].to_numpy()
        zern_act[i,:]=d
    return zern_act

def circle_in_square(Nx, Ny, R):
    x = np.linspace(-1, 1, Nx)
    y = np.linspace(-1, 1, Ny)
    circ = np.zeros((Nx, Ny))
    for i in range(Nx):
        for j in range(Ny):
            if (i - Nx / 2) ** 2 + (j - Nx / 2) ** 2 <= R ** 2:
                circ[i, j] = 100
    return x, y, circ


def zern_fun(wavelen, Nx, Ny, R, zern_coef):
    center = Nx // 2
    phase = np.zeros((Nx, Ny))
    
    for i in range(Nx):
        for k in range(Ny):
            cx = k - center
            cy = i - center
            z0 = zern_coef[0]  # СМЕЩЕНИЕ
            z01 = zern_coef[1] * (cx) / R  # клинХ
            z02 = zern_coef[2] * (cy) / R  # клинУ
            z03 = zern_coef[3] * (((cx) ** 2 + (cy) ** 2) * 2 - 1) / R ** 2  # дефокусdef(op)
            z04 = zern_coef[4] * ((cx) ** 2 - (cy) ** 2) / R ** 2  # астигматизмХ
            z05 = zern_coef[5] * 2 * (cx) * (cy) / R ** 2  # астигматизмУ
            z06 = zern_coef[6] * (cy) * (3 * ((cx) ** 2 + (cy) ** 2) / R ** 2 - 2) / R  # комаХ
            z07 = zern_coef[7] * (cx) * (3 * ((cx) ** 2 + (cy) ** 2) / R ** 2 - 2) / R  # комаУ
            z08 = zern_coef[8] * (
                    6 * ((cx) ** 2 + (cy) ** 2) ** 2 / R ** 4 - 6 * ((cx) ** 2 + (cy) ** 2) / R ** 2 + 1)  # сфера
            z09 = zern_coef[9] * (cx) * ((cx) ** 2 - 3 * (cy) ** 2) / R ** 3  # трилистникХ
            z10 = zern_coef[10] * (cy) * (3 * (cx) ** 2 - (cy) ** 2) / R ** 3  # трилистникУ
            z11 = zern_coef[11] * ((cx) ** 2 - (cy) ** 2) / R ** 2 * (
                    4 * ((cx) ** 2 + (cy) ** 2) / R ** 2 - 3)  # астигмат Х
            z12 = zern_coef[12] * 2 * ((cx) * (cy)) / R ** 2 * (4 * ((cx) ** 2 + (cy) ** 2) / R ** 2 - 3)  # астигмат У
            z13 = zern_coef[13] * (cx) * (10 * ((cx) ** 2 + (cy) ** 2) ** 2 / R ** 4 - 12 * (
                    (cx) ** 2 + (cy) ** 2) / R ** 2 + 3) / R  # кома Х
            z14 = zern_coef[14] * (cy) * (10 * ((cx) ** 2 + (cy) ** 2) ** 2 / R ** 4 - 12 * (
                    (cx) ** 2 + (cy) ** 2) / R ** 2 + 3) / R  # кома У
            z15 = zern_coef[15] * (
                    20 * ((cx) ** 2 + (cy) ** 2) ** 3 / R ** 6 - 30 * ((cx) ** 2 + (cy) ** 2) ** 2 / R ** 4 + 12 * (
                    (cx) ** 2 + (cy) ** 2) / R ** 2 - 1)  # сфера
            z16 = zern_coef[16] * ((cx) ** 4 / R ** 4 - 6 * (cx) ** 2 * (cy) ** 2 / R ** 4 + (
                cy) ** 4 / R ** 4)  # четверолистник Х
            z17 = zern_coef[17] * 4 * (cx) * (cy) * ((cx) ** 2 - (cy) ** 2) / R ** 4  # четверолистник Y
            z18 = zern_coef[18] * (cx) / R * ((cx) ** 2 - 3 * (cy) ** 2) / R ** 2 * (
                    5 * ((cx) ** 2 + (cy) ** 2) / R ** 2 - 4)  # трилистнХ
            z19 = zern_coef[19] * (cy) / R * (3 * (cx) ** 2 - (cy) ** 2) / R ** 2 * (
                    5 * ((cx) ** 2 + (cy) ** 2) / R ** 2 - 4)  # трилистнY
            z20 = zern_coef[20] * ((cx) ** 2 - (cy) ** 2) / R ** 2 * (
                    15 * ((cx) ** 2 + (cy) ** 2) ** 2 / R ** 4 - 20 * (
                    (cx) ** 2 + (cy) ** 2) / R ** 2 + 6)  # астигм Х
            z21 = zern_coef[21] * (2 * (cx) * (cy)) / R ** 2 * (15 * ((cx) ** 2 + (cy) ** 2) ** 2 / R ** 4 - 20 * (
                    (cx) ** 2 + (cy) ** 2) / R ** 2 + 6)  # астигм Y
            z22 = zern_coef[22] * (cx) / R * (
                    35 * ((cx) ** 2 + (cy) ** 2) ** 3 / R ** 6 - 60 * ((cx) ** 2 + (cy) ** 2) ** 2 / R ** 4 + 30 * (
                    (cx) ** 2 + (cy) ** 2) / R ** 2 - 4)  # комаХ
            z23 = zern_coef[23] * (cx) / R * (
                    35 * ((cx) ** 2 + (cy) ** 2) ** 3 / R ** 6 - 60 * ((cx) ** 2 + (cy) ** 2) ** 2 / R ** 4 + 30 * (
                    (cx) ** 2 + (cy) ** 2) / R ** 2 - 4)  # комаY
            z24 = zern_coef[24] * (70 * ((cx) ** 2 + (cy) ** 2) ** 4 / R ** 8 - 140 * (
                    (cx) ** 2 + (cy) ** 2) ** 3 / R ** 6 + 90 * ((cx) ** 2 + (cy) ** 2) ** 2 / R ** 4 - 20 * (
                                           (cx) ** 2 + (cy) ** 2) / R ** 2 + 1)  # сфера
            phase[i][k] = 1e-6 * (
                    z01 + z02 + z03 + z04 + z05 + z06 + z07 + z08 + z09 + z10 + z11 + z12 + z13 + z14 + z15 + z16 + z17 + z18 + z19 + z20 + z21 + z22 + z23 + z24 + z0)
    return phase


def scaled(nf, newN, oldN):
    new_nf = np.zeros((newN, newN))
    center = oldN // 2
    center1 = newN // 2
    new_nf[center1 - center:center1 + center, center1 - center:center1 + center] = nf
    return new_nf


def farfield(new_N, Nx, nf, phase_surf, wavelen):
    # new_N = 1001
    nf_sc = scaled(nf, new_N, Nx)
    phase_surf_sc = scaled(phase_surf, new_N, Nx)
    e_field = nf_sc * np.exp(2 * np.pi * 1j * phase_surf_sc / wavelen)
    ff = np.fft.fft2(e_field, norm=None)
    ff = np.fft.fftshift(ff)
    ff = np.abs(ff)**2
    return ff

def gaus(x,a,x0,sigma):
    return a/sigma*np.exp(-(x-x0)**2/(2*sigma**2))

def size_spot(ff):
    slise_x=np.sum(ff, axis = 0)
    slise_x=slise_x/np.max(slise_x)
    slise_y=np.sum(ff, axis = 1)
    slise_y=slise_y/np.max(slise_y)
    x = range(len(slise_x))
    y = range(len(slise_x))
    poptx,pcovx = curve_fit(gaus,x,slise_x,p0=[1,500,80])
    popty,pcovy = curve_fit(gaus,y,slise_y,p0=[1,500,80])
    sig_x=poptx[2]
    sig_y=popty[2]
    #return np.sqrt(sig_x**2+sig_y**2)
    return (sig_x+sig_y)*np.max([sig_x,sig_y])

    
def ff_coord(lam,F,R_real,R_pix,new_N):
    x_ff=np.linspace(-1/(2*(R_real/R_pix)),1/(2*(R_real/R_pix)),new_N)*lam*F
    return x_ff


def interp_ff(ff,x_ff_old,x_ff_new):
    ff_interp = interpolate.interp2d(x_ff_old, x_ff_old, ff, kind='cubic')
    ff_new=ff_interp(x_ff_new, x_ff_new)
    return ff_new
# Параметры
wavelen = 1e-06 

# Единичная окружность для ближней зоны
new_N = 501
Nx = Ny = 100
R = Nx / 2
[x_nf, y_nf, nf] = circle_in_square(Nx, Ny, R)

zern_coef_0 = np.zeros(25)
zern_coef_ideal = np.zeros(25)
phase_surf_ideal = zern_fun(wavelen, Nx, Ny, R, zern_coef_ideal)
phase_surf_ideal = phase_surf_ideal * (nf > 1)
ff_ideal = farfield(new_N, Nx, nf, phase_surf_ideal, wavelen)


zern_turb = turb(100,r_0=20)
phase_surf_turb = zern_fun(wavelen, Nx, Ny, R, zern_turb[int(np.random.random(1)*100),:])
phase_surf_turb = phase_surf_turb * (nf > 1)
ff_turb = farfield(new_N, Nx, nf, phase_surf_turb, wavelen)
plt.figure()
plt.imshow(phase_surf_turb*1e6)
plt.colorbar()
plt.show()
#%matplotlib qt
plt.figure()
plt.imshow((ff_turb))
plt.show()

print('Штрель', np.max(ff_turb) / np.max(ff_ideal))
print('RMS_0', np.std(phase_surf_turb) * 1e6)
print('PV_0', np.max(phase_surf_turb) * 1e6 - np.min(phase_surf_turb) * 1e6)
a=np.where(ff_ideal==np.max(ff_ideal))
def datagen(batch_size,b=20):
    zern_turb = turb(batch_size,r_0=b)
    l=len(zern_turb)
    x_out=[]
    y_out=[]
    ph=[]
    for i in range(l):
        zern_coef=zern_turb[i]
        zern_coef[0:2]=0
        #zern_coef[6:]=0
        zern_coef[3]=0.2
        phase_surf = zern_fun(wavelen, Nx, Ny, R, zern_coef)
        zern_coef[3]=0
        ph.append(phase_surf)
        phase_surf = phase_surf * (nf > 1)
        ff = farfield(new_N, Nx, nf, phase_surf, wavelen)
        """ #ff=np.log(ff)
        x=ff[int(a[0]-da):int(a[0]+da),int(a[1]-da):int(a[1]+da)]
        #x=cv2.resize(x, (50,50))
        x=x/np.max(x) """
        """ ff = np.fft.fft2(ff, norm=None)
        ff = np.fft.fftshift(ff)
        ff = np.abs(ff)**2/np.max(np.abs(ff)**2) """
        ff=ff[ff.shape[0]//2-20:ff.shape[0]//2+20,ff.shape[0]//2-20:ff.shape[0]//2+20]
        ff/=np.max(ff)
        #ff=(ff>0.05)
        PolarImage, r_grid, theta_grid = abel.tools.polar.reproject_image_into_polar(ff)
        x_out.append(PolarImage)
        y_out.append(zern_coef[3:])
        ''' plt.figure()
        plt.imshow((ff))
        plt.figure()
        plt.imshow(((PolarImage))) '''
        
    x_out=np.array(x_out)
    y_out=np.array(y_out)
    ph=np.array(ph)
    return x_out,y_out,ph
datagen(1)
model = Sequential()
x_train, y_train,c = datagen(1)
x_train = x_train.reshape(len(x_train), x_train.shape[1]*x_train.shape[2])
model.add(Dense(256, activation='relu', input_shape=(x_train.shape[1],)))
model.add(Dense(128, activation='relu'))
model.add(Dense(64, activation='relu'))
model.add(Dense(y_train.shape[1]))
model.compile(optimizer='adam', loss='mse', metrics=['mae'])
b_size=200
for i in range(1000):
    x_train, y_train,c = datagen(b_size,b=8)
    x_train = x_train.reshape(len(x_train), x_train.shape[1]*x_train.shape[2])
    print('start')
    
    model.fit(x_train, y_train, epochs=100, verbose=2)
x_train, y_train,c = datagen(1,b=20)
#x_train=x_out
x_train = x_train.reshape(len(x_train), x_train.shape[1]*x_train.shape[2])
y_pred=model.predict(x_train)
zern_pred=np.zeros(25)
zern_0=np.zeros(25)
zern_pred[3:]=y_pred[0,:]
zern_0[3:]=y_train[0,:]

''' zern_pred[3]=0.2
zern_0[3]=0.2 '''

phase_surf_0 = zern_fun(wavelen, Nx, Ny, R, zern_0)
phase_surf_0 = phase_surf_0 * (nf > 1)
phase_surf_pred = zern_fun(wavelen, Nx, Ny, R, zern_pred)
phase_surf_pred = phase_surf_pred * (nf > 1)
print('RMS_0', np.std(phase_surf_0) * 1e6)
print('RMS_pred', np.std(phase_surf_pred) * 1e6)
print('RMS_0-pred', np.std(phase_surf_0-phase_surf_pred) * 1e6)

fig,[ax1,ax2,ax3]=plt.subplots(1,3)
ax1.imshow(phase_surf_0,clim=(-0.2*1e-6, 0.2*1e-6))
ax1.set_title('Начальное')
ax2.imshow(phase_surf_pred,clim=(-0.2*1e-6, 0.2*1e-6))
ax2.set_title('По предсказанным \n цернике')
ax3.imshow(phase_surf_0-phase_surf_pred,clim=(-0.2*1e-6, 0.2*1e-6))
ax3.set_title('Скорректированное\n по предсказанию')
ff1=farfield(new_N, Nx, nf, phase_surf_0, wavelen)
ff2=farfield(new_N, Nx, nf, phase_surf_pred, wavelen)
ff3=farfield(new_N, Nx, nf, phase_surf_0-phase_surf_pred, wavelen)
ff0=farfield(new_N, Nx, nf, phase_surf_0*0, wavelen)

print('S_0', np.max(ff1)/np.max(ff0))
print('S_res', np.max(ff3)/np.max(ff0))

fig,[ax1,ax2,ax3]=plt.subplots(1,3)
ax1.imshow(ff1[220:280,220:280])
ax2.imshow(ff2[220:280,220:280])

ax3.imshow(ff3[220:280,220:280])

x_train, y_train,c = datagen(1,b=12)
plt.imshow(x_train[0,:,:])
px.imshow(np.sqrt(farfield(new_N, Nx, nf, phase_surf_0, wavelen)))
#fig2=px.imshow(np.sqrt(farfield(new_N, Nx, nf, phase_surf_pred, wavelen)))


model.save('model2')


ff_filename='ff2_002.tif'
FF = np.array(tif.imread(ff_filename))
FF=FF/np.max(FF)
FF[1667:1675,1235:1246]=0
[i_m,j_m]=np.where(FF==np.max(FF))
#print(i_m,j_m)\n",
i_m=int(np.mean(i_m))
j_m=int(np.mean(j_m))
d_cam=300
d_pix_ff=4.54
''' i_m=1056
j_m=1060  '''
FF = FF[i_m-d_cam:i_m+d_cam+1,j_m-d_cam:j_m+d_cam+1]
x_ff_real = np.linspace(-d_cam*d_pix_ff,d_cam*d_pix_ff,2*d_cam+1)
x_ff=ff_coord(0.91,1,0.00834/2,50,501)
d_x_ff=x_ff[1]-x_ff[0]
ff_new=interp_ff(FF,x_ff_real,x_ff)
x_out=[]
ff_new1=ff_new[ff_new.shape[0]//2-20:ff_new.shape[0]//2+20,ff_new.shape[0]//2-20:ff_new.shape[0]//2+20]
ff_new1/=np.max(ff_new1)
#ff=(ff>0.05)
PolarImage, r_grid, theta_grid = abel.tools.polar.reproject_image_into_polar(ff_new1)
x_out.append(PolarImage)
x_out=np.array(x_out)

fig,[ax1,ax2]=plt.subplots(1,2)
ax1.imshow(FF)
ax2.imshow(x_out[0,:,:])
y_pred=model.predict(x_train)

phase_surf_pred = zern_fun(wavelen, Nx, Ny, R, zern_pred)
phase_surf_pred = phase_surf_pred * (nf > 1)

print('RMS_pred', np.std(phase_surf_pred) * 1e6)


plt.figure()
plt.imshow(phase_surf_pred,clim=(-0.2*1e-6, 0.2*1e-6))



ff2=farfield(new_N, Nx, nf, phase_surf_pred, wavelen)

plt.figure()
plt.imshow(ff2[220:280,220:280])<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.2017419628947741, 'fraction_numerical': 0.06677524429967427, 'mean_word_length': 2.7591163162097416, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '46804', 'license': 'MIT', 'n_tokens_mistral': 7288, 'n_tokens_neox': 6804, 'n_words': 1199, 'path': 'nr_zrn.ipynb', 'repo_name': 'kitiv/DeepFaceLive', 'size': '14122'}