%matplotlib inline
import matplotlib.pyplot as plot

import numpy as np
from montepetro.models import Model
from montepetro.properties import RandomProperty, OriginalOilInPlace, ModelOriginalOilInPlace
from montepetro.regions import Region

seed = 300
model = Model("A simple reservoir model", seed)

region_a = Region(name="Region A-Sandstone")
region_b = Region(name="Region B-Carbonate")
region_c = Region(name="Region C-Dolomite")

#For simplicity let's add these to our model
model.add_region(region_a)
model.add_region(region_b)
model.add_region(region_c)

from montepetro.utils import truncated_normal_rvs

#Define number of samples we want to take
n = 100000

#We pass the RandomProperty our functions from which we sample the probability distributions
#Here we make use of the numpy libraries built in random number generators
area = RandomProperty(name="Area", n=n, random_number_function=np.random.uniform)
porosity = RandomProperty(name="Porosity", n=n, random_number_function=truncated_normal_rvs)
sw = RandomProperty(name="Sw", n=n, random_number_function=np.random.triangular)

#Lets add a fictional set of properties.
config = {"Region A-Sandstone": {"Area": {"low": 100.0, "high": 1000.0},
                       "Porosity": {"low": 0.0, "high": 0.4812, "mean": 0.0929, "std":0.0108},
                       "Sw": {"left": 0.00, "right": 0.15, "mode": 0.05}},
          "Region B-Carbonate": {"Area": {"low": 20.0, "high": 200.0},
                       "Porosity": {"low": 0.0, "high": 0.4812, "mean": 0.1, "std":0.1},
                       "Sw": {"left": 0.00, "right": 0.20, "mode": 0.10}},
          "Region C-Dolomite": {"Area": {"low": 200.0, "high": 1000.0},
                       "Porosity": {"low": 0.0, "high": 0.4812, "mean": 0.08, "std":0.02},
                       "Sw": {"left": 0.00, "right": 0.15, "mode": 0.05}}}

#Let's add these to the model.
model.add_property(area)
model.add_property(porosity)
model.add_property(sw)

#Some Model container magic! We add all these properties to the regions.
model.add_defined_properties_to_regions()

#We pass the model our configuration and run the model
#This will generate all the sampled distributions for each region
model.run(config)

model.add_regional_property("ooip", OriginalOilInPlace)

ooip = ModelOriginalOilInPlace(model)
ooip.generate_values()
ooip.calculate_property_statistics()
plot.figure(figsize=(8,8))
plot.hist(ooip.values, bins=100, cumulative=False, normed=True)
plot.xlabel(u"$Original \ Oil \ in \ Place$", fontsize=18)
plot.ylabel(u"$Frequency$", fontsize=18)
plot.title(u"$Model \ Original \ Oil \ in \ Place$", fontsize=18)
plot.show()
print "P90: ", ooip.p10, ", P50: ", ooip.p50, ", P10: ", ooip.p90

areas = {}
porosities = {}
sws = {}
ooips = {}
for region_name, region in model.regions.iteritems():
    areas[region_name] = region.properties["Area"].values
    porosities[region_name] = region.properties["Porosity"].values
    sws[region_name] = region.properties["Sw"].values
    ooips[region_name] = region.properties["ooip"].values

import matplotlib.pyplot as plot

def plot_regional_distributions(prop_name, prop_dic, model, bins, color="white"):
    region_names = model.regions.keys()
    f, axarr = plot.subplots(len(region_names), 1,  figsize=(6,6), sharex=True, sharey=True)
    for region_name, ax in zip(region_names, axarr):
        ax.set_title(u'$'+prop_name+' \ for \ region \ '+region_name+'$')
        ax.hist(prop_dic[region_name], bins=bins)
    plot.setp([a.get_xticklabels() for a in f.axes[:-1]], visible=False)

plot_regional_distributions("Area", areas, model, 50)

plot.show()

plot_regional_distributions("Porosity", porosities, model, 50)
    
plot.show()

plot_regional_distributions("Water Saturation", sws, model, 50)
    
plot.show()

plot_regional_distributions("Oil In Place", ooips, model, 50)
    
plot.show()