diff --git a/Modules/Biophotonics/python/iMC/mc/batches.py b/Modules/Biophotonics/python/iMC/mc/batches.py
index f2bc8040d1..76e0b2afb9 100644
--- a/Modules/Biophotonics/python/iMC/mc/batches.py
+++ b/Modules/Biophotonics/python/iMC/mc/batches.py
@@ -1,321 +1,331 @@
'''
Created on Oct 15, 2015
@author: wirkert
'''
import numpy as np
from pandas import DataFrame
import pandas as pd
+
class AbstractBatch(object):
"""summarizes a batch of simulated mc spectra"""
def __init__(self):
self._nr_layers = 0 # internally keeps track of number of layers
my_index = pd.MultiIndex(levels=[[], []],
labels=[[], []])
self.df = DataFrame(columns=my_index)
def create_parameters(self, nr_samples):
"""create the parameters for the batch, the simulation has
to create the resulting reflectances"""
pass
def nr_elements(self):
return self.df.shape[0]
class GenericBatch(AbstractBatch):
"""generic n-layer batch with each layer having the same oxygenation """
def __init__(self):
super(GenericBatch, self).__init__()
def append_one_layer(self, saO2, nr_samples):
"""helper function to create parameters for one layer"""
# scales data to lie between maxi and mini instead of 0 and 1
scale = lambda x, mini, maxi: x * (maxi - mini) + mini
# shortcut to random generator
gen = np.random.random_sample
gen_n = np.random.normal
# create layer elements
self.df["layer" + str(self._nr_layers), "vhb"] = \
scale(gen(nr_samples), 0, 1.)
self.df["layer" + str(self._nr_layers), "sao2"] = \
saO2
self.df["layer" + str(self._nr_layers), "a_mie"] = \
np.clip(gen_n(loc=18.9, scale=10.2, size=nr_samples),
0.1, np.inf) * 100 # to 1/m
self.df["layer" + str(self._nr_layers), "b_mie"] = \
np.clip(gen_n(loc=1.286, scale=0.521, size=nr_samples), 0, np.inf)
self.df["layer" + str(self._nr_layers), "d"] = \
scale(gen(nr_samples), 0, 1.)
self.df["layer" + str(self._nr_layers), "n"] = \
scale(gen(nr_samples), 1.33, 1.54)
self.df["layer" + str(self._nr_layers), "g"] = \
scale(gen(nr_samples), 0.8, 0.95)
self._nr_layers += 1
def create_parameters(self, nr_samples):
"""Create generic three layer batch with a total diameter of 2mm.
saO2 is the same in all layers, but all other parameters vary randomly
within each layer"""
saO2 = np.random.random_sample(size=nr_samples)
# create three layers with random samples
self.append_one_layer(saO2, nr_samples)
self.append_one_layer(saO2, nr_samples)
self.append_one_layer(saO2, nr_samples)
# "normalize" d to 2mm
# first extract all layers from df
self.df
layers = [l for l in self.df.columns.levels[0] if "layer" in l]
# summarize all ds
sum_d = 0
for l in layers:
sum_d += self.df[l, "d"]
for l in layers:
self.df[l, "d"] = self.df[l, "d"] / sum_d * 2000. * 10 ** -6
self.df[l, "d"] = np.clip(self.df[l, "d"], 25 * 10 ** -6, np.inf)
+ return self.df
class GenericBatch(AbstractBatch):
"""generic n-layer batch with each layer having the same oxygenation """
def __init__(self):
super(GenericBatch, self).__init__()
def append_one_layer(self, saO2, nr_samples):
"""helper function to create parameters for one layer"""
# scales data to lie between maxi and mini instead of 0 and 1
scale = lambda x, mini, maxi: x * (maxi - mini) + mini
# shortcut to random generator
gen = np.random.random_sample
gen_n = np.random.normal
# create layer elements
self.df["layer" + str(self._nr_layers), "vhb"] = \
scale(gen(nr_samples), 0, 1.)
self.df["layer" + str(self._nr_layers), "sao2"] = \
saO2
self.df["layer" + str(self._nr_layers), "a_mie"] = \
np.clip(gen_n(loc=18.9, scale=10.2, size=nr_samples),
0.1, np.inf) * 100 # to 1/m
self.df["layer" + str(self._nr_layers), "b_mie"] = \
np.clip(gen_n(loc=1.286, scale=0.521, size=nr_samples), 0, np.inf)
self.df["layer" + str(self._nr_layers), "d"] = \
scale(gen(nr_samples), 0, 1.)
self.df["layer" + str(self._nr_layers), "n"] = \
scale(gen(nr_samples), 1.33, 1.54)
self.df["layer" + str(self._nr_layers), "g"] = \
scale(gen(nr_samples), 0.8, 0.95)
self._nr_layers += 1
def create_parameters(self, nr_samples):
"""Create generic three layer batch with a total diameter of 2mm.
saO2 is the same in all layers, but all other parameters vary randomly
within each layer"""
saO2 = np.random.random_sample(size=nr_samples)
# create three layers with random samples
self.append_one_layer(saO2, nr_samples)
self.append_one_layer(saO2, nr_samples)
self.append_one_layer(saO2, nr_samples)
# "normalize" d to 2mm
# first extract all layers from df
self.df
layers = [l for l in self.df.columns.levels[0] if "layer" in l]
# summarize all ds
sum_d = 0
for l in layers:
sum_d += self.df[l, "d"]
for l in layers:
self.df[l, "d"] = self.df[l, "d"] / sum_d * 2000. * 10 ** -6
self.df[l, "d"] = np.clip(self.df[l, "d"], 25 * 10 ** -6, np.inf)
+ return self.df
+
class LessGenericBatch(AbstractBatch):
"""less generic three layer batch. This only varies blood volume fraction
w.r.t. the ColonMuscleBatch. Let's see if DA works in this case."""
def __init__(self):
super(LessGenericBatch, self).__init__()
def append_one_layer(self, saO2, n, d_ranges, nr_samples):
"""helper function to create parameters for one layer"""
# scales data to lie between maxi and mini instead of 0 and 1
scale = lambda x, mini, maxi: x * (maxi - mini) + mini
# shortcut to random generator
gen = np.random.random_sample
# create as generic batch
super(LessGenericBatch, self).append_one_layer(saO2, nr_samples)
self._nr_layers -= 1 # we're not finished
# but some changes in specific layer elements
# more specific layer thicknesses
self.df["layer" + str(self._nr_layers), "d"] = \
scale(gen(nr_samples), d_ranges[0], d_ranges[1])
# more specific n
self.df["layer" + str(self._nr_layers), "n"] = \
n
self._nr_layers += 1
def create_parameters(self, nr_samples):
"""Create generic three layer batch with a total diameter of 2mm.
saO2 is the same in all layers, but all other parameters vary randomly
within each layer"""
saO2 = np.random.random_sample(size=nr_samples)
n = np.ones_like(saO2)
# create three layers with random samples
# muscle
self.append_one_layer(saO2, n * 1.36, (600.*10 ** -6, 1010.*10 ** -6),
nr_samples)
# submucosa
self.append_one_layer(saO2, n * 1.36, (415.*10 ** -6, 847.*10 ** -6),
nr_samples)
# mucosa
self.append_one_layer(saO2, n * 1.38, (395.*10 ** -6, 603.*10 ** -6),
nr_samples)
+ return self.df
+
class ColonMuscleBatch(GenericBatch):
"""three layer batch simulating colonic tissue"""
def __init__(self):
super(ColonMuscleBatch, self).__init__()
def append_one_layer(self, saO2, n, d_ranges, nr_samples):
"""helper function to create parameters for one layer"""
# scales data to lie between maxi and mini instead of 0 and 1
scale = lambda x, mini, maxi: x * (maxi - mini) + mini
# shortcut to random generator
gen = np.random.random_sample
# create as generic batch
super(ColonMuscleBatch, self).append_one_layer(saO2, nr_samples)
self._nr_layers -= 1 # we're not finished
# but some changes in specific layer elements
# less blood
self.df["layer" + str(self._nr_layers), "vhb"] = \
scale(gen(nr_samples), 0, 0.1)
# more specific layer thicknesses
self.df["layer" + str(self._nr_layers), "d"] = \
scale(gen(nr_samples), d_ranges[0], d_ranges[1])
# more specific n
self.df["layer" + str(self._nr_layers), "n"] = \
n
self._nr_layers += 1
def create_parameters(self, nr_samples):
"""Create generic three layer batch with a total diameter of 2mm.
saO2 is the same in all layers, but all other parameters vary randomly
within each layer"""
saO2 = np.random.random_sample(size=nr_samples)
n = np.ones_like(saO2)
# create three layers with random samples
# muscle
self.append_one_layer(saO2, n * 1.36, (600.*10 ** -6, 1010.*10 ** -6),
nr_samples)
# submucosa
self.append_one_layer(saO2, n * 1.36, (415.*10 ** -6, 847.*10 ** -6),
nr_samples)
# mucosa
self.append_one_layer(saO2, n * 1.38, (395.*10 ** -6, 603.*10 ** -6),
nr_samples)
+ return self.df
+
class GenericMeanScatteringBatch(GenericBatch):
"""three layer batch simulating colonic tissue"""
def __init__(self):
super(GenericMeanScatteringBatch, self).__init__()
def append_one_layer(self, saO2, nr_samples):
"""helper function to create parameters for one layer"""
# create as generic batch
super(GenericMeanScatteringBatch, self).append_one_layer(saO2,
nr_samples)
self._nr_layers -= 1 # we're not finished
# restrict exponential scattering to mean value for soft tissue.
self.df["layer" + str(self._nr_layers), "b_mie"] = 1.286
self._nr_layers += 1
class ColonMuscleMeanScatteringBatch(ColonMuscleBatch):
"""three layer batch simulating colonic tissue"""
def __init__(self):
super(ColonMuscleMeanScatteringBatch, self).__init__()
def append_one_layer(self, saO2, n, d_ranges, nr_samples):
"""helper function to create parameters for one layer"""
# create as generic batch
super(ColonMuscleMeanScatteringBatch, self).append_one_layer(saO2,
n,
d_ranges,
nr_samples)
self._nr_layers -= 1 # we're not finished
# restrict exponential scattering to mean value for soft tissue.
self.df["layer" + str(self._nr_layers), "b_mie"] = 1.286
self._nr_layers += 1
class VisualizationBatch(AbstractBatch):
"""batch used for visualization of different spectra. Feel free to adapt
for your visualization purposes."""
def __init__(self):
super(VisualizationBatch, self).__init__()
def append_one_layer(self, vhb, sao2, a_mie, b_mie, d, n, g, nr_samples):
"""helper function to create parameters for one layer"""
# create layer elements
self.df["layer" + str(self._nr_layers), "vhb"] = vhb
self.df["layer" + str(self._nr_layers), "sao2"] = sao2
self.df["layer" + str(self._nr_layers), "a_mie"] = a_mie
self.df["layer" + str(self._nr_layers), "b_mie"] = b_mie
self.df["layer" + str(self._nr_layers), "d"] = d
self.df["layer" + str(self._nr_layers), "n"] = n
self.df["layer" + str(self._nr_layers), "g"] = g
self._nr_layers += 1
def create_parameters(self, nr_samples):
# bvf = np.linspace(0.0, .1, nr_samples)
# saO2 = np.linspace(0., 1., nr_samples)
# d = np.linspace(175, 735, nr_samples) * 10 ** -6
# a_mie = np.linspace(5., 30., nr_samples) * 100
# a_ray = np.linspace(0., 60., nr_samples) * 100
# n = np.linspace(1.33, 1.54, nr_samples)
# g = np.linspace(0, 0.95, nr_samples)
# create three layers with random samples
self.append_one_layer([0.1, 0.02], [0.7, 0.1], 18.9*100., 1.286,
500 * 10 ** -6, 1.38, 0.9,
nr_samples)
self.append_one_layer(0.04, 0.7, 18.9*100., 1.286, 500 * 10 ** -6,
1.36, 0.9,
nr_samples)
self.append_one_layer(0.04, 0.7, 18.9*100., 1.286, 500 * 10 ** -6,
1.36, 0.9,
nr_samples)
+
+ return self.df
diff --git a/Modules/Biophotonics/python/iMC/mc/create_spectrum.py b/Modules/Biophotonics/python/iMC/mc/create_spectrum.py
new file mode 100644
index 0000000000..8d42ea197c
--- /dev/null
+++ b/Modules/Biophotonics/python/iMC/mc/create_spectrum.py
@@ -0,0 +1,41 @@
+
+import logging
+import time
+import os
+
+from mc.sim import get_diffuse_reflectance
+
+
+def create_spectrum(tissue_model, sim_wrapper, wavelengths):
+ """
+ Create a whole spectrum from one instance (dataframe_row) using our
+ tissue model at wavelength wavelength.
+
+ Args:
+ tissue_model: the model which should be used to generate the
+ spectrum
+ sim_wrapper: the simulation which should be used to generate the
+ reflectances
+ wavelengths: the wavelengths which shall be simulated
+
+ Returns: the simulated reflectances
+ """
+ start = time.time()
+ # map the _wavelengths array to reflectance list
+
+ def wavelength_to_reflectance(wavelength):
+ # helper function to determine the reflectance for a given
+ # wavelength using the current model and simulation
+ tissue_model.set_wavelength(wavelength)
+ tissue_model.create_mci_file()
+ sim_wrapper.run_simulation()
+ simulation_path = os.path.split(sim_wrapper.mcml_executable)[0]
+ return get_diffuse_reflectance(os.path.join(simulation_path,
+ tissue_model.
+ get_mco_filename()))
+ reflectances = map(wavelength_to_reflectance, wavelengths)
+ end = time.time()
+ # success!
+ logging.info("successfully ran simulation in " +
+ "{:.2f}".format(end - start) + " seconds")
+ return reflectances
diff --git a/Modules/Biophotonics/python/iMC/mc/dfmanipulations.py b/Modules/Biophotonics/python/iMC/mc/dfmanipulations.py
index b2ac14f675..2d94ee8616 100644
--- a/Modules/Biophotonics/python/iMC/mc/dfmanipulations.py
+++ b/Modules/Biophotonics/python/iMC/mc/dfmanipulations.py
@@ -1,43 +1,46 @@
'''
Created on Oct 19, 2015
@author: wirkert
'''
from scipy.interpolate import interp1d
import pandas as pd
def fold_by_sliding_average(df, window_size):
"""take a batch and apply a sliding average with given window size to
the reflectances.
window_size is elements to the left and to the right.
There will be some boundary effect on the edges."""
# next line does the folding.
df.reflectances = pd.rolling_mean(df.reflectances.T, window_size,
- center=True).T
+ center=True).T
# let's get rid of NaN columns which are created at the boundaries
df.dropna(axis="columns", inplace=True)
+ return df
def switch_reflectances(df, new_wavelengths, new_reflectances):
df.drop(df["reflectances"].columns, axis=1, level=1, inplace=True)
for i, nw in enumerate(new_wavelengths):
df["reflectances", nw] = new_reflectances[:, i]
+ return df
def interpolate_wavelengths(df, new_wavelengths):
""" interpolate image data to fit new_wavelengths. Current implementation
performs simple linear interpolation. Neither existing nor new _wavelengths
need to be sorted. """
# build an interpolator using the inormation provided by the dataframes
# reflectance column
interpolator = interp1d(df.reflectances.columns.astype(float),
df.reflectances.as_matrix(), assume_sorted=False,
bounds_error=False)
# use this to create new reflectances
new_reflectances = interpolator(new_wavelengths)
# build a new dataframe out of this information and set the original df
# to the new information. This seems hacky, can't it be done easier?
switch_reflectances(df, new_wavelengths, new_reflectances)
+ return df
diff --git a/Modules/Biophotonics/python/iMC/mc/sim.py b/Modules/Biophotonics/python/iMC/mc/sim.py
index 5cbf268c57..f629a25d99 100644
--- a/Modules/Biophotonics/python/iMC/mc/sim.py
+++ b/Modules/Biophotonics/python/iMC/mc/sim.py
@@ -1,208 +1,211 @@
'''
Created on Sep 8, 2015
This file contains methods which wrap the mcml simulation so it can be
conveniently called from python. One example for a mcml simulation would be
gpumcml:
https://code.google.com/p/gpumcml/
@author: wirkert
'''
import os
import contextlib
import logging
import subprocess32
""" helper method to change to the correct path and back again """
@contextlib.contextmanager
def cd(newPath):
savedPath = os.getcwd()
os.chdir(newPath)
yield
os.chdir(savedPath)
class MciWrapper(object):
'''
this class provides a wrapper to the mcml monte carlo file.
Its purpose is to create a .mci file which the mcml simulation can use to
create the simulation
'''
def set_mci_filename(self, mci_filename):
self.mci_filename = mci_filename
def set_mco_filename(self, mco_filename):
"""path of the mco file.
This can be either a path relative to the mcml executable
or an absolute path.
BUG: it seems that it can only be relative file name
"""
self.mco_filename = mco_filename
def set_nr_photons(self, nr_photons):
self.nr_photons = nr_photons
def add_layer(self, n=None, ua=None, us=None, g=None, d=None):
"""adds a layer below the currently existing ones.
Arguments:
n: Refraction index of medium
ua: absorption coefficient [1/m]
us: scattering coefficient [1/m]
g: anisotropy factor
d: thickness of layer [m]
"""
if n is None:
n = 1.
if ua is None:
ua = 0.
if us is None:
us = 0.
if g is None:
g = 1.
if d is None:
d = 500.*10 ** -6
self.layers.append([n, ua, us, g, d])
def set_layer(self, layer_nr, n, ua, us, g, d):
"""set a layer with a specific layer_nr (stariting with layer_nr 0).
Note that the layer must already exist, otherwise an error will occure
"""
self.layers[layer_nr] = [n, ua, us, g, d]
def set_file_version(self, file_version):
self.file_version = file_version
def set_nr_runs(self, nr_runs):
self.nr_runs = nr_runs
def set_dz_dr(self, dz, dr):
self.dz = dz
self.dr = dr
def set_nr_dz_dr_da(self, nr_dz, nr_dr, nr_da):
self.nr_dz = nr_dz
self.nr_dr = nr_dr
self.nr_da = nr_da
def set_n_medium_above(self, n_above):
self.n_above = n_above
def set_n_medium_below(self, n_below):
self.n_below = n_below
def create_mci_file(self):
"""this method creates the mci file at the location self.mci_filename"""
open(self.mci_filename, 'a').close()
f = open(self.mci_filename, 'w')
# write general information
f.write(str(self.file_version) + " # file version\n")
f.write(str(self.nr_runs) + " # number of runs\n\n")
# write the data for run
f.write(self.mco_filename + " A # output filename, ASCII/Binary\n")
f.write(str(self.nr_photons) + " # No. of photons\n")
f.write(repr(self.dz) + " " + repr(self.dr) + " # dz, dr\n")
f.write(repr(self.nr_dz) + " " +
repr(self.nr_dr) + " " +
repr(self.nr_da) + " # No. of dz, dr & da.\n\n")
# write layer information
f.write(str(len(self.layers)) + " # No. of layers\n")
f.write("# n mua mus g d # One line for each layer\n")
f.write(repr(self.n_above) + " # n for medium above.\n")
for layer in self.layers:
# factors (/100.; *100.) to convert to mcml expected units:
f.write("%.3f" % layer[0] + " " + # n
"%.5f" % (layer[1] / 100.) + " " + # ua
"%.5f" % (layer[2] / 100.) + " " + # us
"%.3f" % layer[3] + " " + # g
"%.3f" % (layer[4] * 100.) + "\n") # d
f.write(repr(self.n_below) + " # n for medium below.\n")
f.close()
+ if not os.path.isfile(self.mci_filename):
+ raise IOError("input file for monte carlo simulation not " +
+ "succesfully created")
def __init__(self):
# set standard parameters
self.file_version = 1.0
self.nr_photons = 10**6
self.nr_runs = 1
self.dz = 0.002
self.dr = 2
self.nr_dz = 500
self.nr_dr = 1
self.nr_da = 1
self.n_above = 1.0
self.n_below = 1.0
# initialize to 0 layers
self.layers = []
class SimWrapper(object):
def set_mci_filename(self, mci_filename):
"""the full path to the input file. E.g. ./data/my.mci
"""
self.mci_filename = mci_filename
def set_mcml_executable(self, mcml_executable):
""" the full path of the excutable. E.g. ./mcml/mcml.exe"""
self.mcml_executable = mcml_executable
def run_simulation(self):
"""this method runs a monte carlo simulation"""
mcml_path, mcml_file = os.path.split(self.mcml_executable)
abs_mci_filename = os.path.abspath(self.mci_filename)
# note: the -A option makes gpumcml much faster, but is not available
# in original mcml. Maybe a switch should be introduced here
args = ("./" + mcml_file, "-A", abs_mci_filename)
# switch to folder where mcml resides in and execute it.
with cd(mcml_path):
try:
popen = subprocess32.Popen(args, stdout=subprocess32.PIPE)
popen.wait(timeout=100)
except:
logging.error("couldn't run simulation")
# popen.kill()
def __init__(self):
pass
def get_diffuse_reflectance(mco_filename):
"""
extract reflectance from mco file.
Attention: mco_filename specifies full path.
Returns: the reflectance
"""
with open(mco_filename) as myFile:
for line in myFile:
if "Diffuse reflectance" in line:
return float(line.split(' ', 1)[0])
def get_specular_reflectance(mco_filename):
"""
extract reflectance from mco file.
Attention: mco_filename specifies full path.
Returns: the reflectance
"""
with open(mco_filename) as myFile:
for line in myFile:
if "Specular reflectance" in line:
return float(line.split(' ', 1)[0])
def get_total_reflectance(mco_filename):
"""
extract reflectance from mco file.
Attention: mco_filename specifies full path.
Returns: the reflectance
"""
return get_diffuse_reflectance(mco_filename) + \
get_specular_reflectance(mco_filename)
diff --git a/Modules/Biophotonics/python/iMC/mc/tissuemodels.py b/Modules/Biophotonics/python/iMC/mc/tissuemodels.py
index 4b35c1a42e..7bcdc7806e 100644
--- a/Modules/Biophotonics/python/iMC/mc/tissuemodels.py
+++ b/Modules/Biophotonics/python/iMC/mc/tissuemodels.py
@@ -1,135 +1,141 @@
'''
Created on Sep 9, 2015
@author: wirkert
'''
import numpy as np
from mc.sim import MciWrapper
from mc.usuag import Ua, UsgJacques
class AbstractTissue(object):
'''
Initializes a abstract tissue model"
'''
def set_nr_photons(self, nr_photons):
self._mci_wrapper.set_nr_photons(nr_photons)
def set_mci_filename(self, mci_filename):
self._mci_wrapper.set_mci_filename(mci_filename)
def set_mco_filename(self, mco_filename):
self._mci_wrapper.set_mco_filename(mco_filename)
+ def get_mco_filename(self):
+ return self._mci_wrapper.mco_filename
+
def set_wavelength(self, wavelength):
self.wavelength = wavelength
def create_mci_file(self):
# set layers
for i, ua in enumerate(self.uas):
self._mci_wrapper.set_layer(i, # layer nr
self.ns[i], # refraction index
self.uas[i](self.wavelength), # ua
self.usgs[i](self.wavelength)[0], # us
self.usgs[i](self.wavelength)[1], # g
self.ds[i]) # d
# now that the layers have been updated: create file
self._mci_wrapper.create_mci_file()
def __str__(self):
""" Overwrite this method!
print the current model"""
model_string = ""
return model_string
def __init__(self, ns, uas, usgs, ds):
self._mci_wrapper = MciWrapper()
self.wavelength = 500.*10**9 # standard wavelength, should be set.
self.uas = uas
self.usgs = usgs
self.ds = ds
self.ns = ns
# initially create layers. these will be overwritten as soon
# as create_mci_file is called.
for i in enumerate(uas):
self._mci_wrapper.add_layer()
class GenericTissue(AbstractTissue):
'''
Initializes a 3-layer generic tissue model
'''
- def set_dataframe_element(self, df, element):
- """take the element element of the batch and set the tissue to
- resemble the structure specified by this"""
- layers = [l for l in df.columns.levels[0] if "layer" in l]
+ def set_dataframe_row(self, df_row):
+ """take one example (one row) of a created batch and set the tissue to
+ resemble the structure specified by this row
+
+ Args:
+ df_row: one row of a dataframe created by a batch."""
+ layers = [l for l in df_row.index.levels[0] if "layer" in l]
for i, l in enumerate(layers):
self.set_layer(i,
- df[l, "vhb"][element],
- df[l, "sao2"][element],
- df[l, "a_mie"][element],
- df[l, "b_mie"][element],
- df[l, "d"][element],
- df[l, "n"][element],
- df[l, "g"][element])
+ df_row[l, "vhb"],
+ df_row[l, "sao2"],
+ df_row[l, "a_mie"],
+ df_row[l, "b_mie"],
+ df_row[l, "d"],
+ df_row[l, "n"],
+ df_row[l, "g"])
def set_layer(self, layer_nr=0,
bvf=None, saO2=None, a_mie=None, b_mie=None, d=None,
n=None, g=None):
"""Helper function to set one layer."""
if bvf is None:
bvf = 0.02
if saO2 is None:
saO2 = 0.7
if a_mie is None:
a_mie = 10. * 100
if d is None:
d = 500. * 10 ** -6
if b_mie is None:
b_mie = 1.286
if n is None:
n = 1.38
if g is None:
g = 0.
# build obejct for absorption coefficient determination
self.uas[layer_nr].bvf = bvf
self.uas[layer_nr].saO2 = saO2
# and one for scattering coefficient
self.usgs[layer_nr].a_mie = a_mie
self.usgs[layer_nr].a_ray = 0.
self.usgs[layer_nr].b_mie = b_mie
self.usgs[layer_nr].g = g
self.ds[layer_nr] = d
self.ns[layer_nr] = n
def __str__(self):
"""print the current model"""
model_string = ""
for i, ua in enumerate(self.uas):
layer_string = "layer " + str(i) + \
" - vhb: " + "%.1f" % (self.uas[i].bvf * 100.) + \
"%; sao2: " + "%.1f" % (self.uas[i].saO2 * 100.) + \
"%; a_mie: " + "%.2f" % (self.usgs[i].a_mie / 100.) + \
"cm^-1; a_ray: " + "%.2f" % (self.usgs[i].a_ray / 100.) + \
"cm^-1; b_mie: " + "%.3f" % self.usgs[i].b_mie + \
"; d: " + "%.0f" % (self.ds[i] * 10 ** 6) + "um" + \
"; n: " + "%.2f" % (self.ns[i]) + \
"; g: " + "%.2f" % self.usgs[i].g + "\n"
model_string += layer_string
return model_string
def __init__(self, nr_layers=3):
uas = []
usgs = []
for i in range(nr_layers):
uas.append(Ua())
usgs.append(UsgJacques())
ds = np.ones(nr_layers, dtype=float) * 500.*10 ** -6
ns = np.ones(nr_layers, dtype=float) * 1.38
super(GenericTissue, self).__init__(ns, uas, usgs, ds)
diff --git a/Modules/Biophotonics/python/iMC/regression/demo.py b/Modules/Biophotonics/python/iMC/regression/demo.py
deleted file mode 100644
index c4717ca45e..0000000000
--- a/Modules/Biophotonics/python/iMC/regression/demo.py
+++ /dev/null
@@ -1,107 +0,0 @@
-'''
-Created on Oct 19, 2015
-
-@author: wirkert
-'''
-
-
-from tabulate import tabulate # do sudo pip install tabulate
-import numpy as np
-from sklearn.linear_model import LinearRegression
-from sklearn.ensemble import RandomForestRegressor
-from sklearn.grid_search import GridSearchCV
-from sklearn.cross_validation import KFold
-from sklearn.preprocessing import Normalizer
-from sklearn.metrics import r2_score
-
-from linear import LinearSaO2Unmixing
-from domain_adaptation import estimate_logistic_regressor, resample
-
-def get_data_from_batch(filenameprefix):
- """returns (X,y): X are absorptions represented by an array of shape
- [nr_samples x nr_wavelengths]. y an array of shape [nr_samples]
- representing the oxygenation values """
- reflectances = np.load(filenameprefix + "_reflectances.npy")
- saO2 = np.load(filenameprefix + "_saO2.npy")
- # add noise to reflectances
- noises = np.random.normal(loc=0., scale=0.15, size=reflectances.shape)
- reflectances += noises * reflectances
- reflectances = np.clip(reflectances, 0.00001, 1.)
- # normalize reflectances
- normalizer = Normalizer(norm='l1')
- reflectances = normalizer.transform(reflectances)
- # transform to absorption
- absorptions = -np.log(reflectances)
- return absorptions, saO2
-
-
-if __name__ == '__main__':
- # first load training and testing data
- X_train, oxy_train = get_data_from_batch("/media/wirkert/data/Data/temp/" +
- "hard_train_2")
- X_gaussian_train, oxy_gaussian_train = get_data_from_batch("/media/wirkert/data/Data/temp/" +
- "gaussian_train_2")
- X_test, oxy_test = get_data_from_batch("/media/wirkert/data/Data/temp/" +
- "gaussian_test")
-
- # do domain adaptation
- lr = estimate_logistic_regressor(X_train, X_gaussian_train)
- weights = lr.predict_proba(X_train)[:, 1] / lr.predict_proba(X_train)[:, 0]
-
- # do standard linear unmixing
- linear_unmixing = LinearSaO2Unmixing()
- predicted_oxy_lu = linear_unmixing.predict(X_test)
-
- # learn linear regression from data
- X_train, oxy_train, weights = resample(X_train, oxy_train, weights)
- # do weighting of samples
- X_train_w = X_train # * np.sqrt(weights)[:, None]
- oxy_train_w = oxy_train # * np.sqrt(weights)
- linear_regression = LinearRegression(fit_intercept=False)
- linear_regression.fit(X_train, oxy_train_w)
- predicted_oxy_lr = linear_regression.predict(X_test)
-
- # use the random forest regressor
- kf = KFold(X_train.shape[0], 5, shuffle=True)
- param_grid = [
- {"n_estimators": np.array([50]),
- "max_depth": np.arange(5, 15, 2),
- "min_samples_leaf": np.array([5, 7, 10, 15])}]
- # train forest
- rf = GridSearchCV(RandomForestRegressor(50, max_depth=10, n_jobs=-1),
- param_grid, cv=kf, n_jobs=-1)
- rf = RandomForestRegressor(50, max_depth=11, min_samples_leaf=5,
- n_jobs=-1)
- rf.fit(X_train, oxy_train)
- predicted_oxy_rf = rf.predict(X_test)
-
- # nicely print each sample, its true value and its predicted values
- np.set_printoptions(precision=2) # why doesn't this work?
- data = np.vstack((oxy_test,
- predicted_oxy_lu, predicted_oxy_lr, predicted_oxy_rf)).T
- print "detailed results:"
- regressor_list = ["lin unmixing", "lin regression", "rf regression"]
- print tabulate(data, headers=["true value"] + regressor_list)
-
- print "\nr2 scores"
- print "lin unmixing, lin regression, rf regression"
- print r2_score(oxy_test, predicted_oxy_lu), \
- linear_regression.score(X_test, oxy_test), \
- rf.score(X_test, oxy_test)
-
- print "\nabs mean errors"
- f = lambda x: np.mean(np.abs(x - oxy_test))
- print "lin unmixing, lin regression, rf regression"
- print f(predicted_oxy_lu), f(predicted_oxy_lr), f(predicted_oxy_rf)
-
- print "\nabs error lower percentile"
- f = lambda x: np.percentile(np.abs(x - oxy_test), 25, axis=0)
- print "lin unmixing, lin regression, rf regression"
- print f(predicted_oxy_lu), f(predicted_oxy_lr), f(predicted_oxy_rf)
-
- print "\nabs error higher percentile"
- f = lambda x: np.percentile(np.abs(x - oxy_test), 75, axis=0)
- print "lin unmixing, lin regression, rf regression"
- print f(predicted_oxy_lu), f(predicted_oxy_lr), f(predicted_oxy_rf)
-
-
diff --git a/Modules/Biophotonics/python/iMC/scripts/ipcai2016/script_calculate_spectra.py b/Modules/Biophotonics/python/iMC/scripts/ipcai2016/script_calculate_spectra.py
index 34abeb1d76..74cce1661f 100644
--- a/Modules/Biophotonics/python/iMC/scripts/ipcai2016/script_calculate_spectra.py
+++ b/Modules/Biophotonics/python/iMC/scripts/ipcai2016/script_calculate_spectra.py
@@ -1,132 +1,118 @@
'''
Created on Sep 9, 2015
@author: wirkert
'''
import logging
import datetime
import os
import time
import numpy as np
import luigi
import scriptpaths as sp
import mc.factories as mcfac
-from mc.sim import SimWrapper, get_diffuse_reflectance
+from mc.sim import SimWrapper
+from mc.create_spectrum import create_spectrum
# parameter setting
NR_BATCHES = 100
NR_ELEMENTS_IN_BATCH = 1000
# the wavelengths to be simulated
WAVELENGHTS = np.arange(450, 720, 2) * 10 ** -9
NR_PHOTONS = 10 ** 6
# experiment configuration
MCI_FILENAME = "./temp.mci"
MCO_FILENAME = "temp.mco"
# this path definitly needs to be adapted by you
PATH_TO_MCML = "/home/wirkert/workspace/monteCarlo/gpumcml/fast-gpumcml/"
EXEC_MCML = "gpumcml.sm_20"
class CreateSpectraTask(luigi.Task):
df_prefix = luigi.Parameter()
batch_nr = luigi.IntParameter()
nr_samples = luigi.IntParameter()
factory = luigi.Parameter()
def output(self):
return luigi.LocalTarget(os.path.join(sp.MC_DATA_FOLDER,
self.df_prefix + "_" +
str(self.batch_nr) + ".txt"))
def run(self):
start = time.time()
# setup simulation wrapper
- self.sim_wrapper = SimWrapper()
- self.sim_wrapper.set_mci_filename(MCI_FILENAME)
- self.sim_wrapper.set_mcml_executable(PATH_TO_MCML + EXEC_MCML)
+ sim_wrapper = SimWrapper()
+ sim_wrapper.set_mci_filename(MCI_FILENAME)
+ sim_wrapper.set_mcml_executable(os.path.join(PATH_TO_MCML, EXEC_MCML))
# setup model
- self.tissue_model = self.factory.create_tissue_model()
- self.tissue_model.set_mci_filename(self.sim_wrapper.mci_filename)
- self.tissue_model.set_mco_filename(MCO_FILENAME)
- self.tissue_model.set_nr_photons(NR_PHOTONS)
+ tissue_model = self.factory.create_tissue_model()
+ tissue_model.set_mci_filename(sim_wrapper.mci_filename)
+ tissue_model.set_mco_filename(MCO_FILENAME)
+ tissue_model.set_nr_photons(NR_PHOTONS)
# setup array in which data shall be stored
batch = self.factory.create_batch_to_simulate()
batch.create_parameters(self.nr_samples)
# dataframe created by batch:
df = batch.df
# add reflectance column to dataframe
for w in WAVELENGHTS:
df["reflectances", w] = np.NAN
# for each instance of our tissue model
for i in range(df.shape[0]):
- self.tissue_model.set_dataframe_element(df, i)
+ # set the desired element in the dataframe to be simulated
+ tissue_model.set_dataframe_row(df.loc[i, :])
logging.info("running simulation " + str(i) + " for\n" +
- str(self.tissue_model))
- start = time.time()
- # map the _wavelengths array to reflectance list
- reflectances = map(self.wavelength_to_reflectance,
- WAVELENGHTS)
- end = time.time()
+ str(tissue_model))
+ reflectances = create_spectrum(tissue_model, sim_wrapper,
+ WAVELENGHTS)
# store in dataframe
for r, w in zip(reflectances, WAVELENGHTS):
df["reflectances", w][i] = r
- # success!
- logging.info("successfully ran simulation in " +
- "{:.2f}".format(end - start) + " seconds")
+
# clean up temporarily created files
os.remove(MCI_FILENAME)
- created_mco_file = PATH_TO_MCML + "/" + MCO_FILENAME
+ created_mco_file = os.path.join(PATH_TO_MCML, MCO_FILENAME)
if os.path.isfile(created_mco_file):
os.remove(created_mco_file)
# save the created output
f = open(self.output().path, 'w')
df.to_csv(f)
end = time.time()
logging.info("time for creating batch of mc data: %.f s" %
(end - start))
- def wavelength_to_reflectance(self, wavelength):
- """helper function to determine the reflectance for a given
- wavelength using the current model and simulation """
- self.tissue_model.set_wavelength(wavelength)
- self.tissue_model.create_mci_file()
- if os.path.isfile(PATH_TO_MCML + EXEC_MCML):
- self.sim_wrapper.run_simulation()
- return get_diffuse_reflectance(PATH_TO_MCML + MCO_FILENAME)
- else:
- raise IOError("path to gpumcml not valid")
-
if __name__ == '__main__':
logging.basicConfig(filename=os.path.join(sp.LOG_FOLDER,
"calculate_spectra" +
str(datetime.datetime.now()) +
'.log'),
level=logging.INFO)
ch = logging.StreamHandler()
ch.setLevel(logging.INFO)
logger = logging.getLogger()
logger.addHandler(ch)
luigi.interface.setup_interface_logging()
# create a folder for the results if necessary
sp.create_folder_if_necessary(sp.INTERMEDIATES_FOLDER)
sch = luigi.scheduler.CentralPlannerScheduler()
w = luigi.worker.Worker(scheduler=sch)
BATCH_NUMBERS = np.arange(0, NR_BATCHES, 1)
for i in BATCH_NUMBERS:
colon_task = CreateSpectraTask("ipcai_revision_generic_mean_scattering",
i,
NR_ELEMENTS_IN_BATCH,
mcfac.GenericMeanScatteringFactory())
w.add(colon_task)
w.run()
diff --git a/Modules/Biophotonics/python/iMC/tutorials/Monte Carlo Spectra Generation - Basic tutorial.ipynb b/Modules/Biophotonics/python/iMC/tutorials/Monte Carlo Spectra Generation - Basic tutorial.ipynb
new file mode 100644
index 0000000000..ad3851bc29
--- /dev/null
+++ b/Modules/Biophotonics/python/iMC/tutorials/Monte Carlo Spectra Generation - Basic tutorial.ipynb
@@ -0,0 +1,1018 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Creating and manipulating monte carlo spectra using MITK-MSI\n",
+ "\n",
+ "In this tutorial we will learn how to\n",
+ "1. create reflectance spectra from examplary tissues\n",
+ "2. how to analyse and visualize the created spectra\n",
+ "3. how to manipulate them\n",
+ "\n",
+ "The MITK-MSI software provides a wrapper to the popular MCML approach to simulate how light travels through tissue. This wrapper can be found in mc/sim.py.\n",
+ "In this tutorial we will utilize our tissue model which uses this wrapper to create the reflectance spectra.\n",
+ "\n",
+ "As a prerequisit, you need a MCML monte carlo simulation which uses the format specified [here](http://omlc.org/software/mc/).\n",
+ "I tested this software with the GPU accelerated version which can be found [here](https://code.google.com/archive/p/gpumcml/)."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "# 1.1 create spectra - setup simulation environment\n",
+ "\n",
+ "# some necessary imports\n",
+ "import logging\n",
+ "import numpy as np\n",
+ "import os\n",
+ "# everything related to the simulation wrapper\n",
+ "from mc import sim\n",
+ "# the factories create batches (tissue samples) and suited tissue models\n",
+ "from mc import factories\n",
+ "# function which runs simulations for each wavelength\n",
+ "from mc.create_spectrum import create_spectrum\n",
+ "\n",
+ "# Where does your monte carlo simulation executable resides in?\n",
+ "MCML_EXECUTABLE = \"/home/wirkert/workspace/monteCarlo/gpumcml/fast-gpumcml/gpumcml.sm_20\"\n",
+ "# The MCML needs a simulation input file, where shall it be created?\n",
+ "MCI_FILENAME = \"./temp.mci\"\n",
+ "# filename of the file with the simulation results. Well, due to a bug in GPUMCML will reside\n",
+ "# in the same folder as the MCML executable\n",
+ "MCO_FILENAME = \"temp.mco\"\n",
+ "# The wavelengths for which we want to run our simulation\n",
+ "WAVELENGTHS = np.arange(450, 720, 2) * 10 ** -9\n",
+ "\n",
+ "# we want to create standard colonic tissue as specified in the IPCAI 2016 publication\n",
+ "# \"Robust Near Real-Time Estimation of Physiological Parameters from Megapixel\n",
+ "# Multispectral Images with Inverse Monte Carlo and Random Forest Regression\"\n",
+ "factory = factories.ColonMuscleMeanScatteringFactory()\n",
+ "\n",
+ "# create a simulation wrapper\n",
+ "sim_wrapper = sim.SimWrapper()\n",
+ "# our simulation needs to know where the input file for the simulation\n",
+ "# shall resign (will be automatically created)\n",
+ "sim_wrapper.set_mci_filename(MCI_FILENAME)\n",
+ "# also it needs to know where the simulation executable shall lie in\n",
+ "sim_wrapper.set_mcml_executable(MCML_EXECUTABLE)\n",
+ "\n",
+ "# create the tissue model\n",
+ "# it is responsible for writing the simulation input file\n",
+ "tissue_model = factory.create_tissue_model()\n",
+ "# tell it where the input file shall lie in\n",
+ "tissue_model.set_mci_filename(sim_wrapper.mci_filename)\n",
+ "# also set the output filename\n",
+ "tissue_model.set_mco_filename(MCO_FILENAME)\n",
+ "# tell it how much photons shall be simulated. Will be set to 10**6 by standard,\n",
+ "# this is just an example\n",
+ "tissue_model.set_nr_photons(10**6)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "text/html": [
+ "\n",
+ "
\n",
+ " \n",
+ " \n",
+ " | \n",
+ " layer0 | \n",
+ " layer1 | \n",
+ " layer2 | \n",
+ "
\n",
+ " \n",
+ " | \n",
+ " vhb | \n",
+ " sao2 | \n",
+ " a_mie | \n",
+ " b_mie | \n",
+ " d | \n",
+ " n | \n",
+ " g | \n",
+ " vhb | \n",
+ " sao2 | \n",
+ " a_mie | \n",
+ " ... | \n",
+ " d | \n",
+ " n | \n",
+ " g | \n",
+ " vhb | \n",
+ " sao2 | \n",
+ " a_mie | \n",
+ " b_mie | \n",
+ " d | \n",
+ " n | \n",
+ " g | \n",
+ "
\n",
+ " \n",
+ " \n",
+ " \n",
+ " | 0 | \n",
+ " 0.060727 | \n",
+ " 0.405108 | \n",
+ " 2732.072167 | \n",
+ " 1.286 | \n",
+ " 0.000643 | \n",
+ " 1.36 | \n",
+ " 0.911375 | \n",
+ " 0.012529 | \n",
+ " 0.405108 | \n",
+ " 2239.618718 | \n",
+ " ... | \n",
+ " 0.000674 | \n",
+ " 1.36 | \n",
+ " 0.871536 | \n",
+ " 0.044635 | \n",
+ " 0.405108 | \n",
+ " 3870.245726 | \n",
+ " 1.286 | \n",
+ " 0.000518 | \n",
+ " 1.38 | \n",
+ " 0.854300 | \n",
+ "
\n",
+ " \n",
+ " | 1 | \n",
+ " 0.008601 | \n",
+ " 0.380418 | \n",
+ " 2715.509993 | \n",
+ " 1.286 | \n",
+ " 0.000959 | \n",
+ " 1.36 | \n",
+ " 0.864808 | \n",
+ " 0.075130 | \n",
+ " 0.380418 | \n",
+ " 1480.405725 | \n",
+ " ... | \n",
+ " 0.000594 | \n",
+ " 1.36 | \n",
+ " 0.904884 | \n",
+ " 0.007109 | \n",
+ " 0.380418 | \n",
+ " 2532.045757 | \n",
+ " 1.286 | \n",
+ " 0.000457 | \n",
+ " 1.38 | \n",
+ " 0.828825 | \n",
+ "
\n",
+ " \n",
+ " | 2 | \n",
+ " 0.092514 | \n",
+ " 0.354746 | \n",
+ " 1522.467374 | \n",
+ " 1.286 | \n",
+ " 0.000626 | \n",
+ " 1.36 | \n",
+ " 0.938954 | \n",
+ " 0.010242 | \n",
+ " 0.354746 | \n",
+ " 1089.628438 | \n",
+ " ... | \n",
+ " 0.000799 | \n",
+ " 1.36 | \n",
+ " 0.921668 | \n",
+ " 0.080251 | \n",
+ " 0.354746 | \n",
+ " 2374.918328 | \n",
+ " 1.286 | \n",
+ " 0.000410 | \n",
+ " 1.38 | \n",
+ " 0.846570 | \n",
+ "
\n",
+ " \n",
+ " | 3 | \n",
+ " 0.014844 | \n",
+ " 0.271838 | \n",
+ " 1817.987790 | \n",
+ " 1.286 | \n",
+ " 0.000947 | \n",
+ " 1.36 | \n",
+ " 0.822284 | \n",
+ " 0.098037 | \n",
+ " 0.271838 | \n",
+ " 1214.562685 | \n",
+ " ... | \n",
+ " 0.000739 | \n",
+ " 1.36 | \n",
+ " 0.816509 | \n",
+ " 0.061391 | \n",
+ " 0.271838 | \n",
+ " 650.190139 | \n",
+ " 1.286 | \n",
+ " 0.000477 | \n",
+ " 1.38 | \n",
+ " 0.858360 | \n",
+ "
\n",
+ " \n",
+ " | 4 | \n",
+ " 0.051730 | \n",
+ " 0.336171 | \n",
+ " 1545.569179 | \n",
+ " 1.286 | \n",
+ " 0.000840 | \n",
+ " 1.36 | \n",
+ " 0.855278 | \n",
+ " 0.054194 | \n",
+ " 0.336171 | \n",
+ " 2574.712829 | \n",
+ " ... | \n",
+ " 0.000518 | \n",
+ " 1.36 | \n",
+ " 0.934114 | \n",
+ " 0.070809 | \n",
+ " 0.336171 | \n",
+ " 3457.117633 | \n",
+ " 1.286 | \n",
+ " 0.000459 | \n",
+ " 1.38 | \n",
+ " 0.924634 | \n",
+ "
\n",
+ " \n",
+ " | 5 | \n",
+ " 0.033220 | \n",
+ " 0.675550 | \n",
+ " 2920.462296 | \n",
+ " 1.286 | \n",
+ " 0.000854 | \n",
+ " 1.36 | \n",
+ " 0.879755 | \n",
+ " 0.052119 | \n",
+ " 0.675550 | \n",
+ " 2121.255538 | \n",
+ " ... | \n",
+ " 0.000716 | \n",
+ " 1.36 | \n",
+ " 0.870148 | \n",
+ " 0.017853 | \n",
+ " 0.675550 | \n",
+ " 2663.702101 | \n",
+ " 1.286 | \n",
+ " 0.000417 | \n",
+ " 1.38 | \n",
+ " 0.901358 | \n",
+ "
\n",
+ " \n",
+ " | 6 | \n",
+ " 0.042772 | \n",
+ " 0.713138 | \n",
+ " 2343.482296 | \n",
+ " 1.286 | \n",
+ " 0.001007 | \n",
+ " 1.36 | \n",
+ " 0.889161 | \n",
+ " 0.052464 | \n",
+ " 0.713138 | \n",
+ " 1494.564335 | \n",
+ " ... | \n",
+ " 0.000795 | \n",
+ " 1.36 | \n",
+ " 0.857640 | \n",
+ " 0.000794 | \n",
+ " 0.713138 | \n",
+ " 1365.344192 | \n",
+ " 1.286 | \n",
+ " 0.000491 | \n",
+ " 1.38 | \n",
+ " 0.902503 | \n",
+ "
\n",
+ " \n",
+ " | 7 | \n",
+ " 0.055083 | \n",
+ " 0.493522 | \n",
+ " 2862.632265 | \n",
+ " 1.286 | \n",
+ " 0.000674 | \n",
+ " 1.36 | \n",
+ " 0.901605 | \n",
+ " 0.011494 | \n",
+ " 0.493522 | \n",
+ " 876.784363 | \n",
+ " ... | \n",
+ " 0.000730 | \n",
+ " 1.36 | \n",
+ " 0.847137 | \n",
+ " 0.040634 | \n",
+ " 0.493522 | \n",
+ " 2460.497781 | \n",
+ " 1.286 | \n",
+ " 0.000494 | \n",
+ " 1.38 | \n",
+ " 0.879252 | \n",
+ "
\n",
+ " \n",
+ " | 8 | \n",
+ " 0.033641 | \n",
+ " 0.121889 | \n",
+ " 10.000000 | \n",
+ " 1.286 | \n",
+ " 0.000716 | \n",
+ " 1.36 | \n",
+ " 0.942503 | \n",
+ " 0.093615 | \n",
+ " 0.121889 | \n",
+ " 2470.880600 | \n",
+ " ... | \n",
+ " 0.000738 | \n",
+ " 1.36 | \n",
+ " 0.824153 | \n",
+ " 0.050557 | \n",
+ " 0.121889 | \n",
+ " 1990.760622 | \n",
+ " 1.286 | \n",
+ " 0.000536 | \n",
+ " 1.38 | \n",
+ " 0.857569 | \n",
+ "
\n",
+ " \n",
+ " | 9 | \n",
+ " 0.090767 | \n",
+ " 0.933581 | \n",
+ " 2049.898620 | \n",
+ " 1.286 | \n",
+ " 0.000914 | \n",
+ " 1.36 | \n",
+ " 0.900553 | \n",
+ " 0.092938 | \n",
+ " 0.933581 | \n",
+ " 2080.923970 | \n",
+ " ... | \n",
+ " 0.000626 | \n",
+ " 1.36 | \n",
+ " 0.859222 | \n",
+ " 0.048168 | \n",
+ " 0.933581 | \n",
+ " 2073.291693 | \n",
+ " 1.286 | \n",
+ " 0.000503 | \n",
+ " 1.38 | \n",
+ " 0.813647 | \n",
+ "
\n",
+ " \n",
+ "
\n",
+ "
10 rows × 21 columns
\n",
+ "
\n",
+ "
\n",
+ " \n",
+ " \n",
+ " | \n",
+ " 4.5e-07 | \n",
+ " 4.52e-07 | \n",
+ " 4.54e-07 | \n",
+ " 4.56e-07 | \n",
+ " 4.58e-07 | \n",
+ " 4.6e-07 | \n",
+ " 4.62e-07 | \n",
+ " 4.64e-07 | \n",
+ " 4.66e-07 | \n",
+ " 4.68e-07 | \n",
+ " ... | \n",
+ " 7e-07 | \n",
+ " 7.02e-07 | \n",
+ " 7.04e-07 | \n",
+ " 7.06e-07 | \n",
+ " 7.08e-07 | \n",
+ " 7.1e-07 | \n",
+ " 7.12e-07 | \n",
+ " 7.14e-07 | \n",
+ " 7.16e-07 | \n",
+ " 7.18e-07 | \n",
+ "
\n",
+ " \n",
+ " \n",
+ " \n",
+ " | 0 | \n",
+ " 0.052139 | \n",
+ " 0.075604 | \n",
+ " 0.105583 | \n",
+ " 0.117613 | \n",
+ " 0.128199 | \n",
+ " 0.137660 | \n",
+ " 0.147927 | \n",
+ " 0.154320 | \n",
+ " 0.162327 | \n",
+ " 0.170435 | \n",
+ " ... | \n",
+ " 0.505471 | \n",
+ " 0.506173 | \n",
+ " 0.506923 | \n",
+ " 0.508379 | \n",
+ " 0.509661 | \n",
+ " 0.509957 | \n",
+ " 0.511044 | \n",
+ " 0.511899 | \n",
+ " 0.512353 | \n",
+ " 0.513397 | \n",
+ "
\n",
+ " \n",
+ " | 1 | \n",
+ " 0.249317 | \n",
+ " 0.294264 | \n",
+ " 0.336788 | \n",
+ " 0.349672 | \n",
+ " 0.360439 | \n",
+ " 0.368425 | \n",
+ " 0.377414 | \n",
+ " 0.382268 | \n",
+ " 0.387558 | \n",
+ " 0.393322 | \n",
+ " ... | \n",
+ " 0.514512 | \n",
+ " 0.514646 | \n",
+ " 0.514554 | \n",
+ " 0.515599 | \n",
+ " 0.515551 | \n",
+ " 0.515706 | \n",
+ " 0.515688 | \n",
+ " 0.515692 | \n",
+ " 0.515558 | \n",
+ " 0.515568 | \n",
+ "
\n",
+ " \n",
+ " | 2 | \n",
+ " 0.015135 | \n",
+ " 0.024211 | \n",
+ " 0.038433 | \n",
+ " 0.044910 | \n",
+ " 0.051130 | \n",
+ " 0.056346 | \n",
+ " 0.062850 | \n",
+ " 0.067001 | \n",
+ " 0.072242 | \n",
+ " 0.077058 | \n",
+ " ... | \n",
+ " 0.369385 | \n",
+ " 0.370703 | \n",
+ " 0.372229 | \n",
+ " 0.373423 | \n",
+ " 0.375172 | \n",
+ " 0.375954 | \n",
+ " 0.376863 | \n",
+ " 0.377677 | \n",
+ " 0.379136 | \n",
+ " 0.380009 | \n",
+ "
\n",
+ " \n",
+ " | 3 | \n",
+ " 0.126252 | \n",
+ " 0.165786 | \n",
+ " 0.208991 | \n",
+ " 0.221679 | \n",
+ " 0.233692 | \n",
+ " 0.241259 | \n",
+ " 0.250605 | \n",
+ " 0.255267 | \n",
+ " 0.261180 | \n",
+ " 0.265604 | \n",
+ " ... | \n",
+ " 0.392249 | \n",
+ " 0.393718 | \n",
+ " 0.394732 | \n",
+ " 0.395270 | \n",
+ " 0.396005 | \n",
+ " 0.396826 | \n",
+ " 0.398016 | \n",
+ " 0.398674 | \n",
+ " 0.398748 | \n",
+ " 0.399876 | \n",
+ "
\n",
+ " \n",
+ " | 4 | \n",
+ " 0.032156 | \n",
+ " 0.048988 | \n",
+ " 0.073223 | \n",
+ " 0.082547 | \n",
+ " 0.092055 | \n",
+ " 0.099348 | \n",
+ " 0.107670 | \n",
+ " 0.112801 | \n",
+ " 0.118988 | \n",
+ " 0.125407 | \n",
+ " ... | \n",
+ " 0.435080 | \n",
+ " 0.436994 | \n",
+ " 0.439287 | \n",
+ " 0.440044 | \n",
+ " 0.441715 | \n",
+ " 0.442364 | \n",
+ " 0.443499 | \n",
+ " 0.444120 | \n",
+ " 0.446017 | \n",
+ " 0.447964 | \n",
+ "
\n",
+ " \n",
+ " | 5 | \n",
+ " 0.114718 | \n",
+ " 0.138649 | \n",
+ " 0.164083 | \n",
+ " 0.175158 | \n",
+ " 0.183507 | \n",
+ " 0.192087 | \n",
+ " 0.201838 | \n",
+ " 0.206713 | \n",
+ " 0.215729 | \n",
+ " 0.222988 | \n",
+ " ... | \n",
+ " 0.532057 | \n",
+ " 0.531468 | \n",
+ " 0.531926 | \n",
+ " 0.531250 | \n",
+ " 0.532234 | \n",
+ " 0.531852 | \n",
+ " 0.531362 | \n",
+ " 0.531587 | \n",
+ " 0.531800 | \n",
+ " 0.531074 | \n",
+ "
\n",
+ " \n",
+ " | 6 | \n",
+ " 0.075428 | \n",
+ " 0.092161 | \n",
+ " 0.111252 | \n",
+ " 0.120155 | \n",
+ " 0.126419 | \n",
+ " 0.133840 | \n",
+ " 0.141624 | \n",
+ " 0.146338 | \n",
+ " 0.153872 | \n",
+ " 0.160242 | \n",
+ " ... | \n",
+ " 0.491204 | \n",
+ " 0.491536 | \n",
+ " 0.491277 | \n",
+ " 0.491265 | \n",
+ " 0.491106 | \n",
+ " 0.490230 | \n",
+ " 0.490610 | \n",
+ " 0.490270 | \n",
+ " 0.490565 | \n",
+ " 0.490454 | \n",
+ "
\n",
+ " \n",
+ " | 7 | \n",
+ " 0.062895 | \n",
+ " 0.086219 | \n",
+ " 0.114160 | \n",
+ " 0.125318 | \n",
+ " 0.134272 | \n",
+ " 0.142476 | \n",
+ " 0.152179 | \n",
+ " 0.157698 | \n",
+ " 0.165529 | \n",
+ " 0.172328 | \n",
+ " ... | \n",
+ " 0.460213 | \n",
+ " 0.461519 | \n",
+ " 0.462544 | \n",
+ " 0.462657 | \n",
+ " 0.463183 | \n",
+ " 0.462922 | \n",
+ " 0.463642 | \n",
+ " 0.463763 | \n",
+ " 0.464488 | \n",
+ " 0.464716 | \n",
+ "
\n",
+ " \n",
+ " | 8 | \n",
+ " 0.002687 | \n",
+ " 0.009663 | \n",
+ " 0.027088 | \n",
+ " 0.034975 | \n",
+ " 0.043935 | \n",
+ " 0.050369 | \n",
+ " 0.058389 | \n",
+ " 0.063868 | \n",
+ " 0.069380 | \n",
+ " 0.074805 | \n",
+ " ... | \n",
+ " 0.335355 | \n",
+ " 0.337836 | \n",
+ " 0.339730 | \n",
+ " 0.342233 | \n",
+ " 0.344490 | \n",
+ " 0.346606 | \n",
+ " 0.348375 | \n",
+ " 0.350101 | \n",
+ " 0.351822 | \n",
+ " 0.353292 | \n",
+ "
\n",
+ " \n",
+ " | 9 | \n",
+ " 0.032634 | \n",
+ " 0.036502 | \n",
+ " 0.041451 | \n",
+ " 0.045510 | \n",
+ " 0.047463 | \n",
+ " 0.050845 | \n",
+ " 0.054574 | \n",
+ " 0.056508 | \n",
+ " 0.060603 | \n",
+ " 0.064730 | \n",
+ " ... | \n",
+ " 0.484122 | \n",
+ " 0.483369 | \n",
+ " 0.482663 | \n",
+ " 0.482151 | \n",
+ " 0.480324 | \n",
+ " 0.479257 | \n",
+ " 0.478852 | \n",
+ " 0.478074 | \n",
+ " 0.477372 | \n",
+ " 0.475579 | \n",
+ "
\n",
+ " \n",
+ "
\n",
+ "
10 rows × 135 columns
\n",
+ "