sample_script.py

#!/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

# sys.path.insert(0, '/scratch/shannon/c/aether/Projects/Flow_Reconstruction_development/Codes/Python/')
from DensityIntegrationUncertaintyQuantification import Density_integration_Poisson_uncertainty

# Load the data: 
# all loaded variables is in standard physical units: X, Y: [m], rho_x, rho_y: [kg/m^4] 
data = loadmat('sample-data.mat', squeeze_me=True)

# set the density uncertainty at the boundary points [kg/m^3]
sigma_rho_dirichlet = 0.01

# calculate density and uncertainty
rho, sigma_rho = Density_integration_Poisson_uncertainty(data['X'], data['Y'], data['mask'],
                                                         data['rho_x'], data['rho_y'],
                                                         data['dirichlet_label'], data['rho_dirichlet'],
                                                         uncertainty_quantification=True,
                                                         sigma_grad_x=data['sigma_rho_x'], sigma_grad_y=data['sigma_rho_y'],
                                                         sigma_dirichlet=sigma_rho_dirichlet)

# save the results to file
savemat(file_name='sample-result.mat', mdict={'X': data['X'], 'Y': data['Y'], 'rho': rho, 'sigma_rho': sigma_rho}, 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.pcolor(data['X'], data['Y'], data['rho_x'])
plt.colorbar()
plt.title('rho_x')

# plot y gradient
plt.axes(ax2)
plt.pcolor(data['X'], data['Y'], data['rho_y'])
plt.colorbar()
plt.title('rho_y')

# plot density
plt.axes(ax3)
plt.pcolor(data['X'], data['Y'], rho)
plt.colorbar()
plt.title('rho')

# plot density uncertainty7
plt.axes(ax4)
plt.pcolor(data['X'], data['Y'], sigma_rho)
plt.colorbar()
plt.title('sigma rho')

plt.tight_layout()

# save plot to file
plt.savefig('sample-result.png')
plt.close()