diff --git a/._DensityIntegrationUncertaintyQuantification.py b/._DensityIntegrationUncertaintyQuantification.py new file mode 100755 index 0000000..e1118fc Binary files /dev/null and b/._DensityIntegrationUncertaintyQuantification.py differ diff --git a/._sample_script.py b/._sample_script.py new file mode 100755 index 0000000..d09e7de Binary files /dev/null and b/._sample_script.py differ diff --git a/DensityIntegrationUncertaintyQuantification.py b/DensityIntegrationUncertaintyQuantification.py new file mode 100755 index 0000000..c2aea57 --- /dev/null +++ b/DensityIntegrationUncertaintyQuantification.py @@ -0,0 +1,507 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Mon Apr 3 15:51:35 2017 + +@author: jiachengzhang +""" + +import numpy as np +import sys +import scipy.sparse as scysparse +import scipy.sparse.linalg as splinalg +import scipy.linalg as linalg + + +def Density_integration_Poisson_uncertainty(Xn,Yn,fluid_mask,grad_x,grad_y,dirichlet_label,dirichlet_value, + uncertainty_quantification=True, sigma_grad_x=None, sigma_grad_y=None, sigma_dirichlet=None): + # Evaluate the density field from the gradient fields by solving the Poisson equation. + # The uncertainty of the density field can be quantified. + """ + Inputs: + Xn,Yn: 2d array of mesh grid. + fluid_mask: 2d array of binary mask of flow field. Boundary points are considered in the flow (mask should be True) + grad_x, grad_y: 2d array of gradient field. + dirichlet_label: 2d array of binary mask indicating the Dirichlet BC locations. Ohterwise Neumann. At least one point should be dirichlet. + dirichlet_value: 2d array of diriclet BC values. + uncertainty_quantification: True to perform the uncertainty quantification, False only perform integration. + sigma_grad_x, sigma_grad_y, sigma_dirichlet: the uncertainty given as std of the input fields. 2d array of fields. + Returns: + Pn: the integrated density field. + sigma_Pn: the uncertainty (std) of the integrated density field. + """ + + Ny, Nx = np.shape(Xn) + dx = Xn[1,1] - Xn[0,0] + dy = Yn[1,1] - Yn[0,0] + invdx = 1.0/dx + invdy = 1.0/dy + invdx2 = invdx**2 + invdy2 = invdy**2 + + fluid_mask_ex = np.zeros((Ny+2,Nx+2)).astype('bool') + fluid_mask_ex[1:-1,1:-1] = fluid_mask + + grad_x_ex = np.zeros((Ny+2,Nx+2)) + grad_x_ex[1:-1,1:-1] = grad_x + grad_y_ex = np.zeros((Ny+2,Nx+2)) + grad_y_ex[1:-1,1:-1] = grad_y + dirichlet_label_ex = np.zeros((Ny+2,Nx+2)).astype('bool') + dirichlet_label_ex[1:-1,1:-1] = dirichlet_label + dirichlet_value_ex = np.zeros((Ny+2,Nx+2)) + dirichlet_value_ex[1:-1,1:-1] = dirichlet_value + + sigma_grad_x_ex = np.zeros((Ny+2,Nx+2)) + sigma_grad_y_ex = np.zeros((Ny+2,Nx+2)) + sigma_dirichlet_ex = np.zeros((Ny+2,Nx+2)) + if uncertainty_quantification==True: + sigma_grad_x_ex[1:-1,1:-1] = sigma_grad_x + sigma_grad_y_ex[1:-1,1:-1] = sigma_grad_y + sigma_dirichlet_ex[1:-1,1:-1] = sigma_dirichlet + + # Generate the linear operator. + j,i = np.where(fluid_mask_ex==True) + Npts = len(j) + fluid_index = -np.ones(fluid_mask_ex.shape).astype('int') + fluid_index[j,i] = range(Npts) + iC = fluid_index[j,i] + iC_label = dirichlet_label_ex[j,i] + iE = fluid_index[j,i+1] + iW = fluid_index[j,i-1] + iN = fluid_index[j+1,i] + iS = fluid_index[j-1,i] + + LaplacianOperator = scysparse.csr_matrix((Npts,Npts),dtype=np.float) + # RHS = np.zeros(Npts) + # Var_RHS = np.zeros(Npts) # variance + + # Also generate the linear operator which maps from the grad_x, grad_y, dirichlet val to the rhs. + Map_grad_x = scysparse.csr_matrix((Npts,Npts),dtype=np.float) + Map_grad_y = scysparse.csr_matrix((Npts,Npts),dtype=np.float) + Map_dirichlet_val = scysparse.csr_matrix((Npts,Npts),dtype=np.float) + + # First, construct the linear operator and RHS as if all they are all Nuemanan Bc + + # if the east and west nodes are inside domain + loc = (iE!=-1)*(iW!=-1) + LaplacianOperator[iC[loc],iC[loc]] += -2.0*invdx2 + LaplacianOperator[iC[loc],iE[loc]] += +1.0*invdx2 + LaplacianOperator[iC[loc],iW[loc]] += +1.0*invdx2 + # RHS[iC[loc]] += (grad_x_ex[j[loc],i[loc]+1] - grad_x_ex[j[loc],i[loc]-1])*(invdx*0.5) + # Var_RHS[iC[loc]] += (invdx*0.5)**2 * (sigma_grad_x_ex[j[loc],i[loc]+1]**2 + sigma_grad_x_ex[j[loc],i[loc]-1]**2) + Map_grad_x[iC[loc],iE[loc]] += invdx*0.5 + Map_grad_x[iC[loc],iW[loc]] += -invdx*0.5 + + # if the east node is ouside domian + loc = (iE==-1) + LaplacianOperator[iC[loc],iC[loc]] += -2.0*invdx2 + LaplacianOperator[iC[loc],iW[loc]] += +2.0*invdx2 + # RHS[iC[loc]] += -(grad_x_ex[j[loc],i[loc]] + grad_x_ex[j[loc],i[loc]-1])*invdx + # Var_RHS[iC[loc]] += invdx**2 * (sigma_grad_x_ex[j[loc],i[loc]]**2 + sigma_grad_x_ex[j[loc],i[loc]-1]**2) + Map_grad_x[iC[loc],iC[loc]] += -invdx + Map_grad_x[iC[loc],iW[loc]] += -invdx + + # if the west node is ouside domian + loc = (iW==-1) + LaplacianOperator[iC[loc],iC[loc]] += -2.0*invdx2 + LaplacianOperator[iC[loc],iE[loc]] += +2.0*invdx2 + # RHS[iC[loc]] += (grad_x_ex[j[loc],i[loc]] + grad_x_ex[j[loc],i[loc]+1])*invdx + # Var_RHS[iC[loc]] += invdx**2 * (sigma_grad_x_ex[j[loc],i[loc]]**2 + sigma_grad_x_ex[j[loc],i[loc]+1]**2) + Map_grad_x[iC[loc],iC[loc]] += invdx + Map_grad_x[iC[loc],iE[loc]] += invdx + + # if the north and south nodes are inside domain + loc = (iN!=-1)*(iS!=-1) + LaplacianOperator[iC[loc],iC[loc]] += -2.0*invdy2 + LaplacianOperator[iC[loc],iN[loc]] += +1.0*invdy2 + LaplacianOperator[iC[loc],iS[loc]] += +1.0*invdy2 + # RHS[iC[loc]] += (grad_y_ex[j[loc]+1,i[loc]] - grad_y_ex[j[loc]-1,i[loc]])*(invdy*0.5) + # Var_RHS[iC[loc]] += (invdy*0.5)**2 * (sigma_grad_y_ex[j[loc]+1,i[loc]]**2 + sigma_grad_y_ex[j[loc]-1,i[loc]]**2) + Map_grad_y[iC[loc],iN[loc]] += invdy*0.5 + Map_grad_y[iC[loc],iS[loc]] += -invdy*0.5 + + # if the north node is ouside domian + loc = (iN==-1) + LaplacianOperator[iC[loc],iC[loc]] += -2.0*invdy2 + LaplacianOperator[iC[loc],iS[loc]] += +2.0*invdy2 + # RHS[iC[loc]] += -(grad_y_ex[j[loc],i[loc]] + grad_y_ex[j[loc]-1,i[loc]])*invdy + # Var_RHS[iC[loc]] += invdy**2 * (sigma_grad_y_ex[j[loc],i[loc]]**2 + sigma_grad_y_ex[j[loc]-1,i[loc]]**2) + Map_grad_y[iC[loc],iC[loc]] += -invdy + Map_grad_y[iC[loc],iS[loc]] += -invdy + + # if the south node is ouside domian + loc = (iS==-1) + LaplacianOperator[iC[loc],iC[loc]] += -2.0*invdy2 + LaplacianOperator[iC[loc],iN[loc]] += +2.0*invdy2 + # RHS[iC[loc]] += (grad_y_ex[j[loc],i[loc]] + grad_y_ex[j[loc]+1,i[loc]])*invdy + # Var_RHS[iC[loc]] += invdy**2 * (sigma_grad_y_ex[j[loc],i[loc]]**2 + sigma_grad_y_ex[j[loc]+1,i[loc]]**2) + Map_grad_y[iC[loc],iC[loc]] += invdy + Map_grad_y[iC[loc],iN[loc]] += invdy + + # Then change the boundary conidtion at locatiosn of Dirichlet. + loc = (iC_label==True) + LaplacianOperator[iC[loc],:] = 0.0 + LaplacianOperator[iC[loc],iC[loc]] = 1.0*invdx2 + # RHS[iC[loc]] = dirichlet_value_ex[j[loc],i[loc]] * invdx2 + # Var_RHS[iC[loc]] = sigma_dirichlet_ex[j[loc],i[loc]]**2 * invdx2**2 + Map_grad_x[iC[loc],:] = 0.0 + Map_grad_y[iC[loc],:] = 0.0 + Map_dirichlet_val[iC[loc],iC[loc]] = 1.0*invdx2 + + LaplacianOperator.eliminate_zeros() + Map_grad_x.eliminate_zeros() + Map_grad_y.eliminate_zeros() + Map_dirichlet_val.eliminate_zeros() + + # Solve for the field. + grad_x_vect = grad_x_ex[j,i] + grad_y_vect = grad_y_ex[j,i] + dirichlet_val_vect = dirichlet_value_ex[j,i] + RHS = Map_grad_x.dot(grad_x_vect) + Map_grad_y.dot(grad_y_vect) + Map_dirichlet_val.dot(dirichlet_val_vect) + # Solve the linear system + Pn_ex = np.zeros((Ny+2,Nx+2)) + p_vect = splinalg.spsolve(LaplacianOperator, RHS) + Pn_ex[j,i] = p_vect + + # Uncertainty propagation + sigma_Pn_ex = np.zeros((Ny+2,Nx+2)) + if uncertainty_quantification==True: + # Propagate to get the covariance matrix for RHS + Cov_grad_x = scysparse.diags(sigma_grad_x_ex[j,i]**2,shape=(Npts,Npts),format='csr') + Cov_grad_y = scysparse.diags(sigma_grad_y_ex[j,i]**2,shape=(Npts,Npts),format='csr') + Cov_dirichlet_val = scysparse.diags(sigma_dirichlet_ex[j,i]**2,shape=(Npts,Npts),format='csr') + Cov_RHS = Map_grad_x*Cov_grad_x*Map_grad_x.transpose() + Map_grad_y*Cov_grad_y*Map_grad_y.transpose() + \ + Map_dirichlet_val*Cov_dirichlet_val*Map_dirichlet_val.transpose() + + Laplacian_inv = linalg.inv(LaplacianOperator.A) + Cov_p = np.matmul(np.matmul(Laplacian_inv, Cov_RHS.A), Laplacian_inv.T) + Var_p_vect = np.diag(Cov_p) + + sigma_Pn_ex[j,i] = Var_p_vect**0.5 + + return Pn_ex[1:-1,1:-1], sigma_Pn_ex[1:-1,1:-1] + + + +def Density_integration_WLS_uncertainty(Xn,Yn,fluid_mask,grad_x,grad_y,dirichlet_label,dirichlet_value, + uncertainty_quantification=True, sigma_grad_x=None, sigma_grad_y=None, sigma_dirichlet=None): + # Evaluate the density field from the gradient fields by solving the WLS system. + # The uncertainty of the density field is also quantified and returned. + """ + Inputs: + Xn,Yn: 2d array of mesh grid. + fluid_mask: 2d array of binary mask of flow field. Boundary points are considered in the flow (mask should be True) + grad_x, grad_y: 2d array of gradient field. + dirichlet_label: 2d array of binary mask indicating the Dirichlet BC locations. Ohterwise Neumann. At least one point should be dirichlet. + dirichlet_value: 2d array of diriclet BC values. + uncertainty_quantification: True to perform the uncertainty quantification, False only perform integration. + sigma_grad_x, sigma_grad_y, sigma_dirichlet: the uncertainty given as std of the input fields. 2d array of fields. + Returns: + Pn: the integrated density field. + sigma_Pn: the uncertainty (std) of the integrated density field. + """ + + Ny, Nx = np.shape(Xn) + dx = Xn[1,1] - Xn[0,0] + dy = Yn[1,1] - Yn[0,0] + invdx = 1.0/dx + invdy = 1.0/dy + invdx2 = invdx**2 + invdy2 = invdy**2 + + fluid_mask_ex = np.zeros((Ny+2,Nx+2)).astype('bool') + fluid_mask_ex[1:-1,1:-1] = fluid_mask + + dirichlet_label_ex = np.zeros((Ny+2,Nx+2)).astype('bool') + dirichlet_label_ex[1:-1,1:-1] = dirichlet_label + dirichlet_value_ex = np.zeros((Ny+2,Nx+2)) + dirichlet_value_ex[1:-1,1:-1] = dirichlet_value + + grad_x_ex = np.zeros((Ny+2,Nx+2)) + grad_x_ex[1:-1,1:-1] = grad_x + grad_y_ex = np.zeros((Ny+2,Nx+2)) + grad_y_ex[1:-1,1:-1] = grad_y + + sigma_grad_x_ex = np.zeros((Ny+2,Nx+2)) + sigma_grad_y_ex = np.zeros((Ny+2,Nx+2)) + sigma_dirichlet_ex = np.zeros((Ny+2,Nx+2)) + sigma_grad_x_ex[1:-1,1:-1] = sigma_grad_x + sigma_grad_y_ex[1:-1,1:-1] = sigma_grad_y + sigma_dirichlet_ex[1:-1,1:-1] = sigma_dirichlet + + # Generate the index for mapping the pressure and pressure gradients. + j,i = np.where(fluid_mask_ex==True) + Npts = len(j) + fluid_index = -np.ones(fluid_mask_ex.shape).astype('int') + fluid_index[j,i] = range(Npts) + + # Generate the mask for the gradients + fluid_mask_Gx_ex = np.logical_and(fluid_mask_ex[:,1:],fluid_mask_ex[:,:-1]) + fluid_mask_Gy_ex = np.logical_and(fluid_mask_ex[1:,:],fluid_mask_ex[:-1,:]) + + # Generate the linear operator and the mapping matrix for generating the rhs. + # For Gx + jx,ix = np.where(fluid_mask_Gx_ex==True) + Npts_x = len(jx) + iC = fluid_index[jx,ix] + iE = fluid_index[jx,ix+1] + Operator_Gx = scysparse.csr_matrix((Npts_x,Npts),dtype=np.float) + Map_Gx = scysparse.csr_matrix((Npts_x,Npts),dtype=np.float) + Operator_Gx[range(Npts_x),iC] += -invdx + Operator_Gx[range(Npts_x),iE] += invdx + Map_Gx[range(Npts_x),iC] += 0.5 + Map_Gx[range(Npts_x),iE] += 0.5 + + # For Gy + jy,iy = np.where(fluid_mask_Gy_ex==True) + Npts_y = len(jy) + iC = fluid_index[jy,iy] + iN = fluid_index[jy+1,iy] + Operator_Gy = scysparse.csr_matrix((Npts_y,Npts),dtype=np.float) + Map_Gy = scysparse.csr_matrix((Npts_y,Npts),dtype=np.float) + Operator_Gy[range(Npts_y),iC] += -invdy + Operator_Gy[range(Npts_y),iN] += invdy + Map_Gy[range(Npts_y),iC] += 0.5 + Map_Gy[range(Npts_y),iN] += 0.5 + + # For Dirichlet BC + j_d, i_d = np.where(dirichlet_label_ex==True) + Npts_d = len(j_d) + iC = fluid_index[j_d,i_d] + Operator_d = scysparse.csr_matrix((Npts_d,Npts),dtype=np.float) + Map_d = scysparse.eye(Npts_d,Npts_d,format='csr') * invdx + Operator_d[range(Npts_d),iC] += 1.0*invdx + # dirichlet value vector and cov + dirichlet_vect = dirichlet_value_ex[j_d,i_d] + dirichlet_sigma_vect = sigma_dirichlet_ex[j_d,i_d] + cov_dirichlet = scysparse.diags(dirichlet_sigma_vect**2, format='csr') + # Generate the vector and cov for pgrad. + pgrad_x_vect = grad_x_ex[j,i] + pgrad_y_vect = grad_y_ex[j,i] + pgrad_vect = np.concatenate((pgrad_x_vect,pgrad_y_vect)) + cov_pgrad_x = sigma_grad_x_ex[j,i]**2 + cov_pgrad_y = sigma_grad_y_ex[j,i]**2 + cov_pgrad_vect = np.concatenate((cov_pgrad_x, cov_pgrad_y)) + cov_pgrad = scysparse.diags(cov_pgrad_vect, format='csr') + + # Construct the full operator. + Operator_GLS = scysparse.bmat([[Operator_Gx],[Operator_Gy],[Operator_d]]) + + # Construct the full mapping matrics and get the rhs. + Map_pgrad = scysparse.bmat([[Map_Gx, None],[None, Map_Gy],[scysparse.csr_matrix((Npts_d,Npts),dtype=np.float), None]]) + Map_dirichlet = scysparse.bmat([[scysparse.csr_matrix((Npts_x+Npts_y,Npts_d),dtype=np.float)],[Map_d]]) + rhs = Map_pgrad.dot(pgrad_vect) + Map_dirichlet.dot(dirichlet_vect) + + # Evaluate the covriance matrix for the rhs + cov_rhs = Map_pgrad * cov_pgrad * Map_pgrad.transpose() + Map_dirichlet * cov_dirichlet * Map_dirichlet.transpose() + + # Solve for the WLS solution + weights_vect = cov_rhs.diagonal()**(-1) + weights_matrix = scysparse.diags(weights_vect,format='csr') + # Operator_WLS = weights_matrix * Operator_GLS + # rhs_WLS = weights_matrix.dot(rhs) + sys_LHS = Operator_GLS.transpose() * weights_matrix * Operator_GLS + sys_rhs = (Operator_GLS.transpose() * weights_matrix).dot(rhs) + + # Get the solution from lsqr + # p_vect_wls = splinalg.lsqr(Operator_WLS,rhs_WLS)[0] + # Pn_WLS = np.zeros(fluid_mask_ex.shape) + # Pn_WLS[j,i] = p_vect_wls + + # Solve for the WLS solution + p_vect_wls = splinalg.spsolve(sys_LHS, sys_rhs) + Pn_WLS_ex = np.zeros(fluid_mask_ex.shape) + Pn_WLS_ex[j,i] = p_vect_wls + + # Perform the uncertainty propagation + sigma_Pn_ex = np.zeros((Ny+2,Nx+2)) + if uncertainty_quantification == True: + cov_sys_rhs = (Operator_GLS.transpose() * weights_matrix) * cov_rhs * (Operator_GLS.transpose() * weights_matrix).transpose() + sys_LHS_inv = linalg.inv(sys_LHS.A) + Cov_p = np.matmul(np.matmul(sys_LHS_inv, cov_sys_rhs.A), sys_LHS_inv.T) + Var_p_vect = np.diag(Cov_p) + sigma_Pn_ex[j,i] = Var_p_vect**0.5 + + return Pn_WLS_ex[1:-1,1:-1], sigma_Pn_ex[1:-1,1:-1] + + + +def Density_integration_WLS_uncertainty_weighted_average(Xn,Yn,fluid_mask,grad_x,grad_y,dirichlet_label,dirichlet_value, + uncertainty_quantification=True, sigma_grad_x=None, sigma_grad_y=None, sigma_dirichlet=None): + # Evaluate the density field from the gradient fields by solving the WLS system. + # The uncertainty of the density field is also quantified and returned. + # The gradient interpolation (from grid points to staggered location) is done by a weighted average approach + # which minimizes the sum of squared bias error and random error. + """ + Inputs: + Xn,Yn: 2d array of mesh grid. + fluid_mask: 2d array of binary mask of flow field. Boundary points are considered in the flow (mask should be True) + grad_x, grad_y: 2d array of gradient field. + dirichlet_label: 2d array of binary mask indicating the Dirichlet BC locations. Ohterwise Neumann. At least one point should be dirichlet. + dirichlet_value: 2d array of diriclet BC values. + uncertainty_quantification: True to perform the uncertainty quantification, False only perform integration. + sigma_grad_x, sigma_grad_y, sigma_dirichlet: the uncertainty given as std of the input fields. 2d array of fields. + Returns: + Pn: the integrated density field. + sigma_Pn: the uncertainty (std) of the integrated density field. + """ + + Ny, Nx = np.shape(Xn) + dx = Xn[1,1] - Xn[0,0] + dy = Yn[1,1] - Yn[0,0] + invdx = 1.0/dx + invdy = 1.0/dy + invdx2 = invdx**2 + invdy2 = invdy**2 + + fluid_mask_ex = np.zeros((Ny+2,Nx+2)).astype('bool') + fluid_mask_ex[1:-1,1:-1] = fluid_mask + + dirichlet_label_ex = np.zeros((Ny+2,Nx+2)).astype('bool') + dirichlet_label_ex[1:-1,1:-1] = dirichlet_label + dirichlet_value_ex = np.zeros((Ny+2,Nx+2)) + dirichlet_value_ex[1:-1,1:-1] = dirichlet_value + + grad_x_ex = np.zeros((Ny+2,Nx+2)) + grad_x_ex[1:-1,1:-1] = grad_x + grad_y_ex = np.zeros((Ny+2,Nx+2)) + grad_y_ex[1:-1,1:-1] = grad_y + + sigma_grad_x_ex = np.zeros((Ny+2,Nx+2)) + sigma_grad_y_ex = np.zeros((Ny+2,Nx+2)) + sigma_dirichlet_ex = np.zeros((Ny+2,Nx+2)) + sigma_grad_x_ex[1:-1,1:-1] = sigma_grad_x + sigma_grad_y_ex[1:-1,1:-1] = sigma_grad_y + sigma_dirichlet_ex[1:-1,1:-1] = sigma_dirichlet + + # Generate the index for mapping the pressure and pressure gradients. + j,i = np.where(fluid_mask_ex==True) + Npts = len(j) + fluid_index = -np.ones(fluid_mask_ex.shape).astype('int') + fluid_index[j,i] = range(Npts) + + # Generate the mask for the gradients + fluid_mask_Gx_ex = np.logical_and(fluid_mask_ex[:,1:],fluid_mask_ex[:,:-1]) + fluid_mask_Gy_ex = np.logical_and(fluid_mask_ex[1:,:],fluid_mask_ex[:-1,:]) + + # Generate the linear operator and the mapping matrix for generating the rhs. + # For Gx + jx,ix = np.where(fluid_mask_Gx_ex==True) + Npts_x = len(jx) + iC = fluid_index[jx,ix] + iE = fluid_index[jx,ix+1] + Operator_Gx = scysparse.csr_matrix((Npts_x,Npts),dtype=np.float) + Operator_Gx[range(Npts_x),iC] += -invdx + Operator_Gx[range(Npts_x),iE] += invdx + # The Mapping Gx that maps grided gradients to staggered location is by weighted average. + Map_Gx = scysparse.csr_matrix((Npts_x,Npts),dtype=np.float) + V_C = grad_x_ex[jx,ix] + V_E = grad_x_ex[jx,ix+1] + sigma_C = sigma_grad_x_ex[jx,ix] + sigma_E = sigma_grad_x_ex[jx,ix+1] + weight_C = ((V_C-V_E)**2 + 2*sigma_E**2) / (2*(V_C-V_E)**2 + 2*sigma_C**2 + 2*sigma_E**2) + weight_E = 1.0 - weight_C + Map_Gx[range(Npts_x),iC] += weight_C + Map_Gx[range(Npts_x),iE] += weight_E + + # For Gy + jy,iy = np.where(fluid_mask_Gy_ex==True) + Npts_y = len(jy) + iC = fluid_index[jy,iy] + iN = fluid_index[jy+1,iy] + Operator_Gy = scysparse.csr_matrix((Npts_y,Npts),dtype=np.float) + Operator_Gy[range(Npts_y),iC] += -invdy + Operator_Gy[range(Npts_y),iN] += invdy + # The Mapping Gy that maps grided gradients to staggered location is by weighted average. + Map_Gy = scysparse.csr_matrix((Npts_y,Npts),dtype=np.float) + V_C = grad_y_ex[jy,iy] + V_N = grad_y_ex[jy+1,iy] + sigma_C = sigma_grad_y_ex[jy,iy] + sigma_N = sigma_grad_y_ex[jy+1,iy] + weight_C = ((V_C-V_N)**2 + 2*sigma_N**2) / (2*(V_C-V_N)**2 + 2*sigma_C**2 + 2*sigma_N**2) + weight_N = 1.0 - weight_C + Map_Gy[range(Npts_y),iC] += weight_C + Map_Gy[range(Npts_y),iN] += weight_N + + # For Dirichlet BC + j_d, i_d = np.where(dirichlet_label_ex==True) + Npts_d = len(j_d) + iC = fluid_index[j_d,i_d] + Operator_d = scysparse.csr_matrix((Npts_d,Npts),dtype=np.float) + Map_d = scysparse.eye(Npts_d,Npts_d,format='csr') * invdx + Operator_d[range(Npts_d),iC] += 1.0*invdx + # dirichlet value vector and cov + dirichlet_vect = dirichlet_value_ex[j_d,i_d] + dirichlet_sigma_vect = sigma_dirichlet_ex[j_d,i_d] + cov_dirichlet = scysparse.diags(dirichlet_sigma_vect**2, format='csr') + # Generate the vector and cov for pgrad. + pgrad_x_vect = grad_x_ex[j,i] + pgrad_y_vect = grad_y_ex[j,i] + pgrad_vect = np.concatenate((pgrad_x_vect,pgrad_y_vect)) + cov_pgrad_x = sigma_grad_x_ex[j,i]**2 + cov_pgrad_y = sigma_grad_y_ex[j,i]**2 + cov_pgrad_vect = np.concatenate((cov_pgrad_x, cov_pgrad_y)) + cov_pgrad = scysparse.diags(cov_pgrad_vect, format='csr') + + # Construct the full operator. + Operator_GLS = scysparse.bmat([[Operator_Gx],[Operator_Gy],[Operator_d]]) + + # Construct the full mapping matrics and get the rhs. + Map_pgrad = scysparse.bmat([[Map_Gx, None],[None, Map_Gy],[scysparse.csr_matrix((Npts_d,Npts),dtype=np.float), None]]) + Map_dirichlet = scysparse.bmat([[scysparse.csr_matrix((Npts_x+Npts_y,Npts_d),dtype=np.float)],[Map_d]]) + rhs = Map_pgrad.dot(pgrad_vect) + Map_dirichlet.dot(dirichlet_vect) + + # Evaluate the covriance matrix for the rhs + cov_rhs = Map_pgrad * cov_pgrad * Map_pgrad.transpose() + Map_dirichlet * cov_dirichlet * Map_dirichlet.transpose() + + # Solve for the WLS solution + weights_vect = cov_rhs.diagonal()**(-1) + weights_matrix = scysparse.diags(weights_vect,format='csr') + # Operator_WLS = weights_matrix * Operator_GLS + # rhs_WLS = weights_matrix.dot(rhs) + sys_LHS = Operator_GLS.transpose() * weights_matrix * Operator_GLS + sys_rhs = (Operator_GLS.transpose() * weights_matrix).dot(rhs) + + # Get the solution from lsqr + # p_vect_wls = splinalg.lsqr(Operator_WLS,rhs_WLS)[0] + # Pn_WLS = np.zeros(fluid_mask_ex.shape) + # Pn_WLS[j,i] = p_vect_wls + + # Solve for the WLS solution + p_vect_wls = splinalg.spsolve(sys_LHS, sys_rhs) + Pn_WLS_ex = np.zeros(fluid_mask_ex.shape) + Pn_WLS_ex[j,i] = p_vect_wls + + # Perform the uncertainty propagation + sigma_Pn_ex = np.zeros((Ny+2,Nx+2)) + if uncertainty_quantification == True: + cov_sys_rhs = (Operator_GLS.transpose() * weights_matrix) * cov_rhs * (Operator_GLS.transpose() * weights_matrix).transpose() + sys_LHS_inv = linalg.inv(sys_LHS.A) + Cov_p = np.matmul(np.matmul(sys_LHS_inv, cov_sys_rhs.A), sys_LHS_inv.T) + Var_p_vect = np.diag(Cov_p) + sigma_Pn_ex[j,i] = Var_p_vect**0.5 + + return Pn_WLS_ex[1:-1,1:-1], sigma_Pn_ex[1:-1,1:-1] + + + + + + + + + + + + + + + + + + + + + + + diff --git a/__pycache__/DensityIntegrationUncertaintyQuantification.cpython-38.pyc b/__pycache__/DensityIntegrationUncertaintyQuantification.cpython-38.pyc new file mode 100755 index 0000000..39e0b06 Binary files /dev/null and b/__pycache__/DensityIntegrationUncertaintyQuantification.cpython-38.pyc differ diff --git a/__pycache__/helper_functions.cpython-38.pyc b/__pycache__/helper_functions.cpython-38.pyc new file mode 100755 index 0000000..24413b3 Binary files /dev/null and b/__pycache__/helper_functions.cpython-38.pyc differ diff --git a/__pycache__/loadmat_functions.cpython-38.pyc b/__pycache__/loadmat_functions.cpython-38.pyc new file mode 100755 index 0000000..a1448e4 Binary files /dev/null and b/__pycache__/loadmat_functions.cpython-38.pyc differ diff --git a/loadmat_functions.py b/loadmat_functions.py new file mode 100755 index 0000000..721cb44 --- /dev/null +++ b/loadmat_functions.py @@ -0,0 +1,35 @@ +# The following functions convert an object to a struct so that it can be saved to a mat file +import scipy.io as sio + +def loadmat(filename): + ''' + this function should be called instead of direct spio.loadmat + as it cures the problem of not properly recovering python dictionaries + from mat files. It calls the function check keys to cure all entries + which are still mat-objects + ''' + data = sio.loadmat(filename, struct_as_record=False, squeeze_me=True) + return _check_keys(data) + +def _check_keys(dict): + ''' + checks if entries in dictionary are mat-objects. If yes + todict is called to change them to nested dictionaries + ''' + for key in dict: + if isinstance(dict[key], sio.matlab.mio5_params.mat_struct): + dict[key] = _todict(dict[key]) + return dict + +def _todict(matobj): + ''' + A recursive function which constructs from matobjects nested dictionaries + ''' + dict = {} + for strg in matobj._fieldnames: + elem = matobj.__dict__[strg] + if isinstance(elem, sio.matlab.mio5_params.mat_struct): + dict[strg] = _todict(elem) + else: + dict[strg] = elem + return dict \ No newline at end of file diff --git a/sample-displacements.mat b/sample-displacements.mat new file mode 100755 index 0000000..7d57dcb Binary files /dev/null and b/sample-displacements.mat differ diff --git a/sample-result.mat b/sample-result.mat new file mode 100755 index 0000000..01d9fed Binary files /dev/null and b/sample-result.mat differ diff --git a/sample-result.png b/sample-result.png new file mode 100755 index 0000000..7100a48 Binary files /dev/null and b/sample-result.png differ diff --git a/sample_script.py b/sample_script.py new file mode 100755 index 0000000..4ed1e75 --- /dev/null +++ b/sample_script.py @@ -0,0 +1,285 @@ +#!/usr/bin/env python2 +# -*- coding: utf-8 -*- + +import numpy as np +import matplotlib +# matplotlib.use('Agg') +import matplotlib.pyplot as plt +import sys +import os +from scipy.io import savemat +from scipy.io import loadmat +import timeit + +# import density integration functions +from DensityIntegrationUncertaintyQuantification import Density_integration_Poisson_uncertainty +from DensityIntegrationUncertaintyQuantification import Density_integration_WLS_uncertainty + +# import helper functions +import loadmat_functions + +def set_bc(X, experimental_parameters): + + # define dirichlet boundary points (minimum one point) - here defined to be all boundaries + dirichlet_label = create_border_array(X.shape[0], X.shape[1], fill_val=1) + + # set density and density uncertainty at these boundary points + rho_dirichlet = experimental_parameters['rho_0'] * dirichlet_label + sigma_rho_dirichlet = experimental_parameters['sigma_rho_0'] * dirichlet_label + + # convert to boolean array + dirichlet_label = dirichlet_label.astype('bool') + + return dirichlet_label, rho_dirichlet, sigma_rho_dirichlet + + +def convert_displacements_to_physical_units(X_pix, Y_pix, U, V, sigma_U, sigma_V, experimental_parameters, mask): + + # calculate co-ordinates in the density gradient plane (m) + X = (X_pix - experimental_parameters['x0']) * experimental_parameters['pixel_pitch'] * 1 / experimental_parameters['magnification_grad'] + Y = (Y_pix - experimental_parameters['y0']) * experimental_parameters['pixel_pitch'] * 1 / experimental_parameters['magnification_grad'] + + # calculate density gradients from displacements (kg/m^4) + rho_x = U * experimental_parameters['pixel_pitch'] * 1 / ( + experimental_parameters['Z_D'] * experimental_parameters['magnification'] * + experimental_parameters[ + 'delta_z']) * 1 / experimental_parameters['gladstone_dale'] * experimental_parameters['n_0'] + rho_y = V * experimental_parameters['pixel_pitch'] * 1 / ( + experimental_parameters['Z_D'] * experimental_parameters['magnification'] * + experimental_parameters[ + 'delta_z']) * 1 / experimental_parameters['gladstone_dale'] * experimental_parameters['n_0'] + + # calculate density gradient uncertainty from displacement uncertainty (kg/m^4) + sigma_rho_x = sigma_U * experimental_parameters['pixel_pitch'] * \ + 1 / ( + experimental_parameters['Z_D'] * experimental_parameters['magnification'] * experimental_parameters[ + 'delta_z']) * \ + 1 / experimental_parameters['gladstone_dale'] * experimental_parameters['n_0'] + sigma_rho_y = sigma_V * experimental_parameters['pixel_pitch'] * \ + 1 / ( + experimental_parameters['Z_D'] * experimental_parameters['magnification'] * experimental_parameters[ + 'delta_z']) * \ + 1 / experimental_parameters['gladstone_dale'] * experimental_parameters['n_0'] + + # -------------------------- + # enforce mask and ensure valid values + # -------------------------- + # set uncertainties in masked regions to be zero (e.g. for image matching) + sigma_rho_x[mask == False] = 0.0 + sigma_rho_y[mask == False] = 0.0 + + # set nan values to be zero + sigma_rho_x[~np.isfinite(sigma_rho_x)] = 0.0 + sigma_rho_y[~np.isfinite(sigma_rho_y)] = 0.0 + + # set zero values to a small non-zero number + sigma_rho_x[sigma_rho_x == 0] = 1e-8 + sigma_rho_y[sigma_rho_y == 0] = 1e-8 + + return X, Y, rho_x, rho_y, sigma_rho_x, sigma_rho_y + + +def create_border_array(nr, nc, fill_val=1): + # create array with fill val + border_array = np.ones((nr, nc)) * fill_val + border_array[1:-1, 1:-1] = 0 + + return border_array + + +def create_mask(X, Eval): + # -------------------------- + # set mask + # -------------------------- + # create binary mask. True for flow, False for blocked region. + if len(Eval.shape) == 1: + mask = np.reshape(a=Eval, newshape=(X.shape[1], X.shape[0])) + mask = np.transpose(mask) + else: + mask = Eval + + return mask + + +def load_displacements(filename, displacement_uncertainty_method): + # load the data: + data = loadmat_functions.loadmat(filename=filename) # , squeeze_me=True) + + # extract co-ordinates + X_pix = data['X'].astype('float') + Y_pix = data['Y'].astype('float') + + # extract displacements (sign conversion for the ray tracing code) + U = -data['U'] + V = -data['V'] + + # extract displacement uncertainties + sigma_U = data['uncertainty2D'][displacement_uncertainty_method + 'x'] + sigma_V = data['uncertainty2D'][displacement_uncertainty_method + 'y'] + + # bias correction for MC + if displacement_uncertainty_method == 'MC': + sigma_U = np.sqrt(sigma_U ** 2 - data['uncertainty2D']['biasx'] ** 2) + sigma_V = np.sqrt(sigma_V ** 2 - data['uncertainty2D']['biasy'] ** 2) + + # extract mask + Eval = data['Eval'] + 1 + + return X_pix, Y_pix, U, V, sigma_U, sigma_V, Eval + + +def main(): + + # file containing displacements and uncertainties + filename = 'sample-displacements.mat' + + # set integration method ('p' for poisson or 'w' for wls) + density_integration_method = 'p' + + # displacement uncertainty method + displacement_uncertainty_method = 'MC' + + # ------------------------------------------------- + # experimental parameters for density integration + # ------------------------------------------------- + experimental_parameters = dict() + + # ambient density (kg/m^3) + experimental_parameters['rho_0'] = 1.225 + + # uncertainty in the free-stream density (kg/m^3) + experimental_parameters['sigma_rho_0'] = 0.0 #0.05 * peak_density_offset #5% of peak difference + + # gladstone dale constant (m^3/kg) + experimental_parameters['gladstone_dale'] = 0.225e-3 + + # ambient refractive index + experimental_parameters['n_0'] = 1.0 + experimental_parameters['gladstone_dale']*experimental_parameters['rho_0'] + + # thickness of the density gradient field (m) + experimental_parameters['delta_z'] = 0.01 + + # distance between lens and dot target (object / working distance) (m) + experimental_parameters['object_distance'] = 1.0 + + # distance between the mid-point of the density gradient field and the dot pattern (m) + experimental_parameters['Z_D'] = 0.25 + + # distance between the mid-point of the density gradient field and the camera lens (m) + experimental_parameters['Z_A'] = experimental_parameters['object_distance'] - experimental_parameters['Z_D'] + + # distance between the dot pattern and the camera lens (m) + experimental_parameters['Z_B'] = experimental_parameters['object_distance'] + + # pixels corresponding to center of the FOV + experimental_parameters['x0'] = 256 + experimental_parameters['y0'] = 256 + + # size of a pixel on the camera sensor (m) + experimental_parameters['pixel_pitch'] = 10e-6 + + # focal length of camera lens (m) + experimental_parameters['focal_length'] = 105e-3 + + # non-dimensional magnification of the dot pattern (can also set it directly) + experimental_parameters['magnification'] = experimental_parameters['focal_length'] / ( + experimental_parameters['object_distance'] - experimental_parameters['focal_length']) + + # non-dimensional magnification of the midpoint of the density gradient field + experimental_parameters['magnification_grad'] = experimental_parameters['magnification'] \ + * experimental_parameters['Z_B'] / experimental_parameters['Z_A'] + + # -------------------------- + # processing + # -------------------------- + # load displacements from file + X_pix, Y_pix, U, V, sigma_U, sigma_V, Eval = load_displacements(filename, displacement_uncertainty_method) + + # set mask + mask = create_mask(X_pix, Eval) + + # convert displacements to density gradients and co-ordinates to physical units + X, Y, rho_x, rho_y, sigma_rho_x, sigma_rho_y = convert_displacements_to_physical_units(X_pix, Y_pix, U, V, sigma_U, sigma_V, experimental_parameters, mask) + + # define dirichlet boundary points (minimum one point) - here defined to be all boundaries + dirichlet_label, rho_dirichlet, sigma_rho_dirichlet = set_bc(X, experimental_parameters) + + # calculate density and uncertainty + if density_integration_method == 'p': + rho, sigma_rho = Density_integration_Poisson_uncertainty(X, Y, mask,rho_x, rho_y, + dirichlet_label, rho_dirichlet, + uncertainty_quantification=True, + sigma_grad_x=sigma_rho_x, sigma_grad_y=sigma_rho_y, + sigma_dirichlet=sigma_rho_dirichlet) + elif density_integration_method == 'w': + rho, sigma_rho = Density_integration_WLS_uncertainty(X, Y, mask,rho_x, rho_y, + dirichlet_label, rho_dirichlet, + uncertainty_quantification=True, + sigma_grad_x=sigma_rho_x, sigma_grad_y=sigma_rho_y, + sigma_dirichlet=sigma_rho_dirichlet) + + # # save the results to file + # savemat(file_name='sample-result.mat', mdict={'X': X, 'Y': Y, 'rho': rho, 'sigma_rho': sigma_rho, + # 'dirichlet_label': dirichlet_label, 'rho_dirichlet':rho_dirichlet, 'sigma_rho_dirichlet':sigma_rho_dirichlet + # }, long_field_names=True) + + # Plot the results + fig1 = plt.figure(1, figsize=(12,8)) + plt.figure(1) + ax1 = fig1.add_subplot(2,2,1) + ax2 = fig1.add_subplot(2,2,2) + ax3 = fig1.add_subplot(2,2,3) + ax4 = fig1.add_subplot(2,2,4) + + # plot x gradient + plt.axes(ax1) + plt.pcolormesh(X, Y, rho_x) + plt.colorbar() + ax1.set_aspect('equal') + plt.title('rho_x') + + # plot y gradient + plt.axes(ax2) + plt.pcolormesh(X, Y, rho_y) + plt.colorbar() + ax2.set_aspect('equal') + plt.title('rho_y') + + # plot density + plt.axes(ax3) + plt.pcolormesh(X, Y, rho) + plt.colorbar() + ax3.set_aspect('equal') + plt.title('rho') + + # plot density uncertainty7 + plt.axes(ax4) + plt.pcolormesh(X, Y, sigma_rho) + plt.colorbar() + ax4.set_aspect('equal') + plt.title('sigma rho') + + # plt.tight_layout() + + # save plot to file + plt.savefig('sample-result.png') + plt.close() + +if __name__ == '__main__': + main() + + + + + + + + + + + + + + + +