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

from DensityIntegrationUncertaintyQuantification import Density_integration_Poisson_uncertainty
from DensityIntegrationUncertaintyQuantification import Density_integration_WLS_uncertainty
from DensityIntegrationUncertaintyQuantification import Density_integration_WLS_uncertainty_weighted_average

# 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 (Poisson)
rho_poisson, sigma_rho_poisson = 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)

# calculate density and uncertainty (WLS)
rho_wls, sigma_rho_wls = Density_integration_WLS_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_poisson': rho_poisson, 'sigma_rho_poisson': sigma_rho_poisson,
                                    'rho_wls': rho_wls, 'sigma_rho_wls': sigma_rho_wls}, long_field_names=True)

# Plot the results
fig1 = plt.figure(1, figsize=(12,8))
plt.figure(1)
ax1 = fig1.add_subplot(3,2,1)
ax2 = fig1.add_subplot(3,2,2)
ax3 = fig1.add_subplot(3,2,3)
ax4 = fig1.add_subplot(3,2,4)
ax5 = fig1.add_subplot(3,2,5)
ax6 = fig1.add_subplot(3,2,6)

# plot x gradient
plt.axes(ax1)
plt.pcolor(data['X'], data['Y'], data['rho_x'], vmin=-60, vmax=60)
plt.colorbar()
plt.title('rho_x')

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

# plot density (poisson)
plt.axes(ax3)
plt.pcolor(data['X'], data['Y'], rho_poisson, vmin=1.2, vmax=1.5)
plt.colorbar()
plt.title('rho, Poisson')

# plot density uncertainty (poisson)
plt.axes(ax4)
plt.pcolor(data['X'], data['Y'], sigma_rho_poisson, vmin=0.0, vmax=0.01)
plt.colorbar()
plt.title('sigma rho, Poisson')

# plot density (wls)
plt.axes(ax5)
plt.pcolor(data['X'], data['Y'], rho_wls, vmin=1.2, vmax=1.5)
plt.colorbar()
plt.title('rho, WLS')

# plot density uncertainty (wls)
plt.axes(ax6)
plt.pcolor(data['X'], data['Y'], sigma_rho_wls, vmin=0.0, vmax=0.01)
plt.colorbar()
plt.title('sigma rho, WLS')

plt.tight_layout()

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