From b7d7695faf261e5c5ad0e0a7977b2f40b64ef643 Mon Sep 17 00:00:00 2001 From: jeremias Date: Mon, 29 Jun 2020 12:38:24 +0200 Subject: [PATCH] asd --- codes/CS.py | 940 ++++++++++------------ codes/Graphics.py | 610 ++++++-------- codes/MRI.py | 191 +---- codes/PostCheck.py | 567 +++++++++---- codes/__pycache__/CS.cpython-36.pyc | Bin 15426 -> 13607 bytes codes/__pycache__/Graphics.cpython-36.pyc | Bin 36990 -> 32550 bytes codes/monolithic.py | 101 ++- input/Graphics_input.yaml | 59 +- input/PostCheck_input.yaml | 101 ++- input/cfd_input.yaml | 50 -- input/mri_input.yaml | 116 ++- 11 files changed, 1297 insertions(+), 1438 deletions(-) delete mode 100755 input/cfd_input.yaml diff --git a/codes/CS.py b/codes/CS.py index 1d0d2bd..85a1a1c 100644 --- a/codes/CS.py +++ b/codes/CS.py @@ -14,487 +14,403 @@ rank = comm.Get_rank() # to allow MR imaging scientists to design undersampled # acquisitions and reconstruct the resulting data with CS without # needing to be a CS expert. -# +# # The Cartesian reconstruction is based on the split Bregman # code written by Tom Goldstein, originally available here: # -def pdf(k,kw,klo,q): - - p = (np.abs(k)/kw)**(-q) - p[np.where(k==0)] = 0 - p[np.where(np.abs(k)<=kw)] = 1 - p[np.where(k=norm: + P = pdf(ks, k+1, klo, q) + if np.sum(P) >= norm: break - - P = np.fft.fftshift(P) - - - #if np.mod(n,2)!=0: - # P = np.concatenate(([1],P),axis=None) - - + P = np.fft.fftshift(P) + return P -def mask_pdf_2d(dims,norm,q,pf): - - nz = dims[1] - ny = dims[0] - yc = round(ny/2) - zc = round(nz/2) - rmax = np.sqrt((ny-yc)**2 + (nz-zc)**2) - [Z,Y] = np.meshgrid(np.arange(0,nz),np.arange(0,ny)) - RR = np.sqrt( (Y-yc)**2 + (Z-zc)**2) - Z = np.abs(Z - nz/2 - 0.5) - Y = np.abs(Y - ny/2 - 0.5) - - - for rw in range(1,int(rmax)+1): - P = np.ones([ny,nz])/pf - C = np.logical_and( Z <= rw , Y <= rw) - W = np.logical_or( Z > rw , Y > rw) - P[W] = (RR[W]/rw)**(-q) - if np.sum(P) >= norm: - break - - - - return [P,C] +def mask_pdf_2d(dims, norm, q, pf): + + nz = dims[1] + ny = dims[0] + yc = round(ny/2) + zc = round(nz/2) + rmax = np.sqrt((ny-yc)**2 + (nz-zc)**2) + [Z, Y] = np.meshgrid(np.arange(0, nz), np.arange(0, ny)) + RR = np.sqrt((Y-yc)**2 + (Z-zc)**2) + Z = np.abs(Z - nz/2 - 0.5) + Y = np.abs(Y - ny/2 - 0.5) + + for rw in range(1, int(rmax)+1): + P = np.ones([ny, nz])/pf + C = np.logical_and(Z <= rw, Y <= rw) + W = np.logical_or(Z > rw, Y > rw) + P[W] = (RR[W]/rw)**(-q) + if np.sum(P) >= norm: + break + + return [P, C] + +def GeneratePattern(dim, R): -def GeneratePattern(dim,R): - # 3D CASE - if np.size(dim)==3: - - nro = dim[0] - npe = dim[1] - nacq = round(npe/R) - q = 1 - pf = 1 - P = mask_pdf_1d(npe, nacq, q, pf) + if np.size(dim) == 3: + + nro = dim[0] + npe = dim[1] + nacq = round(npe/R) + q = 1 + pf = 1 + P = mask_pdf_1d(npe, nacq, q, pf) while True: - M = np.random.rand(npe) - M = 1*(M<=P) - if np.sum(M)==nacq: + M = np.random.rand(npe) + M = 1*(M <= P) + if np.sum(M) == nacq: break # remove partial Fourier plane and compensate sampling density - M = M!=0 - M = np.tile(M,[nro,1]); + M = M != 0 + M = np.tile(M, [nro, 1]) #M = M.T # 4D CASE - if np.size(dim)==4: - - nro = dim[0] - npe1 = dim[1] - npe2 = dim[2] - nacq = round(npe1*npe2/R) - q = 1 - pf = 1 - [P,C] = mask_pdf_2d([npe1,npe2], nacq, q, pf) - RR = np.random.rand(npe1,npe2) - M = (RR <= P) - nchosen = np.sum(M) + if np.size(dim) == 4: + + nro = dim[0] + npe1 = dim[1] + npe2 = dim[2] + nacq = round(npe1*npe2/R) + q = 1 + pf = 1 + [P, C] = mask_pdf_2d([npe1, npe2], nacq, q, pf) + RR = np.random.rand(npe1, npe2) + M = (RR <= P) + nchosen = np.sum(M) if nchosen > nacq: # Correct for inexact number chosen #outerOn = np.logical_and( M , P!=1 ) - outerOn = np.where((M)*(P!=1)) - numToFlip = nchosen-nacq - idxs = np.random.permutation(outerOn[0].size) - idxx = outerOn[0][idxs[0:numToFlip]] - idxy = outerOn[1][idxs[0:numToFlip]] - M[idxx,idxy] = False - + outerOn = np.where((M)*(P != 1)) + numToFlip = nchosen-nacq + idxs = np.random.permutation(outerOn[0].size) + idxx = outerOn[0][idxs[0:numToFlip]] + idxy = outerOn[1][idxs[0:numToFlip]] + M[idxx, idxy] = False + elif nchosen < nacq: - outerOff = np.where(~M) - idxs = np.random.permutation(outerOff[0].size) - numToFlip = nacq - nchosen - idxx = outerOff[0][idxs[0:numToFlip]] - idxy = outerOff[1][idxs[0:numToFlip]] - M[idxx,idxy] = True - - - - - M = np.rollaxis(np.tile(np.rollaxis(M,1),[nro,1,1]),2) - M = np.fft.ifftshift(M) - M = M.transpose((1,0,2)) - + outerOff = np.where(~M) + idxs = np.random.permutation(outerOff[0].size) + numToFlip = nacq - nchosen + idxx = outerOff[0][idxs[0:numToFlip]] + idxy = outerOff[1][idxs[0:numToFlip]] + M[idxx, idxy] = True + + M = np.rollaxis(np.tile(np.rollaxis(M, 1), [nro, 1, 1]), 2) + M = np.fft.ifftshift(M) + M = M.transpose((1, 0, 2)) return M -def get_norm_factor(MASK,uu): - UM = MASK==1 - return UM.shape[0]/LA.norm(uu) +def get_norm_factor(MASK, uu): + UM = MASK == 1 + return UM.shape[0]/LA.norm(uu) def Dxyzt(X): - - if np.ndim(X)==3: - dd0 = X[:,:,0] - dd1 = X[:,:,1] - DA = dd0 - np.vstack((dd0[1::,:],dd0[0,:])) - DB = dd1 - np.hstack((dd1[:,1::],dd1[:,0:1])) - + + if np.ndim(X) == 3: + dd0 = X[:, :, 0] + dd1 = X[:, :, 1] + DA = dd0 - np.vstack((dd0[1::, :], dd0[0, :])) + DB = dd1 - np.hstack((dd1[:, 1::], dd1[:, 0:1])) + return DA + DB - if np.ndim(X)==4: - dd0 = X[:,:,:,0] - dd1 = X[:,:,:,1] - dd2 = X[:,:,:,2] - - DA = dd0 - np.vstack((dd0[1::,:,:],dd0[0,:,:][np.newaxis,:,:])) - DB = dd1 - np.hstack((dd1[:,1::,:],dd1[:,0,:][:,np.newaxis,:])) - DC = dd2 - np.dstack((dd2[:,:,1::],dd2[:,:,0][:,:,np.newaxis])) - - return DA + DB + DC + if np.ndim(X) == 4: + dd0 = X[:, :, :, 0] + dd1 = X[:, :, :, 1] + dd2 = X[:, :, :, 2] + + DA = dd0 - np.vstack((dd0[1::, :, :], dd0[0, :, :][np.newaxis, :, :])) + DB = dd1 - np.hstack((dd1[:, 1::, :], dd1[:, 0, :][:, np.newaxis, :])) + DC = dd2 - np.dstack((dd2[:, :, 1::], dd2[:, :, 0][:, :, np.newaxis])) + + return DA + DB + DC def Dxyz(u): - - if np.ndim(u)==2: - - dx = u[:,:]- np.vstack((u[-1,:],u[0:-1,:])) - dy = u[:,:]- np.hstack((u[:,-1:],u[:,0:-1])) - D = np.zeros([dx.shape[0],dx.shape[1],2],dtype=complex) - D[:,:,0] = dx - D[:,:,1] = dy - + + if np.ndim(u) == 2: + + dx = u[:, :] - np.vstack((u[-1, :], u[0:-1, :])) + dy = u[:, :] - np.hstack((u[:, -1:], u[:, 0:-1])) + D = np.zeros([dx.shape[0], dx.shape[1], 2], dtype=complex) + D[:, :, 0] = dx + D[:, :, 1] = dy + return D - if np.ndim(u)==3: - - dx = u[:,:,:]- np.vstack((u[-1,:,:][np.newaxis,:,:],u[0:-1,:,:])) - dy = u[:,:,:]- np.hstack((u[:,-1,:][:,np.newaxis,:],u[:,0:-1,:])) - dz = u[:,:,:]- np.dstack((u[:,:,-1][:,:,np.newaxis],u[:,:,0:-1])) - - D = np.zeros([dx.shape[0],dx.shape[1],dx.shape[2],3],dtype=complex) - - D[:,:,:,0] = dx - D[:,:,:,1] = dy - D[:,:,:,2] = dz - + if np.ndim(u) == 3: + + dx = u[:, :, :] - \ + np.vstack((u[-1, :, :][np.newaxis, :, :], u[0:-1, :, :])) + dy = u[:, :, :] - \ + np.hstack((u[:, -1, :][:, np.newaxis, :], u[:, 0:-1, :])) + dz = u[:, :, :] - \ + np.dstack((u[:, :, -1][:, :, np.newaxis], u[:, :, 0:-1])) + + D = np.zeros([dx.shape[0], dx.shape[1], dx.shape[2], 3], dtype=complex) + + D[:, :, :, 0] = dx + D[:, :, :, 1] = dy + D[:, :, :, 2] = dz + return D -def shrink(X,pgam): - p = 1 - s = np.abs(X) - tt = pgam/(s)**(1-p) +def shrink(X, pgam): + p = 1 + s = np.abs(X) + tt = pgam/(s)**(1-p) # t = pgam/np.sqrt(s) - ss = s-tt - ss = ss*(ss>0) - s = s + 1*(s 0) + s = s + 1*(s < tt) + ss = ss/s return ss*X -def CSMETHOD(ITOT,R): - - ''' Compressed Function. +def CSMETHOD(ITOT, R): + ''' Compressed Sensing Function. Args: ITOT: a numpy matrix with the full sampled (3D or 4D) dynamical data R: the acceleration factor ''' # Method parameters - ninner = 5 - nbreg = 10 - lmbda = 4 - mu = 20 - gam = 1 - - - if np.ndim(ITOT)==3: - [row,col,numt2] = ITOT.shape - elif np.ndim(ITOT)==4: - [row,col,dep,numt2] = ITOT.shape + ninner = 5 + nbreg = 10 + lmbda = 4 + mu = 20 + gam = 1 + + if np.ndim(ITOT) == 3: + [row, col, numt2] = ITOT.shape + elif np.ndim(ITOT) == 4: + [row, col, dep, numt2] = ITOT.shape else: raise Exception('Dynamical data is requested') - - - MASK = GeneratePattern(ITOT.shape,R) - CS1 = np.zeros(ITOT.shape,dtype=complex) + MASK = GeneratePattern(ITOT.shape, R) + CS1 = np.zeros(ITOT.shape, dtype=complex) + + nit = 0 + nit_tot = (numt2-1)/20 + + if np.ndim(ITOT) == 3: - nit = 0 - nit_tot = (numt2-1)/20 - - - if np.ndim(ITOT)==3: - for t in range(numt2): - - if rank==0: - print('{3D COMPRESSED SENSING} t = ',t) - - Kdata = np.fft.fft2(ITOT[:,:,t])*MASK - - data_ndims = Kdata.ndim - mask = Kdata!=0 # not perfect, but good enough - # normalize the data so that standard parameter values work - norm_factor = get_norm_factor(mask, Kdata) - Kdata = Kdata*norm_factor - # Reserve memory for the auxillary variables - Kdata0 = Kdata - img = np.zeros([row,col],dtype=complex) - X = np.zeros([row,col, data_ndims]) - B = np.zeros([row,col, data_ndims]) - # Build Kernels - scale = np.sqrt(row*col) - murf = np.fft.ifft2(mu*mask*Kdata)*scale - uker = np.zeros([row,col]) - uker[0,0] = 4 - uker[0,1] = -1 ; uker[1,0] = -1 - uker[-1,0] = -1 ; uker[0,-1] = -1 - uker = 1/(mu*mask + lmbda*np.fft.fftn(uker) + gam) + if rank == 0: + print('{3D COMPRESSED SENSING} t = ', t) + + Kdata = np.fft.fft2(ITOT[:, :, t])*MASK + + data_ndims = Kdata.ndim + mask = Kdata != 0 # not perfect, but good enough + # normalize the data so that standard parameter values work + norm_factor = get_norm_factor(mask, Kdata) + Kdata = Kdata*norm_factor + # Reserve memory for the auxillary variables + Kdata0 = Kdata + img = np.zeros([row, col], dtype=complex) + X = np.zeros([row, col, data_ndims]) + B = np.zeros([row, col, data_ndims]) + # Build Kernels + scale = np.sqrt(row*col) + murf = np.fft.ifft2(mu*mask*Kdata)*scale + uker = np.zeros([row, col]) + uker[0, 0] = 4 + uker[0, 1] = -1 + uker[1, 0] = -1 + uker[-1, 0] = -1 + uker[0, -1] = -1 + + uker = 1/(mu*mask + lmbda*np.fft.fftn(uker) + gam) # Do the reconstruction for outer in range(nbreg): for inner in range(ninner): # update u - rhs = murf + lmbda*Dxyzt(X-B) + gam*img - img = np.fft.ifft2(np.fft.fft2(rhs)*uker) + rhs = murf + lmbda*Dxyzt(X-B) + gam*img + img = np.fft.ifft2(np.fft.fft2(rhs)*uker) # update x and y - A = Dxyz(img) + B - X = shrink(A, 1/lmbda) + A = Dxyz(img) + B + X = shrink(A, 1/lmbda) # update bregman parameters - B = A - X - Kdata = Kdata + Kdata0 - mask*np.fft.fftn(img)/scale - murf = np.fft.ifftn(mu*mask*Kdata)*scale - + B = A - X + Kdata = Kdata + Kdata0 - mask*np.fft.fftn(img)/scale + murf = np.fft.ifftn(mu*mask*Kdata)*scale + # undo the normalization so that results are scaled properly img = img / norm_factor / scale - - - CS1[:,:,t] = img - - if np.ndim(ITOT)==4: - + + CS1[:, :, t] = img + + if np.ndim(ITOT) == 4: + for t in range(numt2): - - - if rank==0: - print('[4D CS] R = {re} t = {te}/{tef}'.format(re=R,te=t,tef=numt2)) - #if np.mod(t,nit_tot)<1: - # sys.stdout.write('\r') - # # Progress bar - # if numt2==3: - # sys.stdout.write("{4d-CS} [%-6s] %d%%" % ('=='*nit, 100*t/numt2)) - # else: - # sys.stdout.write("{4d-CS} [%-40s] %d%%" % ('=='*nit, 100*t/numt2)) - # sys.stdout.flush() - # nit = nit +1 - + if rank == 0: + print( + '[4D CS] R = {re} t = {te}/{tef}'.format(re=R, te=t, tef=numt2)) - Kdata_0 = np.fft.fftn(ITOT[:,:,:,t]) - Kdata = Kdata_0*MASK - data_ndims = Kdata.ndim - mask = Kdata!=0 # not perfect, but good enough + Kdata_0 = np.fft.fftn(ITOT[:, :, :, t]) + Kdata = Kdata_0*MASK + data_ndims = Kdata.ndim + mask = Kdata != 0 # not perfect, but good enough # normalize the data so that standard parameter values work - norm_factor = get_norm_factor(mask, Kdata) - Kdata = Kdata*norm_factor + norm_factor = get_norm_factor(mask, Kdata) + Kdata = Kdata*norm_factor # Reserve memory for the auxillary variables - Kdata0 = Kdata - img = np.zeros([row,col,dep],dtype=complex) - X = np.zeros([row,col,dep, data_ndims]) - B = np.zeros([row,col,dep, data_ndims]) + Kdata0 = Kdata + img = np.zeros([row, col, dep], dtype=complex) + X = np.zeros([row, col, dep, data_ndims]) + B = np.zeros([row, col, dep, data_ndims]) # Build Kernels - scale = np.sqrt(row*col*dep) - murf = np.fft.ifftn(mu*mask*Kdata)*scale - uker = np.zeros([row,col,dep]) - uker[0,0,0] = 8 - uker[1,0,0] = -1 ; uker[0,1,0] = -1 ; uker[0,0,1] = -1 - uker[-1,0,0] = -1 ; uker[0,-1,0] = -1 ; uker[0,0,-1] = -1 - uker = 1/(mu*mask + lmbda*np.fft.fftn(uker) + gam) + scale = np.sqrt(row*col*dep) + murf = np.fft.ifftn(mu*mask*Kdata)*scale + uker = np.zeros([row, col, dep]) + uker[0, 0, 0] = 8 + uker[1, 0, 0] = -1 + uker[0, 1, 0] = -1 + uker[0, 0, 1] = -1 + uker[-1, 0, 0] = -1 + uker[0, -1, 0] = -1 + uker[0, 0, -1] = -1 + uker = 1/(mu*mask + lmbda*np.fft.fftn(uker) + gam) # Do the reconstruction for outer in range(nbreg): for inner in range(ninner): # update u - rhs = murf + lmbda*Dxyzt(X-B) + gam*img - img = np.fft.ifft2(np.fft.fft2(rhs)*uker) + rhs = murf + lmbda*Dxyzt(X-B) + gam*img + img = np.fft.ifft2(np.fft.fft2(rhs)*uker) # update x and y - A = Dxyz(img) + B - X = shrink(A, 1/lmbda) + A = Dxyz(img) + B + X = shrink(A, 1/lmbda) # update bregman parameters - B = A - X - Kdata = Kdata + Kdata0 - mask*np.fft.fftn(img)/scale - murf = np.fft.ifftn(mu*mask*Kdata)*scale - + B = A - X + Kdata = Kdata + Kdata0 - mask*np.fft.fftn(img)/scale + murf = np.fft.ifftn(mu*mask*Kdata)*scale + # undo the normalization so that results are scaled properly img = img / norm_factor / scale - - - CS1[:,:,:,t] = img - - - - + + CS1[:, :, :, t] = img + return CS1 -def CSMETHOD_SENSE(ITOT,R,R_SENSE): - ''' Compressed Function. +def CSMETHOD_SENSE(ITOT, R, R_SENSE): + ''' Compressed sense algorith with SENSE... in contruction!. Args: ITOT: a numpy matrix with the full sampled (3D or 4D) dynamical data R: the acceleration factor ''' # Method parameters - ninner = 5 - nbreg = 10 - lmbda = 4 - mu = 20 - gam = 1 + ninner = 5 + nbreg = 10 + lmbda = 4 + mu = 20 + gam = 1 - [row,col,dep,numt2] = ITOT.shape - MASK = {} - ITOTCS = {} - MASK[0] = GeneratePattern([row,int(np.ceil(col/2)),dep,numt2],R) - MASK[1] = GeneratePattern([row,int(np.ceil(col/2)),dep,numt2],R) - SenseMAP = {} - [SenseMAP[0],SenseMAP[1]] = Sensitivity_Map([row,col,dep]) + [row, col, dep, numt2] = ITOT.shape + MASK = {} + ITOTCS = {} + MASK[0] = GeneratePattern([row, int(np.ceil(col/2)), dep, numt2], R) + MASK[1] = GeneratePattern([row, int(np.ceil(col/2)), dep, numt2], R) + SenseMAP = {} + [SenseMAP[0], SenseMAP[1]] = Sensitivity_Map([row, col, dep]) + + col = int(np.ceil(col/2)) + ITOTCS[0] = np.zeros([row, col, dep, numt2], dtype=complex) + ITOTCS[1] = np.zeros([row, col, dep, numt2], dtype=complex) - col = int(np.ceil(col/2)) - ITOTCS[0] = np.zeros([row,col,dep,numt2],dtype=complex) - ITOTCS[1] = np.zeros([row,col,dep,numt2],dtype=complex) - for rs in range(R_SENSE): for t in range(numt2): - if rank==0: - print('[4D CS] R = {re} t = {te}/{tef}'.format(re=R,te=t,tef=numt2)) + if rank == 0: + print( + '[4D CS] R = {re} t = {te}/{tef}'.format(re=R, te=t, tef=numt2)) - Kdata_0 = np.fft.fftn(ITOT[:,:,:,t]) - Kdata_0 = Kdata_0*SenseMAP[rs] - Kdata_0 = Kdata_0[:,0::R_SENSE,:] - - Kdata = Kdata_0*MASK[rs] - data_ndims = Kdata.ndim - mask = Kdata!=0 # not perfect, but good enough + Kdata_0 = np.fft.fftn(ITOT[:, :, :, t]) + Kdata_0 = Kdata_0*SenseMAP[rs] + Kdata_0 = Kdata_0[:, 0::R_SENSE, :] + + Kdata = Kdata_0*MASK[rs] + data_ndims = Kdata.ndim + mask = Kdata != 0 # not perfect, but good enough # normalize the data so that standard parameter values work - norm_factor = get_norm_factor(mask, Kdata) - Kdata = Kdata*norm_factor + norm_factor = get_norm_factor(mask, Kdata) + Kdata = Kdata*norm_factor # Reserve memory for the auxillary variables - Kdata0 = Kdata - img = np.zeros([row,col,dep],dtype=complex) - X = np.zeros([row,col,dep, data_ndims]) - B = np.zeros([row,col,dep, data_ndims]) + Kdata0 = Kdata + img = np.zeros([row, col, dep], dtype=complex) + X = np.zeros([row, col, dep, data_ndims]) + B = np.zeros([row, col, dep, data_ndims]) # Build Kernels - scale = np.sqrt(row*col*dep) - murf = np.fft.ifftn(mu*mask*Kdata)*scale - uker = np.zeros([row,col,dep]) - uker[0,0,0] = 8 - uker[1,0,0] = -1 ; uker[0,1,0] = -1 ; uker[0,0,1] = -1 - uker[-1,0,0] = -1 ; uker[0,-1,0] = -1 ; uker[0,0,-1] = -1 - uker = 1/(mu*mask + lmbda*np.fft.fftn(uker) + gam) + scale = np.sqrt(row*col*dep) + murf = np.fft.ifftn(mu*mask*Kdata)*scale + uker = np.zeros([row, col, dep]) + uker[0, 0, 0] = 8 + uker[1, 0, 0] = -1 + uker[0, 1, 0] = -1 + uker[0, 0, 1] = -1 + uker[-1, 0, 0] = -1 + uker[0, -1, 0] = -1 + uker[0, 0, -1] = -1 + uker = 1/(mu*mask + lmbda*np.fft.fftn(uker) + gam) # Do the reconstruction for outer in range(nbreg): for inner in range(ninner): # update u - rhs = murf + lmbda*Dxyzt(X-B) + gam*img - img = np.fft.ifft2(np.fft.fft2(rhs)*uker) + rhs = murf + lmbda*Dxyzt(X-B) + gam*img + img = np.fft.ifft2(np.fft.fft2(rhs)*uker) # update x and y - A = Dxyz(img) + B - X = shrink(A, 1/lmbda) + A = Dxyz(img) + B + X = shrink(A, 1/lmbda) # update bregman parameters - B = A - X - Kdata = Kdata + Kdata0 - mask*np.fft.fftn(img)/scale - murf = np.fft.ifftn(mu*mask*Kdata)*scale + B = A - X + Kdata = Kdata + Kdata0 - mask*np.fft.fftn(img)/scale + murf = np.fft.ifftn(mu*mask*Kdata)*scale # undo the normalization so that results are scaled properly img = img / norm_factor / scale - ITOTCS[rs][:,:,:,t] = img + ITOTCS[rs][:, :, :, t] = img - - return [ITOTCS[0],ITOTCS[1]] - -def CSMETHOD_peaksystole(ITOT,R,tstar): - - ''' Compressed Function. + return [ITOTCS[0], ITOTCS[1]] - Args: - ITOT: a numpy matrix with the full sampled (3D or 4D) dynamical data - R: the acceleration factor - tstar: the time when the flux in the inlet it's maximum - ''' - # Method parameters - ninner = 5 - nbreg = 10 - lmbda = 4 - mu = 20 - gam = 1 - - - - [row,col,dep,numt2] = ITOT.shape - - MASK = GeneratePattern(ITOT.shape,R) - CS1 = np.zeros([row,col,dep],dtype=complex) - - for t in range(tstar,tstar+1): - Kdata = np.fft.fftn(ITOT[:,:,:,t])*MASK - data_ndims = Kdata.ndim - mask = Kdata!=0 # not perfect, but good enough - # normalize the data so that standard parameter values work - norm_factor = get_norm_factor(mask, Kdata) - Kdata = Kdata*norm_factor - # Reserve memory for the auxillary variables - Kdata0 = Kdata - img = np.zeros([row,col,dep],dtype=complex) - X = np.zeros([row,col,dep, data_ndims]) - B = np.zeros([row,col,dep, data_ndims]) - # Build Kernels - scale = np.sqrt(row*col*dep) - murf = np.fft.ifftn(mu*mask*Kdata)*scale - uker = np.zeros([row,col,dep]) - uker[0,0,0] = 8 - uker[1,0,0] = -1 ; uker[0,1,0] = -1 ; uker[0,0,1] = -1 - uker[-1,0,0] = -1 ; uker[0,-1,0] = -1 ; uker[0,0,-1] = -1 - uker = 1/(mu*mask + lmbda*np.fft.fftn(uker) + gam) - - # Do the reconstruction - for outer in range(nbreg): - for inner in range(ninner): - # update u - rhs = murf + lmbda*Dxyzt(X-B) + gam*img - img = np.fft.ifft2(np.fft.fft2(rhs)*uker) - # update x and y - A = Dxyz(img) + B - X = shrink(A, 1/lmbda) - # update bregman parameters - B = A - X - Kdata = Kdata + Kdata0 - mask*np.fft.fftn(img)/scale - murf = np.fft.ifftn(mu*mask*Kdata)*scale - - # undo the normalization so that results are scaled properly - img = img / norm_factor / scale - CS1[:,:,:] = img - - - return CS1 - -def phase_contrast(M1,M0,VENC,scantype='0G'): - param = 1 - if scantype=='-G+G': - param = 0.5 +def phase_contrast(M1, M0, VENC, scantype='0G'): + param = 1 + if scantype == '-G+G': + param = 0.5 return VENC*param*(np.angle(M1) - np.angle(M0))/np.pi - -def GenerateMagnetization(Sq,VENC,noise,scantype='0G'): + +def GenerateMagnetization(Sq, VENC, noise, scantype='0G'): + ''' Simulation of a typical magnetization. A x-dependent plane is added into the + reference phase. + ''' # MRI PARAMETERS gamma = 267.513e6 # rad/Tesla/sec Gyromagnetic ratio for H nuclei B0 = 1.5 # Tesla Magnetic Field Strenght @@ -504,202 +420,178 @@ def GenerateMagnetization(Sq,VENC,noise,scantype='0G'): RHO0 = np.zeros(Sq.shape, dtype=complex) RHO1 = np.zeros(Sq.shape, dtype=complex) - if np.ndim(Sq)==3: - [row,col,numt2] = Sq.shape - [X,Y] = np.meshgrid(np.linspace(0,col,col),np.linspace(0,row,row)) + if np.ndim(Sq) == 3: + [row, col, numt2] = Sq.shape + [X, Y] = np.meshgrid(np.linspace(0, col, col), + np.linspace(0, row, row)) for k in range(numt2): - if noise: - Drho = np.random.normal(0,0.2,[row,col]) - Drho2 = np.random.normal(0,0.2,[row,col]) - else: - Drho = np.zeros([row,col]) - Drho2 = np.zeros([row,col]) - - varPHASE0 = np.random.randint(-10,11,size=(row,col))*np.pi/180*(np.abs(Sq[:,:,k])<0.001) #Hugo's observation - modulus = 0.5 + 0.5*(np.abs(Sq[:,:,k])>0.001) - - if scantype=='0G': - PHASE0[:,:,k] = (gamma*B0*TE+0.01*X)*(np.abs(Sq[:,:,k])>0.001) + 10*varPHASE0 - PHASE1[:,:,k] = (gamma*B0*TE+0.01*X)*(np.abs(Sq[:,:,k])>0.001) + 10*varPHASE0 + np.pi*Sq[:,:,k]/VENC - - if scantype=='-G+G': - PHASE0[:,:,k] = gamma*B0*TE*np.ones([row,col]) + 10*varPHASE0 - np.pi*Sq[:,:,k]/VENC - PHASE1[:,:,k] = gamma*B0*TE*np.ones([row,col]) + 10*varPHASE0 + np.pi*Sq[:,:,k]/VENC - - RHO0[:,:,k] = modulus*np.cos(PHASE0[:,:,k]) + Drho + 1j*modulus*np.sin(PHASE0[:,:,k]) + 1j*Drho2 - RHO1[:,:,k] = modulus*np.cos(PHASE1[:,:,k]) + Drho + 1j*modulus*np.sin(PHASE1[:,:,k]) + 1j*Drho2 + if noise: + Drho = np.random.normal(0, 0.2, [row, col]) + Drho2 = np.random.normal(0, 0.2, [row, col]) + else: + Drho = np.zeros([row, col]) + Drho2 = np.zeros([row, col]) - - if np.ndim(Sq)==4: + varPHASE0 = np.random.randint(-10, 11, size=(row, col))*np.pi/180*( + np.abs(Sq[:, :, k]) < 0.001) # Hugo's observation + modulus = 0.5 + 0.5*(np.abs(Sq[:, :, k]) > 0.001) + + if scantype == '0G': + PHASE0[:, :, k] = (gamma*B0*TE+0.01*X) * \ + (np.abs(Sq[:, :, k]) > 0.001) + 10*varPHASE0 + PHASE1[:, :, k] = (gamma*B0*TE+0.01*X)*(np.abs(Sq[:, :, k]) + > 0.001) + 10*varPHASE0 + np.pi*Sq[:, :, k]/VENC + + if scantype == '-G+G': + PHASE0[:, :, k] = gamma*B0*TE * \ + np.ones([row, col]) + 10*varPHASE0 - np.pi*Sq[:, :, k]/VENC + PHASE1[:, :, k] = gamma*B0*TE * \ + np.ones([row, col]) + 10*varPHASE0 + np.pi*Sq[:, :, k]/VENC + + RHO0[:, :, k] = modulus*np.cos(PHASE0[:, :, k]) + \ + Drho + 1j*modulus*np.sin(PHASE0[:, :, k]) + 1j*Drho2 + RHO1[:, :, k] = modulus*np.cos(PHASE1[:, :, k]) + \ + Drho + 1j*modulus*np.sin(PHASE1[:, :, k]) + 1j*Drho2 + + if np.ndim(Sq) == 4: [row, col, dep, numt2] = Sq.shape [X, Y, Z] = np.meshgrid(np.linspace(0, col, col), np.linspace( 0, row, row), np.linspace(0, dep, dep)) for k in range(numt2): - if noise: - Drho = np.random.normal(0, 0.2, [row, col, dep]) - Drho2 = np.random.normal(0, 0.2, [row, col, dep]) - else: - Drho = np.zeros([row, col, dep]) - Drho2 = np.zeros([row, col, dep]) + if noise: + Drho = np.random.normal(0, 0.2, [row, col, dep]) + Drho2 = np.random.normal(0, 0.2, [row, col, dep]) + else: + Drho = np.zeros([row, col, dep]) + Drho2 = np.zeros([row, col, dep]) - varPHASE0 = np.random.randint(-10, 11, size=(row, col, dep)) * \ - np.pi/180*(np.abs(Sq[:, :, :, k]) < 0.001) - modulus = 0.5 + 0.5*(np.abs(Sq[:, :, :, k]) > 0.001) + varPHASE0 = np.random.randint(-10, 11, size=(row, col, dep)) * \ + np.pi/180*(np.abs(Sq[:, :, :, k]) < 0.001) + modulus = 0.5 + 0.5*(np.abs(Sq[:, :, :, k]) > 0.001) - if scantype == '0G': - PHASE0[:, :, :, k] = (gamma*B0*TE+0.01*X) * \ - (np.abs(Sq[:, :, :, k]) > 0.001) + 10*varPHASE0 - PHASE1[:, :, :, k] = (gamma*B0*TE+0.01*X)*(np.abs(Sq[:, :, :, k]) - > 0.001) + 10*varPHASE0 + np.pi*Sq[:, :, :, k]/VENC + if scantype == '0G': + PHASE0[:, :, :, k] = (gamma*B0*TE+0.01*X) * \ + (np.abs(Sq[:, :, :, k]) > 0.001) + 10*varPHASE0 + PHASE1[:, :, :, k] = (gamma*B0*TE+0.01*X)*(np.abs(Sq[:, :, :, k]) + > 0.001) + 10*varPHASE0 + np.pi*Sq[:, :, :, k]/VENC - if scantype == '-G+G': - PHASE0[:, :, :, k] = gamma*B0*TE * \ - np.ones([row, col, dep]) + varPHASE0 - np.pi*Sq[:, :, :, k]/VENC - PHASE1[:, :, :, k] = gamma*B0*TE * \ - np.ones([row, col, dep]) + varPHASE0 + np.pi*Sq[:, :, :, k]/VENC + if scantype == '-G+G': + PHASE0[:, :, :, k] = gamma*B0*TE * \ + np.ones([row, col, dep]) + varPHASE0 - \ + np.pi*Sq[:, :, :, k]/VENC + PHASE1[:, :, :, k] = gamma*B0*TE * \ + np.ones([row, col, dep]) + varPHASE0 + \ + np.pi*Sq[:, :, :, k]/VENC - RHO0[:, :, :, k] = modulus*np.cos(PHASE0[:, :, :, k]) + \ - Drho + 1j*modulus*np.sin(PHASE0[:, :, :, k]) + 1j*Drho2 - RHO1[:, :, :, k] = modulus*np.cos(PHASE1[:, :, :, k]) + \ - Drho + 1j*modulus*np.sin(PHASE1[:, :, :, k]) + 1j*Drho2 + RHO0[:, :, :, k] = modulus*np.cos(PHASE0[:, :, :, k]) + \ + Drho + 1j*modulus*np.sin(PHASE0[:, :, :, k]) + 1j*Drho2 + RHO1[:, :, :, k] = modulus*np.cos(PHASE1[:, :, :, k]) + \ + Drho + 1j*modulus*np.sin(PHASE1[:, :, :, k]) + 1j*Drho2 + return [RHO0, RHO1] +def undersampling(Sqx, Sqy, Sqz, options, savepath): - return [RHO0,RHO1] + R = options['cs']['R'] -def undersampling(Sqx,Sqy,Sqz,options,savepath): - - R = options['cs']['R'] - for r in R: - - if rank==0: + + if rank == 0: print('Using Acceleration Factor R = ' + str(r)) print('Component x of M0') - - [M0,M1] = GenerateMagnetization(Sqx,options['cs']['VENC'],options['cs']['noise']) - + + [M0, M1] = GenerateMagnetization( + Sqx, options['cs']['VENC'], options['cs']['noise']) print('\n Component x of M0') - M0_cs = CSMETHOD(M0,r) + M0_cs = CSMETHOD(M0, r) print('\n Component x of M1') - M1_cs = CSMETHOD(M1,r) + M1_cs = CSMETHOD(M1, r) - Sqx_cs = phase_contrast(M1_cs,M0_cs,options['cs']['VENC']) - del M0,M1 + Sqx_cs = phase_contrast(M1_cs, M0_cs, options['cs']['VENC']) + del M0, M1 del M0_cs, M1_cs - - [M0,M1] = GenerateMagnetization(Sqy,options['cs']['VENC'],options['cs']['noise']) - + + [M0, M1] = GenerateMagnetization( + Sqy, options['cs']['VENC'], options['cs']['noise']) print('\n Component y of M0') - M0_cs = CSMETHOD(M0,r) + M0_cs = CSMETHOD(M0, r) print('\n Component y of M1') - M1_cs = CSMETHOD(M1,r) - + M1_cs = CSMETHOD(M1, r) - Sqy_cs = phase_contrast(M1_cs,M0_cs,options['cs']['VENC']) - - del M0,M1 + Sqy_cs = phase_contrast(M1_cs, M0_cs, options['cs']['VENC']) + + del M0, M1 del M0_cs, M1_cs - - [M0,M1] = GenerateMagnetization(Sqz,options['cs']['VENC'],options['cs']['noise']) - - if rank==0: + + [M0, M1] = GenerateMagnetization( + Sqz, options['cs']['VENC'], options['cs']['noise']) + + if rank == 0: print('\n Component z of M0') - M0_cs = CSMETHOD(M0,r) - if rank==0: + M0_cs = CSMETHOD(M0, r) + if rank == 0: print('\n Component z of M1') - M1_cs = CSMETHOD(M1,r) - if rank==0: + M1_cs = CSMETHOD(M1, r) + if rank == 0: print(' ') - - Sqz_cs = phase_contrast(M1_cs,M0_cs,options['cs']['VENC']) - - - if rank==0: + + Sqz_cs = phase_contrast(M1_cs, M0_cs, options['cs']['VENC']) + + if rank == 0: print('saving the sequences in ' + savepath) - seqname = options['cs']['name'] +'_R' + str(r) + '.npz' - print('sequence name: ' + seqname) - np.savez_compressed( savepath + seqname, x=Sqx_cs, y=Sqy_cs,z=Sqz_cs) - - del Sqx_cs,Sqy_cs,Sqz_cs - -def undersampling_peakpv(Sqx,Sqy,Sqz,options,R): - - Sqx_cs = {} - Sqy_cs = {} - Sqz_cs = {} - [Mx0,Mx1] = GenerateMagnetization(Sqx,options['cs']['VENC'],options['cs']['noise'],scantype='0G') - [My0,My1] = GenerateMagnetization(Sqy,options['cs']['VENC'],options['cs']['noise'],scantype='0G') - [Mz0,Mz1] = GenerateMagnetization(Sqz,options['cs']['VENC'],options['cs']['noise'],scantype='0G') - + seqname = options['cs']['name'] + '_R' + str(r) + '.npz' + print('sequence name: ' + seqname) + np.savez_compressed(savepath + seqname, + x=Sqx_cs, y=Sqy_cs, z=Sqz_cs) - Mx0_cs = CSMETHOD(Mx0,R) - Mx1_cs = CSMETHOD(Mx1,R) - My0_cs = CSMETHOD(My0,R) - My1_cs = CSMETHOD(My1,R) - Mz0_cs = CSMETHOD(Mz0,R) - Mz1_cs = CSMETHOD(Mz1,R) + del Sqx_cs, Sqy_cs, Sqz_cs - Sqx_cs = phase_contrast(Mx1_cs,Mx0_cs,options['cs']['VENC'],scantype='0G') - Sqy_cs = phase_contrast(My1_cs,My0_cs,options['cs']['VENC'],scantype='0G') - Sqz_cs = phase_contrast(Mz1_cs,Mz0_cs,options['cs']['VENC'],scantype='0G') - - - - return [Sqx_cs,Sqy_cs,Sqz_cs] - -def undersampling_short(Mx,My,Mz,options): - - R = options['cs']['R'] - savepath = options['cs']['savepath'] +def undersampling_short(Mx, My, Mz, options): + R = options['cs']['R'] + savepath = options['cs']['savepath'] - R_SENSE = 1 + R_SENSE = 1 if 'R_SENSE' in options['cs']: - R_SENSE = options['cs']['R_SENSE'][0] + R_SENSE = options['cs']['R_SENSE'][0] - for r in R: - if rank==0: + for r in R: + if rank == 0: print('Using Acceleration Factor R = ' + str(r)) - - if R_SENSE==2: - [MxS0_cs,MxS1_cs] = CSMETHOD_SENSE(Mx,r,2) - [MyS0_cs,MyS1_cs] = CSMETHOD_SENSE(My,r,2) - [MzS0_cs,MzS1_cs] = CSMETHOD_SENSE(Mz,r,2) - if rank==0: + if R_SENSE == 2: + [MxS0_cs, MxS1_cs] = CSMETHOD_SENSE(Mx, r, 2) + [MyS0_cs, MyS1_cs] = CSMETHOD_SENSE(My, r, 2) + [MzS0_cs, MzS1_cs] = CSMETHOD_SENSE(Mz, r, 2) + if rank == 0: print('saving the sequences in ' + savepath) - seqname_s0 = options['cs']['name'] +'S0_R' + str(r) + '.npz' - seqname_s1 = options['cs']['name'] +'S1_R' + str(r) + '.npz' - print('sequence name: ' + seqname_s0) - np.savez_compressed( savepath + seqname_s0, x=MxS0_cs, y=MyS0_cs,z=MzS0_cs) - print('sequence name: ' + seqname_s1) - np.savez_compressed( savepath + seqname_s1, x=MxS1_cs, y=MyS1_cs,z=MzS1_cs) + seqname_s0 = options['cs']['name'] + 'S0_R' + str(r) + '.npz' + seqname_s1 = options['cs']['name'] + 'S1_R' + str(r) + '.npz' + print('sequence name: ' + seqname_s0) + np.savez_compressed(savepath + seqname_s0, + x=MxS0_cs, y=MyS0_cs, z=MzS0_cs) + print('sequence name: ' + seqname_s1) + np.savez_compressed(savepath + seqname_s1, + x=MxS1_cs, y=MyS1_cs, z=MzS1_cs) del MxS0_cs, MyS0_cs, MzS0_cs del MxS1_cs, MyS1_cs, MzS1_cs - elif R_SENSE==1: - Mx_cs = CSMETHOD(Mx,r) - My_cs = CSMETHOD(My,r) - Mz_cs = CSMETHOD(Mz,r) - if rank==0: + elif R_SENSE == 1: + Mx_cs = CSMETHOD(Mx, r) + My_cs = CSMETHOD(My, r) + Mz_cs = CSMETHOD(Mz, r) + if rank == 0: print('saving the sequences in ' + savepath) - seqname = options['cs']['name'] +'_R' + str(r) + '.npz' - print('sequence name: ' + seqname) - np.savez_compressed( savepath + seqname, x=Mx_cs, y=My_cs,z=Mz_cs) - del Mx_cs,My_cs,Mz_cs + seqname = options['cs']['name'] + '_R' + str(r) + '.npz' + print('sequence name: ' + seqname) + np.savez_compressed(savepath + seqname, + x=Mx_cs, y=My_cs, z=Mz_cs) + del Mx_cs, My_cs, Mz_cs else: raise Exception('Only implemented for 2-fold SENSE!!') - - - - # THE END - - - diff --git a/codes/Graphics.py b/codes/Graphics.py index 91d2685..89bf6ed 100644 --- a/codes/Graphics.py +++ b/codes/Graphics.py @@ -586,25 +586,15 @@ def CenterComparison(A,B,C,R,center): def VelocityChannel(M,repeat,recorder): [row,col,numt2] = M.shape - - [X,Y] = np.meshgrid(np.linspace(0,col,col),np.linspace(0,row,row)) - plt.ion() - rown = row*6/row coln = col*6/row - fig = plt.figure()#figsize=(4, 6) , dpi=200) ax = fig.add_subplot(111, projection='3d') - - for rr in range(-1,repeat): - for t in range(numt2): - V = M[:,:,t] - #ax.plot_surface(X, Y, V, cmap=plt.cm.magma, vmin=-30, vmax=50, linewidth=0, antialiased=False) #ax.plot_wireframe(X, Y, V,rcount=20,ccount=20,linewidth=0.5) ax.plot_surface(X, Y, V,cmap='magma',vmin=-30, vmax=100, linewidth=0, antialiased=False) @@ -615,7 +605,6 @@ def VelocityChannel(M,repeat,recorder): ax.set_xlabel('$x$') ax.set_ylabel('$y$') ax.set_zlabel('$v(r)$') - ax.set_title('phase ' + str(t)) plt.pause(0.001) plt.draw() @@ -634,116 +623,7 @@ def PlotTri(tri,pos,p): #ax.plot_trisurf(pos[:,0], pos[:,1], pos[:,2], triangles=tri.simplices, cmap=plt.cm.Spectral,linewidth=0.1,edgecolors='k') #ax.tripcolor(pos[:,0], pos[:,1], pos[:,2], triangles=tri.simplices, facecolors=p2.T, edgecolor='black') plt.show() - -def PlotPressureDrop(mode): - - barye2mmHg = 1/1333.22387415 - - CT = np.loadtxt('Pressure/DROPS/pd_NS_coarse.txt') - CT2 = np.loadtxt('Pressure/DROPS/pd_NS_coarse2.txt') - - ref_ao = np.loadtxt('Pressure/DROPS2/refPPE1_coarse.txt') - ref = np.loadtxt('Pressure/DROPS2/refPPE1_coarse_leo1.txt') - CS2 = np.loadtxt('Pressure/DROPS2/CSPPE2_coarse_leo1.txt') - - - PPE_CT = np.loadtxt('Pressure/DROPS/pd_PPE_NS_coarse.txt') - PPE_CT_leo = np.loadtxt('Pressure/DROPS/pd_PPE_NS_coarse_leo1.txt') - - - STE_CT = np.loadtxt('Pressure/DROPS/test_STE_R1_coarse.txt') - test_STE_CT = np.loadtxt('Pressure/DROPS/test2_STE_R1_coarse_leo1.txt') - - STEint_CT = np.loadtxt('Pressure/DROPS/test_STEint_R1_coarse.txt') - test_STEint_CT = np.loadtxt('Pressure/DROPS/test2_STEint_R1_coarse_leo1.txt') - - - # UNDERSAMPLING - PPE_R1_leo1 = np.loadtxt('Pressure/DROPS/test_PPE_R1_coarse_leo1.txt') - - PPE_KT_R2_leo1 = np.loadtxt('Pressure/DROPS/KT_PPE_R2_coarse_leo1.txt') - PPE_KT_R4_leo1 = np.loadtxt('Pressure/DROPS/KT_PPE_R4_coarse_leo1.txt') - PPE_KT_R8_leo1 = np.loadtxt('Pressure/DROPS/KT_PPE_R8_coarse_leo1.txt') - PPE_KT_R12_leo1 = np.loadtxt('Pressure/DROPS/KT_PPE_R12_coarse_leo1.txt') - PPE_KT_R20_leo1 = np.loadtxt('Pressure/DROPS/KT_PPE_R20_coarse_leo1.txt') - - PPE_CS_R2_leo1 = np.loadtxt('Pressure/DROPS/CS_PPE_R2_coarse_leo1.txt') - PPE_CS_R4_leo1 = np.loadtxt('Pressure/DROPS/CS_PPE_R4_coarse_leo1.txt') - PPE_CS_R8_leo1 = np.loadtxt('Pressure/DROPS/CS_PPE_R8_coarse_leo1.txt') - PPE_CS_R12_leo1 = np.loadtxt('Pressure/DROPS/CS_PPE_R12_coarse_leo1.txt') - PPE_CS_R20_leo1 = np.loadtxt('Pressure/DROPS/CS_PPE_R20_coarse_leo1.txt') - #CT_fine = np.loadtxt('Pressure/DROPS/pd_NS_fine.txt') - #CT_fine_T = np.loadtxt('Pressure/DROPS/pd_NS_fine_T.txt') - #PPE_CT_fine = np.loadtxt('Pressure/DROPS/pd_PPE_NS_fine.txt') - #PPE_CT_fine_leo3 = np.loadtxt('Pressure/DROPS/pd_PPE_NS_fine_leo3.txt') - - - tvec2 = np.linspace(0, 2.5 ,PPE_CT_leo.size) - tvec = np.linspace(tvec2[1]*0.5 ,2.5+tvec2[1]*0.5 , CT.size) - #tvec_T = np.linspace(0,2.5,CT_fine_T.size) - - fig = plt.figure() - - if mode=='KT': - plt.plot(tvec,CT,'-k',linewidth=2,label='$ref$') - plt.plot(tvec2,PPE_R1_leo1,'xkcd:red',linewidth=2,linestyle='-' , label='$R = 1$') - plt.plot(tvec2,PPE_KT_R2_leo1,'xkcd:blue',linewidth=2,linestyle='-' ,label='$R = 2$') - plt.plot(tvec2,PPE_KT_R4_leo1,'xkcd:green',linewidth=2,linestyle='-',label='$R = 4$') - plt.plot(tvec2,PPE_KT_R8_leo1,'xkcd:orange',linewidth=2,linestyle='-',label='$R = 8$') - plt.plot(tvec2,PPE_KT_R20_leo1,'xkcd:magenta',linewidth=2,linestyle='-', label='$R = 20$') - plt.title('$kt-BLAST$',fontsize=20) - - if mode=='CS': - plt.plot(tvec,CT,'-k',linewidth=2,label='$ref$') - plt.plot(tvec2,ref_ao,'xkcd:red',linewidth=2,linestyle='--' , label='$aorta$') - plt.plot(tvec2,ref,'xkcd:red',linewidth=2,linestyle='-' , label='$leo$') - plt.plot(tvec2,CS2,'xkcd:blue',linewidth=2,linestyle='-' ,label='$R = 2$') - - #plt.plot(tvec2,PPE_R1_leo1,'xkcd:red',linewidth=2,linestyle='-' , label='$R = 1$') - #plt.plot(tvec2,PPE_CS_R2_leo1,'xkcd:blue',linewidth=2,linestyle='-' ,label='$R = 2$') - #plt.plot(tvec2,PPE_CS_R4_leo1,'xkcd:green',linewidth=2,linestyle='-',label='$R = 4$') - #plt.plot(tvec2,PPE_CS_R8_leo1,'xkcd:orange',linewidth=2,linestyle='-',label='$R = 8$') - #plt.plot(tvec2,PPE_CS_R20_leo1,'xkcd:magenta',linewidth=2,linestyle='-', label='$R = 20$') - plt.title('$Compressed \ \ Sensing$',fontsize=20) - - - - if mode=='STE': - plt.plot(tvec,CT,'-k',linewidth=2,label='$ref$') - plt.plot(tvec2,STE_CT , 'xkcd:purple' , label='STE-CT aorta') - plt.plot(tvec2,test_STE_CT , 'xkcd:purple' , linestyle='--', marker='o', label='STE-CT leo') - - - if mode=='STEint': - plt.plot(tvec,CT,'-k',linewidth=2,label='$ref$') - plt.plot(tvec2,STEint_CT, 'xkcd:aquamarine', label='STEint-CT aorta') - plt.plot(tvec2,test_STEint_CT, 'xkcd:aquamarine', linestyle='--', marker='o',label='STEint-CT aorta') - - - plt.ylim([-2,7]) - plt.xlabel(r'$time \ \ \ (s)$',fontsize=20) - plt.ylabel(r'$\delta p \ \ \ (mmHg) $',fontsize=20) - plt.legend(fontsize=16) - ############################################################################################################################## - - - #fig = plt.figure() - #plt.plot(tvec2,CT_fine,'-ok',label='CT') - ##plt.plot(tvec_T,CT_fine_T,'--k',label='CT T') - #plt.plot(tvec2 ,PPE_CT_fine ,'-om',label='PPE-CT aorta') - #plt.plot(tvec2 ,PPE_CT_fine_leo3,'-oc',label='PPE-CT leo3') - #plt.title('aorta fine') - #plt.xlabel(r'$time \ \ (s) $',fontsize=20) - #plt.ylabel(r'$\delta p \ \ \ (mmHg) $',fontsize=20) - #plt.legend(fontsize=14) - - - - - - plt.show() - def Plot_flux(masterpath,meshpath,options,mode,R): mesh = Mesh() @@ -1090,8 +970,7 @@ def Plot_dP(masterpath,options,mode,R): tcat[l-2] = float(row[0]) tcat = tcat+shift_t ax.plot(tcat[cstar[0]:tcat.size-cstar[1]],catheter[cstar[0]:tcat.size-cstar[1]],'white',linewidth=linesize,linestyle='--',label='$catheter$') - - + tm = np.linspace(0,Dt*tend,tend) ax.set_xlim([-0.05,0.81]) @@ -1131,226 +1010,7 @@ def Plot_dP(masterpath,options,mode,R): plt.show() - -def Plot_peaksystole(datapath,options,meshes,dt,R): - import pickle - barye2mmHg = 1/1333.22387415 - - - for mesh_size in meshes: - t_star = 6 - - - PPE_MEAN = np.zeros([len(R)]) - PPE_STD = np.zeros([len(R)]) - STE_MEAN = np.zeros([len(R)]) - STE_STD = np.zeros([len(R)]) - V_MEAN = np.zeros([len(R)]) - V_STD = np.zeros([len(R)]) - - ref_P = 0 - ref_V = 0 - - - - - if mesh_size=='Ucoarse': - ref_V = 326.95828118309191 - if mesh_size=='Ufine': - ref_V = 232.95021682714497 - if mesh_size=='Uffine': - ref_V = 234.66445211879045 - - - - - for l in range(len(R)): - if R[l]==0: - ref = np.loadtxt('/home/yeye/N_MRI/codes/pressure_drop/'+mesh_size+'/dt' + str(dt) + '/ref_'+mesh_size+'.txt') - PPE0_raw = open('/home/yeye/N_MRI/codes/pressure_drop/'+mesh_size+'/dt' + str(dt) + '/R0/pdrop_PPE_impl_stan.dat','rb') - STE0_raw = open('/home/yeye/N_MRI/codes/pressure_drop/'+mesh_size+'/dt' + str(dt) + '/R0/pdrop_STE_impl_stan.dat','rb') - PPE0 = pickle.load(PPE0_raw)['pdrop']*(-barye2mmHg) - STE0 = pickle.load(STE0_raw)['pdrop']*(-barye2mmHg) - - curpath = datapath + 'sequences/aorta_'+mesh_size+'.npz' - p = np.load(curpath) - px = p['x'] - py = p['y'] - pz = p['z'] - v = np.sqrt(px[:,:,:,t_star]**2 + py[:,:,:,t_star]**2 + pz[:,:,:,t_star]**2) - max0 = np.where(v==np.max(v)) - - ref_P = ref[6] - V_MEAN[l] = ref_V - V_STD[l] = 0 - PPE_MEAN[l] = PPE0[6] - PPE_STD[l] = 0 - STE_MEAN[l] = STE0[6] - STE_STD[l] = 0 - else: - PPE_MEAN[l] = np.loadtxt(datapath + 'maxpv/'+mesh_size + '/ppemean_R'+ str(R[l]) +'.txt') - PPE_STD[l] = np.loadtxt(datapath + 'maxpv/'+mesh_size + '/ppestd_R'+ str(R[l]) +'.txt') - STE_MEAN[l] = np.loadtxt(datapath + 'maxpv/'+mesh_size + '/stemean_R'+ str(R[l]) +'.txt') - STE_STD[l] = np.loadtxt(datapath + 'maxpv/'+mesh_size + '/stestd_R'+ str(R[l]) +'.txt') - V_MEAN[l] = np.loadtxt(datapath + 'maxpv/'+mesh_size + '/vmean_R'+ str(R[l]) +'.txt') - V_STD[l] = np.loadtxt(datapath + 'maxpv/'+mesh_size + '/vstd_R'+ str(R[l]) +'.txt') - - - - plt.figure(figsize=(10, 6), dpi=100) - Rvec = np.linspace(-10,R[-1]+5,100) - hline = Rvec*0+ref_P - plt.subplot(1,2,1) - plt.plot([0],[ref_P],color='k',marker='o',label= '$ref$') - plt.plot(Rvec,hline,color='k',linestyle='--') - plt.errorbar(R,PPE_MEAN, yerr=PPE_STD,color=options['ppecol'],marker='o',label= '$PPE$') - plt.errorbar(R,STE_MEAN, yerr=STE_STD,color=options['stecol'],marker='o',label= '$STE$') - plt.xlim([-2,R[-1]+5]) - plt.ylim([-7,18]) - plt.xlabel(r'$R$',fontsize=20) - plt.title('$Peak \ Systole \ pressure$',fontsize=18) - plt.ylabel(r'$ max \ \big ( \delta P \big ) \ \ mmHg$',fontsize=16) - - - plt.subplot(1,2,2) - #plt.plot([0],[ref_V],color='k',marker='o',label= '$ref$') - Rvec = np.linspace(-10,R[-1]+5,100) - hline = Rvec*0+ref_V - plt.plot(Rvec,hline,color='k',linestyle='--') - plt.errorbar(R,V_MEAN, yerr=V_STD,color='royalblue',marker='o',label= '$' + mesh_size + '$') - plt.xlim([-2,R[-1]+5]) - plt.ylim([0,350]) - plt.xlabel(r'$R$',fontsize=20) - plt.ylabel(r'$ max \ \big ( v \big ) \ \ cm/s$',fontsize=16) - plt.title('$Peak \ Systole \ velocity$',fontsize=18) - #plt.legend(fontsize=15) - - plt.annotate('$'+mesh_size+'$', xy=(-0.2, 1.1), xycoords='axes fraction',fontsize=15) - - - plt.show() - -def Plot_peaksystole_flux(datapath,options,meshes,dt,R): - import pickle - barye2mmHg = 1/1333.22387415 - - - for mesh_size in meshes: - t_star = 6 - PPE_MEAN = np.zeros([len(R)]) - PPE_STD = np.zeros([len(R)]) - STE_MEAN = np.zeros([len(R)]) - STE_STD = np.zeros([len(R)]) - V_MEAN = np.zeros([len(R)]) - V_STD = np.zeros([len(R)]) - Q_MEAN = np.zeros([len(R)]) - Q_STD = np.zeros([len(R)]) - - ref_V = 0 - - if mesh_size=='Ucoarse': - ref_V = 326.95828118309191 - if mesh_size=='Ufine': - ref_V = 232.95021682714497 - if 'Uffine' in mesh_size: - ref_V = 226.93523675462458 - maxv1 = 184.05091675316763 - ref_Q = 454.77517472437495 - maxq1 = 429.01393994253556 - - - for l in range(len(R)): - if R[l]==0: - ref = np.loadtxt('/home/yeye/N_MRI/codes/pressure_drop/'+mesh_size+'/dt' + str(dt) + '/ref_'+mesh_size+'.txt') - PPE0_raw = open('/home/yeye/N_MRI/codes/pressure_drop/'+mesh_size+'/dt' + str(dt) + '/R0/pdrop_PPE_impl_stan.dat','rb') - STE0_raw = open('/home/yeye/N_MRI/codes/pressure_drop/'+mesh_size+'/dt' + str(dt) + '/R0/pdrop_STE_impl_stan.dat','rb') - PPE0 = pickle.load(PPE0_raw)['pdrop']*(-barye2mmHg) - STE0 = pickle.load(STE0_raw)['pdrop']*(-barye2mmHg) - - curpath = datapath + 'sequences/aorta_'+mesh_size+'.npz' - p = np.load(curpath) - px = p['x'] - py = p['y'] - pz = p['z'] - v = np.sqrt(px[:,:,:,t_star]**2 + py[:,:,:,t_star]**2 + pz[:,:,:,t_star]**2) - max0 = np.where(v==np.max(v)) - V_MEAN[l] = ref_V - V_STD[l] = 0 - PPE_MEAN[l] = PPE0[6] - PPE_STD[l] = 0 - STE_MEAN[l] = STE0[6] - STE_STD[l] = 0 - elif R[l]==1: - PPE_MEAN[l] = 0 - PPE_STD[l] = 0 - STE_MEAN[l] = 0 - STE_STD[l] = 0 - V_MEAN[l] = maxv1 - V_STD[l] = 0 - Q_MEAN[l] = maxq1 - Q_STD[l] = 0 - - else: - PPE_MEAN[l] = np.loadtxt(datapath + 'maxpv/'+mesh_size + '/ppemean_R'+ str(R[l]) +'.txt') - PPE_STD[l] = np.loadtxt(datapath + 'maxpv/'+mesh_size + '/ppestd_R'+ str(R[l]) +'.txt') - STE_MEAN[l] = np.loadtxt(datapath + 'maxpv/'+mesh_size + '/stemean_R'+ str(R[l]) +'.txt') - STE_STD[l] = np.loadtxt(datapath + 'maxpv/'+mesh_size + '/stestd_R'+ str(R[l]) +'.txt') - V_MEAN[l] = np.loadtxt(datapath + 'maxpv/'+mesh_size + '/vmean_R'+ str(R[l]) +'.txt') - V_STD[l] = np.loadtxt(datapath + 'maxpv/'+mesh_size + '/vstd_R'+ str(R[l]) +'.txt') - Q_MEAN[l] = -np.loadtxt(datapath + 'maxpv/'+mesh_size + '/Q2_R'+ str(R[l]) +'.txt') - Q_STD[l] = -np.loadtxt(datapath + 'maxpv/'+mesh_size + '/Qstd2_R'+ str(R[l]) +'.txt') - - - - plt.figure(figsize=(15, 5), dpi=100) - plt.annotate('$'+mesh_size+'$', xy=(-0.2, 1.1), xycoords='axes fraction',fontsize=15) - - - Rvec = np.linspace(-10,R[-1]+5,100) - plt.subplot(1,3,1) - #plt.plot(Rvec,hline,color='k',linestyle='--') - plt.errorbar(R,PPE_MEAN, yerr=PPE_STD,color=options['ppecol'],marker='o',label= '$PPE$') - plt.errorbar(R,STE_MEAN, yerr=STE_STD,color=options['stecol'],marker='o',label= '$STE$') - plt.xlim([-2,R[-1]+5]) - plt.ylim([-0.1,2.0]) - plt.xlabel(r'$R$',fontsize=20) - plt.legend(fontsize=12,loc=2) - plt.title('$ l_2 \ error \ in \ Peak \ pressure$',fontsize=16) - plt.xticks(R) - plt.ylabel(r'$ || p - p^{ref} ||/|| p^{ref} ||_{l2} $',fontsize=16) - - - plt.subplot(1,3,2) - #plt.plot([0],[ref_V],color='k',marker='o',label= '$ref$') - Rvec = np.linspace(-10,R[-1]+7,100) - hline = Rvec*0+ref_V - plt.plot(Rvec,hline,color='royalblue',linestyle='--') - plt.errorbar(R,V_MEAN, yerr=V_STD,color='royalblue',marker='o',label= '$' + mesh_size + '$') - plt.xlim([-2,R[-1]+5]) - plt.ylim([0,350]) - plt.xlabel(r'$R$',fontsize=20) - plt.ylabel(r'$ v \ \ cm/s$',fontsize=16) - plt.xticks(R) - plt.title('$Peak \ velocity$',fontsize=16) - #plt.legend(fontsize=15) - - plt.subplot(1,3,3) - Rvec = np.linspace(-10,R[-1]+7,100) - hline = Rvec*0+ref_Q - plt.plot(Rvec,hline,color='mediumvioletred',linestyle='--') - plt.errorbar(R,Q_MEAN, yerr=Q_STD,color='mediumvioletred',marker='o',label= '$' + mesh_size + '$') - plt.xlim([-2,R[-1]+5]) - plt.xlabel(r'$R$',fontsize=20) - plt.ylabel(r'$ Q \ \ ml/s$',fontsize=16) - plt.xticks(R) - plt.ylim([150,550]) - plt.title('$Peak \ Flux $',fontsize=16) - - - - - plt.show() - + def CLOCK(t1, t2): tot_time = np.round(t2 - t1, 2) if tot_time < 60: @@ -1411,6 +1071,9 @@ def ROUTINE(options): if options['Error-curves']['apply']: print('--- Error-curves analysis ---') ratio = False + fac=1 + if '11mm' in options['Error-curves']['subfolders']: + fac=100 for types in options['Error-curves']['type']: if types=='mean_ratio': types = 'mean' @@ -1418,6 +1081,9 @@ def ROUTINE(options): if types=='max_ratio': types = 'max' ratio = True + if types=='norm2_m': + types='norm2' + meas=True nc = 0 if len(options['Error-curves']['subfolders'])==0: ucomp = [] @@ -1445,46 +1111,262 @@ def ROUTINE(options): for subf in options['Error-curves']['subfolders']: ucomp = [] wcomp = [] - colorset = options['Error-curves']['colors'] + if 'colors' in options['Error-curves']: + print('colors setted...') + colorset = options['Error-curves']['colors'] + colorsetted = True + else: + colorsetted = False + styles = options['Error-curves']['styles'] labelset = options['Error-curves']['labels'] path = options['Error-curves']['folder'] + subf + '/' - try: - ucomp = np.loadtxt(path + '/u'+types+'.txt') - wcomp = np.loadtxt(path + '/w'+types+'.txt') - times = np.loadtxt(path + '/times.txt') - except IOError: - print('no cError-curves for ' + subf) + wcomp = np.loadtxt(path + 'w'+types+'.txt') + times = np.loadtxt(path + 'times.txt') + if types != 'grad': + ucomp = np.loadtxt(path + 'u'+types+'.txt') - if not ratio: - plt.plot( - times, ucomp, color=colorset[nc], linestyle='-', label= '$u'+ subf +'$' ) - plt.plot( - times, wcomp, color=colorset[nc], linestyle='--', label='$w'+subf+'$') + if colorsetted: + if not ratio: + if types not in ['grad']: + if meas: + plt.plot(times, fac*ucomp, color=colorset[nc], linestyle='--', label= '$ u \ '+ labelset[nc] +'$',linewidth=2.5 ) + plt.plot(times, fac*wcomp, color=colorset[nc], linestyle=styles[nc], label= '$ w \ '+ labelset[nc] +'$',linewidth=2.5) + else: + wu = wcomp/ucomp + plt.plot( + times, wu, color=colorset[nc], linestyle=styles[nc], label= '$'+ labelset[nc] +'$' ) + nc +=1 else: - wu = wcomp/ucomp - plt.plot( - times, wu, color=colorset[nc], linestyle='-', label= '$'+ labelset[nc] +'$' ) - nc +=1 + if not ratio: + if types not in ['grad']: + if meas: + plt.plot(times, fac*ucomp, linestyle='--', label= '$'+ labelset[nc] +'$' ) + else: + plt.plot(times, fac*ucomp, linestyle='--', label= '$'+ labelset[nc] +'$' ) + plt.plot( + times, fac*wcomp, linestyle=styles[nc], label= '$'+ labelset[nc] +'$') + else: + wu = wcomp/ucomp + plt.plot( + times, wu, linestyle=styles[nc], label= '$'+ labelset[nc] +'$' ) + nc +=1 + - #plt.ylim([0, 170]) plt.xlabel('$time \ \ (s)$', fontsize=18) - plt.legend(fontsize=16) + plt.legend(fontsize=14) if options['Error-curves']['title']: plt.title(options['Error-curves']['title'], fontsize=18) if not ratio: - plt.ylabel('$velocity \ \ (cm/s)$', fontsize=18) + if types == 'grad': + plt.ylim([0, 230]) + plt.ylabel('$ |grad(w)|_2 \ \ (1/s)$', fontsize=18) + elif types == 'norm2': + plt.ylim([0, 52]) + plt.ylabel('$ |w|_2 \ \ (cm/s)$', fontsize=18) + else: + plt.ylabel('$velocity \ \ (cm/s)$', fontsize=18) plt.savefig(options['Error-curves']['outpath'] + types + '.png', dpi=500, bbox_inches='tight') else: - plt.ylabel('$w/u$', fontsize=18) - if 'max' in types: - plt.ylim([0, 1.8]) + plt.ylabel('$|w/u|$', fontsize=18) + #if 'max' in types: + #plt.ylim([0, 1.8]) plt.savefig(options['Error-curves']['outpath'] + types + '_ratio.png', dpi=500, bbox_inches='tight') plt.show() + if 'Histograms_checkpoint' in options: + if options['Histograms_checkpoint']['apply']: + print('--- Histograms ---') + path = options['Histograms_checkpoint']['path'] + freq = np.loadtxt(path + 'hist_freq.txt') + edges = np.loadtxt(path + 'hist_edges.txt') + fig, ax = plt.subplots() + #ax.bar(edges[:-1], freq, width=np.diff(edges), edgecolor="black", align="edge") + ax.bar(edges[:-1], freq, width=np.diff(edges), align="edge") + plt.title(options['Histograms_checkpoint']['title'], fontsize=18) + plt.xlim([0 , 50]) + plt.ylim([0 , 1.8]) + plt.savefig(path + 'hist.png', dpi=500, bbox_inches='tight') + plt.show() + if 'Pressure_drops' in options: + if options['Pressure_drops']['apply']: + print('--- Pressure_drops ---') + import pickle + nc = 0 + tommhg = options['Pressure_drops']['convertor'] + + for subf in options['Pressure_drops']['subfolders']: + ucomp = [] + wcomp = [] + if 'colors' in options['Pressure_drops']: + colorset = options['Pressure_drops']['colors'] + colorsetted = True + else: + colorsetted = False + styles = options['Pressure_drops']['styles'] + labelset = options['Pressure_drops']['labels'] + path = options['Pressure_drops']['folder'] + subf + '/' + dataname = 'pdrop_COR_impl_stan.dat' + if 'STE' in path: + dataname = 'pdrop_STE_impl_stan.dat' + if labelset[nc]=='ref': + dataname = 'pdrop.dat' + data = open(path+dataname, 'rb') + p_drop = pickle.load(data)['pdrop']/tommhg + data = open(path+dataname, 'rb') + times = pickle.load(data)['time'] + + if labelset[nc] == 'ref': + plt.plot(times, p_drop,color='black', linestyle='-', label= '$ref$' ) + else: + + if colorsetted: + plt.plot( + times, p_drop, color=colorset[nc], linestyle=styles[nc], label= '$'+ labelset[nc] +'$' ) + else: + plt.plot(times, p_drop, linestyle=styles[nc], label= '$'+ labelset[nc] +'$' ) + + if options['Pressure_drops']['catheter']: + + c_path = '/home/yeye/Desktop/PhD/MEDICAL_DATA/Study_David/catheter_data/catheter_'+ labelset[nc]+'_rest_stats.csv' + + + with open(c_path, 'r') as csvfile: + mylist = [row[0] for row in csv.reader(csvfile, delimiter=';')] + + Values = np.array(mylist) + catheter = np.zeros([len(Values)-2]) + tcat = np.zeros([len(Values)-2]) + for l in range(2,len(Values)): + row = Values[l].split(',') + catheter[l-2] = float(row[5]) + tcat[l-2] = float(row[0]) + + + if '11mm' in subf: + tdelay = 0.015 + elif '9mm' in subf: + tdelay = -0.12 + elif '13mm' in subf: + tdelay = 0.1 + elif 'Normal' in subf: + tdelay = -0.01 + + plt.plot(tcat+tdelay,catheter,color=colorset[nc],linestyle='--')# ,label='$cat' + subf + '$') + + nc +=1 + + + #plt.ylim([0, 170]) + plt.xlabel('$time \ \ (s)$', fontsize=18) + plt.legend(fontsize=14) + if options['Pressure_drops']['title']: + plt.title(options['Pressure_drops']['title'], fontsize=18) + plt.ylabel('$ \delta P \ \ (mmHg)$', fontsize=18) + plt.savefig(options['Pressure_drops']['outpath'] + 'pressure_drops.png', dpi=500, bbox_inches='tight') + plt.show() + + if 'l2_comp' in options: + if options['l2_comp']['apply']: + print('--- L2 component analysis ---') + fig, ax = plt.subplots() + for subf in options['l2_comp']['subfolder']: + path = options['l2_comp']['folder'] + subf + '/' + colors = options['l2_comp']['colors'] + mode = options['l2_comp']['mode']['type'] + + if mode in ['gain','gain_compressed']: + gain = True + path_to_comp = options['l2_comp']['mode']['comp'] + wx = np.loadtxt(path_to_comp + '/wx.txt') + wy = np.loadtxt(path_to_comp + '/wy.txt') + wz = np.loadtxt(path_to_comp + '/wz.txt') + else: + gain = False + wx = np.loadtxt(path + '/wx.txt') + wy = np.loadtxt(path + '/wy.txt') + wz = np.loadtxt(path + '/wz.txt') + + times = np.loadtxt(path + '/times.txt') + if subf != 'SNRinfV120' and gain: + varux = np.loadtxt(path + '/varux.txt') + varuy = np.loadtxt(path + '/varuy.txt') + varuz = np.loadtxt(path + '/varuz.txt') + if 'SNRinfV120' in subf: + lsty = '--' + labels = ['','','','',''] + else: + lsty = '-' + lsty2 = '--' + labels = ['$wx$','$wy$','$wz$','$div$'] + labels2 = ['$\delta u_x$','$\delta u_y$','$\delta u_z$'] + labels3 = ['$x$','$y$','$z$'] + + if mode == 'gain': + plt.plot(times, varux, color = colors[0], linestyle=lsty2 , label= labels2[0] ) + plt.plot(times, varuy, color = colors[1], linestyle=lsty2 , label= labels2[1] ) + plt.plot(times, varuz, color = colors[2], linestyle=lsty2 , label= labels2[2] ) + plt.plot(times, wx, color = colors[0], linestyle=lsty, label= labels[0] ) + plt.plot(times, wy, color = colors[1], linestyle=lsty, label= labels[1] ) + plt.plot(times, wz, color = colors[2], linestyle=lsty, label= labels[2] ) + elif mode == 'gain_compressed': + plt.plot(times, varux-wx, color = colors[0], linestyle=lsty, label= labels3[0] ) + plt.plot(times, varuy-wy, color = colors[1], linestyle=lsty, label= labels3[1] ) + plt.plot(times, varuz-wz, color = colors[2], linestyle=lsty, label= labels3[2] ) + else: + plt.plot(times, wx, color = colors[0], linestyle=lsty, label= labels[0] ) + plt.plot(times, wy, color = colors[1], linestyle=lsty, label= labels[1] ) + plt.plot(times, wz, color = colors[2], linestyle=lsty, label= labels[2] ) + + plt.xlabel('$time \ \ (s)$', fontsize=18) + + if options['l2_comp']['div']: + div = np.loadtxt(path + 'div.txt') + div_rescaled = div*0.01 + div_rescaled = div_rescaled + 0.1 + plt.plot(times, div_rescaled, color = 'indigo', linestyle=lsty, label= labels[3] ) + + if 'title' in options['l2_comp']: + if options['l2_comp']['title']: + plt.title(options['l2_comp']['title'], fontsize=18) + + + + if options['l2_comp']['aliasing']: + if 'Poiseuille' in options['l2_comp']['folder']: + print('adding alaising color in Poiseuille') + import matplotlib.transforms as mtransforms + trans = mtransforms.blended_transform_factory(ax.transData, ax.transAxes) + if 'SNR10' in subf: + time_al = [0.15,0.78] + elif 'SNRinf' in subf: + time_al = [0.21,0.78] + pm_al = 0.5*(time_al[1] + time_al[0]) + r_al = 0.5*(time_al[1] - time_al[0]) + ax.fill_between(times, 0, 1, where=np.abs(times -pm_al)<= r_al ,facecolor='gold', alpha=0.4, transform=trans) + if mode == 'gain_compressed': + ax.text(pm_al/times[-1], 0.55, '$aliasing$', horizontalalignment='center',verticalalignment='center', transform=ax.transAxes,fontsize=17) + else: + ax.text(pm_al/times[-1], 0.82, '$aliasing$', horizontalalignment='center',verticalalignment='center', transform=ax.transAxes,fontsize=17) + + leg = plt.legend(fancybox=True,fontsize=16) + leg.get_frame().set_linewidth(0.0) + ax.tick_params(axis='both', which='major', labelsize=14) + #plt.ylim([-0.005 , 0.235]) + + if mode == 'gain_compressed': + plt.ylim([-0.25 , 1.75]) + plt.ylabel('$|| \delta u ||_{L2} - || w ||_{L2}$', fontsize=18) + if not gain: + plt.ylim([-0.005 , 0.5]) + plt.ylabel('$ \sqrt{\int w^2 dx /| \Omega|}/ venc$', fontsize=18) + + plt.savefig(options['l2_comp']['folder'] + options['l2_comp']['name'] + '.png', dpi=500, bbox_inches='tight') + plt.show() + diff --git a/codes/MRI.py b/codes/MRI.py index a1eab4b..9707027 100644 --- a/codes/MRI.py +++ b/codes/MRI.py @@ -3,13 +3,10 @@ # Workspace for MRI analysis of the 4Dflow data # # written by Jeremias Garay L: j.e.garay.labra@rug.nl -# Fernanda te amo # for autoreload in ipython3 # %load_ext autoreload # %autoreload 2 ################################################################ -import h5py -from dPdirectestim.dPdirectestim import * from dolfin import * import dolfin import numpy as np @@ -999,14 +996,8 @@ def CLOCK(rank,t1,t2): print('Total time: ' + str(time_hour) + ' hrs, ' + str(time_min) + ' min, ' + str(time_sec) + ' s') - def SCANNER(options): - ######################################## - # - # Basic Tools - # - ######################################## if 'kspace_cib' in options: if options['kspace_cib']['apply']: print('--- kspace from CIB data ---') @@ -1416,24 +1407,18 @@ def SCANNER(options): else: CIBtoH5(path_to_cib,times,dt,outpath,flip=flip) - ######################################## - # - # Undersampling - # - ######################################## if 'cs' in options: if options['cs']['apply']: if rank==0: print('Applying Compressed Sensing') - [Sqx, Sqy, Sqz] = LOADsequences(options['cs']['seqpath']) - import CS - if options['cs']['short']: - [Mx,My,Mz] = LOADsequences(options['cs']['Mpath']) - CS.undersampling_short(Mx,My,Mz,options) + if 'short' in options['cs']: + if options['cs']['short']: + [Mx,My,Mz] = LOADsequences(options['cs']['Mpath']) + CS.undersampling_short(Mx,My,Mz,options) else: CS.undersampling(Sqx,Sqy,Sqz,options,options['cs']['savepath']) @@ -1457,12 +1442,6 @@ def SCANNER(options): print('saving the sequences' + options['SENSE']['savepath']) np.savez_compressed(options['SENSE']['savepath'], x=MxS,y=MyS,z=MzS) - ######################################## - # - # Writing Checkpoint from Sequence - # - ######################################## - if 'create_checkpoint' in options: if options['create_checkpoint']['apply']: print('--- Create Checkpoint ---') @@ -1529,7 +1508,7 @@ def SCANNER(options): comm = MESH['mesh'].mpi_comm() dt = options['create_checkpoint']['dt'] if options['create_checkpoint']['xdmf']: - xdmf_u = XDMFFile(options['create_checkpoint']['savepath']+'u.xdmf') + xdmf_u = XDMFFile(options['create_checkpoint']['savepath']+ 'R' + str(r) + '/u.xdmf') for l in range(len(vel_seq)): if rank == 0: @@ -1545,11 +1524,6 @@ def SCANNER(options): inout.write_HDF5_data( comm, path + '/u.h5', vel_seq[l], '/u', t=l*dt) - ######################################## - # - # Relative Pressure Estimators - # - ######################################## if 'reference' in options: if options['reference']['apply']: if rank == 0: @@ -1688,157 +1662,6 @@ def SCANNER(options): if rank==0: print(' ') - if 'peak_pv' in options: - - if options['peak_pv']['apply']: - import CS - import pickle - import sys - import logging - # DPestim - logging.getLogger().setLevel(logging.INFO) - parameters['form_compiler']['optimize'] = True - parameters['form_compiler']['cpp_optimize'] = True - parameters['form_compiler']['cpp_optimize_flags'] = '-O3 -ffast-math -march=native' - infile_dp = options['peak_pv']['infile_dp'] - estimator = DPDirectEstim(infile_dp) - - barye2mmHg = 1/1333.22387415 - t_star = 0.185 - if rank==0: - print('Computing Velocity and Pressure at Peak Systole') - print('The inlet max occurs at ' + str(t_star) + ' sec') - - - [BOX,AORTA,LEO] = CREATEmeshes(options) - - if options['peak_pv']['mesh_into']=='leo': - MESH = LEO - del AORTA - if options['peak_pv']['mesh_into']=='aorta': - MESH = AORTA - del LEO - - if options['peak_pv']['p_comp']=='error': - if rank==0: - print('Reading Pressure Reference') - P0_PPE = READcheckpoint(MESH,'p',options,options['checkpoint_path'],'p_PPE_impl_stan') - P0_STE = READcheckpoint(MESH,'p',options,options['checkpoint_path'],'p_STE_impl_stan') - - [Sqx,Sqy,Sqz] = LOADsequences(options['peak_pv']['orig_seq']) - Nite = options['peak_pv']['N_CS'] - [row,col,dep,numt] = Sqz.shape - frms_num = 3 - peakslice = options['peak_pv']['peak_slice'] - R = options['peak_pv']['R'] - - v0 = np.sqrt(Sqx[:,:,:,peakslice]**2 + Sqy[:,:,:,peakslice]**2 + Sqz[:,:,:,peakslice]**2) - max0 = np.where(v0==np.max(v0)) - - qqvec = {} - for ss in options['peak_pv']['flux_bnd']: - qqvec[ss] = np.zeros([Nite]) - - ppemax = np.zeros([Nite]) - stemax = np.zeros([Nite]) - vmax = np.zeros([Nite]) - slrdmax = peakslice - # Selecting around the max only - slrdmax = 1 - Sqx = Sqx[:,:,:,peakslice-1:peakslice+2] - Sqy = Sqy[:,:,:,peakslice-1:peakslice+2] - Sqz = Sqz[:,:,:,peakslice-1:peakslice+2] - - for l in range(len(R)): - if rank==0: - print('Peak Systole velocity and pressure at R = ' + str(R[l])) - - sx_cs = np.zeros([row,col,dep,frms_num,Nite]) - sy_cs = np.zeros([row,col,dep,frms_num,Nite]) - sz_cs = np.zeros([row,col,dep,frms_num,Nite]) - - for k in range(Nite): - if rank==0: - print('CS iteration number ' + str(k+1)) - [xcs,ycs,zcs] = CS.undersampling_peakpv(Sqx,Sqy,Sqz,options,R[l]) - sx_cs[:,:,:,:,k] = xcs - sy_cs[:,:,:,:,k] = ycs - sz_cs[:,:,:,:,k] = zcs - vk = np.sqrt(sx_cs[:,:,:,1,k]**2 + sy_cs[:,:,:,1,k]**2 + sz_cs[:,:,:,1,k]**2) - vmax[k] = vk[max0] - - if rank==0: - print('\n CS done') - - # To write the checkpoints - vel_seq = SqtoH5(BOX,MESH,sx_cs[:,:,:,:,k],sy_cs[:,:,:,:,k],sz_cs[:,:,:,:,k]) - comm = MESH['mesh'].mpi_comm() - - # Computing the Fluxes - if rank==0: - print('\n Computing the Flux') - QQ = Fluxes(MESH,vel_seq,options,options['peak_pv']['flux_bnd']) - - for ss in options['peak_pv']['flux_bnd']: - qqvec[ss][k] = QQ[ss][slrdmax] - - if rank==0: - print('\n Writing checkpoints') - - for ns in range(len(vel_seq)): - pathss = options['peak_pv']['savepath'] + 'H5/checkpoint/{i}/'.format(i=ns) - if l<10 and l>0: - pathss = options['peak_pv']['savepath'] + 'H5/checkpoint/0{i}/'.format(i=ns) - inout.write_HDF5_data(comm, pathss + '/u.h5', vel_seq[ns], '/u', t=0) - if rank==0: - print('\n The checkpoints were wrote') - - # Computing the Pressure Drop - estimator.estimate() - # Reading the results - if options['peak_pv']['p_comp']=='peak': - ppe_raw = open(options['peak_pv']['savepath'] + '/H5/pdrop_PPE_impl_stan.dat','rb') - ste_raw = open(options['peak_pv']['savepath'] + '/H5/pdrop_STE_impl_stan.dat','rb') - ppe = pickle.load(ppe_raw)['pdrop']*(-barye2mmHg) - ste = pickle.load(ste_raw)['pdrop']*(-barye2mmHg) - p1max[k] = ppe[slrdmax] - p2max[k] = ste[slrdmax] - elif options['peak_pv']['p_comp']=='error': - PPE = READcheckpoint(MESH,'p',options,options['peak_pv']['savepath']+'H5/checkpoint/','p_PPE_impl_stan') - STE = READcheckpoint(MESH,'p',options,options['peak_pv']['savepath']+'H5/checkpoint/','p_STE_impl_stan') - ppe_vec_0 = P0_PPE[peakslice].vector().get_local() - P0_PPE[peakslice].vector().get_local()[0] - ppe_vec = PPE[slrdmax].vector().get_local() - PPE[slrdmax].vector().get_local()[0] - ste_vec_0 = P0_STE[peakslice].vector().get_local() - P0_STE[peakslice].vector().get_local()[0] - ste_vec = STE[slrdmax].vector().get_local() - STE[slrdmax].vector().get_local()[0] - ppemax[k] = np.linalg.norm(ppe_vec_0 - ppe_vec)/np.linalg.norm(ppe_vec_0) - stemax[k] = np.linalg.norm(ste_vec_0 - ste_vec)/np.linalg.norm(ste_vec_0) - else: - raise Exception('Pressure computation not recognize!') - - - # VELOCITIES - vmean = np.mean(vmax) - vstd = np.std(vmax) - # PRESSURES - ppemean = np.mean(ppemax) - stemean = np.mean(stemax) - ppestd = np.std(ppemax) - stestd = np.std(stemax) - - - if options['peak_pv']['save']: - if rank==0: - print('\n saving the files in ' + options['peak_pv']['savepath']) - for ss in options['peak_pv']['flux_bnd']: - np.savetxt( options['peak_pv']['savepath'] + 'Q'+str(ss) +'_R'+str(R[l])+'.txt', [np.mean(qqvec[ss])]) - np.savetxt( options['peak_pv']['savepath'] + 'Qstd'+str(ss) +'_R'+str(R[l])+'.txt', [np.std(qqvec[ss])]) - np.savetxt( options['peak_pv']['savepath'] + 'ppemean_R'+str(R[l])+'.txt', [ppemean] ) - np.savetxt( options['peak_pv']['savepath'] + 'stemean_R'+str(R[l])+'.txt', [stemean] ) - np.savetxt( options['peak_pv']['savepath'] + 'vmean_R'+str(R[l])+'.txt', [vmean] ) - np.savetxt( options['peak_pv']['savepath'] + 'ppestd_R'+str(R[l])+'.txt', [ppestd] ) - np.savetxt( options['peak_pv']['savepath'] + 'stestd_R'+str(R[l])+'.txt', [stestd] ) - np.savetxt( options['peak_pv']['savepath'] + 'vstd_R'+str(R[l])+'.txt', [vstd] ) - if 'change_mesh' in options: if options['change_mesh']['apply']: if rank == 0: @@ -1861,9 +1684,9 @@ def SCANNER(options): changed = {} #W = LEO['FEM'].sub(0).collapse() comm = MESH_out['mesh'].mpi_comm() - v2 = Function(MESH_out['FEM']) for k in range(len(list(origin))): + v2 = Function(MESH_out['FEM']) if rank == 0: print('CHANGING: index', k) LagrangeInterpolator.interpolate(v2, origin[k]) @@ -1878,7 +1701,7 @@ def SCANNER(options): if rank == 0: print('saving checkpoint', l) path = options['change_mesh']['savepath'] + \ - 'R1/checkpoint/{i}/'.format(i=l) + 'checkpoint/{i}/'.format(i=l) writepath = path + '/'+options['change_mesh']['mode']+'.h5' inout.write_HDF5_data( comm, writepath, changed[l] ,'/'+options['change_mesh']['mode'], t=l*dt) diff --git a/codes/PostCheck.py b/codes/PostCheck.py index f3d75e4..0187b37 100644 --- a/codes/PostCheck.py +++ b/codes/PostCheck.py @@ -109,14 +109,10 @@ def WORKcheck(MESH, mode, output_path, checkpoint_path, filename, outname, optio if mode == 'p' or mode == 'p_cib': xdmf_p = XDMFFile(output_path+outname+'.xdmf') - for k in range(0, len(indexes), 1): + dt = options['Pressure']['dt'] + for k in range(0, len(indexes), options['Pressure']['undersampling']): + te = k*dt path = checkpoint_path + str(indexes[k]) + '/'+filename+'.h5' - - if filename == 'p': - if k < 10 and k > 0: - path = checkpoint_path + '0' + \ - str(indexes[k]) + '/'+filename+'.h5' - p = Function(W) if mode == 'p': barye2mmHg = 1/1333.22387415 @@ -126,36 +122,10 @@ def WORKcheck(MESH, mode, output_path, checkpoint_path, filename, outname, optio hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') hdf.read(p, 'p/vector_0') hdf.close() - p1vec = p.vector().get_local() - p1vec = p1vec - np.mean(p1vec) + #p1vec = p1vec - np.mean(p1vec) p.vector()[:] = p1vec*barye2mmHg xdmf_p.write(p, te) - te = te + dt - numt = numt + 1 - - if mode == 'divu': - xdmf_u = XDMFFile(output_path+outname+'.xdmf') - for k in range(0, len(indexes), 1): - path = checkpoint_path + str(indexes[k]) + '/'+filename+'.h5' - - if indexes[k] < 10 and indexes[k] > 0: - path = checkpoint_path + '0' + \ - str(indexes[k]) + '/'+filename+'.h5' - - v = Function(V) - dv = Function(W) - - dv.rename('div', outname) - comm = MPI.COMM_WORLD - hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') - hdf.read(v, 'u/vector_0') - - hdf.close() - dv.assign(project(div(v), W)) - xdmf_u.write(dv, te) - te = te + dt - numt = numt + 1 def READcheckpoint(MESH, mode, output_path, checkpoint_path, filename, outname, options, flow=False, bnds=None): @@ -203,7 +173,7 @@ def READcheckpoint(MESH, mode, output_path, checkpoint_path, filename, outname, freq,edges = np.histogram(ValuesPeak, bins=80, density=True) #Saving the histogram print('Saving at ' + output_path) - np.savetxt(output_path + 'hist_frew.txt', freq) + np.savetxt(output_path + 'hist_freq.txt', freq) np.savetxt(output_path + 'hist_edges.txt', edges) if mode == 'perturbation': @@ -235,6 +205,20 @@ def READcheckpoint(MESH, mode, output_path, checkpoint_path, filename, outname, noise_in_coil = options['Perturbation']['type']['coil'] + umaxima = [] + for k in indexes: + path = checkpoint_path + str(k) + '/'+filename+'.h5' + hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') + hdf.read(u, 'u/vector_0') + time = hdf.attributes('u/vector_0').to_dict()['timestamp'] + hdf.close() + uvec = u.vector().get_local() + umaxima.append(np.max(uvec)) + + ufactor = options['Perturbation']['type']['phase_contrast']/100 + VENC = np.floor(np.ceil(np.max(umaxima))*ufactor) + print('VENC chosen = ',VENC) + for k in indexes: path = checkpoint_path + str(k) + '/'+filename+'.h5' hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') @@ -244,9 +228,7 @@ def READcheckpoint(MESH, mode, output_path, checkpoint_path, filename, outname, uvec = u.vector().get_local() if Phase_Contrast: - ufactor = options['Perturbation']['type']['phase_contrast']/100 - VENC = np.max(np.abs(uvec))*ufactor - gamma = 267.513e6 # rad/Tesla/sec Gyromagnetic ratio for H nuclei + #gamma = 267.513e6 # rad/Tesla/sec Gyromagnetic ratio for H nuclei B0 = 1.5 # Tesla Magnetic Field Strenght TE = 5e-3 # Echo-time Phi1 = 1*B0*TE + 0*uvec @@ -474,7 +456,11 @@ def ERRORmap(MESH, mode, outpath, reference_path, checkpoint_path, refname,comna from dolfin import HDF5File V = MESH['FEM'] + V1 = MESH['FEM'].sub(1).collapse() W = MESH['FEM'].sub(0).collapse() + DGs = FunctionSpace( MESH['mesh'], 'DG',0) + cellv = CellVolume(MESH['mesh']) + h = CellDiameter(MESH['mesh']) unsort_indexes = os.listdir(checkpoint_path) indexes = [int(x) for x in unsort_indexes] indexes.sort() @@ -486,6 +472,80 @@ def ERRORmap(MESH, mode, outpath, reference_path, checkpoint_path, refname,comna if len(indexes)!=len(indexes0): raise Exception('The lengh of the checkpoints are not the same!') + if mode == 'algebra': + outname = options['Algebra']['outname'] + output = XDMFFile(outpath+outname+'.xdmf') + checkpoint = options['Algebra']['checkpoint'] + v1 = Function(W) + v2 = Function(W) + vout = Function(W) + vout.rename('velocity','velocity') + vout_vec = np.zeros(vout.vector().get_local().size) + times_vec = np.zeros(len(indexes0)) + vdict1 = {} + vdict2 = {} + if options['Algebra']['mode'] == 'aliasing': + print('Assuming the corrector in the second path entered') + + # Reading all the timestamps first + for k in range(len(indexes)): + path_v1 = reference_path + str(indexes0[k]) + '/'+refname+'.h5' + path_v2 = checkpoint_path + str(indexes[k]) + '/'+comname+'.h5' + hdf_v1 = HDF5File(MESH['mesh'].mpi_comm(), path_v1, 'r') + + if 'w' in refname: + hdf_v1.read(v1, 'w/vector_0') + time = hdf_v1.attributes('w/vector_0').to_dict()['timestamp'] + else: + hdf_v1.read(v1, 'u/vector_0') + time = hdf_v1.attributes('u/vector_0').to_dict()['timestamp'] + times_vec[k] = time + hdf_v2 = HDF5File(MESH['mesh'].mpi_comm(), path_v2, 'r') + if 'w' in comname: + hdf_v2.read(v2, 'w/vector_0') + else: + hdf_v2.read(v2, 'u/vector_0') + + print('computing algebra for the time',time) + hdf_v1.close() + hdf_v2.close() + vdict1[k] = v1.vector().get_local() + vdict2[k] = v2.vector().get_local() + + + for k in range(len(indexes)): + if options['Algebra']['mode'] == '+': + vout.vector()[:] = vdict1[k] + vdict2[k] + elif options['Algebra']['mode'] == '-': + vout.vector()[:] = vdict1[k] - vdict2[k] + elif options['Algebra']['mode'] == 'aliasing': + VENC = options['Algebra']['VENC'] + aliasing = False + for l in range(len(vout_vec)): + if k>0: + #mean1 = abs(np.mean(vdict2[0][:])) + #mean2 = abs(np.mean(vdict2[1][:])) + #mean3 = abs(np.mean(vdict2[2][:])) + #treshold = 10*max(mean1,mean2,mean3) + #if abs(np.mean(vdict2[k][:]))>treshold: + # vout_vec[l] = 0 + if vdict1[k][l]-vdict1[k-1][l] < -VENC: + vdict1[k][l] = vdict1[k][l] + 2*VENC + vout_vec[l] = vdict1[k][l] + else: + vout_vec[l] = vdict1[k][l] + else: + vout_vec[l] = vdict1[k][l] # first case is suppouse to be aliasing-free + + vout.vector()[:] = vout_vec + else: + raise Exception('Not supported operation between vectors!') + + output.write(vout, times_vec[k]) + if checkpoint: + path = outpath + 'checkpoint/' + str(indexes[k]) + '/' + 'u.h5' + inout.write_HDF5_data(MESH['mesh'].mpi_comm(), path , vout, '/u', t=times_vec[k]) + if mode =='curves': if options['Error-curves']['undersampling']>1: @@ -495,27 +555,54 @@ def ERRORmap(MESH, mode, outpath, reference_path, checkpoint_path, refname,comna u = Function(W) w = Function(W) + ones = interpolate(Constant(1), V.sub(1).collapse()) + L_sys = assemble(ones*dx) + VENC = options['Error-curves']['VENC'] + for typs in options['Error-curves']['type']: ucomp = [] wcomp = [] + div_array = [] + varu = [] times = [] - dt = options['Error-curves']['dt'] - for k in range(1,len(indexes)): - path_w = checkpoint_path + str(indexes[k]) + '/'+comname+'.h5' - path_u = reference_path + str(indexes0[k]) + '/'+refname+'.h5' + + if typs == 'l2_comp': + wy = Function(V1) + wz = Function(V1) + comname2 = comname + wx_array = [] + wy_array = [] + wz_array = [] + if typs == 'utrue-uobs': + u_path = options['Error-curves']['true_check'] + 'checkpoint/' + w_path = options['Error-curves']['ref_check'] + 'checkpoint/' + comname2 = 'u' + varux_array = [] + varuy_array = [] + varuz_array = [] + else: + u_path = reference_path + w_path = checkpoint_path + + dt = options['Error-curves']['dt'] + + for k in range(1,len(indexes)): + path_w = w_path + str(indexes[k]) + '/'+comname2+'.h5' + path_u = u_path + str(indexes0[k]) + '/'+refname+'.h5' u.rename('meas', 'meas') w.rename('w', 'w') hdf_w = HDF5File(MESH['mesh'].mpi_comm(),path_w,'r') - if 'w' in comname: + if 'w' in comname2: hdf_w.read(w, 'w/vector_0') else: hdf_w.read(w, 'u/vector_0') - hdf_u = HDF5File(MESH['mesh'].mpi_comm(), path_u, 'r') - hdf_u.read(u, 'u/vector_0') - hdf_u.close() + if typs != 'l2_comp': + hdf_u = HDF5File(MESH['mesh'].mpi_comm(), path_u, 'r') + hdf_u.read(u, 'u/vector_0') + hdf_u.close() + u_vec = u.vector().get_local() hdf_w.close() - u_vec = u.vector().get_local() w_vec = w.vector().get_local() print('computing error in timestep numer',k) @@ -525,59 +612,53 @@ def ERRORmap(MESH, mode, outpath, reference_path, checkpoint_path, refname,comna elif typs == 'max': ucomp.append(np.max(abs(u_vec))) wcomp.append(np.max(abs(w_vec))) + elif typs == 'norm2': + ucomp.append(np.sqrt(assemble(dot(u,u)*dx)/L_sys)) + wcomp.append(np.sqrt(assemble(dot(w,w)*dx)/L_sys)) + elif typs == 'grad': + wcomp.append(np.sqrt(assemble(inner(grad(w),grad(w))*dx)/L_sys)) + elif typs == 'l2_comp': + wx = w.sub(0,deepcopy=True) + wy = w.sub(1,deepcopy=True) + wz = w.sub(2,deepcopy=True) + wx_array.append(np.sqrt(assemble(wx*wx*dx)/L_sys)/VENC) + wy_array.append(np.sqrt(assemble(wy*wy*dx)/L_sys)/VENC) + wz_array.append(np.sqrt(assemble(wz*wz*dx)/L_sys)/VENC) + elif typs == 'div': + div_array.append(np.sqrt(assemble(div(u)**2/h*dx)/L_sys)/VENC) + elif typs == 'utrue-uobs': + wx = w.sub(0,deepcopy=True) + wy = w.sub(1,deepcopy=True) + wz = w.sub(2,deepcopy=True) + ux = u.sub(0,deepcopy=True) + uy = u.sub(1,deepcopy=True) + uz = u.sub(2,deepcopy=True) + varux_array.append(np.sqrt(assemble((ux-wx)**2*dx)/L_sys)/VENC) + varuy_array.append(np.sqrt(assemble((uy-wy)**2*dx)/L_sys)/VENC) + varuz_array.append(np.sqrt(assemble((uz-wz)**2*dx)/L_sys)/VENC) else: raise Exception('No defined type for curve printing!') times.append(k*dt) - - np.savetxt(outpath+'u' +typs+'.txt',ucomp) - np.savetxt(outpath+'w' +typs+'.txt',wcomp) + + if typs == 'grad': + np.savetxt(outpath+'w' +typs+'.txt',wcomp) + elif typs == 'l2_comp': + np.savetxt(outpath+'wx.txt',wx_array) + np.savetxt(outpath+'wy.txt',wy_array) + np.savetxt(outpath+'wz.txt',wz_array) + elif typs == 'div': + np.savetxt(outpath+'div.txt',div_array) + elif typs == 'utrue-uobs': + np.savetxt(outpath+'varux.txt',varux_array) + np.savetxt(outpath+'varuy.txt',varuy_array) + np.savetxt(outpath+'varuz.txt',varuz_array) + else: + np.savetxt(outpath+'u' +typs+'.txt',ucomp) + np.savetxt(outpath+'w' +typs+'.txt',wcomp) + np.savetxt(outpath+'times.txt',times) - if mode == 'histogram': - from pathlib import Path - import pickle - u = Function(V) - w = Function(V) - errors = {} - - for k in range(len(indexes)): - path_w = checkpoint_path + str(indexes[k]) + '/'+comname+'.h5' - path_u = reference_path + str(indexes0[k]) + '/'+refname+'.h5' - - if indexes0[k] < 10: - path_u = reference_path + '0' + \ - str(indexes0[k]) + '/'+refname+'.h5' - - u.rename('meas', 'meas') - w.rename('w', 'w') - hdf_w = HDF5File(MESH['mesh'].mpi_comm(),path_w,'r') - hdf_w.read(w, 'w/vector_0') - hdf_u = HDF5File(MESH['mesh'].mpi_comm(), path_u, 'r') - hdf_u.read(u, 'u/vector_0') - hdf_u.close() - hdf_w.close() - - u_vec = u.vector().get_local() - w_vec = w.vector().get_local() - errors[k] = np.zeros(u_vec.size) - - for l in range(len(errors[k])): - #errors[k][l] = np.nan - if u_vec[l]<1e-9: - errors[k][l] = -1 - else: - eta = np.abs(w_vec[l]/u_vec[l]) - if np.abs(eta)>50: - errors[k][l] = -1 - else: - errors[k][l] = eta - - - write_path = Path(outpath) - fpath = write_path.joinpath('errors.dat') - pickle.dump(errors, fpath.open('wb')) - if mode == 'h5': xdmf_u = XDMFFile(outpath+'meas.xdmf') #xdmf_ucor = XDMFFile(output_path+'ucor.xdmf') @@ -592,10 +673,6 @@ def ERRORmap(MESH, mode, outpath, reference_path, checkpoint_path, refname,comna path_w = checkpoint_path + str(indexes[k]) + '/'+comname+'.h5' path_u = reference_path + str(indexes0[k]) + '/'+refname+'.h5' - if indexes0[k] < 10: - path_u = reference_path + '0' + \ - str(indexes0[k]) + '/'+refname+'.h5' - u.rename('meas', 'meas') w.rename('w', 'w') #ucor.rename('ucor', 'ucor') @@ -612,34 +689,173 @@ def ERRORmap(MESH, mode, outpath, reference_path, checkpoint_path, refname,comna #ucor.vector()[:] = u_vec + w_vec xdmf_u.write(u, k) - #xdmf_ucor.write(ucor, k) xdmf_w.write(w, k) - if mode == 'colormap': - colormap = XDMFFile(outpath+'colormap.xdmf') + if mode == 'u': + colormap = XDMFFile(outpath+'u.xdmf') + #ds = Measure("ds", subdomain_data=MESH['boundaries']) + u = Function(W) + #P2 = FunctionSpace( MESH['mesh'], 'P',1) + cm = Function(DGs) + vs = TestFunction(DGs) + + for k in range(len(indexes)): + path_u = reference_path + str(indexes0[k]) + '/'+refname+'.h5' + u.rename('meas', 'meas') + cm.rename('vel','vel') + hdf_u = HDF5File(MESH['mesh'].mpi_comm(), path_u, 'r') + hdf_u.read(u, 'u/vector_0') + time = hdf_u.attributes('u/vector_0').to_dict()['timestamp'] + print('making the colormap in the time',time) + hdf_u.close() + r = assemble(dot(u,u)*vs/cellv*dx).get_local() + cm.vector()[:] = np.sqrt(np.abs(r)) + #cm.vector()[:] = np.sqrt(r) + colormap.write(cm, time) + + if mode == 'w': + colormap = XDMFFile(outpath+'w.xdmf') + #ds = Measure("ds", subdomain_data=MESH['boundaries']) + w = Function(W) + #P2 = FunctionSpace( MESH['mesh'], 'P',2) + cm = Function(DGs) + vs = TestFunction(DGs) + + for k in range(len(indexes)): + path_w = checkpoint_path + str(indexes0[k]) + '/'+comname+'.h5' + w.rename('cor', 'cor') + cm.rename('w','w') + hdf_w = HDF5File(MESH['mesh'].mpi_comm(), path_w, 'r') + hdf_w.read(w, 'w/vector_0') + time = hdf_w.attributes('w/vector_0').to_dict()['timestamp'] + print('making the colormap in the time',time) + hdf_w.close() + r = assemble(dot(w,w)*vs/cellv*dx).get_local() + cm.vector()[:] = np.sqrt(np.abs(r)) + colormap.write(cm, time) + + if mode == 'w/u': + colormap = XDMFFile(outpath+'w_u.xdmf') #ds = Measure("ds", subdomain_data=MESH['boundaries']) u = Function(W) w = Function(W) - cm = Function(W) - dt = options['Corrector']['dt'] + cm = Function(DGs) + vs = TestFunction(DGs) for k in range(len(indexes)): - print('making the colormap in the time',np.round(k*dt,2)) path_w = checkpoint_path + str(indexes[k]) + '/'+comname+'.h5' path_u = reference_path + str(indexes0[k]) + '/'+refname+'.h5' u.rename('meas', 'meas') w.rename('w', 'w') - cm.rename('color','color') + cm.rename('w_u','w_u') hdf_w = HDF5File(MESH['mesh'].mpi_comm(),path_w,'r') hdf_w.read(w, 'w/vector_0') hdf_u = HDF5File(MESH['mesh'].mpi_comm(), path_u, 'r') hdf_u.read(u, 'u/vector_0') + time = hdf_u.attributes('u/vector_0').to_dict()['timestamp'] + print('making the colormap in the time',time) hdf_u.close() hdf_w.close() - uvec = u.vector().get_local() - wvec = w.vector().get_local() - cm.vector()[:] = np.sqrt((uvec - wvec)**2) - colormap.write(cm, k*dt) + rw = assemble(dot(w,w)*vs/cellv*dx).get_local() + ru = assemble(dot(u,u)*vs/cellv*dx).get_local() + cm.vector()[:] = np.sqrt(rw)/np.sqrt(ru) + colormap.write(cm, time) + + if mode == 'dot': + colormap = XDMFFile(outpath+'dot.xdmf') + #ds = Measure("ds", subdomain_data=MESH['boundaries']) + u = Function(W) + w = Function(W) + cm = Function(V1) + u.rename('meas', 'meas') + w.rename('w', 'w') + + for k in range(len(indexes)): + path_w = checkpoint_path + str(indexes[k]) + '/'+comname+'.h5' + path_u = reference_path + str(indexes0[k]) + '/'+refname+'.h5' + hdf_w = HDF5File(MESH['mesh'].mpi_comm(),path_w,'r') + hdf_w.read(w, 'w/vector_0') + hdf_u = HDF5File(MESH['mesh'].mpi_comm(), path_u, 'r') + hdf_u.read(u, 'u/vector_0') + time = hdf_u.attributes('u/vector_0').to_dict()['timestamp'] + print('making the colormap in the time',time) + hdf_u.close() + hdf_w.close() + cm = project(inner(u,w),V1) + cm.rename('dot','dot') + colormap.write(cm, time) + + if mode == 'div(u)': + + colormap = XDMFFile(outpath+'divu.xdmf') + #ds = Measure("ds", subdomain_data=MESH['boundaries']) + u = Function(W) + cm = Function(DGs) + vs = TestFunction(DGs) + + for k in range(len(indexes)): + path_u = reference_path + str(indexes0[k]) + '/'+refname+'.h5' + u.rename('meas', 'meas') + cm.rename('divu','divu') + hdf_u = HDF5File(MESH['mesh'].mpi_comm(), path_u, 'r') + hdf_u.read(u, 'u/vector_0') + time = hdf_u.attributes('u/vector_0').to_dict()['timestamp'] + print('making the colormap in the time',time) + hdf_u.close() + r = assemble(div(u)**2*vs/cellv*dx).get_local() + cm.vector()[:] = np.sqrt(r) + colormap.write(cm, time) + + if mode == 'div(u)/u': + + colormap = XDMFFile(outpath+'divu_u.xdmf') + #ds = Measure("ds", subdomain_data=MESH['boundaries']) + u = Function(W) + cm = Function(DGs) + vs = TestFunction(DGs) + + surface = np.zeros(cm.vector().get_local().size) + volume = np.zeros(cm.vector().get_local().size) + Cells = MESH['mesh'].cells() + Points = MESH['mesh'].coordinates() + for k in range(MESH['mesh'].num_cells()): + A = Points[Cells[k]][0] + B = Points[Cells[k]][1] + C = Points[Cells[k]][2] + D = Points[Cells[k]][3] + # volume + Vol = 1/6*np.abs(np.inner(A-D, np.cross(B-D, C-D))) + volume[k] = Vol + # surface area + SA = 0 + for l in range(4): + if l == 0: + [P1, P2, P3, P4] = [A, B, C, D] + if l == 1: + [P1, P2, P3, P4] = [D, A, B, C] + if l == 2: + [P1, P2, P3, P4] = [C, D, A, B] + if l == 3: + [P1, P2, P3, P4] = [B, C, D, A] + + SA += 0.5*np.linalg.norm(np.cross(P2-P1, P3-P1)) + surface[k] = SA + + + + for k in range(len(indexes)): + path_u = reference_path + str(indexes0[k]) + '/'+refname+'.h5' + u.rename('meas', 'meas') + cm.rename('divu','divu') + hdf_u = HDF5File(MESH['mesh'].mpi_comm(), path_u, 'r') + hdf_u.read(u, 'u/vector_0') + time = hdf_u.attributes('u/vector_0').to_dict()['timestamp'] + print('making the colormap in the time',time) + hdf_u.close() + rdiv = assemble(div(u)**2*vs*dx).get_local() + ru = assemble(dot(u,u)*vs*dx).get_local() + cm.vector()[:] = np.sqrt(rdiv)/np.sqrt(ru)*volume/surface/np.sqrt(120) + colormap.write(cm, time) if mode == 'error_u': xdmf_u = XDMFFile(output_path+'error_u.xdmf') @@ -874,6 +1090,11 @@ def ESTIMpressure(MESH, outpath, checkpoint_path, filename, options): #from dolfin import HDF5File V = MESH['FEM'] W = MESH['FEM'].sub(1).collapse() + p = Function(W) + one_mesh = interpolate(Constant(1), W) + dt = options['Estim_Pressure']['dt'] + p_drop_lst = [] + time_ = [] # Checkpoints folders unsort_indexes = os.listdir(checkpoint_path) indexes = [int(x) for x in unsort_indexes] @@ -897,22 +1118,25 @@ def ESTIMpressure(MESH, outpath, checkpoint_path, filename, options): sphere1 = mshr.Sphere(x1, radius1) mesh_s1 = mshr.generate_mesh(sphere0, Npts0) mesh_s2 = mshr.generate_mesh(sphere1, Npts1) + # New elements + VS1 = FunctionSpace(mesh_s1, FiniteElement('P', mesh_s1.ufl_cell(), 1)) + VS2 = FunctionSpace(mesh_s2, FiniteElement('P', mesh_s2.ufl_cell(), 1)) + s1 = Function(VS1) + s2 = Function(VS2) + LagrangeInterpolator.interpolate(s1, one_mesh) + LagrangeInterpolator.interpolate(s2, one_mesh) + vol1 = assemble(s1*dx) + vol2 = assemble(s1*dx) - VS1 = FunctionSpace(mesh_s1, FiniteElement('P', mesh_s1.ufl_cell(), 1)) - VS2 = FunctionSpace(mesh_s2, FiniteElement('P', mesh_s2.ufl_cell(), 1)) - s1 = Function(VS1) - s2 = Function(VS2) - p = Function(W) - - one_mesh = interpolate(Constant(1), W) - LagrangeInterpolator.interpolate(s1, one_mesh) - LagrangeInterpolator.interpolate(s2, one_mesh) - vol1 = assemble(s1*dx) - vol2 = assemble(s1*dx) - - dt = options['Estim_Pressure']['dt'] - p_drop_lst = [] - time_ = [] + if options['Estim_Pressure']['method'] == 'boundaries': + boundaries = MESH['boundaries'] + ds = Measure("ds", subdomain_data=boundaries) + bnd1 = options['Estim_Pressure']['boundaries'][0] + bnd2 = options['Estim_Pressure']['boundaries'][1] + area1 = assemble(one_mesh*ds(bnd1)) + area2 = assemble(one_mesh*ds(bnd2)) + + for k in range(0, len(indexes)): path = checkpoint_path + str(indexes[k]) + '/'+filename+'.h5' @@ -920,11 +1144,14 @@ def ESTIMpressure(MESH, outpath, checkpoint_path, filename, options): hdf.read(p, 'p/vector_0') hdf.close() - LagrangeInterpolator.interpolate(s1, p) - LagrangeInterpolator.interpolate(s2, p) - - P1 = assemble(s1*dx)/vol1 - P2 = assemble(s2*dx)/vol2 + if options['Estim_Pressure']['method'] == 'spheres': + LagrangeInterpolator.interpolate(s1, p) + LagrangeInterpolator.interpolate(s2, p) + P1 = assemble(s1*dx)/vol1 + P2 = assemble(s2*dx)/vol2 + if options['Estim_Pressure']['method'] == 'boundaries': + P1 = assemble(p*ds(bnd1))/area1 + P2 = assemble(p*ds(bnd2))/area2 p_drop_lst.append(P2-P1) time_.append(k*dt) @@ -1012,41 +1239,60 @@ def OUTLETwind(MESH, output_path, checkpoint_path, filename, bnds): def ROUTINE(options): + if 'Algebra' in options: + if options['Algebra']['apply']: + + if rank == 0: + print('Applying Algebra') + print('Choosen mode: ' + options['Algebra']['mode']) + + MESH = LOADmesh(options['Algebra']['meshpath']) + v1path = options['Algebra']['v1']['path'] + 'checkpoint/' + v1name = options['Algebra']['v1']['name'] + v2path = options['Algebra']['v2']['path'] + 'checkpoint/' + v2name = options['Algebra']['v2']['name'] + outpath = options['Algebra']['outpath'] + + ERRORmap(MESH, 'algebra', outpath, v1path, + v2path, v1name, v2name, options) + if 'Outlet_Wind' in options: if options['Outlet_Wind']['apply']: if rank == 0: print('--- Outlet Windkessel ---') - mesh_path = options['Outlet_Wind']['mesh_path'] + meshpath = options['Outlet_Wind']['meshpath'] out_path = options['Outlet_Wind']['out_path'] filename = options['Outlet_Wind']['filename'] checkpoint = options['Outlet_Wind']['checkpoint'] + 'checkpoint/' bnds = options['Outlet_Wind']['bnds'] - MESH = LOADmesh(mesh_path) + MESH = LOADmesh(meshpath) OUTLETwind(MESH, out_path, checkpoint, filename, bnds) - if 'Corrector' in options: - if options['Corrector']['apply']: - if rank == 0: - print('Applying Corrector') + if 'Colormap' in options: + if options['Colormap']['apply']: - MESH = LOADmesh(options['Corrector']['mesh_path']) - u_path = options['Corrector']['u_path'] + 'checkpoint/' - w_path = options['Corrector']['w_path'] + 'checkpoint/' - wname = options['Corrector']['wname'] - uname = options['Corrector']['uname'] - mode = options['Corrector']['mode'] - outpath = options['Corrector']['outpath'] + for mode in options['Colormap']['mode']: + if rank == 0: + print('Applying Colormap') + print('Choosen mode: ' + mode) - ERRORmap(MESH, mode, outpath, u_path, - w_path, uname, wname, options) + MESH = LOADmesh(options['Colormap']['meshpath']) + upath = options['Colormap']['upath'] + 'checkpoint/' + wpath = options['Colormap']['wpath'] + 'checkpoint/' + wname = options['Colormap']['wname'] + uname = options['Colormap']['uname'] + outpath = options['Colormap']['outpath'] + + ERRORmap(MESH, mode, outpath, upath, + wpath, uname, wname, options) if 'Velocity' in options: if options['Velocity']['apply']: if rank == 0: print('--- Reading Velocity ---') - MESH = LOADmesh(options['Velocity']['mesh_path']) + MESH = LOADmesh(options['Velocity']['meshpath']) filename = options['Velocity']['filename'] checkpoint_path = options['Velocity']['checkpoint'] + 'checkpoint/' outpath = options['Velocity']['checkpoint'] @@ -1067,7 +1313,7 @@ def ROUTINE(options): if options['W_map']['apply']: if rank == 0: print('Applying W--MAP') - MESH = LOADmesh(options['mesh_path']) + MESH = LOADmesh(options['meshpath']) filename = options['W_map']['filename'] outname = options['W_map']['out_name'] mode = 'w' @@ -1078,7 +1324,7 @@ def ROUTINE(options): if options['GradW_map']['apply']: if rank == 0: print('Applying Grad W--MAP') - MESH = LOADmesh(options['mesh_path']) + MESH = LOADmesh(options['meshpath']) filename = options['GradW_map']['filename'] outname = options['GradW_map']['out_name'] mode = 'gradw' @@ -1088,21 +1334,21 @@ def ROUTINE(options): if 'Pressure' in options: if options['Pressure']['apply']: if rank == 0: - print('Applying Pressure-MAP') + print('---Applying Pressure---') - MESH = LOADmesh(options['mesh_path']) + MESH = LOADmesh(options['meshpath']) + outpath = options['Pressure']['checkpoint'] + checkpoint_path = options['Pressure']['checkpoint'] + 'checkpoint/' filename = options['Pressure']['filename'] - outname = options['Pressure']['out_name'] - mode = 'p' - WORKcheck(MESH, mode, output_path, checkpoint_path, - filename, outname, options) + WORKcheck(MESH, 'p', outpath, checkpoint_path, + filename, filename, options) if 'Error_P' in options: if options['Error_P']['apply']: if rank == 0: print('Applying L2 error to Pressure') - MESH = LOADmesh(options['mesh_path']) + MESH = LOADmesh(options['meshpath']) refname = options['Error_P']['refname'] reference_path = options['Error_P']['refpath'] filename = options['Error_P']['filename'] @@ -1115,10 +1361,13 @@ def ROUTINE(options): if options['Estim_Pressure']['apply']: if rank == 0: print('Applying Pressure Estimator') + print('Method choosed: '+options['Estim_Pressure']['method']) - MESH = LOADmesh(options['mesh_path']) + MESH = LOADmesh(options['Estim_Pressure']['meshpath']) filename = options['Estim_Pressure']['filename'] - outpath = options['Estim_Pressure']['outpath'] + outpath = options['Estim_Pressure']['checkpath'] + checkpoint_path = options['Estim_Pressure']['checkpath'] + 'checkpoint/' + ESTIMpressure(MESH, outpath, checkpoint_path, filename, options) if 'Error-curves' in options: @@ -1162,7 +1411,11 @@ def ROUTINE(options): MESH = LOADmesh(options['Perturbation']['meshpath']) checkpath = options['Perturbation']['checkpath'] + 'checkpoint/' - name_folder = 'SNR' +options['Perturbation']['type']['SNR'] + 'V' + str(options['Perturbation']['type']['phase_contrast']) + if type(options['Perturbation']['type']['SNR'])=='str': + name_folder = 'SNR' +options['Perturbation']['type']['SNR'] + 'V' + str(options['Perturbation']['type']['phase_contrast']) + else: + name_folder = 'SNR' +str(options['Perturbation']['type']['SNR']) + 'V' + str(options['Perturbation']['type']['phase_contrast']) + out_check = options['Perturbation']['checkpath'] + name_folder + '/' READcheckpoint(MESH,'perturbation', out_check,checkpath,'u','u',options) diff --git a/codes/__pycache__/CS.cpython-36.pyc b/codes/__pycache__/CS.cpython-36.pyc index c0eede00f50d427370bdc2ec230ca0d976fe68db..a3acd0f37170f8e2e15ed7bb03472c6826ef32cc 100644 GIT binary patch delta 2485 zcmb7FTWl0n7(VCB>~?m$4BPFrWn0=Fw%D>QEf6rKP~;Ynh!Av9OqC`}&vd(`oh`FV zD0QG3KnMf{j>bkpf>&-qV+IY0qQpdfF+TVvPZ)jCM2(3LMt$)AXRy1yaB)t~H)qcI z=fBOD`92>xKhU&v)~xWr_#ZttHWTs#DY+C_uZEfb=v33`I zM7?M@M7ramQE<4|3YsFxrVoDYX5}nKV^pU>7GZ8mna|o5qoV`>VH(wF1V#hpFso@b zn*+uMjB4){_szf#g>E<<+Fr&IU>9bxCv=RCxj%(2GkVZ{HF8R69!G!8oS9jQqIiHTQCL52M{O96!@? z1HE2CubN7452iy_I*l@oMF3SAVe;#EQyEkNgUQdDtD>wH;>K8*G8$v^sdjo^ZOA_d zl-HYU0xZGoBkt3+1*P-L(nRu*h{@A~nVg+1^z*$4sr6e@>pMC+c-G>*Ijdmr=`Cb) z);#z7`T{H83)$H(B82H9UJbfFi>K*S;b z`7|rlAT?ya;V5GEE_RaHx}!P)gE&Ehh^V9E>Al{6RA^2m8aicDgLIMq=T26%SB56a zk?GdTx{0U542uJgim32Xy&;#hwmdN9JPu_$w2sV~l#NWbE}HCjYSX~40fp+*YruAb z&x#r%Vvb>8J7(X>%Cm`iy$H?F3qOK}CtId|v1uVZv*}hPA;QHFa!p$mt(4o_Tv{c| z7dDob=UYMZ)O#~yn(btM1nlwq`Jf`wkCZG432tU)NipF4aiwVQ7lVvhS;B$RbDs&o z%So*)TKHU&>GrlIcl9utOa3R3*W26g)a3h!I|J8>dEJ|FzU|n54MxHNIcFeV(dh)i z@6e1HgWlD)W~>;rgo-7URlCSZW%rNUi4$_d83v+LJ(}UE?9iUUH2x@gu8*gAVc&4J zH$BLQ(izJvWJkTD4!(wuv8L)EMylJ13jf>B&e;t3q1F%g3zTluVMnPzz-rWFcR8bmHhFAkl z&;=kfG1|o}7G~SnHf1O4WAg9C3B7@?27}c~fsHczoZGUbi-jpMf5)@% z?GTv(X-4Q=)NfEs`>OlwL%-0fH!#{tg!d3sgrl;jb0d9Ap6^UX-@%P{5zZsLEpK%` zxaI;bk0AsQE+Tjp^uj7>++>?aVD5uCYG1*_Vt}Lg_`ZBK3?$sml=p4dTglqEW70t9&R&*Khs8+3YYN}SH)obP27*Vyy!@=;A+VXz^ Dr5*)B delta 3908 zcmbVPTWnlM8J?MQcK7UhZ(g6h`LfyAn>bEN5<5ZA;*z)yq7d0NG$@kQxOIKjj<0df zZi;Q@upl?RKuw$s5lGON6q+{S(y|ng3Q~mxNJYHB6L~;Xl|YCWQUnhb4}AYw@5Nmw zEu5qI=bt(M{Fj;UpZ)jx+%H#pj)X#iYd`Cncz24iw^_rZfj^E``1R$UEBB_X^pk8n zY3;E(tj;I-c*@#qby*p3X=|62#kY2Ax7Cfl!|Jhe=sPX`Fzf4m^-trXMT_wm4|A^w z8s%XvCk}js6nBF#VQ_;7@n(1yEx{973c@Dd1TE3aS3GAlKBoR*oc4LTAsEoqsQ)W` zqx?tzO~EghF9t7b*%wJ+)>Rl5ALkE18B_^fM1{ris=?@Q_#0)t^<6!9EjxT}aoNrn z3VEv>OZ-dlTgpg(tM|90<0aDJ*ZJ0~h2cb)M}?sUc~sMRP*ZQFdrAy$GeczQyXfR8 z5kw#6ojjujwRTvIX#sAAatO}L4G}=o5fb)6#CVEFF^}~XEucH=>Gs&qz?&&L3VL>< z`O$`1d8vIu+ZKq5CB9KsnZNM96-d}m6Ims)M&t`bBoT+m7m1uy*R#F+DfMQy+k1v1 z^ab@!wvU@Cv%A}0^Gyx5r*`er?5|=;&5;ZS2HG%MVN9q(d#+UElFc#6*M$=-iaAXh z(k}xuT2Z&W8^Q@~iWP6sBZU-}C&|{u3AV}4yyJu+>z(sSX8A1NMbSuQcv?@hbuNP@ zb0RV_XE;$Av5a}mX)XF?tBk-_Oh#l3-%Aj1+IBE(*5_n&MqlSnylEzZnE*L#G8r&o z(n(^*JJ%w;vK4a?nv+SGN!8cRXpl==ewnrcgREq=fC&yVD}+9bK7u}qz7>5NdJ}yD zeG+{dy=P0g$2H}X>2+?hB$K`*>tIfLvI`D)I&ct9Tz1KLiaCi*e*T%$ zh*q}2jn*x|l_HxN|$~$|fH=zlNozjCHuejtM^0c$yV{S{r2NPeC$-%b|Vk zuMMe42tJ<0Q!ErG?fZxWe>R^j&MxM&pFf*la!1n(7gn;f6@GSUAz#eS7JCcX#mSY~ z#S4q-Yne=_kK3a_=BN!+Gj%(2qqEf>0&5>9;vq6daATLi)HM28r%;u?r|MZ zhqD91zMx7-k9xk8o+7bOjfBD7GKU`?`NXLQA3E{qa(;5Yuv#ddTgcy_l=>__vq0@O z(ZjeWOp)UI{@+`{hifRK+0Uu`?w_~)olH<3GF(>%ntio=v|H!AP2JnmRbm!HLb+5% zr!`6?D>%=9;x;Q}g^Fx0DDibiDc!}ADK#rRFE)i8mb&x+2VRqa9j34Ex4vWK$PxPeRq5{$!c@Z)d?$nNA)RaZ088l`szLS3Q;I#c5H^OL&pAm%->`&alS z^|$?*&=pMG@+_U3E#z&b%mdFKyG?iJV;z9 zQN9Fvik}jXbELt)4|E46#DrE9Yr1`{y!W0#;U_Eh^C0R=2TfHv_^Xy{G($Cd#BQo@9~vtS3OMj-jF||EFv?)mgs3MCj#`ligh@`-$C^oQiY!%2%qRtGW z?hL|N`zO#gMZa(Km}p)!w@C9M+*#&$(7b5AVu?w*6?vLjB3~@}!jhG@3zLh>3$sfzZn3fa#1n+0ER0uyoCADwKt=6Ey@6M^e+ro?oh0#; zKX>>o9)AH`C6HH$qmTBp>g2$K{JZMw1AWo!H1QIVmx;WnemC&pp&P`~|A9&=b(1(3 z<5S@J%v)k)1#J~=tzoBJ1J9D|t7`bj41Z1CJktBvCG!3O%pwc$--^;Oiljp}AOfT# z-U_~5gcHiCa;nEx#XfS!ex1zo>QTBU+NC~v^vgNdIoILm$o^F#?h#YIxZ7gCuKsXz s@4!Kf>E9?dGGjyyzmYW3#{GtF#Eo{N#VD~hW5gH=cmpSm`;7E|0XeA$n*aa+ diff --git a/codes/__pycache__/Graphics.cpython-36.pyc b/codes/__pycache__/Graphics.cpython-36.pyc index d0cf857bfc5254147e0ff22b58e628c6daf3b463..1cad2a76d272eeb965ae18a099ff72fb719d1251 100644 GIT binary patch delta 11996 zcmb7K33waFb;d3Z0w4*V;w6fLNJt_i@z!ZeqGU?44$GD$S+>@aArP|!LE_Q?6iImD zbxgZ)4%&1!&0(iGoWqHon{i?%j*~BEo1|@;G}k6glQ!-bC25-`jobQ5)BbOE0UmOa z7Wnq<%$qlFX5PH_X7ekXMeRag68`rdEudgSjMj(>D4`?-)lfG79EN4h>< zxr;TjHTO$O8{7TJPfqFTMz(cQUXZxSwy8A}n^@1BG_iTkG3fzrFY8shRp&&Q^|AhW zNpi4(YiQ20LC|bbYbUnONo;7+g?jG+)wYQq!m*)N$D9P*LMf;8s(oy`+S?4gbjFcx zRQuVEiGfn_oHQrhC#Ai>?xd93IMF!gPzKe8i6OOjPExk3A$2>dI_#1iYL_~w?wFI& zJ@5;QL+(m@) z>>wv@78o!B)qA9)i$|+>`oK?oS5+r!c;_Ou3%&;Z#Mp;c~^~D8OvliOQfTDHOE_weJwXijm85l zJEdXcxt4AS>xGu!T8y!uZ=ve)%m6ZY3InE-u^H0{w%#ReHT2d)p#4JYNw+D-hD@WW z?Na4(p}eAPqa+zmwC(W8wUYHWUTAw(dVKM1YbT`*BoZ+KegY{B=_W$Z@VG=^$k@7W zz$2#e4r6rPK!-J<7Hzq5*e6EKW5!3;?X4N6(i%lR)Q@MtG=95oIBYtj7y!Va0x-0E z-au`7PbD+?vuZBiZm^DHo~NnlCyd38SLx+8uIXT*Ds-}YCepuccH3Au#GfEO6;~2G z!4f(CF6x8K-%X6(LqLp%#sSSNV@=nU6jTf_Is@u*TVanq9*2&q@1y$fC-4CRVt6i7 zJ~?3ebFpYrjifS6eq>GG+32KtYCMzA97ZQ{wZJ*#fb5dJQoY>n_Dg;@ z{}gaXVs0wJk7Vq_GvaV8UKjkj)Lz!JYlOh7jZK}WnCSVfsB*OufZsu>=Lx))Ko{Xo z5g;M)cM>20#UFh${r+uuxlH&csQ4KG)0NJo)v^=~v-7ZC?#b2v!7tomX}dF#N%I)h zdoO_|0%Xga6PN-p>r&Bt7L%MzjP?q1%JJt2?;GQpR3w?mH*Z`#uHB=PvD$gOQc+O@r$=ngs`1{^=eLvL{8dcUHp+NjAR6vHz zze?b10B|9sIi_@nDV;YZg`cI8LjqN>to!SPDXS9Pu$qNmPQWiG6z6PW(CY5r5Y;yT zlq2pKO$%u9pvfz8A#XZo^Mm{yMDlYYskWymA}Bejq3J3Jr&NHFqH4}GKHIl;{VLrI ztkBKD1>>c@^Rg^0-qwGwR4bH`RMeV{2IGZ+9w}lp40d|Y01M5F8v6!Y&o&duWzl35 zjg@>Jkivq+gtSnx6Lw&UTS#<<>dg`;5D*Tr*0Q^K6-C zv}z=i%0{_*DuWh&YcvlXlmpVDGW1@lhn!#>4}L!y7MlrMU%xlZozODw0w~BmlNe6K z=-bhum&S$~SQ~;@)MO@xB|I47#TY zrPJJtYP#k9KSKDA5?Cz-!ug25pFkymDUXcsD#}$8csU8xj=%<0ZpkT;%cM8TSVXP^ zFNf15RY_G?MyjL$9^=H$t#X~TII}Y*Nym*J>{`23C<3|rs~yxhsl*Z6yqDZOzYS>9 z=-hpL_=`v-{sGpSc#urCamkaL-jtdfk7IG=LRdnIUm~VUrrc&cwflrFOzW$}(#3dm z+VaqF&Czg;DV^khMRop~z%u=y?5_y?8bERJ%ye3DVU+yJHyyYs1cN!nwbE7U~ z)`01NUNzFowahNW@+S1%jXfz(n9gh_XG+7UZjn5e`l~mtm0oSQkM8pKfxwiec#UKnIl5IR zV^ONIjAycu(PT6>IVko3d?%6BuqdC5i~<|ynW;3hwymTL@ZC4zCH8x!Gcz`3x~Inz zdG#f!r{E58bu1KiTI%7Wt!pkno5Wf~Jf#feLW;wCFc7clNudpu8^1o44l&w|@A1Os-EIBokBej@76^HB*+Svk&HjbT$NQV~Rcj8yo z_>y_c`8H{<@rm=@QiJiM^FtM;^UT@9g=V8R-E0gf`!Nq)f|HhWAk+>5@`07LDaIJZG7= z%(7QMM*Q8xJERzOml`C+*m$X1hO*yssa@(el9$#ktz?f}+S2v=RnBy)oa5tdhw;Np z*Y7nwu}qZb)U)7z>eS&-xWaVJf&iPHrEfgbGnK{GGiQ2|8G74YHExR@k}esKM#I}i z(Z8#$0AkmpR4y|-9>%!vp_5GVEj04O#=xEJ#=eWq#@-9{#`e*j1D3G)tu%h&16N9$ zpGLvwoKOp>hQotZ!FX(RxcyFK;JV8uXLYKdm_&>hM>n^NHBRPH(3wuGe2O<0y`o04 zS(OhMdt#j)OHYhwtf`y)JSPBy+(GH%09StppMUkg$?=!W6vS)vaG&YPCSsE?eKJF{ zBAJP@{7l|5zXMeLmxOm^vMRo%sMPY9d(b&-%5zst4_Bl3MB-gUOgjqkNi5b`)18B` z^L&U#Lbhyrb80?9dw(mX&(deWlxQKPA`HMO>5Ree%Gm^`mC#z)SM5rKONgOV@sS&i z2n!XP9eh17(<{_sG|H_F;Xz{a7X+pW?4sW6CUE{V<=ZH^mVj{n`}lQ~zMjBd0>ayj zA*J}ia#kHorH;n&^_%B6Q1K1mVaoXwAEDfh#(%K`jb&%S4-oPofN7jmn~#|l5ZPEF zslusZ?76(@ITuY%Axdz9Xk2-i1%HHkMcSh{M-HClbCi96z}pD?jKD7A18RS(aAm$M zS97@tmScP?Q4+q%cu^fd{GfZRw{@ignqtX{uv1@>7c*mTkR;jgjPH<&#%J$&!9xoT zf&x$JJ@rBv+}4UBJoQ(MpXoO=uVUsNx#^aX6E8ZAk$Xe!_hA4rJmN=^ti>~eCwHi3 zaqiw%N&0@x!8%(VOm1~7c=-OjgLP2=iiL_kR$uh9jduB_Rx57y9BU~07hE&<*!az@ zjzf+G*TrP9QmbU$BHs2B7A{t@#-b0HEdui*+giqK14iav)}B9?Q}2Fb+h~Po{0vnj`cM=w8~}&>!-k5F<@f{R$yyv?4Y2j zzul4UbUE^RCKqvAiIT~ns@3DWUr{u{s$`D zx>P}Qt9}l~x*%Qt+U*W@S|i6)!N~$n+l~x5*oD)!3&Nfue_^XEKEvT9T0~B-W4DQH zD>~6K=eQsj+lp(owG*vcD|4`qKLEd|5ow#&3deMg-A5#woxpUB&JxM$%f zMpyX>a1+wGV$%s9y39Z5R)Zx6L5-Ij1f}6-m;3{}tGJFuv^v4?8Fs0-&N3_(g;rGX zudvbLy1YxP)7D{>b#rpDLvv5KwGPOw5`5ZWpzBwc>o!3Q&tn4o0@QOk7ILOB5nNc?r=a;I$OT&y#tY>Ld+nf(-n^}eqo3xO&iP-IR2waW8Z7ksye3$*` z%0{gLQ{>|c%ND!Sz3j4`pJ1OrzJulLe1YYQdUra&rfkH6?3rS>5bZ{-8-6`h^3v?e z3fEI{`RE^UJ)}G2dVKs}z>jRj#w;`PS~ktr>1EGFPJ)QwB-qSqGZ9uYE3Bl39{}Yn zsA$WA9IT{ZyQ{JRFF0q8*i?55%Wz)YUTm;2%avNEW{DXYhSN@K$YWnq3~OQbEJW$& zKY+<>Syrl}bSq}FLEFmi75-zJoqL7I_1L-lM6MT;5GrmFxVjbebvj#Kp^n)VY*Kr# zqBXENEV`Baw6=wyU3$M@-e>0?5V?Lk_a`DZz~;o%R`M*{RUBlmEdAr|;!wJqZ7=GB zDA>*(WUpecX0M66ibL$RY@R(7cNTZBhuH#q1gV{nLlHJ!#c#3=uTOI_D4(`n+d=OY zyPZKJl?Ob4~0Ynb-fOsQ-&(_Wh?l?kTAGyPd~`X15g zXUd)4*5W9)u~W36FW<-YYdwO1dOxV`yr!W6tyh%O2&jCiv9op)|L-B_cEy(+1Im`I z?zZ~94q9K0eV?`SHN-sc(sr>mjgI1OZ8yuqY3;VT?q;t?M~B#-mZpInwYUhSzJRn& zx=)ts=NC|kEY-%8EOo(u`NHkwBMJWoPzLy0*c~Z!^ z_Ad4=O-Jw8yV-k)INez$`jjAg!X{cq0~s~pmce`3`{w21er^9mc+N4gg*_>@8f$SU z;3@WgY{5P-kBE7>aq@uy_8j|Ac`s1QK5WTX0oj&?np~|Q`80d$GLny$rT}FnH$}e7 znu`w4ptSc+~w z{on!pD5XVQmV%aN<99cH@Wvk(?jN>doD_^QeHk82$5oEgdSvnY=`ap)FCf7$5@;iE z7lCsGb`iLP8mm*J_wNJ?>w1fr`J1LK#kUaL2zGb71eI zP@&d3cMT%Gf`*OR4yW7x>-7H;v0!`^*G8RV_dyVqb&qPmacr*qo*VyqweVh%a5DxXfhE`UmZlv z_+Nevze`{H?Qp?Ml|)N8l|yg^_dX&NiO&_9NJCcQNHr8+ER(*Xas)jK%~tGiOsG zKKWcPB1v3?cy$WPKtg{o>iBI}IqFzgd-|5!&K|q@aPMBk5~Io7ek3qV5d`LE&>`5w zS-h+pKDSYdqWbY?Iv)7>tl;;0GCD*C2pgc2YzTS^otq4>|-d^(XHJ2x=cU-0)`iSj9n(n?#zR$9>W z8d5D+XnH0@a28MxO&3C#GM2a!!r`Pn)Tsz4mgPXXEQbot5D1(h;GLmD2q(^oD}@RG z>Xqq9q*)@K;S}Xpyim<3?n}iLKMqrK>Qo|`R0~~EhAT$#U{p^=6I2RO;VwWhkq$0V zDQ*mJpg%q&exDzXzu0}@;m`f%r^AZJ8h8B11@gFG8Mcl7M>u4O!=;rEWWz2xID+@Z zbq-wIvO?+PV@)4U>eIQgjLu7w|Knkt8k9O45+|ou>p%(~4dWS}C}h(4Xi^wj3Nvcf zk)b4@f@KtcjI=F^I`Mt}@4Pd+?`OlDBJ!|P?^rY)o5ig!s~AY}GhO;`zus3O<*aUZXK=>E)Hg9#j#t&q9>yD)(rr1DJLC4$oAEvrnN!JQcByif}DYpbkazQ9F5XM$MU3 z2%9G(qiTLyRny{T1NWk9rZ2AJOkoTwfNECL@dBOQSuVX2xl-Hsqwht#fpG^M@kj@CdHE# zYY(5Iu3jN9O+cR^NCyafmcY{afhP~Q8=bYzF-X&Jj*K(CNIoN6q$y9&pc}Jz3W{5} zKI=7Wy$XseiCqML71eZ8au6p(VmmRY*oh(jhKOJc&qb>fU7BNNWvV=_IY@-2Sg9bv zM50N>olVip#Q|}grN&nup;p-IlaJP{plqaGJEBW2{Oj>F;aH#2%OPP2M^^!4q?BC!N8BhTfd}xq>$W6 z2!Za90xr6jfg1;+m~J2V=?(&}xOB<^XHW{@=7H5n0PiX}D4GbMDY}zDCB%NIJ5|zB zXRN^mP#(Y_P+tz{qHiu~lMME9t+NWv_`!_+#Iy~U9;}{QZl1R1p78S>Vtms!n5S5vg`arye=KkPF?$ z@QSW(LJD#eG_~k5)$)jjyCwd-@yhv-w9|P1{8k)|e|P@w-X~y=*n`r`y6yx`!!4Jc zJczgAVVUF@)_iw*DwUnZ-P>0@v<6oQ-uTcuX^rvpLp?eY_2O{q^2_$C2FxE(w|ppHN%flUMs5fC2uK1w|gV0t4Fve8Hc2XJ&w$`fE= zx~7cF4>yG@CmiMRE0%k}%{y$ha*8jJ&Q9g&C=hp_pL%%vARRd3G365`Np zfLMC>;`+l>Qr!%8kk_WvF;(dp4VWIPX E0qg(!(EtDd delta 16743 zcmdUWdvF}ddEd_L1B=CC0RkYNB)9~Zk0l8#9wd)EJ^&JY9*-nQfR7tKmP_4vOgEnBe=5f|Khec>W;D&w(+Bq(#ekf;`hdMwukNgs6*M#dRZUtJ6MGE(VK%fPP48ZD%y__?VEdFks&jfT+s_WXBuNhT)Eb#2I|!LxwQ;&{ zMPi3$TzKwlkcv$A6OSG0aI8oVb%}qWkW~iMLH4vd(1rrZk!n$g*x~7Yh0+yiMf!-8 z@`8JWQfkX|%ZfwUuQpE~PzP2d8BmmE}gsQcA}D>9lVjXT1_kH{t$!;&fO z2JvK;GRd^jDfdds>RH+E+0=ky-biEvNZ7^wlx`y8Tg`YkNI^wTW)yEKHmgP#V%+q^ zIFBu5O=-%M7G87kKI8MXuSxx@=Y5;xR(jZU$C8=*G1GNV&BZQ;Yk8}&?B6Uk8!!7$ zbQ~bX`${7+UhS_xP%HV3cp>PcvkVs{DhY@UR_l%SDe&bjhJ9-|;u*}a8 z=PXG4s`&e(;qX>|7CB*pk#@N~H~oe6kFXI`dW+A~7ERp4(qBh9<)1>G#L_}|>9~hoRuIXBsO{9p%79m}%E9|lR z6)FFP4as?pgrQ~9q^*{>5RGe5Fs?L=%w9lM{-rOz6F zzokQZebpEIYiUz)?5yb+_ZxwpF{|enq2W;-X}wCRuMznBY%GUWB9~0u>zhm`(>%+6oA|X;JUtstCbBte zY6#nMDse3)6SMqZlZMX`AzO#}E1;$m(!%U~fwJyQY(A@+&dKCdxS6kz@N34u+19fF4VnpE!?$3o|& zKUn>jp|8r_YvzWS5HzI15;Wf3xkdV&)er6bflJzN%Q_OVEW|Dk~9wUPZgDdr1jl*Ork510UGT{cH%5>vt zo-@6(36@Et>82~2OK0$gjV--<4oy%EQrwj8aRW)ECZ#^0r>nc&kiDOn%BgIW-qM9w z=5=GGH+uNnRQEd|_IP@8I4}1-5F_l-7}#OAXoKQW7jtSVXU#Edlo^}*+Vnr9=ZRn# zimCS>iSb87z7CT35)G#<)-3x9aqkkT8ct*+F`TS{^#)~yQGS3=r zJeIVER1f$%Y>~*TJcYF4nMkLT+0%+QuI4g{MK!53o=zlT&FhKW{V}p!mLDOFRTG`Y z9)FVvO?mznk>3NsVnc(@lx~>PbEc&5%T#hgu=3VApubJ5;xIIwOwC{kPhbg8D9)wC zpw*9mPg4H@L~&2jyHVuD9F}}s`5eS@gZx)W;!jAT-hNL}A<1b3E3Q1oG||{1v}f^X z79R+FR|lpQ=VBIX?GlLV_*5cU{n`T+qkYgYz8iT?{=Q@NX#a<#Mlrl;&|B}W$@u0# zuXM}!Ujse5c=-~)gF-C8!38y$o=oJHf}ux*0x&eiLXbq-bj8rF!a@pNVx180w2g|L zFue7~gz>ZB@f?wPB4U+ov}X55!BE_(G@p%4nx09FHq?K68{>JD!cp{@%|h{aq;v6Of7u7QkDg2Mk_Msd;|$o zx>h}g8m`$>Vf23Na4MxHFUs&7_CnYp`JExzEBhtC?3Dt}MmZq2x;kaQ(~&>3CW?8kT%qv{{&<|>;X7n@lpm=&k*Y_^FH(UOU0SXe$p$b3BGrgggGg;a zs!^nxklG+p%}6zgR0~qgBDE2z7LjU2YNJTCA=N5U?MSsDC9@9Jsd=ZJE6(M1%4_XR z(mNcitHpr=WUO?XLz9J2N6t0vUU4jUY7VwZIC^h@Yzqs%E>rZ7lSXr^pSk?Si|bg4@7$3hpkZmkT$tuwd`55N?A1Ns9Ltxi#d8vdBwt+t9Rr zG)?BACp66^nl``&iyAyyyXGdY4#h(iY|ok&*RXwM)$JG6O+TT!T6Tav#SXsYDq6t@ z&Hy{ao@R&1HfZz_q2)_C{qY8^V@KJs3N?OkNX;=z4VvI(C)sdZCYKhuS~gN>C~`h_ zyewD8PQYT8oPWguht~mz_mqCoTIs*|-O=H{eCL-(@>|9Lj|S+y^>n9B#BATTH-dMYIbBV3&t8#^^nu4Sh!qucdV3q_M^dO&E~QLSxs9Lff5 zWhJ?7m2%W@$W=9dSGAB{B~{Y5y-F0$Yc4uIaVmPHq9c3O5gV+OE9tDRKBwxIQUHW1 zy?u|(cyt;g2D3_R3MQd&acV)w&d`IXOQV=wwEHNub+igR>gQ1M^@L}O{wi90ys;LzfSVxM)V$uGK2yt?9AHmXl3wN|&OAW$fv zMTP|FVJslD7%KtFV5(jjtjh1J${)1y4F!430HL7Bz%Vap)X~WCN|lr^S4sJDm6Weg zDHPZ$-Q)=kTBu*R1684qM#6`M_2(o{>zNM|Qk;g>9kL71)JK=(Bfv@*CdYEig{-=p*-5k^Ec z{_xmsY1r609NsA0nr7l|CL)~1HX}N`y`N?c1%$AAnNEt1d9(BPC;*q`*GN?OOJY7P zQ1+AmVzJZ$f9x*-)qPXs}sR6}TW?R=OIQX3Uz&{}*E}R-aD?A{+(P$glB~2J3 zBQ15pV{9d<=Z)w{+deVILj z_Zu#CN*j#5mv%`_#??zhZtjO7sTxg^!Nmxwd=TJ9uFIVQISYaY%DMWQ<=sj6H^5xYiaz zybslaIgZB^&wK_k+pOYArs>WDqWQ*uzIIZ&Z8TgD@1%wQkqBc@C{U>@n}xg$!lpBp zO!9Ukdi_HK){Ef#NQdysE8hg)kAkBvtV&p!93DhRjK&+on{FV3U_#l1qvtz`+%jfv z?A&yb1Y{AfGM$K3Dc)>sL5;#}yvz9Fjh>Um@l=R(zHKL7dz+KZ1%ZYkN>iMJH$eNN z|3hD~!o)RNo(t+?Ha1etVzM10QTn2pt~3I>ri(sxmKVPbLZ)*vyI|rA6~o>E?;+tf zB4T2b$uOH?-5lRc8foUDHIzFe?mv@++`zjGmWGl6%hlp5KJ|Nyti5#US zhL7?Pm2D?dT3k+3wh09LK{+m#gb_%aa&DG~sjSynQbt-ztHw@Z?xF{bZz^r)lp5Fy zF-#Sh1-;JZOwY9#U|E*;lZ-2eIl(_pt)PK}Fa7ZezD(H{h!{lvE0O(1@6G;>2_%(T zfv2K~4JmGs#3E%Z+#JA;+81v2bySW$Q^ZfAZ0xvm_3v(eQj(m;$8PV}ucLD7dj*X_ zT(!6mJ#l2c_=+)!Ig*ZLH#@%UdFWbvPm@nN9=dM-&2kNMAs`}6%Q>fpa7rB8vj~e+ z=6zN9x~jY%dA(dgeU(6fAzXs1QLAC@8JYhQ^N4f}(yuVDNY^6$I`b_zK9mgG9%4FWyjBDXno+Ax)wZ-b5XO=4 zAcc0?5wF?U;aKj#y%l$g-?ZW01&La%gTcLjQ7)&Ncglr6wl_Pp&T_t$E|sFXuIh6A zMY&u6)o(&IK84Y zpL+}X+~{+MZRt*7ewU`>?FjP<(_4zFM7wO!N=Z>wmE@*%B}GH3BsUi%5hN?CUbMAJ zw7Vc$At|g-CAno?NnxkCx2`*K{^@!|YPPZ-if-ZA02|eG7dwe}w+*9UD}s03S~rtg z92y4O@)kVMz+~3W8FXmbhr5a>Iq6hxz zjR>Fh7jvz29;BGFn%*v&PJKNK-RIVHRd;BOc3;=;MZFe4#BDn}C?@DZ*s0XqXQ1-@ zn&#rT$zz(^iYQ&Ffg9Oj)~^@a+g4V^1;pWo*@$h-c8aOZwOqI2Y8NkB`$-htju&EU z;+tyLT(yE}aqq2HyRW8mBc}W66x`4FGoS44z27}MVtTGl3C6cRXdM{6Iy;g-x;=6~ zJ*!5R)Fm}?F?!+3xyU43LF`}KvAQUMph(^sVL7-F=%P zDH}0#-mXU|+nBYtfU^5E{~w#Y`>VtXSCHnS$OAXpZ<(5oE};SyN;G2~ekAfr>jFB@*5ErdsbMF^8w81%3u z3|BJ*Z)60yswemnO9D=ANI3NTHXhdV&d`;LphUhaG^WO8DCV$)&E+I{`~~xd^5IYr z{Rb;}>s}%r+!eGV0Ap4@EaCvWO^5u#JQWSK!nAhgoW>P z3*aT*2ez}YA3;5lO8g^u>lN6Qe&cPZh?4XzkZ=bVm`6auN?6AQLKL{{42hp55<>&w zU+VW&92pXaD8$y*U9bh1G;3^c&e{&VOwW`snZS1f^u|)DbPmUgtUc7A=>_zU zW#Fs`Tm}pU)CcS&kllM!MFF3o|CgZJ^yJ88;5e0o*v5rZ$1ef^AvHd6()41)*epsB z%X7_QxUW%GannQ_GsW$ig~|x=31$5RxP%TX<$~IN(*;faX3gZh&P6M0NLwsw$Kwq8QN(#%tCZteo71pFv$1GwxnP;?SXc1bQb z_*6!hyH!UVVu#d?{}!o5ZbU@GEB>dDvr>NK{fK|;mJ0uUt^nd3emvQZs75FLTSYtr zF%2iS`B2gXeI;I}lP{uwMcf!d>k7CLVQ*~C6@lZk0FIu;y8<|RZeRO<1{@!Y8-eJs z+g6E08>;g7kOCT7pz@0Vm6Rq>`4&Ls23D)p%y{^>0V+4Jy5%P3Uv6gg%PkL`iw7$J zrt|g=U^|xvux10yDb%C8SwMWN@C_$~Dc}XZae@b!h~~8dLavo*!nds%0JCk~QeZ?; zip?zxClW@q0cMBbVKb0;y3M4my;i&0Hhf|UK6(FT^d1Ls=6(&G83 z4Ct14b+(+J?sMO@Uw|wS{u|i#1n?<)2uQW!0j;U3q86>emT$BaqG5zo0Vc&OD?mk< zR=i5_3avV-BrT+VE!`r@JFj-HjQRO-fcXuPqe8p8-1(VT10C*7%#66J$v3-HUYqDT{eorW<2ZNf1oxk2- z=Fal~fXm26h*Qi(pHB`h2eBDQkzjt4>ZtNsEO%f-D2TrWYY))^|4o(;g8GY?-)xvw zSz>;6ApxHx$5ox>1EJt$>#5l!Tn7HXQBPtAg2gDY+kG1Q48c%2=qJxWysB~-|2)}v z31mIIU$Exe6@@yU_%c<+5 z1k$W~9o-QHw`$r!vbVp7s3`Ju+HNP%DdL!H7oa0W@EVtPAhryUYo?Zy5WkxyZ8^SzT>&=c>n0B)pzl3*Yw@R_g3`X#gA0<-Nm0- zr|$+>M9esDrc|_j>tbvUbfl?ztWU22zn(GsH-mcV?|7hkbi)+qRS zjDn9@qujQ{LU9zh`MwHExcSKnOJFPXiIymihpy6y&}rh?)5RVVWES^a5LR)|1&UT_ zDznPTf3{UlJbxjYliUsVLy_P(Sep7(rf^gxRZSKVqm%=VBHv5tgsw-RB zYWUN|MP{rxdr?JcWwYmevpvTz!%d@CO0#?1-shk-p~WtrfJaqceiN(y+*y0xU;CNg zE4?I|_iN@lv<+UjJc;YZAE&|t5Aea^i$s=)*fI=jb7YNC%OdL_`>TKc$XR9W^2C zL1bPFi3m4N#N9FeOS)A*kMg|#6diErot)C(hUbZ$U5ZS7cIg zai{}NQ%}WNBhx*JgBw|KBt;(@>K_s69td@cN0dbsHwDeJ`3aHq*AF_-w-2j z=G-MQbS$tKKP5^DNcsd zTof-hcqK=r|#lVLY59`Ib^R9${*7y0FHn9 zSxT#oXq)7ANd~V6NV!v{RJWbm;cP-An|`X|u`d^*@h!M39+pmGXVG-0=4UfY*11$)>F2B5NsrNf3{K=s zCTL3{OQ(n^*_TM0aw5%JjJf4r{T$Td=V|HLS?oE{uWT>{uwd{SQqWK2G!Zc%U!~OR zME-!t_lf+7$ooX_`E`gN=ZG_rd^csq5T>mIt`Q+r&mV%A-e^?(u->SDp*F1IanpqZ zfZW7BW-gkc-?XMPv<;ROXN&NwJRV=*A!GQ3HvHiH`V0F8$sNS8pA?IV=9Y$!NJU7a z_Xtke&L`Dl5XHGr}&Py_a?nz?+)*7Z>@LMJMPu-