Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

added a realization explorer notebook, patch for arm mac issue with LD flags #40

Open
wants to merge 6 commits into
base: develop
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
304 changes: 304 additions & 0 deletions examples/understanding_variograms_and_realizations.ipynb
Original file line number Diff line number Diff line change
@@ -0,0 +1,304 @@
{
"cells": [
{
"cell_type": "markdown",
"id": "90d365d6-ad4e-4109-954c-4b1cb46f733d",
"metadata": {},
"source": [
"# A very simple demo to explain how the variogram effects stochastic realizations\n",
"\n",
"In this simple notebook, we exploit the speed of `PyPestUtils` to help build understanding of how the variogram parameters effect the resulting stochastic field"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "78a88aba-80be-48d1-b2eb-48bd6d67551f",
"metadata": {},
"outputs": [],
"source": [
"import os\n",
"import sys\n",
"import numpy as np\n",
"import pandas as pd\n",
"import matplotlib.pyplot as plt\n",
"from ipywidgets import interact\n",
"sys.path.append(\"..\")\n",
"from pypestutils.pestutilslib import PestUtilsLib\n",
"lib = PestUtilsLib(logger_level=0)"
]
},
{
"cell_type": "markdown",
"id": "da9d0eee-ebfe-44f8-b230-e2df7ac8fc1e",
"metadata": {},
"source": [
"Just some definitions. We will use a 100 X 100 grid of nodes with delx/dely (and an area) of 1."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "0a7ec106-37ef-42f2-a9d9-5d0cea29080b",
"metadata": {},
"outputs": [],
"source": [
"nrow = ncol = 100\n",
"x = np.cumsum(np.ones((nrow,ncol)),axis=1)\n",
"y = np.cumsum(np.ones((nrow,ncol)),axis=0)\n",
"area = 1\n",
"active = 1\n",
"mean = 1 #log10 HK mean\n",
"var = 0.1 #log10 HK variance\n",
"ec = x.flatten()\n",
"nc = y.flatten()\n",
"transtype = 0 #not log transformed\n",
"vtype = 1 #exponential variogram\n",
"power = 0\n",
"nreal = 1"
]
},
{
"cell_type": "markdown",
"id": "d2dca882-8b41-4827-b27f-87c5138a0608",
"metadata": {},
"source": [
"This is the function that will generate one realization and plot it for us, given a correlation length (`aa`), an anisotropy ratio (`anis`), and a `bearing`) "
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "615acaee-212d-4a08-b3ef-1c690b795307",
"metadata": {},
"outputs": [],
"source": [
"\n",
"def plot_real(corrlen,anisotropy,bearing):\n",
" fig, axes = plt.subplots(1,1,figsize=(7,6))\n",
" #axes.clear()\n",
" #reset the random seed so that the underlying pseudo random numbers dont change\n",
" lib.initialize_randgen(1123455)\n",
" # generate one realization\n",
" arrs = lib.fieldgen2d_sva(ec,nc,area,int(active),mean,var,corrlen,anisotropy,bearing,transtype,vtype,power,nreal)\n",
" # plot\n",
" \n",
" cb = axes.imshow(arrs[:,0].reshape((nrow,ncol)),vmin=0.0,vmax=2,cmap=\"magma\")\n",
" #$plt.colorbar(cb,ax=axes,label=\"$log_{10}$ HK\")\n",
" axes.set_xlabel(\"column\")\n",
" axes.set_ylabel(\"row\")\n",
" axes.set_title(\"a random realization of HK\",loc=\"left\")\n",
" plt.show()\n",
" \n",
" "
]
},
{
"cell_type": "markdown",
"id": "cb1483d1-e062-4567-bd85-715e20f382f9",
"metadata": {},
"source": [
"Move these slowly so it doesnt flicker...\n",
"\n",
" - corrlen = the correlation length of the variogram\n",
" - anisotropy = the anisotropy ratio of the primary to second axes of the anisotropy ellipse\n",
" - bearing = the angle from north of the primary axis of the anisotropy ellipse"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ffdee08f-9202-42e0-b503-fc135c099337",
"metadata": {},
"outputs": [],
"source": [
"_ = interact(plot_real,corrlen=(1,50,1),anisotropy=(1,10,0.5),bearing=(0,90,1))"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "dfb9eade-f0b0-43ac-9213-6978f95d0e72",
"metadata": {},
"outputs": [],
"source": [
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "eba1e58f-48fa-4b48-b268-f3da3d484cf8",
"metadata": {},
"source": [
"# These go to 11...\n",
"\n",
"Now let's explore what happens when we introduce hyper parameteres on the property variogram. In essense, we are now going to treat the bearing, anisotropy, and correlation length of the above variogram as themselves being described by variograms - wat?! Its like geostatistical inception..."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "066d8f57-6c0c-4d1f-8112-baff57fc7c9c",
"metadata": {},
"outputs": [],
"source": [
"nrow = 80\n",
"ncol = 50\n",
"x = np.cumsum(np.ones((nrow,ncol)),axis=1)\n",
"y = np.cumsum(np.ones((nrow,ncol)),axis=0)\n",
"area = 1\n",
"active = 1\n",
"mean = 1 #log10 HK mean\n",
"var = 0.1 #log10 HK variance\n",
"ec = x.flatten()\n",
"nc = y.flatten()\n",
"transtype = 0 #not log transformed\n",
"vtype = 1 #exponential variogram\n",
"power = 0\n",
"nreal = 1"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "4efe47ec-95d4-4644-893b-f240c6385e9a",
"metadata": {},
"outputs": [],
"source": [
"def plot_real(b_mean, b_var,b_corrlen,b_anisotropy,b_bearing,\n",
" a_mean, a_var, a_corrlen,a_anisotropy,a_bearing,\n",
" c_mean, c_var, c_corrlen,c_anisotropy,c_bearing ):\n",
" fig, axes = plt.subplots(1,4,figsize=(10,3))\n",
" #axes.clear()\n",
" #reset the random seed so that the underlying pseudo random numbers dont change\n",
" lib.initialize_randgen(1123455)\n",
" # generate one realization\n",
" barr = lib.fieldgen2d_sva(ec,nc,area,int(active),b_mean,b_var,b_corrlen,b_anisotropy,b_bearing,transtype,vtype,power,nreal)\n",
" \n",
" cb = axes[0].imshow(barr[:,0].reshape((nrow,ncol)),cmap=\"magma\")\n",
" plt.colorbar(cb,ax=axes[0])\n",
" axes[0].set_title(\"bearing\")\n",
" \n",
" aarr = lib.fieldgen2d_sva(ec,nc,area,int(active),a_mean,a_var,a_corrlen,a_anisotropy,a_bearing,transtype,vtype,power,nreal)\n",
" \n",
" cb = axes[1].imshow(aarr[:,0].reshape((nrow,ncol)),cmap=\"magma\")\n",
" axes[1].set_title(\"aniso\")\n",
" plt.colorbar(cb,ax=axes[1])\n",
" \n",
" carr = lib.fieldgen2d_sva(ec,nc,area,int(active),c_mean,c_var,c_corrlen,c_anisotropy,c_bearing,transtype,vtype,power,nreal)\n",
" \n",
" cb = axes[2].imshow(carr[:,0].reshape((nrow,ncol)),cmap=\"magma\")\n",
" axes[2].set_title(\"corrlen\")\n",
" plt.colorbar(cb,ax=axes[2])\n",
" \n",
" arrs = lib.fieldgen2d_sva(ec,nc,area,int(active),mean,var,carr[:,0],aarr[:,0],barr[:,0],transtype,vtype,power,nreal)\n",
" # plot\n",
" \n",
" cb = axes[3].imshow(arrs[:,0].reshape((nrow,ncol)),vmin=0.0,vmax=2,cmap=\"magma\")\n",
" axes[3].set_title(\"HK\")\n",
" plt.colorbar(cb,ax=axes[3])\n",
" #$plt.colorbar(cb,ax=axes,label=\"$log_{10}$ HK\")\n",
" for ax in axes:\n",
" #ax.set_xlabel(\"column\")\n",
" #ax.set_ylabel(\"row\")\n",
" #ax.set_title(\"a random realization of HK\",loc=\"left\")\n",
" ax.set_xticks([])\n",
" ax.set_yticks([])\n",
" \n",
" plt.show()\n",
" \n",
" "
]
},
{
"cell_type": "markdown",
"id": "77c2cc90-803f-4c6d-b857-045ce0788676",
"metadata": {},
"source": [
"The 'b_' prefix is the for bearing variogram, 'a_' is the anisotropy variogram, and, you guessed it, 'c_' is for the correlation length variogram...This is complex stuff, but we found that these sliders really helped us start to build some intution regarding how these hyper parameters interact...."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2c844dc6-c255-4361-addd-c8ef8e0796de",
"metadata": {},
"outputs": [],
"source": [
"_ = interact(plot_real,b_mean=(0,90,1), b_var=(0.1,50,1),b_corrlen=(1,100,10),b_anisotropy=(1,10,0.5),b_bearing=(0,90,1),\n",
" a_mean=(1,10,1), a_var=(0.1,10.0,0.1),a_corrlen=(10,500,10),a_anisotropy=(1,10,0.5),a_bearing=(0,90,1),\n",
" c_mean=(1,50,1), c_var=(0.1,50.0,0.1),c_corrlen=(1,1000,10),c_anisotropy=(1,50,1),c_bearing=(0,90,1),)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "9d33a4ee-4706-402c-a65b-5d833da48ebb",
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"id": "d85e0325-afb4-4aa1-a08c-5643017e6513",
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"language_info": {
"name": "python"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
4 changes: 3 additions & 1 deletion scripts/build_lib.sh
Original file line number Diff line number Diff line change
Expand Up @@ -4,9 +4,11 @@ set -e
# always run from top of repo
cd $(dirname $0)/..



# this needs bash
case "$OSTYPE" in
darwin*) libname=lib/libpestutils.dylib ;;
darwin*) libname=lib/libpestutils.dylib && export LDFLAGS="$LDFLAGS -Wl,-ld_classic";;
linux*) libname=lib/libpestutils.so ;;
msys* ) libname=bin/pestutils.dll ;;
*) echo "unknown \$OSTYPE: $OSTYPE" && exit 1 ;;
Expand Down
Loading