Added preprocessing.

Added plot for ML result for each predicted dataset.

main.py

@@ -37,17 +37,23 @@ else:
     # Settings
     number_of_tests = 5
     rho_spread_factor = 1.5
-    rho_max = 20
+    rho_max = 10
     layers_min = 1
-    layers_max = 3
-    min_depth = 2
-    max_depth = 10
+    layers_max = 1
+    min_depth = 4
+    max_depth = 8
 
     # Generate parameters for tests
-    tests_horizontal = slopestabilitytools.model_params(number_of_tests,
-                                                        rho_spread_factor, rho_max,
-                                                        layers_min, layers_max,
-                                                        min_depth, max_depth)
+    # tests_horizontal = slopestabilitytools.model_params(number_of_tests,
+    #                                                     rho_spread_factor, rho_max,
+    #                                                     layers_min, layers_max,
+    #                                                     min_depth, max_depth)
+
+    tests_horizontal = {'hor_1': {'layer_n': 1, 'rho_values': [[1, 5], [2, 15]], 'layers_pos': np.array([-5])},
+                        'hor_2': {'layer_n': 1, 'rho_values': [[1, 5], [2, 50]], 'layers_pos': np.array([-5])},
+                        'hor_3': {'layer_n': 1, 'rho_values': [[1, 15], [2, 20]], 'layers_pos': np.array([-8])},
+                        'hor_4': {'layer_n': 1, 'rho_values': [[1, 5], [2, 10]], 'layers_pos': np.array([-3])},
+                        'hor_5': {'layer_n': 1, 'rho_values': [[1, 5], [2, 25]], 'layers_pos': np.array([-3])}}
 
     #  Create models and invert them
     test_results = {}
@@ -61,8 +67,8 @@ else:
         slopestabilitytools.plot_and_save(test_name, test_results[test_name], 'Test: ' + test_name)
 
 if not create_new_data_only:
-    for test_name in test_results.keys():
-        slopestabilitytools.plot_and_save(test_name, test_results[test_name], 'Test: ' + test_name)
+    #for test_name in test_results.keys():
+        #slopestabilitytools.plot_and_save(test_name, test_results[test_name], 'Test: ' + test_name)
 
     ml_results = slopestabilityML.run_all_tests(test_results)
     # svm_accuracy_score, svm_accuracy_labels = slopestabilityML.svm_run(test_results)

requirements.txt

@@ -2,4 +2,5 @@ numpy~=1.18.5
 matplotlib~=3.3.2
 pandas~=1.1.5
 pygimli~=1.1.0
-scikit-learn~=0.23.2
\ No newline at end of file
+scikit-learn~=0.23.2
+scipy~=1.5.2
\ No newline at end of file

slopestabilityML/preprocess_data.py

@@ -6,15 +6,13 @@ Created on 19.01.2021
 @author: Feliks Kiszkurno
 """
 
-from sklearn.preprocessing import StandardScaler, OneHotEncoder
-from sklearn.compose import ColumnTransformer
+import pandas as pd
 
 
 def preprocess_data(data_set):
 
-    x_train = data_set.drop(['Z', 'INM', 'CLASS'], axis='columns')
-    y_train = data_set['CLASS']
-
-    num_feat = []
+    x_train = data_set.drop(['X', 'Y', 'Z', 'INM'], axis='columns')
+    #y_train = data_set.drop(['X', 'Y', 'Z', 'INM'], axis='columns')
+    y_train = pd.DataFrame(data_set['CLASS'])
 
     return x_train, y_train

slopestabilityML/run_classification.py

@@ -8,6 +8,14 @@ Created on 19.01.2021
 
 import slopestabilityML
 import numpy as np
+import pandas as pd
+from sklearn.preprocessing import StandardScaler, OneHotEncoder
+from sklearn.compose import ColumnTransformer
+from sklearn.pipeline import make_pipeline
+import matplotlib.pyplot as plt
+from matplotlib import ticker
+from scipy import interpolate
+import slopestabilitytools
 
 
 def run_classification(test_training, test_prediction, test_results, clf):
@@ -15,22 +23,98 @@ def run_classification(test_training, test_prediction, test_results, clf):
     accuracy_score = []
     accuracy_labels = []
 
+    num_feat = ['RES', 'SEN']
+    cat_feat = ['CLASS']
+
+    cat_lab = [0, 1]
+
+    num_trans = StandardScaler()
+    cat_trans = OneHotEncoder(categories=[cat_lab])
+
+    preprocessor = ColumnTransformer(transformers=[('num', num_trans, num_feat),
+                                                   ('cat', cat_trans, cat_feat)])
+
+    clf_pipeline_UM = make_pipeline(preprocessor, clf)
+
     for test_name in test_training:
         # Prepare data
         x_train, y_train = slopestabilityML.preprocess_data(test_results[test_name])
 
         # Train classifier
-        clf.fit(x_train, y_train)
+        # print(type(x_train))
+        # print(type(y_train))
+        clf_pipeline_UM.fit(x_train, y_train)
 
     # Predict with classifier
     for test_name_pred in test_prediction:
         # Prepare data
         x_question, y_answer = slopestabilityML.preprocess_data(test_results[test_name_pred])
 
-        y_pred = clf.predict(x_question)
+        # y_pred = clf_pipeline_UM.score(x_question, y_answer)
+        y_pred = clf_pipeline_UM.predict(x_question)
+        # print(y_pred)
+        score = clf_pipeline_UM.score(x_question, y_answer)
+        # print('score: '+str(score))
+
+        x = test_results[test_name_pred]['X']
+        y = test_results[test_name_pred]['Y']
+        class_in = test_results[test_name]['CLASS']
+        class_out = y_pred
+        x_min = np.min(x)
+        x_max = np.max(x)
+        x_n = len(x)
+
+        y_min = np.min(y)
+        y_max = np.max(y)
+        y_start = y_max
+        y_end = y_min
+        y_n = len(y)
+
+        xi = np.linspace(x_min, x_max, x_n)
+        yi = np.linspace(y_start, y_end, y_n)
+        xx, yy = np.meshgrid(xi, yi)
+        class_in_i = interpolate.griddata((x, y), class_in, (xx, yy), method='nearest')
+        class_out_i = interpolate.griddata((x, y), class_out, (xx, yy), method='nearest')
+        class_diff = np.zeros_like(class_out_i)
+        class_diff[np.where(class_in_i == class_out_i)] = 1
+        cb = []
+
+        fig, _ax = plt.subplots(nrows=3, ncols=1)
+        ax = _ax.flatten()
+
+        fig.suptitle(test_name_pred)
+        fig.subplots_adjust(hspace=0.8)
+
+        im0 = ax[0].contourf(xi, yi, class_in_i)
+        ax[0].set_title('Input classes')
+        ax[0] = slopestabilitytools.set_labels(ax[0])
+        cb.append(plt.colorbar(im0, ax=ax[0], label='Class'))  # , shrink=0.9)
+        tick_locator = ticker.MaxNLocator(nbins=4)
+        cb[0].locator = tick_locator
+        cb[0].update_ticks()
+
+        im1 = ax[1].contourf(xi, yi, class_out_i)
+        ax[1].set_title('Result of classification')
+        ax[1] = slopestabilitytools.set_labels(ax[1])
+        cb.append(plt.colorbar(im1, ax=ax[1], label='Class'))  # , shrink=0.9)
+        tick_locator = ticker.MaxNLocator(nbins=4)
+        cb[1].locator = tick_locator
+        cb[1].update_ticks()
+
+        im2 = ax[2].contourf(xi, yi, class_diff)
+        ax[2].set_title('Difference')
+        ax[2] = slopestabilitytools.set_labels(ax[2])
+        cb.append(plt.colorbar(im2, ax=ax[2], label='Resistivity [om]'))  # , shrink=0.9)
+        tick_locator = ticker.MaxNLocator(nbins=4)
+        cb[2].locator = tick_locator
+        cb[2].update_ticks()
+
+        fig.savefig('results/figures/eps/{}_ML_class_res.eps'.format(test_name_pred))
+        fig.savefig('results/figures/png/{}_ML_class_res.png'.format(test_name_pred))
+        fig.savefig('results/figures/pdf/{}_ML_class_res.pdf'.format(test_name_pred))
 
         # Evaluate result
-        accuracy_score.append(len(np.where(y_pred == y_answer)) / len(y_answer) * 100)
+        accuracy_score.append(len(np.where(y_pred == y_answer.to_numpy())) / len(y_answer.to_numpy()) * 100)
         accuracy_labels.append(test_name_pred)
 
     return accuracy_labels, accuracy_score

slopestabilitytools/datamanagement/import_tests.py

@@ -17,7 +17,7 @@ def import_tests():
     #print(test_names)
 
     for test_name in test_names:
-        test_result_curr = pd.read_csv('results/results/' + test_name + '.csv')
+        test_result_curr = pd.read_csv('results/results/' + test_name + '.csv', index_col=0)
         test_results.update({test_name: test_result_curr})
 
     return test_results

slopestabilitytools/model_params.py

@@ -31,7 +31,10 @@ def model_params(n_of_tests, rho_spread, rho_max, layers_n_min, layers_n_max, de
     for test_id in range(n_of_tests):
 
         test_names[test_id] = 'hor_{}'.format(str(test_id + 1))
-        test_n_layers[test_names[test_id]] = np.random.randint(layers_n_min, layers_n_max)
+        if layers_n_max == layers_n_min:
+            test_n_layers[test_names[test_id]] = layers_n_max
+        else:
+            test_n_layers[test_names[test_id]] = np.random.randint(layers_n_min, layers_n_max)
         rho_temp = []
         rho_used = []
         layer_pos_temp = []

slopestabilitytools/plot_and_save.py

@@ -22,10 +22,22 @@ def plot_and_save(test_name, test_result, plot_title):
     res = test_result['RES']
     cov = test_result['SEN']
 
-    xx, yy = np.meshgrid(x, y)
-    inm_i = interpolate.griddata((x, y), inm, (xx, yy), method='linear')
-    res_i = interpolate.griddata((x, y), res, (xx, yy), method='linear')
-    cov_i = interpolate.griddata((x, y), cov, (xx, yy), method='linear')
+    x_min = np.min(x)
+    x_max = np.max(x)
+    x_n = len(x)
+
+    y_min = np.min(y)
+    y_max = np.max(y)
+    y_start = y_max
+    y_end = y_min
+    y_n = len(y)
+
+    xi = np.linspace(x_min, x_max, x_n)
+    yi = np.linspace(y_start, y_end, y_n)
+    xx, yy = np.meshgrid(xi, yi)
+    inm_i = interpolate.griddata((x, y), inm, (xx, yy), method='nearest')
+    res_i = interpolate.griddata((x, y), res, (xx, yy), method='nearest')
+    cov_i = interpolate.griddata((x, y), cov, (xx, yy), method='nearest')
 
     print('plot_and_save')
 
@@ -38,7 +50,7 @@ def plot_and_save(test_name, test_result, plot_title):
     fig.suptitle(plot_title)
     fig.subplots_adjust(hspace=0.8)
 
-    im0 = ax[0].contourf(xx, yy, inm_i)
+    im0 = ax[0].contourf(xi, yi, inm_i)
     ax[0].set_title('Input model')
     ax[0] = slopestabilitytools.set_labels(ax[0])
     cb.append(plt.colorbar(im0, ax=ax[0], label='Resistivity [om]'))#, shrink=0.9)
@@ -46,7 +58,7 @@ def plot_and_save(test_name, test_result, plot_title):
     cb[0].locator = tick_locator
     cb[0].update_ticks()
 
-    im1 = ax[1].contourf(xx, yy, res_i)
+    im1 = ax[1].contourf(xi, yi, res_i)
     ax[1].set_title('Result')
     ax[1] = slopestabilitytools.set_labels(ax[1])
     cb.append(plt.colorbar(im1, ax=ax[1], label='Resistivity [om]'))#, shrink=0.9)
@@ -54,7 +66,7 @@ def plot_and_save(test_name, test_result, plot_title):
     cb[1].locator = tick_locator
     cb[1].update_ticks()
 
-    im2 = ax[2].contourf(xx, yy, inm_i-res_i)
+    im2 = ax[2].contourf(xi, yi, inm_i-res_i)
     ax[2].set_title('Difference')
     ax[2] = slopestabilitytools.set_labels(ax[2])
     cb.append(plt.colorbar(im2, ax=ax[2], label='Resistivity [om]'))#, shrink=0.9)
@@ -72,7 +84,7 @@ def plot_and_save(test_name, test_result, plot_title):
 
     fig.suptitle(plot_title+' coverage')
 
-    im0 = ax_cov.contourf(xx, yy, cov_i)
+    im0 = ax_cov.contourf(xi, yi, cov_i)
     ax_cov.set_title(plot_title+' coverage')
     ax_cov = slopestabilitytools.set_labels(ax_cov)
     cb_cov = plt.colorbar(im0, ax=ax_cov, label='Coverage')  # , shrink=0.9)

slostabcreatedata/create_data.py

@@ -16,8 +16,10 @@ import pygimli.physics.ert as ert
 
 
 def create_data(test_name, test_config, max_depth):
-    world_boundary_v = [-20 * max_depth, 0]  # [right, left border] relatively to the middle
-    world_boundary_h = [20 * max_depth, -4 * max_depth]  # [top, bottom border]
+    world_boundary_v = [-20 * max_depth, 0]  # [NW edge] relatively to the middle
+    world_boundary_h = [20 * max_depth, -5 * max_depth]  # [SE edge]
+    # world_boundary_v = [-500, 0]  # [right, left border] relatively to the middle
+    # world_boundary_h = [500, -100]  # [top, bottom border]
 
     test_results = {}
 
@@ -29,7 +31,12 @@ def create_data(test_name, test_config, max_depth):
 
     geometry = world  # +block
 
-    measurement_scheme = ert.createERTData(elecs=np.linspace(start=-5*max_depth, stop=5*max_depth, num=6*max_depth+1),
+    fig_geometry, ax_geometry = plt.subplots(1)
+    pg.show(geometry, ax=ax_geometry)
+    ax_geometry.set_title('1 Geometry of the model')
+    fig_geometry.savefig('results/figures/png/'+test_name+'_1_geometry.png')
+
+    measurement_scheme = ert.createERTData(elecs=np.linspace(start=-8*max_depth, stop=8*max_depth, num=16*max_depth+1),
                                            schemeName='dd')
 
     for electrode in measurement_scheme.sensors():
@@ -40,6 +47,11 @@ def create_data(test_name, test_config, max_depth):
 
     resistivity_map = test_config['rho_values']  # [0]
 
+    fig_model, ax_model = plt.subplots(1)
+    pg.show(mesh, data=resistivity_map, label=pg.unit('res'), showMesh=True, ax=ax_model)
+    ax_model.set_title('2 Mesh and resistivity distribution')
+    fig_model.savefig('results/figures/png/' + test_name + '_2_meshdist.png')
+
     input_model = pg.solver.parseMapToCellArray(resistivity_map, mesh)  # rename to input_mesh
 
     # INPUT MODEL - SUBSURFACE MODEL END ###
@@ -53,16 +65,26 @@ def create_data(test_name, test_config, max_depth):
 
     # RUN INVERSION #
     k0 = pg.physics.ert.createGeometricFactors(data)
-    model_inverted = ert_manager.invert(data=data, lam=20, paraDX=0.1, paraMaxCellSize=2, paraDepth=max_depth,
+    model_inverted = ert_manager.invert(data=data, lam=20, paraDX=0.25, paraMaxCellSize=2, paraDepth=2*max_depth,
                                         quality=34, zPower=0.4)
 
     result = ert_manager.inv.model
     result_array = result.array()
 
+    fig_result, ax_result = plt.subplots(1)
+    pg.show(ert_manager.paraDomain, result, label=pg.unit('res'), showMesh=True, ax=ax_result)
+    ax_result.set_title('3 Result')
+    fig_result.savefig('results/figures/png/' + test_name + '_3_result.png')
+
     input_model2 = pg.interpolate(srcMesh=mesh, inVec=input_model, destPos=ert_manager.paraDomain.cellCenters())
 
     input_model2_array = input_model2.array()
 
+    fig_input, ax_input = plt.subplots(1)
+    pg.show(ert_manager.paraDomain, input_model2, label=pg.unit('res'), showMesh=True, ax=ax_input)
+    ax_input.set_title('4 Model on inv mesh')
+    fig_input.savefig('results/figures/png/' + test_name + '_4_modelinv.png')
+
     # Create classes labels
     classes = []
     resistivity_values = []