%matplotlib inline
from montepetro.models import Model
from montepetro.regions import Region
from montepetro.properties import RandomProperty
from montepetro.seed_generators import SeedGenerator

import numpy

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

region_a = Region(name="Inner Bone")
region_b = Region(name="Outer Bone")

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

#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
porosity = RandomProperty(name="Porosity", n=n, \
                          random_number_function=numpy.random.uniform)
density_calcite = RandomProperty(name="Rho_Calc", n=n, \
                                 random_number_function=numpy.random.triangular)
density_bone_fluid = RandomProperty(name="Rho_Fluid", n=n, \
                                    random_number_function=numpy.random.triangular)

config = {"Inner Bone": {"Porosity": {"low": 0.0, "high": 0.2},
                         "Rho_Calc": {"left": 2.54, "right": 2.7, "mode": 2.58},
                         "Rho_Fluid": {"left": 1.01, "right": 1.5, "mode": 1.2}},
          "Outer Bone": {"Porosity": {"low": 0.1, "high": 0.2},
                       "Rho_Calc": {"left": 2.50, "right": 3.0, "mode": 2.8},
                       "Rho_Fluid": {"left": 1.1, "right": 1.3, "mode": 1.2}}}

#Let's add these to the model.
model.add_property(porosity)
model.add_property(density_calcite)
model.add_property(density_bone_fluid)

#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)

import matplotlib.pyplot as plot
densities = []
for region_name, region in model.regions.iteritems():
    porosity = region.properties["Porosity"].values
    rho_calc = region.properties["Rho_Calc"].values
    rho_bone_fluid = region.properties["Rho_Fluid"].values
    ensemble_density = rho_calc*(1-porosity)+rho_bone_fluid*porosity
    densities.append(ensemble_density)
    #And we can plot these values of course for our example.
    #plot.hist(ensemble_density, bins=50)

total_density = numpy.add(densities[0], densities[1])
plot.hist(total_density, bins=500)
plot.xlabel(u"$Bone \ Density$", fontsize=14)
plot.ylabel(u"$Frequency$",fontsize=14)
plot.show()