diff --git a/Modules/Biophotonics/python/optimization.py b/Modules/Biophotonics/python/optimization.py index 86efedb004..38a22738c7 100644 --- a/Modules/Biophotonics/python/optimization.py +++ b/Modules/Biophotonics/python/optimization.py @@ -1,152 +1,150 @@ # -*- coding: utf-8 -*- """ Created on Sat Feb 07 23:42:41 2015 @author: Seb """ -from scipy.interpolate import RectBivariateSpline import numpy as np import matplotlib.pyplot as plt -from scipy.optimize import minimize -from reflectanceError import ReflectanceError +from scipy.interpolate import RectBivariateSpline +from scipy.optimize import minimize +from reflectanceError import ReflectanceError from mpl_toolkits.mplot3d import Axes3D -from matplotlib import cm +from matplotlib import cm +import setupData dataFolder = "data/output/" # load data -trainingParameters = np.load(dataFolder + "2015February0511:02PMparamters2D.npy") -trainingReflectances = np.load(dataFolder + "2015February0511:02PMreflectances2D.npy") - -testParameters = np.load(dataFolder + "2015February1107:43PMparamtersRandom2D.npy") -testReflectances = np.load(dataFolder + "2015February1107:43PMreflectancesRandom2D.npy") +trainingParameters, trainingReflectances, testParameters, testReflectances = \ + setupData.setupTwoDimensionalData(dataFolder) BVFs = np.unique(trainingParameters[:,0]) Vss = np.unique(trainingParameters[:,1]) reflectanceGrid3D = np.reshape(trainingReflectances, (len(BVFs), len(Vss), trainingReflectances.shape[1])) functionToMinimize = ReflectanceError(BVFs, Vss, reflectanceGrid3D) absErrors = np.zeros_like(testParameters) for idx, (testParameter, testReflectance) in enumerate(zip(testParameters, testReflectances)): functionToMinimize.setReflectanceToMatch(testReflectance) minimization = minimize(functionToMinimize.f, [0.05, 0.3], method="Nelder-Mead") # interpolation extrapolates with constant values. Since the minimization function is # covex, we can just crop it to the bounds estimatedBVF= np.clip(minimization.x[0], min(BVFs), max(BVFs)) estimatedVs = np.clip(minimization.x[1], min(Vss), max(Vss)) absErrors[idx,:] = np.abs([estimatedBVF, estimatedVs] - testParameter) # plot print("error distribution BVF, Volume fraction") print("median: " + str(np.median(absErrors, axis=0))) print("lower quartile: " + str(np.percentile(absErrors, 25, axis=0))) print("higher quartile: " + str(np.percentile(absErrors, 75, axis=0))) if __name__ == "__main__": rbs = RectBivariateSpline(BVFs, Vss, np.reshape(trainingReflectances[:,0],(100,100))) rbsFrom3DGrid = RectBivariateSpline(BVFs, Vss, reflectanceGrid3D[:,:,0]) print("check if reflectance grid is build correctly. Reflectance at " + \ "(0.04258016, 0.29442218): " + str(rbs( 0.04258016 , 0.29442218, dx=0)) + "; expected " + str(654.222)) print("check if reflectance grid is build correctly. Reflectance at " + \ "(0.1, 0.60): " + str(rbs( 0.1 , 0.6, dx=0)) + "; expected " + str(818.945)) print("check if reflectance grid is build correctly. Reflectance at " + \ "(0.01, 0.60): " + str(rbs(0.01 , 0.6, dx=0)) + "; expected " + str(1120.118)) print("check if reflectance grid is build correctly. Reflectance at " + \ "(0.0112, 0.041): " + str(rbsFrom3DGrid(0.0112, 0.041, dx=0)) + "; expected " + str(rbs(0.0112, 0.041, dx=0))) reflectanceError = ReflectanceError(BVFs, Vss, reflectanceGrid3D) # reflectance at point (0.062, 0.149) exampleReflectance = [ 312.09769118, 327.67866117, 336.20970583, 343.10626114, 344.77202411, 322.68062559, 280.01722363, 265.45441892, 292.59760112, 294.58676191, 256.88475134, 296.10442177, 388.44056814, 422.99479049, 425.6602738 , 426.88353542, 427.09604971, 413.01675144, 424.50008822, 423.78125335, 422.71219033, 421.51283858, 421.32106797] reflectanceError.setReflectanceToMatch(exampleReflectance) print("check if reflectance error value is correctly computed. Expected: 0, actual: " + str(reflectanceError.f([6.191087469731736126e-02, 1.149252955326078018e-01]))) print("check if derivate of reflectance error value is correctly computed. Expected: close to [0,0], actual: " + str(reflectanceError.df([6.191087469731736126e-02, 1.149252955326078018e-01]))) # reflectance at point (0.01, 0.6) exampleReflectance2 = [1.120118468811987896e+03, 1.131318509608477598e+03, 1.137543168078165081e+03, 1.141159082538952589e+03, 1.139712500764461311e+03, 1.120915174518419235e+03, 1.081868496849364647e+03, 1.064766763752600582e+03, 1.084980179526470238e+03, 1.081898683840950525e+03, 1.038063983415024040e+03, 1.073340106470938281e+03, 1.117816020857113244e+03, 1.122940442828165033e+03, 1.116782763468411758e+03, 1.108377361568318520e+03, 1.104745713217586172e+03, 1.098370113333694690e+03, 1.093525893516002043e+03, 1.089713462981120756e+03, 1.085745981406935471e+03, 1.087101893762737973e+03, 1.083097058199290814e+03] reflectanceError.setReflectanceToMatch(exampleReflectance2) print("check if reflectance error value is correctly computed (second one). Expected: 0, actual: " + str(reflectanceError.f([0.01, 0.6]))) print("check if derivate of reflectance error value is correctly computed (second one). Expected: close to [0,0], actual: " + str(reflectanceError.df([0.01, 0.6]))) functionToMinimize.setReflectanceToMatch(trainingReflectances[0,:]) absError = minimize(functionToMinimize.f, [0.05, 0.3], method="Nelder-Mead") #absError = minimize(functionToMinimize.f, [0.05, 0.3], method="BFGS", jac=functionToMinimize.df) print(" function evaluation yields: " + str(absError.x) + " (unclipped) with absolute error: " + str(functionToMinimize.f(absError.x))) print(" expected evaluation result: " + str([0.01, 0.04]) + " real result error: " + str(functionToMinimize.f([0.01, 0.04]))) # check if grid interpolation looks good. #%% - grid_x, grid_y = np.mgrid[0.04:0.6:100j, 0.01:0.1:100j] + grid_x, grid_y = np.mgrid[0.01:0.1:100j, 0.04:0.6:100j] grid_z = rbs.ev(grid_x, grid_y) fig = plt.figure(0) ax = fig.gca(projection='3d') surf = ax.plot_surface(grid_x, grid_y, grid_z, cmap=cm.jet, linewidth=1, antialiased=True) - ax.scatter(testParameters[:,1], testParameters[:,0], testReflectances[:,0]) + ax.scatter(testParameters[:,0], testParameters[:,1], testReflectances[:,0]) ax.set_zlim3d(np.min(grid_z), np.max(grid_z)) fig.colorbar(surf) plt.show() # plt.imshow(grid_z.T, extent=(0.04,0.6,0.01,0.1)) diff --git a/Modules/Biophotonics/python/randomForest.py b/Modules/Biophotonics/python/randomForest.py index 318d3b786e..5e84957ab3 100644 --- a/Modules/Biophotonics/python/randomForest.py +++ b/Modules/Biophotonics/python/randomForest.py @@ -1,68 +1,52 @@ # -*- coding: utf-8 -*- """ Created on Fri Feb 6 10:49:45 2015 @author: wirkert """ import numpy as np -import monteCarloHelper as mch import matplotlib.pyplot as plt -from sklearn import tree +from sklearn import tree from sklearn.ensemble import RandomForestRegressor +import setupData -# todo we: -# 2. optimization -# what to find out in this study: + +# additional things that could be checked in this study: # 1. band selection results # 2. effect of normalizations # 3. parameter study # 4. optimal image quotient #%% initialize # the folder with the reflectance spectra dataFolder = 'data/output/' # load data -trainingParameters = np.load(dataFolder + "2015February0511:02PMparamters2D.npy") -trainingReflectances = np.load(dataFolder + "2015February0511:02PMreflectances2D.npy") - -testParameters = np.load(dataFolder + "2015February1107:43PMparamtersRandom2D.npy") -testReflectances = np.load(dataFolder + "2015February1107:43PMreflectancesRandom2D.npy") - +trainingParameters, trainingReflectances, testParameters, testReflectances = \ + setupData.setupTwoDimensionalData(dataFolder) -# normalize data -#trainingReflectances = trainingReflectances[:,[0, 1, 4, 9, 11, 13, 18, 22]] -#trainingReflectances = trainingReflectances[:,[ 0, 3, 6, 9, 12, 15, 18, 22]] -#trainingReflectances = trainingReflectances[:,[ 0, 1, 2, 3, 4, 5, 6, 7]] -trainingReflectances = mch.normalizeImageQuotient(trainingReflectances, iqBand=4) - #%% train forest rf = RandomForestRegressor(50, max_depth=8, max_features="log2", n_jobs=5) rf.fit(trainingReflectances, trainingParameters) #with open("iris.dot", 'w') as f: # f = tree.export_graphviz(rf, out_file=f) #%% test -#testReflectances = testReflectances[:,[0, 1, 4, 9, 11, 13, 18, 22]] -#testReflectances = testReflectances[:,[ 0, 3, 6, 9, 12, 15, 18, 22]] -#testReflectances = testReflectances[:,[0, 1, 2, 3, 4, 5, 6, 7]] -testReflectances = mch.normalizeImageQuotient(testReflectances, iqBand=4) - # predict test reflectances and get absolute errors. absErrors = np.abs(rf.predict(testReflectances) - testParameters) print("error distribution BVF, Volume fraction") print("median: " + str(np.median(absErrors, axis=0))) print("lower quartile: " + str(np.percentile(absErrors, 25, axis=0))) print("higher quartile: " + str(np.percentile(absErrors, 75, axis=0)))